use strict;
use warnings;
my $F = [qw(F D)];
pp_add_exported('', 'random_cluster', 'iv_cluster');
pp_addpm({At=>'Top'}, <<'EOD');
use strict;
use warnings;
use Carp;
use PDL::LiteF;
use PDL::Stats::Basic;
=head1 NAME
PDL::Stats::Kmeans -- classic k-means cluster analysis
=head1 DESCRIPTION
Assumes that we have data pdl dim [observation, variable] and the goal is to put observations into clusters based on their values on the variables. The terms "observation" and "variable" are quite arbitrary but serve as a reminder for "that which is being clustered" and "that which is used to cluster".
The terms FUNCTIONS and METHODS are arbitrarily used to refer to methods that are broadcastable and methods that are non-broadcastable, respectively.
=head1 SYNOPSIS
Implement a basic k-means procedure,
use PDL::LiteF;
use PDL::Stats;
my ($data, $idv, $ido) = rtable( $file );
# or generate random data:
$data = grandom(200, 2); # two vars as below
my ($cluster, $centroid, $ss_centroid, $cluster_last);
# start out with 8 random clusters
$cluster = random_cluster( $data->dim(0), 8 );
# iterate to minimize total ss
# stop when no more changes in cluster membership
do {
$cluster_last = $cluster;
($centroid, $ss_centroid) = $data->centroid( $cluster );
$cluster = $data->assign( $centroid );
} while sum(abs($cluster - $cluster_last)) > 0;
or, use the B<kmeans> function provided here,
my %k = $data->kmeans( \%opt );
print "$_\t$k{$_}\n" for sort keys %k;
plot the clusters if there are only 2 vars in $data,
use PDL::Graphics::Simple;
my ($win, $c);
$win = pgswin();
$win->plot(map +(with=>'points', $data->dice_axis(0,which($k{cluster}->(,$_)))->dog), 0 .. $k{cluster}->dim(1)-1);
=cut
EOD
pp_addhdr('
#include <math.h>
#include <stdlib.h>
#include <time.h>
'
);
pp_addpm pp_line_numbers(__LINE__, <<'EOD' );
=head2 random_cluster
=for ref
Creates masks for random mutually exclusive clusters. Accepts two
parameters, num_obs and num_cluster. Extra parameter turns into extra
dim in mask. May loop a long time if num_cluster approaches num_obs
because empty cluster is not allowed.
=for usage
my $cluster = random_cluster( $num_obs, $num_cluster );
=cut
# can't be called on pdl
sub random_cluster {
my ($obs, $clu, @extra) = @_;
# extra param in @_ made into extra dim
my $cluster = zeroes short(), @_;
do {
(random($obs, @extra) * $obs)->_random_cluster($cluster);
} while (PDL::any $cluster->sumover == 0 );
$cluster;
}
EOD
pp_def('_random_cluster',
Pars => 'rand_data(o); [o]b(o,c)',
GenericTypes => ['S'],
RedoDimsCode => 'if ($SIZE(c) > $SIZE(o)) $CROAK("more cluster than obs!");',
Code => '
PDL_Indx nc = $SIZE(c);
broadcastloop %{
loop (o) %{
PDL_Indx cl = $rand_data() % nc;
loop (c) %{
$b() = (c == cl) ? 1 : 0;
%}
%}
%}
',
Doc => undef,
);
pp_def('which_cluster',
Pars => 'short a(o,c); indx [o]b(o)',
GenericTypes => [qw(U L)],
HandleBad => 1,
Code => '
loop(o) %{
PDL_Indx cl=-1;
loop(c) %{
PDL_IF_BAD(if ($ISBAD($a())) continue;,)
if (!$a()) continue;
cl = c;
break;
%}
PDL_IF_BAD(if (cl==-1) { $SETBAD(b()); continue; },)
$b() = cl;
%}
',
Doc => 'Given cluster mask dim [obs x clu], returns the cluster index to which an obs belong.
Does not support overlapping clusters. If an obs has TRUE value for multiple clusters, the returned index is the first cluster the obs belongs to. If an obs has no TRUE value for any cluster, the return val is set to -1 or BAD if the input mask has badflag set.
Usage:
# create a cluster mask dim [obs x clu]
pdl> p $c_mask = iv_cluster [qw(a a b b c c)]
[
[1 1 0 0 0 0]
[0 0 1 1 0 0]
[0 0 0 0 1 1]
]
# get cluster membership list dim [obs]
pdl> p $ic = $c_mask->which_cluster
[0 0 1 1 2 2]
',
);
pp_def('assign',
Pars => 'data(o,v); centroid(c,v); short [o]cluster(o,c)',
GenericTypes => $F,
HandleBad => 1,
Code => '
PDL_Indx cl = 0;
loop (o) %{
$GENERIC(centroid) ssmin = -1;
loop (c) %{
$GENERIC(centroid) ssc = 0;
PDL_Indx nvc = 0;
loop (v) %{
PDL_IF_BAD(if ($ISBAD($data()) || $ISBAD($centroid())) continue;,)
$GENERIC() diff = $data() - $centroid();
ssc += diff*diff;
nvc ++;
%}
if (nvc) {
ssc /= nvc;
}
else {
/* taking advantage of the fact that 1st valid ssmin takes precedence */
/* so ssc has no effect if there is already ssmin. or it is -1 */
ssc = ssmin;
}
/* notice that if multiple ssc == ssmin the 1st is taken as cluster */
if (ssmin < 0 || ssmin > ssc) {
cl = c;
ssmin = ssc;
}
%}
loop (c) %{
PDL_IF_BAD(if (ssmin < 0) { $SETBAD($cluster()); continue; },)
$cluster() = (c == cl)? 1 : 0;
%}
%}
',
Doc => '
=for ref
Takes data pdl dim [obs x var] and centroid pdl dim [cluster x var] and returns mask dim [obs x cluster] to cluster membership. An obs is assigned to the first cluster with the smallest distance (ie sum squared error) to cluster centroid. With bad value, obs is assigned by smallest mean squared error across variables.
=for usage
pdl> p $centroid = xvals 2, 3
[
[0 1]
[0 1]
[0 1]
]
pdl> p $b = qsort( random 4, 3 )
[
[0.022774068 0.032513883 0.13890034 0.30942479]
[ 0.16943853 0.50262636 0.56251531 0.7152271]
[ 0.23964483 0.59932745 0.60967495 0.78452117]
]
# notice that 1st 3 obs in $b are on average closer to 0
# and last obs closer to 1
pdl> p $b->assign( $centroid )
[
[1 1 1 0] # cluster 0 membership
[0 0 0 1] # cluster 1 membership
]
',
);
pp_def('centroid',
Pars => 'data(o,v); cluster(o,c); [o]m(c,v); [o]ss(c,v)',
GenericTypes => $F,
HandleBad => 1,
Code => '
$GENERIC(m) s[ $SIZE(c) ][ $SIZE(v) ], s2[ $SIZE(c) ][ $SIZE(v) ];
PDL_Indx n[ $SIZE(c) ]PDL_IF_BAD([ $SIZE(v) ],);
loop (c) %{
loop (v) %{
s[c][v] = 0.0;
s2[c][v] = 0.0;
PDL_IF_BAD(n[c][v] = 0;,)
%}
PDL_IF_BAD(,n[c] = 0;)
loop (o) %{
PDL_IF_BAD(if ($ISBAD($cluster())) continue;,)
if (!$cluster()) continue;
PDL_IF_BAD(,n[c]++;)
loop (v) %{
PDL_IF_BAD(if ($ISBAD($data())) continue;,)
s[c][v] += $data();
s2[c][v] += $data()*$data();
PDL_IF_BAD(n[c][v]++;,)
%}
%}
loop (v) %{
if (n[c]PDL_IF_BAD([v],)) {
$m() = s[c][v] / n[c]PDL_IF_BAD([v],);
$ss() = s2[c][v] PDL_IF_BAD(/ n[c][v],) - pow(s[c][v] / n[c]PDL_IF_BAD([v],), 2) PDL_IF_BAD(,* n[c]);
}
else {
$m() = 0;
$ss() = 0;
/* $CROAK("please make sure there is no empty cluster!"); */
}
%}
%}
',
Doc => '
=for ref
Takes data dim [obs x var] and mask dim [obs x cluster], returns mean and ss (ms when data contains bad values) dim [cluster x var], using data where mask == 1. Multiple cluster membership for an obs is okay. If a cluster is empty all means and ss are set to zero for that cluster.
=for usage
# data is 10 obs x 3 var
pdl> p $d = sequence 10, 3
[
[ 0 1 2 3 4 5 6 7 8 9]
[10 11 12 13 14 15 16 17 18 19]
[20 21 22 23 24 25 26 27 28 29]
]
# create two clusters by value on 1st var
pdl> p $a = $d( ,(0)) <= 5
[1 1 1 1 1 1 0 0 0 0]
pdl> p $b = $d( ,(0)) > 5
[0 0 0 0 0 0 1 1 1 1]
pdl> p $c = cat $a, $b
[
[1 1 1 1 1 1 0 0 0 0]
[0 0 0 0 0 0 1 1 1 1]
]
pdl> p $d->centroid($c)
# mean for 2 cluster x 3 var
[
[ 2.5 7.5]
[12.5 17.5]
[22.5 27.5]
]
# ss for 2 cluster x 3 var
[
[17.5 5]
[17.5 5]
[17.5 5]
]
',
);
pp_addpm pp_line_numbers(__LINE__, <<'EOD');
sub _scree_ind {
# use as scree cutoff the point with max distance to the line formed
# by the 1st and last points in $self
# it's a heuristic--whether we can get "good" results depends on
# the number of components in $self.
my ($self) = @_;
$self = $self->squeeze;
$self->ndims > 1 and
croak "1D pdl only please";
my $a = zeroes 2, $self->nelem;
$a->slice('(0)') .= sequence $self->nelem;
$a->slice('(1)') .= $self;
my $d = _d_point2line( $a, $a->slice(':,(0)'), $a->slice(':,(-1)') );
return $d->maximum_ind;
}
sub _d_point2line {
my ($self, $p1, $p2) = @_;
for ($self, $p1, $p2) {
$_->dim(0) != 2 and
carp "point pdl dim(0) != 2";
}
return _d_p2l( $self->mv(0,-1)->dog, $p1->mv(0,-1)->dog, $p2->mv(0,-1)->dog );
}
EOD
pp_def('_d_p2l',
Pars => 'xc(); yc(); xa(); ya(); xb(); yb(); [o]d()',
GenericTypes => $F,
HandleBad => 1,
Code => '
PDL_IF_BAD(if ($ISBAD(xc()) || $ISBAD(yc()) || $ISBAD(xa()) || $ISBAD(ya()) || $ISBAD(xb()) || $ISBAD(yb()) ) { $SETBAD(d()); continue; },)
$GENERIC(d) xba, yba;
xba = $xb() - $xa();
yba = $yb() - $ya();
$d() = fabs( xba * ($ya() - $yc()) - ($xa() - $xc()) * yba )
/ hypot( xba, yba );
',
Doc => undef,
);
pp_addpm pp_line_numbers(__LINE__, <<'EOD');
=head2 kmeans
=for ref
Implements classic k-means cluster analysis.
=for example
$data = grandom(200, 2); # two rows = two dimensions
%k = $data->kmeans; # use default of 3 clusters
print "$_\t$k{$_}\n" for sort keys %k;
$w->plot(
(map +(with=>'points', style=>$_+1, ke=>"Cluster ".($_+1),
$data->dice_axis(0,which($k{cluster}->slice(",$_")))->dog),
0 .. $k{cluster}->dim(1)-1),
(map +(with=>'circles', style=>$_+1, ke=>"Centroid ".($_+1), $k{centroid}->slice($_)->dog, 0.1),
0 .. $k{centroid}->dim(0)-1),
{le=>'tr'},
);
Given a number of observations with values on a set of variables,
kmeans puts the observations into clusters that maximizes within-cluster
similarity with respect to the variables. Tries several different random
seeding and clustering in parallel. Stops when cluster assignment of the
observations no longer changes. Returns the best result in terms of R2
from the random-seeding trials.
Instead of random seeding, kmeans also accepts manual seeding. This is
done by providing a centroid to the function, in which case clustering
will proceed from the centroid and there is no multiple tries.
There are two distinct advantages from seeding with a centroid compared to
seeding with predefined cluster membership of a subset of the observations
ie "seeds":
=over
=item *
a centroid could come from a previous study with a different set of observations;
=item *
a centroid could even be "fictional", or in more proper parlance,
an idealized prototype with respect to the actual data. For example,
if there are 10 person's ratings of 1 to 5 on 4 movies, ie a ratings
pdl of dim [10 obs x 4 var], providing a centroid like
[
[5 0 0 0]
[0 5 0 0]
[0 0 5 0]
[0 0 0 5]
]
will produce 4 clusters of people with each cluster favoring a different
one of the 4 movies. Clusters from an idealized centroid may not give the
best result in terms of R2, but they sure are a lot more interpretable.
=back
If clustering has to be done from predefined clusters of seeds, simply
calculate the centroid using the B<centroid> function and feed it
to kmeans,
my ($centroid, $ss) = $rating($iseeds, )->centroid( $seeds_cluster );
my %k = $rating->kmeans( { CNTRD=>$centroid } );
kmeans supports bad value*.
=for options
Default options (case insensitive):
V => 1, # prints simple status
FULL => 0, # returns results for all seeding trials
CNTRD => PDL->null, # optional. pdl [clu x var]. disables next 3 opts
NTRY => 5, # num of random seeding trials
NSEED => 1000, # num of initial seeds, use NSEED up to max obs
NCLUS => 3, # num of clusters
=for usage
Usage:
# suppose we have 4 person's ratings on 5 movies
pdl> p $rating = ceil( random(4, 5) * 5 )
[
[3 2 2 3]
[2 4 5 4]
[5 3 2 3]
[3 3 1 5]
[4 3 3 2]
]
# we want to put the 4 persons into 2 groups
pdl> %k = $rating->kmeans( {NCLUS=>2} )
# by default prints back options used
# as well as info for all tries and iterations
CNTRD => Null
FULL => 0
NCLUS => 3
NSEED => 4
NTRY => 5
V => 1
ss total: 20.5
iter 0 R2 [0.024390244 0.024390244 0.26829268 0.4796748 0.4796748]
iter 1 R2 [0.46341463 0.46341463 0.4796748 0.4796748 0.4796748]
pdl> p "$_\t$k{$_}\n" for sort keys %k
R2 0.479674796747968
centroid # mean ratings for 2 group x 5 movies
[
[ 3 2.3333333]
[ 2 4.3333333]
[ 5 2.6666667]
[ 3 3]
[ 4 2.6666667]
]
cluster # 4 persons' membership in two groups
[
[1 0 0 0]
[0 1 1 1]
]
n [1 3] # cluster size
ss
[
[ 0 0.66666667]
[ 0 0.66666667]
[ 0 0.66666667]
[ 0 8]
[ 0 0.66666667]
]
Now, for the valiant, kmeans is broadcastable. Say you gathered 10
persons' ratings on 5 movies from 2 countries, so the data is dim
[10,5,2], and you want to put the 10 persons from each country into
3 clusters, just specify NCLUS => [3,1], and there you have it. The
key is for NCLUS to include $data->ndims - 1 numbers. The 1 in [3,1]
turns into a dummy dim, so the 3-cluster operation is repeated on both
countries. Similarly, when seeding, CNTRD needs to have ndims that at
least match the data ndims. Extra dims in CNTRD will lead to broadcasting
(convenient if you want to try out different centroid locations,
for example, but you will have to hand pick the best result). See
F<t/kmeans.t> for examples w 3D and 4D data.
*With bad value, R2 is based on average of variances instead of sum squared error.
=cut
*kmeans = \&PDL::kmeans;
sub PDL::kmeans {
my ($self, $opt) = @_;
my %opt = (
V => 1, # prints simple status
FULL => 0, # returns results for all seeding trials
CNTRD => PDL->null, # optional. pdl [clu x var]. disables next 3 opts
NTRY => 5, # num of random seeding trials
NSEED => 1000, # num of initial seeds, use NSEED up to max obs
NCLUS => 3, # num of clusters
);
if ($opt) { $opt{uc $_} = $opt->{$_} for keys %$opt; }
if (defined($opt{CNTRD}) and $opt{CNTRD}->nelem) {
$opt{NTRY} = 1;
$opt{NSEED} = $self->dim(0);
$opt{NCLUS} = $opt{CNTRD}->dim(0);
}
else {
$opt{NSEED} = pdl($self->dim(0), $opt{NSEED})->min->sclr;
}
$opt{V} and print "$_\t=> $opt{$_}\n" for sort keys %opt;
my $ss_ms = $self->badflag? 'ms' : 'ss';
my $ss_total
= $self->badflag? $self->var->average : $self->ss->sumover;
$opt{V} and print "overall $ss_ms:\t$ss_total\n";
my ($centroid, $ss_cv, $R2, $clus_this, $clus_last);
# NTRY made into extra dim in $cluster for broadcasting
my @nclus = (ref $opt{NCLUS} eq 'ARRAY')? @{$opt{NCLUS}} : ($opt{NCLUS});
$clus_this
= (defined($opt{CNTRD}) and $opt{CNTRD}->nelem) ?
$self->assign( $opt{CNTRD}->dummy(-1) ) # put dummy(-1) to match NTRY
: random_cluster($opt{NSEED}, @nclus, $opt{NTRY} )
;
($centroid, $ss_cv) = $self->slice([0,$opt{NSEED} - 1])->centroid( $clus_this );
# now obs in $clus_this matches $self
$clus_this = $self->assign( $centroid );
($centroid, $ss_cv) = $self->centroid( $clus_this );
my $iter = 0;
do {
$R2 = $self->badflag? 1 - $ss_cv->average->average / $ss_total
: 1 - $ss_cv->sumover->sumover / $ss_total
;
$opt{V} and print join(' ',('iter', $iter++, 'R2', $R2)) . "\n";
$clus_last = $clus_this;
$clus_this = $self->assign( $centroid );
($centroid, $ss_cv) = $self->centroid( $clus_this );
}
while ( any long(abs($clus_this - $clus_last))->sumover->sumover > 0 );
$opt{FULL} and
return (
centroid => PDL::squeeze( $centroid ),
cluster => PDL::squeeze( $clus_this ),
n => PDL::squeeze( $clus_this )->sumover,
R2 => PDL::squeeze( $R2 ),
$ss_ms => PDL::squeeze( $ss_cv ),
);
# xchg/mv(-1,0) leaves it as was if single dim--unlike transpose
my $i_best = $R2->mv(-1,0)->maximum_ind;
$R2->getndims == 1 and
return (
centroid => $centroid->dice_axis(-1,$i_best)->sever->squeeze,
cluster => $clus_this->dice_axis(-1,$i_best)->sever->squeeze,
n => $clus_this->dice_axis(-1,$i_best)->sever->squeeze->sumover,
R2 => $R2->dice_axis(-1,$i_best)->sever->squeeze,
$ss_ms => $ss_cv->dice_axis(-1,$i_best)->sever->squeeze,
);
# now for broadcasting beyond 2D data
# can't believe i'm using a perl loop :P
$i_best = $i_best->flat->sever;
my @i_best = map { $opt{NTRY} * $_ + $i_best->slice("($_)") }
0 .. $i_best->nelem - 1;
my @shapes;
for ($centroid, $clus_this, $R2) {
my @dims = $_->dims;
pop @dims;
push @shapes, \@dims;
}
$clus_this = $clus_this->mv(-1,2)->clump(2..$clus_this->ndims-1)->dice_axis(2,\@i_best)->sever->reshape( @{ $shapes[1] } )->sever,
return (
centroid =>
$centroid->mv(-1,2)->clump(2..$centroid->ndims-1)->dice_axis(2,\@i_best)->sever->reshape( @{ $shapes[0] } )->sever,
cluster => $clus_this,
n => $clus_this->sumover,
R2 =>
$R2->mv(-1,0)->clump(0..$R2->ndims-1)->dice_axis(0,\@i_best)->sever->reshape( @{ $shapes[2] } )->sever,
$ss_ms =>
$ss_cv->mv(-1,2)->clump(2..$ss_cv->ndims-1)->dice_axis(2,\@i_best)->sever->reshape( @{ $shapes[0] } )->sever,
);
}
=head1 METHODS
=head2 iv_cluster
=for ref
Turns an independent variable into a cluster pdl. Returns cluster pdl and level-to-pdl_index mapping in list context and cluster pdl only in scalar context.
This is the method used for mean and var in anova. The difference between iv_cluster and dummy_code is that iv_cluster returns pdl dim [obs x level] whereas dummy_code returns pdl dim [obs x (level - 1)].
=for usage
Usage:
pdl> @bake = qw( y y y n n n )
# accepts @ ref or 1d pdl
pdl> p $bake = iv_cluster( \@bake )
[
[1 1 1 0 0 0]
[0 0 0 1 1 1]
]
pdl> p $rating = sequence 6
[0 1 2 3 4 5]
pdl> p $rating->centroid( $bake )
# mean for each iv level
[
[1 4]
]
# ss
[
[2 2]
]
=cut
*iv_cluster = \&PDL::iv_cluster;
sub PDL::iv_cluster {
my ($var_ref) = @_;
my ($var, $map_ref) = PDL::Stats::Basic::code_ivs( $var_ref );
my $var_a = yvals( short, $var->nelem, $var->max->sclr + 1 ) == $var;
$var_a = $var_a->setbadif( $var->isbad ) if $var->badflag;
wantarray ? ($var_a, $map_ref) : $var_a;
}
=head2 pca_cluster
Assign variables to components ie clusters based on pca loadings or
scores. One way to seed kmeans (see Ding & He, 2004, and Su & Dy, 2004
for other ways of using pca with kmeans). Variables are assigned to
their most associated component. Note that some components may not have
any variable that is most associated with them, so the returned number
of clusters may be smaller than NCOMP.
Default options (case insensitive):
V => 1,
ABS => 1, # high pos and neg loadings on a comp in same cluster
NCOMP => undef, # max number of components to consider. determined by
# scree plot black magic if not specified
PLOT => 0, # pca scree plot with cutoff at NCOMP
WIN => undef, # pass pgswin object for more plotting control
Usage:
# say we need to cluster a group of documents
# $data is pdl dim [word x doc]
($data, $idd, $idw) = get_data 'doc_word_info.txt';
pdl> %p = $data->pca;
# $cluster is pdl mask dim [doc x ncomp]
pdl> $cluster = $p{loading}->pca_cluster;
# pca clusters var while kmeans clusters obs. hence transpose
pdl> ($m, $ss) = $data->transpose->centroid( $cluster );
pdl> %k = $data->transpose->kmeans( { cntrd=>$m } );
# take a look at cluster 0 doc ids
pdl> p join("\n", @$idd[ list which $k{cluster}->( ,0) ]);
=cut
*pca_cluster = \&PDL::pca_cluster;
sub PDL::pca_cluster {
my ($self, $opt) = @_;
my %opt = (
V => 1,
ABS => 1, # high pos and neg loadings on a comp in same cluster
NCOMP => undef, # max number of components to consider. determined by
# scree plot black magic if not specified
PLOT => 0, # pca scree plot with cutoff at NCOMP
WIN => undef, # pass pgswin object for more plotting control
);
if ($opt) { $opt{uc $_} = $opt->{$_} for keys %$opt; }
my $var = sumover($self ** 2) / $self->dim(0);
if (!$opt{NCOMP}) {
# here's the black magic part
my $comps = ($self->dim(1) > 300)? int($self->dim(1) * .1)
: pdl($self->dim(1), 30)->min
;
$var = $var->slice([0,$comps-1])->sever;
$opt{NCOMP} = _scree_ind( $var );
}
$opt{PLOT} and do {
require PDL::Stats::GLM;
$var->plot_screes({NCOMP=>$var->dim(0), CUT=>$opt{NCOMP}, WIN=>$opt{WIN}});
};
my $c = $self->slice(':',[0,$opt{NCOMP}-1])->transpose->abs->maximum_ind;
if ($opt{ABS}) {
$c = $c->iv_cluster;
}
else {
my @c = map { ($self->slice($_,$c->slice($_)) >= 0)? $c->slice($_)*2 : $c->slice($_)*2 + 1 }
( 0 .. $c->dim(0)-1 );
$c = iv_cluster( \@c );
}
$opt{V} and print "cluster membership mask as " . $c->info . "\n";
return $c;
}
=head1 REFERENCES
Ding, C., & He, X. (2004). K-means clustering via principal component analysis. Proceedings of the 21st International Conference on Machine Learning, 69, 29.
Su, T., & Dy, J. (2004). A deterministic method for initializing K-means clustering. 16th IEEE International Conference on Tools with Artificial Intelligence, 784-786.
Romesburg, H.C. (1984). Cluster Analysis for Researchers. NC: Lulu Press.
Wikipedia (retrieved June, 2009). K-means clustering. http://en.wikipedia.org/wiki/K-means_algorithm
=head1 AUTHOR
Copyright (C) 2009 Maggie J. Xiong <maggiexyz users.sourceforge.net>
All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation as described in the file COPYING in the PDL distribution.
=cut
EOD
pp_done();