Algorithm-KMeans
view release on metacpan or search on metacpan
lib/Algorithm/KMeans.pm view on Meta::CPAN
or die "Unable to open a temp file in this directory: $!\n";
foreach my $cluster (@{$self->{_clusters}}) {
foreach my $item (@$cluster) {
print OUTPUT "@{$visualization_data{$item}}";
print OUTPUT "\n";
}
print OUTPUT "\n\n";
}
close OUTPUT;
my $plot;
my $hardcopy_plot;
if (!defined $pause_time) {
$plot = Graphics::GnuplotIF->new( persist => 1 );
$hardcopy_plot = Graphics::GnuplotIF->new();
$hardcopy_plot->gnuplot_cmd('set terminal png', "set output \"clustering_results.png\"");
} else {
$plot = Graphics::GnuplotIF->new();
}
$plot->gnuplot_cmd( "set noclip" );
$plot->gnuplot_cmd( "set pointsize 2" );
my $arg_string = "";
if ($visualization_data_field_width > 2) {
foreach my $i (0..$K-1) {
my $j = $i + 1;
$arg_string .= "\"$temp_file\" index $i using 1:2:3 title \"Cluster $i\" with points lt $j pt $j, ";
}
} elsif ($visualization_data_field_width == 2) {
foreach my $i (0..$K-1) {
my $j = $i + 1;
$arg_string .= "\"$temp_file\" index $i using 1:2 title \"Cluster $i\" with points lt $j pt $j, ";
}
} elsif ($visualization_data_field_width == 1 ) {
foreach my $i (0..$K-1) {
my $j = $i + 1;
$arg_string .= "\"$temp_file\" index $i using 1 title \"Cluster $i\" with points lt $j pt $j, ";
}
}
$arg_string = $arg_string =~ /^(.*),[ ]+$/;
$arg_string = $1;
if ($visualization_data_field_width > 2) {
$plot->gnuplot_cmd( "splot $arg_string" );
$hardcopy_plot->gnuplot_cmd( "splot $arg_string" ) unless defined $pause_time;
$plot->gnuplot_pause( $pause_time ) if defined $pause_time;
} elsif ($visualization_data_field_width == 2) {
$plot->gnuplot_cmd( "plot $arg_string" );
$hardcopy_plot->gnuplot_cmd( "plot $arg_string" ) unless defined $pause_time;
$plot->gnuplot_pause( $pause_time ) if defined $pause_time;
} elsif ($visualization_data_field_width == 1) {
croak "No provision for plotting 1-D data\n";
}
}
# It makes sense to call visualize_data() only AFTER you have called the method
# read_data_from_file().
#
# The visualize_data() is meant for the visualization of the original data in its
# various 2D or 3D subspaces. The method can also be used to visualize the normed
# data in a similar manner. Recall the normed data is the original data after each
# data dimension is normalized by the standard-deviation along that dimension.
#
# Whether you see the original data or the normed data depends on the second
# argument supplied in the method call. It must be either the string 'original' or
# the string 'normed'.
sub visualize_data {
my $self = shift;
my $v_mask = shift || croak "visualization mask missing";
my $datatype = shift; # must be either 'original' or 'normed'
croak "\n\nABORTED: You called visualize_data() for normed data " .
"but without first turning on data normalization in the " .
"in the KMeans constructor"
if ($datatype eq 'normed') && ! $self->{_var_normalize};
my $master_datafile = $self->{_datafile};
my @v_mask = split //, $v_mask;
my $visualization_mask_width = @v_mask;
my $original_data_mask = $self->{_mask};
my @mask = split //, $original_data_mask;
my $data_field_width = scalar grep {$_ eq '1'} @mask;
croak "\n\nABORTED: The width of the visualization mask (including " .
"all its 1s and 0s) must equal the width of the original mask " .
"used for reading the data file (counting only the 1's)"
if $visualization_mask_width != $data_field_width;
my $visualization_data_field_width = scalar grep {$_ eq '1'} @v_mask;
my %visualization_data;
my $data_source;
if ($datatype eq 'original') {
$data_source = $self->{_original_data};
} elsif ($datatype eq 'normed') {
$data_source = $self->{_data};
} else {
croak "\n\nABORTED: improper call to visualize_data()";
}
while ( my ($record_id, $data) = each %{$data_source} ) {
my @fields = @$data;
croak "\nABORTED: Visualization mask size exceeds data record size\n"
if $#v_mask > $#fields;
my @data_fields;
foreach my $i (0..@fields-1) {
if ($v_mask[$i] eq '0') {
next;
} elsif ($v_mask[$i] eq '1') {
push @data_fields, $fields[$i];
} else {
croak "Misformed visualization mask. It can only have 1s and 0s\n";
}
}
$visualization_data{ $record_id } = \@data_fields;
}
my $filename = basename($master_datafile);
my $temp_file;
if ($datatype eq 'original') {
$temp_file = "__temp_data_" . $filename;
} elsif ($datatype eq 'normed') {
$temp_file = "__temp_normed_data_" . $filename;
} else {
croak "ABORTED: Improper call to visualize_data()";
}
unlink $temp_file if -e $temp_file;
open OUTPUT, ">$temp_file"
or die "Unable to open a temp file in this directory: $!\n";
foreach my $datapoint (values %visualization_data) {
print OUTPUT "@$datapoint";
lib/Algorithm/KMeans.pm view on Meta::CPAN
die "matrix multiplication called with non-matching matrix arguments"
unless $nrows1 == $ncols2 && $ncols1 == $nrows2;
if ($nrows1 == 1) {
my @row = $matrix1->row(0)->as_list;
my @col = $matrix2->col(0)->as_list;
my $result;
foreach my $j (0..$ncols1-1) {
$result += $row[$j] * $col[$j];
}
return $result;
}
my $product = Math::GSL::Matrix->new($nrows1, $nrows1);
foreach my $i (0..$nrows1-1) {
my $row = $matrix1->row($i);
my @product_row;
foreach my $j (0..$ncols2-1) {
my $col = $matrix2->col($j);
my $row_times_col = matrix_multiply($row, $col);
push @product_row, $row_times_col;
}
$product->set_row($i, \@product_row);
}
return $product;
}
1;
=pod
=head1 NAME
Algorithm::KMeans - for clustering multidimensional data
=head1 SYNOPSIS
# You now have four different choices for clustering your data with this module:
#
# 1) With Euclidean distances and with random cluster seeding
#
# 2) With Mahalanobis distances and with random cluster seeding
#
# 3) With Euclidean distances and with smart cluster seeding
#
# 4) With Mahalanobis distances and with smart cluster seeding
#
# Despite the qualifier 'smart' in 'smart cluster seeding', it may not always
# produce results that are superior to those obtained with random seeding. (If you
# also factor in the choice regarding variance normalization, you actually have
# eight different choices for data clustering with this module.)
#
# In all cases, you'd obviously begin with
use Algorithm::KMeans;
# You'd then name the data file:
my $datafile = "mydatafile.csv";
# Next, set the mask to indicate which columns of the datafile to use for
# clustering and which column contains a symbolic ID for each data record. For
# example, if the symbolic name is in the first column, you want the second column
# to be ignored, and you want the next three columns to be used for 3D clustering,
# you'd set the mask to:
my $mask = "N0111";
# Now construct an instance of the clusterer. The parameter K controls the number
# of clusters. If you know how many clusters you want (let's say 3), call
my $clusterer = Algorithm::KMeans->new( datafile => $datafile,
mask => $mask,
K => 3,
cluster_seeding => 'random',
terminal_output => 1,
write_clusters_to_files => 1,
);
# By default, this constructor call will set you up for clustering based on
# Euclidean distances. If you want the module to use Mahalanobis distances, your
# constructor call will look like:
my $clusterer = Algorithm::KMeans->new( datafile => $datafile,
mask => $mask,
K => 3,
cluster_seeding => 'random',
use_mahalanobis_metric => 1,
terminal_output => 1,
write_clusters_to_files => 1,
);
# For both constructor calls shown above, you can use smart seeding of the clusters
# by changing 'random' to 'smart' for the cluster_seeding option. See the
# explanation of smart seeding in the Methods section of this documentation.
# If your data is such that its variability along the different dimensions of the
# data space is significantly different, you may get better clustering if you first
# normalize your data by setting the constructor parameter
# do_variance_normalization as shown below:
my $clusterer = Algorithm::KMeans->new( datafile => $datafile,
mask => $mask,
K => 3,
cluster_seeding => 'smart', # or 'random'
terminal_output => 1,
do_variance_normalization => 1,
write_clusters_to_files => 1,
);
# But bear in mind that such data normalization may actually decrease the
# performance of the clusterer if the variability in the data is more a result of
# the separation between the means than a consequence of intra-cluster variance.
# Set K to 0 if you want the module to figure out the optimum number of clusters
# from the data. (It is best to run this option with the terminal_output set to 1
# so that you can see the different value of QoC for the different K):
my $clusterer = Algorithm::KMeans->new( datafile => $datafile,
mask => $mask,
K => 0,
cluster_seeding => 'random', # or 'smart'
terminal_output => 1,
write_clusters_to_files => 1,
);
# Although not shown above, you can obviously set the 'do_variance_normalization'
# flag here also if you wish.
# For very large data files, setting K to 0 will result in searching through too
# many values for K. For such cases, you can range limit the values of K to search
# through by
my $clusterer = Algorithm::KMeans->new( datafile => $datafile,
mask => "N111",
Kmin => 3,
Kmax => 10,
cluster_seeding => 'random', # or 'smart'
terminal_output => 1,
write_clusters_to_files => 1,
);
# FOR ALL CASES ABOVE, YOU'D NEED TO MAKE THE FOLLOWING CALLS ON THE CLUSTERER
# INSTANCE TO ACTUALLY CLUSTER THE DATA:
$clusterer->read_data_from_file();
$clusterer->kmeans();
# If you want to directly access the clusters and the cluster centers in your own
# top-level script, replace the above two statements with:
$clusterer->read_data_from_file();
my ($clusters_hash, $cluster_centers_hash) = $clusterer->kmeans();
# You can subsequently access the clusters directly in your own code, as in:
lib/Algorithm/KMeans.pm view on Meta::CPAN
Version 1.30 includes a bug fix for the case when the datafile contains empty lines,
that is, lines with no data records. Another bug fix in Version 1.30 deals with the
case when you want the module to figure out how many clusters to form (this is the
C<K=0> option in the constructor call) and the number of data records is close to the
minimum.
Version 1.21 includes fixes to handle the possibility that, when clustering the data
for a fixed number of clusters, a cluster may become empty during iterative
calculation of cluster assignments of the data elements and the updating of the
cluster centers. The code changes are in the C<assign_data_to_clusters()> and
C<update_cluster_centers()> subroutines.
Version 1.20 includes an option to normalize the data with respect to its variability
along the different coordinates before clustering is carried out.
Version 1.1.1 allows for range limiting the values of C<K> to search through. C<K>
stands for the number of clusters to form. This version also declares the module
dependencies in the C<Makefile.PL> file.
Version 1.1 is a an object-oriented version of the implementation presented in
version 1.0. The current version should lend itself more easily to code extension.
You could, for example, create your own class by subclassing from the class presented
here and, in your subclass, use your own criteria for the similarity distance between
the data points and for the QoC (Quality of Clustering) metric, and, possibly a
different rule to stop the iterations. Version 1.1 also allows you to directly
access the clusters formed and the cluster centers in your calling script.
=head1 SPECIAL USAGE NOTE
If you were directly accessing in your own scripts the clusters produced by the older
versions of this module, you'd need to make changes to your code if you wish to use
Version 2.0 or higher. Instead of returning arrays of clusters and cluster centers,
Versions 2.0 and higher return hashes. This change was made necessary by the logic
required for implementing the two new C<which_cluster> methods that were introduced
in Version 2.0. These methods return the best cluster for a new data sample from the
clusters you created using the existing data. One of the C<which_cluster> methods is
based on the Euclidean metric for finding the cluster that is closest to the new data
sample, and the other on the Mahalanobis metric. Another point of incompatibility
with the previous versions is that you must now explicitly set the C<cluster_seeding>
parameter in the call to the constructor to either C<random> or C<smart>. This
parameter does not have a default associated with it starting with Version 2.0.
=head1 DESCRIPTION
Clustering with K-Means takes place iteratively and involves two steps: 1) assignment
of data samples to clusters on the basis of how far the data samples are from the
cluster centers; and 2) Recalculation of the cluster centers (and cluster covariances
if you are using the Mahalanobis distance metric for clustering).
Obviously, before the two-step approach can proceed, we need to initialize the the
cluster centers. How this initialization is carried out is important. The module
gives you two very different ways for carrying out this initialization. One option,
called the C<smart> option, consists of subjecting the data to principal components
analysis to discover the direction of maximum variance in the data space. The data
points are then projected on to this direction and a histogram constructed from the
projections. Centers of the smoothed histogram are used to seed the clustering
operation. The other option is to choose the cluster centers purely randomly. You
get the first option if you set C<cluster_seeding> to C<smart> in the constructor,
and you get the second option if you set it to C<random>.
How to specify the number of clusters, C<K>, is one of the most vexing issues in any
approach to clustering. In some case, we can set C<K> on the basis of prior
knowledge. But, more often than not, no such prior knowledge is available. When the
programmer does not explicitly specify a value for C<K>, the approach taken in the
current implementation is to try all possible values between 2 and some largest
possible value that makes statistical sense. We then choose that value for C<K>
which yields the best value for the QoC (Quality of Clustering) metric. It is
generally believed that the largest value for C<K> should not exceed C<sqrt(N/2)>
where C<N> is the number of data samples to be clustered.
What to use for the QoC metric is obviously a critical issue unto itself. In the
current implementation, the value of QoC is the ratio of the average radius of the
clusters and the average distance between the cluster centers.
Every iterative algorithm requires a stopping criterion. The criterion implemented
here is that we stop iterations when there is no re-assignment of the data points
during the assignment step.
Ordinarily, the output produced by a K-Means clusterer will correspond to a local
minimum for the QoC values, as opposed to a global minimum. The current
implementation protects against that when the module constructor is called with the
C<random> option for C<cluster_seeding> by trying different randomly selected initial
cluster centers and then selecting the one that gives the best overall QoC value.
A K-Means clusterer will generally produce good results if the overlap between the
clusters is minimal and if each cluster exhibits variability that is uniform in all
directions. When the data variability is different along the different directions in
the data space, the results you obtain with a K-Means clusterer may be improved by
first normalizing the data appropriately, as can be done in this module when you set
the C<do_variance_normalization> option in the constructor. However, as pointed out
elsewhere in this documentation, such normalization may actually decrease the
performance of the clusterer if the overall data variability along any dimension is
more a result of separation between the means than a consequence of intra-cluster
variability.
=head1 METHODS
The module provides the following methods for clustering, for cluster visualization,
for data visualization, for the generation of data for testing a clustering
algorithm, and for determining the cluster membership of a new data sample:
=over 4
=item B<new():>
my $clusterer = Algorithm::KMeans->new(datafile => $datafile,
mask => $mask,
K => $K,
cluster_seeding => 'random', # also try 'smart'
use_mahalanobis_metric => 1, # also try '0'
terminal_output => 1,
write_clusters_to_files => 1,
);
A call to C<new()> constructs a new instance of the C<Algorithm::KMeans> class. When
C<$K> is a non-zero positive integer, the module will construct exactly that many
clusters. However, when C<$K> is 0, the module will find the best number of clusters
to partition the data into. As explained in the Description, setting
C<cluster_seeding> to C<smart> causes PCA (principal components analysis) to be used
for discovering the best choices for the initial cluster centers. If you want purely
random decisions to be made for the initial choices for the cluster centers, set
C<cluster_seeding> to C<random>.
The data file is expected to contain entries in the following format
c20 0 10.7087017086940 9.63528386251712 10.9512155258108 ...
c7 0 12.8025925026787 10.6126270065785 10.5228482095349 ...
b9 0 7.60118206283120 5.05889245193079 5.82841781759102 ...
....
....
where the first column contains the symbolic ID tag for each data record and the rest
of the columns the numerical information. As to which columns are actually used for
clustering is decided by the string value of the mask. For example, if we wanted to
cluster on the basis of the entries in just the 3rd, the 4th, and the 5th columns
above, the mask value would be C<N0111> where the character C<N> indicates that the
ID tag is in the first column, the character C<0> that the second column is to be
ignored, and the C<1>'s that follow that the 3rd, the 4th, and the 5th columns are to
be used for clustering.
If you wish for the clusterer to search through a C<(Kmin,Kmax)> range of values for
C<K>, the constructor should be called in the following fashion:
my $clusterer = Algorithm::KMeans->new(datafile => $datafile,
mask => $mask,
Kmin => 3,
Kmax => 10,
cluster_seeding => 'smart', # try 'random' also
terminal_output => 1,
);
where obviously you can choose any reasonable values for C<Kmin> and C<Kmax>. If you
choose a value for C<Kmax> that is statistically too large, the module will let you
know. Again, you may choose C<random> for C<cluster_seeding>, the default value being
C<smart>.
If you believe that the variability of the data is very different along the different
dimensions of the data space, you may get better clustering by first normalizing the
data coordinates by the standard-deviations along those directions. When you set the
constructor option C<do_variance_normalization> as shown below, the module uses the
overall data standard-deviation along a direction for the normalization in that
direction. (As mentioned elsewhere in the documentation, such a normalization could
backfire on you if the data variability along a dimension is more a result of the
separation between the means than a consequence of the intra-cluster variability.):
my $clusterer = Algorithm::KMeans->new( datafile => $datafile,
mask => "N111",
K => 2,
cluster_seeding => 'smart', # try 'random' also
terminal_output => 1,
do_variance_normalization => 1,
);
=back
=head2 Constructor Parameters
=over 8
=item C<datafile>:
This parameter names the data file that contains the multidimensional data records
you want the module to cluster.
=item C<mask>:
This parameter supplies the mask to be applied to the columns of your data file. See
the explanation in Synopsis for what this mask looks like.
=item C<K>:
This parameter supplies the number of clusters you are looking for. If you set this
option to 0, that means that you want the module to search for the best value for
C<K>. (Keep in mind the fact that searching for the best C<K> may take a long time
for large data files.)
=item C<Kmin>:
lib/Algorithm/KMeans.pm view on Meta::CPAN
produces bizarre results, try C<random>.
=item C<use_mahalanobis_metric>:
When set to 1, this option causes Mahalanobis distances to be used for clustering.
The default is 0 for this parameter. By default, the module uses the Euclidean
distances for clustering. In general, Mahalanobis distance based clustering will
fail if your data resides on a lower-dimensional hyperplane in the data space, if you
seek too many clusters, and if you do not have a sufficient number of samples in your
data file. A necessary requirement for the module to be able to compute Mahalanobis
distances is that the cluster covariance matrices be non-singular. (Let's say your
data dimensionality is C<D> and the module is considering a cluster that has only
C<d> samples in it where C<d> is less than C<D>. In this case, the covariance matrix
will be singular since its rank will not exceed C<d>. For the covariance matrix to
be non-singular, it must be of full rank, that is, its rank must be C<D>.)
=item C<do_variance_normalization>:
When set, the module will first normalize the data variance along the different
dimensions of the data space before attempting clustering. Depending on your data,
this option may or may not result in better clustering.
=item C<terminal_output>:
This boolean parameter, when not supplied in the call to C<new()>, defaults to 0.
When set, you will see in your terminal window the different clusters as lists of the
symbolic IDs and their cluster centers. You will also see the QoC (Quality of
Clustering) values for the clusters displayed.
=item C<write_clusters_to_files>:
This parameter is also boolean. When set to 1, the clusters are written out to files
that are named in the following manner:
cluster0.txt
cluster1.txt
cluster2.txt
...
...
Before the clusters are written to these files, the module destroys all files with
such names in the directory in which you call the module.
=back
=over
=item B<read_data_from_file()>
$clusterer->read_data_from_file()
=item B<kmeans()>
$clusterer->kmeans();
or
my ($clusters_hash, $cluster_centers_hash) = $clusterer->kmeans();
The first call above works solely by side-effect. The second call also returns the
clusters and the cluster centers. See the C<cluster_and_visualize.pl> script in the
C<examples> directory for how you can in your own code extract the clusters and the
cluster centers from the variables C<$clusters_hash> and C<$cluster_centers_hash>.
=item B<get_K_best()>
$clusterer->get_K_best();
This call makes sense only if you supply either the C<K=0> option to the constructor,
or if you specify values for the C<Kmin> and C<Kmax> options. The C<K=0> and the
C<(Kmin,Kmax)> options cause the module to determine the best value for C<K>.
Remember, C<K> is the number of clusters the data is partitioned into.
=item B<show_QoC_values()>
$clusterer->show_QoC_values();
presents a table with C<K> values in the left column and the corresponding QoC
(Quality-of-Clustering) values in the right column. Note that this call makes sense
only if you either supply the C<K=0> option to the constructor, or if you specify
values for the C<Kmin> and C<Kmax> options.
=item B<visualize_clusters()>
$clusterer->visualize_clusters( $visualization_mask )
The visualization mask here does not have to be identical to the one used for
clustering, but must be a subset of that mask. This is convenient for visualizing
the clusters in two- or three-dimensional subspaces of the original space.
=item B<visualize_data()>
$clusterer->visualize_data($visualization_mask, 'original');
$clusterer->visualize_data($visualization_mask, 'normed');
This method requires a second argument and, as shown, it must be either the string
C<original> or the string C<normed>, the former for the visualization of the raw data
and the latter for the visualization of the data after its different dimensions are
normalized by the standard-deviations along those directions. If you call the method
with the second argument set to C<normed>, but do so without turning on the
C<do_variance_normalization> option in the KMeans constructor, it will let you know.
=item B<which_cluster_for_new_data_element()>
If you wish to determine the cluster membership of a new data sample after you have
created the clusters with the existing data samples, you would need to call this
method. The C<which_cluster_for_new_data.pl> script in the C<examples> directory
shows how to use this method.
=item B<which_cluster_for_new_data_element_mahalanobis()>
This does the same thing as the previous method, except that it determines the
cluster membership using the Mahalanobis distance metric. As for the previous
method, see the C<which_cluster_for_new_data.pl> script in the C<examples> directory
for how to use this method.
=item B<cluster_data_generator()>
Algorithm::KMeans->cluster_data_generator(
input_parameter_file => $parameter_file,
output_datafile => $out_datafile,
number_data_points_per_cluster => 20 );
for generating multivariate data for clustering if you wish to play with synthetic
data for clustering. The input parameter file contains the means and the variances
for the different Gaussians you wish to use for the synthetic data. See the file
C<param.txt> provided in the examples directory. It will be easiest for you to just
edit this file for your data generation needs. In addition to the format of the
parameter file, the main constraint you need to observe in specifying the parameters
is that the dimensionality of the covariance matrix must correspond to the
dimensionality of the mean vectors. The multivariate random numbers are generated by
calling the C<Math::Random> module. As you would expect, this module requires that
the covariance matrices you specify in your parameter file be symmetric and positive
definite. Should the covariances in your parameter file not obey this condition, the
C<Math::Random> module will let you know.
=back
=head1 HOW THE CLUSTERS ARE OUTPUT
When the option C<terminal_output> is set in the call to the constructor, the
clusters are displayed on the terminal screen.
When the option C<write_clusters_to_files> is set in the call to the constructor, the
module dumps the clusters in files named
cluster0.txt
cluster1.txt
cluster2.txt
...
...
in the directory in which you execute the module. The number of such files will
equal the number of clusters formed. All such existing files in the directory are
destroyed before any fresh ones are created. Each cluster file contains the symbolic
ID tags of the data samples in that cluster.
( run in 0.668 second using v1.01-cache-2.11-cpan-39bf76dae61 )