#!/usr/bin/perl -w
#
#$Id: runkmeans,v 1.36 2004/01/21 01:48:59 dpelleg Exp $
#
# runkmeans --- run the kmeans implementation
# Dan Pelleg, November 1998
# Copyright Dan Pelleg

use IPC::Open3;
use Symbol;
use Getopt::Std;

use vars qw($cache_hit $cache_miss $cache_write);

# Euclidean distance, squared
# subd_sqd(@$v1, @$v2)

#given point and list of candidates, returns index of nearest
# candidate to point
# sub find_nearest_neighbor(@$point, @$candidates);

#debug directives
# some debug directives
local($debug_level, 
	  $D_CENTERS, $D_EXEC, $D_SHOW_SUBPROCESS, $D_EXTRA) 
	= (0, 0x01, 0x02, 0x04, 0x08);

my %opts;

getopt('bufFmdksnlcirpBvLCrRSt', \%opts);

local @orig_dist;				# each line is original center, stored as
                                # ref to list
local @centers;					# each line is claimed center, stored as
                                # ref to list

local @closest_orig;            # for center i stores the neareset orig neighbor
                                # (we want this to be a permutation to mean anything)

local $bad_grouping = 0;		# if couldn't get good centers

local $tempdir = $ENV{"TEMP"} || "/tmp";
local $bindir = ".";
local $tmpclust = "$tempdir/clust.$$";
local $ANALYSE_CACHE = 0;		# BUGBUG: setting to 1 increases runtime significantly

$debug_level = $opts{d} || 0;
$num_classes = $opts{k} || 5;
$sigma = $opts{s} || 0.02;
$num_points = $opts{n} || 500;
$num_dims = $opts{m} || 2;
$min_bwidth = $opts{b} || "0.03";
$method = $opts{f} || "blacklist";
$leaf_size = $opts{l} || 40;
$cutoff_factor = $opts{c} || "0.5";
$max_iter = defined($opts{i}) ? $opts{i} : 200;
$num_splits = defined($opts{p}) ? $opts{p} : 1;
$kmeans_binary_default = "kmeans";
$kmeans_binary = $opts{B} || $kmeans_binary_default ;
$save_file = $opts{v};
$load_file = $opts{L};
$clusterfile = $opts{C};
$del_steps_ratio = $opts{R};
$split_stat = $opts{S} || "BIC";
$max_ctrs_overshoot = $opts{t} || 1;
$forced_split_fraction = $opts{F};

if($opts{u}) {
	$create_datafile = ($opts{u} eq "create");
	$destroy_datafile = 0;
	$run_kmeans = ($opts{u} ne "create");
	$datafile = $create_datafile ? "$tempdir/kmeans.$$.dat" : $opts{u};
} else {
	$create_datafile = 1;
	$destroy_datafile = 1;
	$run_kmeans = 1;
	$datafile = "$tempdir/kmeans.$$.dat";
}

$unifile = "$datafile.universe"; # contains universe dimensions

if($opts{h}) {
	print <<EOF;
Usage: $0 [-d debug_level] [-k classes] [-s sigma] [-n num_points]
	[-m num_dims] [-f method] [-u datafile] [-b min_box_width]
	[-l max_leaf_size] [-c cutoff_factor] [-i max-iter] [-p num_splits]
	[-B binary] [-v savefile] [-L loadfile] [-C clusterfile]
    [-R del_steps_ratio]
  Generate a random data file, run kmeans on it and measure how
  close the reported centers are to the original distribution.
  The "-f" switch stands for method: slow, hplane, BL-RS, BL-nocache,
  BL-search, or blacklist.
  The "-u" switch specifies which data file to use.
  If it is omitted, a random data file is created and deleted after
  usage. If it is the word "create", a random data file is created,
  its name is printed, and no other action is taken. If it is any
  other value, it is taken to be the name of a data file to use.
  If num_splits is negative, runs kmeans with just two starting 
  centers and let it produce at most <number of classes> centers. 
  The "-B" specifies a non-standard binary to execute instead of "kmeans".
  If max-iter is zero, will try and read the initial centers and report
  the scores on them.
  With "-v", will save the resulting centers to savefile.
  With "-L", will read the initial centers from loadfile.
  With "-C", will store point membership into clusterfile.
  The "-r" switch is just passed to kmeans and is the inverse of the ratio
  of deletion steps to split steps. So if this is <r>, then one of every <r>
  splitting steps will not be performed, and a delete-worst-center step will
  take its place (I know, the exact ratio is not 1/r, but its close enough).
EOF
	exit(1);
}

# pipes for child process write, read & error
$WRITER = gensym();
$READER = gensym();
$ERROR = gensym();

# first run kmeans to generate a data file
# keep the original distribution in orig_dist
if($create_datafile) {
	$prog = "$bindir/$kmeans_binary gen -num_classes $num_classes -out $datafile"
		." -num_rows $num_points -draw false -sigma $sigma -num_cols $num_dims"
			." -D_SHOW_START_CENTERS"
				." -seed " . time();

	dbg_print($D_EXEC, "Execing $prog\n");

	open3($WRITER,$READER,$ERROR, $prog)
		or die("cant do process!");
	while(<$ERROR>) {
		dbg_print($D_SHOW_SUBPROCESS, $_);
		if(/^.*\[([^\]]*)]\s*$/) {
		    my $str = $1;
			$str =~ s/^\s+//;
			push(@orig_dist, [ split(/\s+/,$str) ]);
		}
    }
    close($WRITER);
    close($READER);
    close($ERROR);

    die "Got wrong number of rows" if($num_classes != @orig_dist);

    dbg_print($D_CENTERS, "Original distribution:\n");
    for($i=0; $i<@orig_dist; $i++) {
		for($j=0; $j<=$#{$orig_dist[$i]}; $j++) {
			dbg_print($D_CENTERS, "$orig_dist[$i][$j] ");
		}
		dbg_print($D_CENTERS, "\n");
	}
}

# read the datafile to get the parameters
open(DATA, "$datafile") or die("Can't open datafile $datafile");
while(<DATA>) {
	$num_points = $1 if /^\# num_rows = (\d+)\s*$/;
    $num_dims = $1 if /^\# num_cols = (\d+)\s*$/;
    ($sigma = $1, last) if /^\# sigma = ([0-9.]+)\s*$/;
}
close(DATA);

if($run_kmeans) {
	# now run kmeans to find clusters
	# keep the result in centers
	
	#accumulators for program counters
	local ($sum_cost, $num_iter, $sum_fail, $etime, $cache_hit,
		   $cache_miss, $cache_write, $sum_leaves) = (0,0,0,0,0,0,0,0);

    local (%hits, %unhits);
    # see if a universal node has to be created
    if(! -s $unifile) {
		my $WRITER = gensym();
		my $READER = gensym();
		my $ERROR = gensym();
		open3($WRITER,$READER,$ERROR, "$bindir/$kmeans_binary makeuni -in $datafile") or die("cant do process!");
		while(<$READER>) {
		}
	}
    if($max_iter == 0) {
		$num_splits = 1;
	}
    $prog = "$bindir/$kmeans_binary kmeans -draw false -D_SHOW_END_CENTERS"
	." -D_SHOWCOUNTERS -in $datafile -method $method"
	." -max_leaf_size $leaf_size -min_box_width $min_bwidth"
	." -cutoff_factor $cutoff_factor -max_iter $max_iter";
    $prog .= " -D_KDCACHE" if($ANALYSE_CACHE);
    $prog .= " -seed $opts{r}" if($opts{r});
    $prog .= " -printclusters /dev/stderr" if($clusterfile);
    $prog .= " -del_steps_ratio $del_steps_ratio" if($del_steps_ratio);
    $prog .= " -split_stat $split_stat " if($split_stat);
    $prog .= " -forced_split_fraction $forced_split_fraction " if(defined($forced_split_fraction));
    # are we cheating by giving the program just 2 ctrs to start with?
    my $max_ctrs_allowed = $num_classes * $max_ctrs_overshoot;
    if(0 && $num_splits > 0) {	# was: just $num_splits > 0
	  $prog .= " -max_ctrs $max_ctrs_allowed -k 2 -num_splits $num_splits";
	} else {
	  $prog .= " -num_splits $num_splits -max_ctrs $max_ctrs_allowed";
	}
    if($max_iter == 0) {
	  $prog .= " -init_ctrs $datafile.ctrs";
	  $method = "true";
	}
	if($save_file) {
	  $prog .= " -save_ctrs $save_file";
	}
	if($load_file) {
	  $prog .= " -init_ctrs $load_file";
	}

    dbg_print($D_EXEC, "Execing $prog\n");

    if($clusterfile) {
      open(CLUSTERS, ">$tmpclust") or die ("Can't open $tmpclust");
    }
    open3($WRITER,$READER,$ERROR, $prog) or die("cant do process!");

    while(<$ERROR>) {
		dbg_print($D_SHOW_SUBPROCESS, $_);
        if(/^\+/) {             # a line with cluster information
          print CLUSTERS;
          next;
        }
		chop;
		if(/^\#\#END CENTERS:$/) {	# reading centers state
			my $state = "BEFORE_CENTER";
			@centers = ();
			while(<$ERROR>) {
				dbg_print($D_SHOW_SUBPROCESS, $_);
				chop;
				if(/^[^\(]*\(\s*(-?[0-9.]+)\s*\)\s*$/) {	# another coordinate
					if($state eq "BEFORE_CENTER") {
						push(@centers, [ ]); # create a new center
						$state = "IN_CENTER";
					}
					push(@{$centers[$#centers]}, $1); # put in new coordinate
				}
				if(/^$/) {			# empty line = end of center
					$state = "BEFORE_CENTER";
				}
				last if (/^\#\#END CENTERS END.$/);
			}
		}

		if(/^\#\#cost: (\-?\d+)$/) {
			$sum_cost += $1;
			$sum_cost = 0 if ($1 < 0); # counter overflowed
			$num_iter++;
		}
		if(/^\#\#avg. len total: ([0-9.]+)$/) {
			$sum_fail += $1;
		}
		if(/^\#\#S?TIME ([0-9.]+)$/) {
			$etime += $1;
		}
		if(/^\#\#DISTORTION ([0-9.]+)$/) {
			$distortion = $1;
		}
		if(/^\#\#(BICSCORE) (\-?[0-9.]+)$/) {
			$bicscore = $2;
		}
		if(/^\#\#NPOINTS (\d+)$/) {
			$num_points = $1;
		}
		if(/^\#\#NCOLS (\d+)$/) {
			$num_dims = $1;
		}
		if(/^\#leaves.*: (\d+)$/) {
			$sum_leaves += $1;
		}
 		if(/^Cache miss/) {
 			$cache_miss++;
 		}
 		if(/^Cache hit/) {
 			$cache_hit++;
 		}
 		if(/^Cache write/) {
 			$cache_write++;
 		}
		if(/^Hit count (-?\d+)/) {
			$hits{$1}++;
		}
		if(/^Unhit count (-?\d+)/) {
			$unhits{$1}++;
		}
		if(/^\#(NCENTERS): (\d+)/) {
			print STDERR sprintf("%s %4d BICSCORE %6.2e\n", $1, $2, $bicscore);
		}
		if(/^\#(ITERS): (\d+)/) {
			print STDERR sprintf("%s %4d ", $1, $2);
		}
		if(/^\#(INFO)(.*)$/) {
			print STDERR sprintf("INFO %s", $2);
		}
#		if(/Wctrs/) {			# BUGBUG
#			print STDERR "$_\n";
#		}
#		if(/^\#\#([BS]TIME) ([0-9.]+)$/) {
#			print STDERR sprintf("%s %6.4f ", $1, $2);
#		}
	}

    close($WRITER);
    close($READER);
    close($ERROR);
    close(CLUSTERS);

    dbg_print($D_CENTERS, "Centers:\n");
    for($i=0; $i<@centers; $i++) {
		for($j=0; $j<=$#{$centers[$i]}; $j++) {
			dbg_print($D_CENTERS, "$centers[$i][$j] ");
		}
		dbg_print($D_CENTERS, "\n");
	}

    if($create_datafile) {
        # find closest neighbors for centers
    	for($i=0; $i<@centers; $i++) {
			$closest_orig[$i] = &find_nearest_neighbor($centers[$i], \@orig_dist);
		}

		# check to see if closest_orig is a permutation
		local @occupied;
		for($i=0; $i<@closest_orig; $i++) {
			$occupied[$closest_orig[$i]] = 1;
		}
		for($i=0; $i<@closest_orig; $i++) {
			if(!defined($occupied[$i])) {
				$bad_grouping = 1;
				dbg_print($D_CENTERS, "No clear groupings!\n");
				last;
			}
		}
	} else {					# no knowledge about orig centers
		$bad_grouping = 1;		# hack -- classify as bad grouping.
	}

    my $sum_error;
    if(!$bad_grouping) {
		$sum_error = 0;
		#calculate sum-of-error
		for($i=0; $i<@centers; $i++) {
			$sum_error += d_sqd($centers[$i], $orig_dist[$closest_orig[$i]]);
		}
	}

#show results
my $sum_err_str = $bad_grouping ? "N/A" : sprintf("%e", $sum_error);
$num_iter = 1 if ($num_iter == 0);
dbg_print($D_EXTRA, "nc=num classes(actual), sig = sigma, np = num points\n"
		  ."nd = num dims, it =iterations, cost = avg. cost per iteration\n"
		  ."comp = avg. tests to find competitor, err = sum distances from resp. centers.\n"
		  ."time = elapsed time, mtd = method\n, lsz = leaf size"
		  ."dst = avg distortion/point\n"
		  ."cut = cutoff factor\n"
		  ."nspl = num-splits\n"
		  ."bic = avg. BIC score/point\n"
		  ."lvs = avg. #leaves scanned"
		  );
    printf("%2d(%d) nc %6.4f sig %d np %d nd %d it %9.2f cost %7.4f comp "
	   ."%e dst %6.2f time %6s mtd %3d lsz %f bwidth %g cut %d nspl %e bic %s splitstat"
	   ."\n",
	   $num_classes, scalar(@centers), $sigma, $num_points, $num_dims,
	   $num_iter, $sum_cost/$num_iter, $sum_fail/$num_iter,
	   $distortion/$num_points, $etime, 
	   ($kmeans_binary eq $kmeans_binary_default) ? $method : "$kmeans_binary#$method",
	   $leaf_size, $min_bwidth,
	   $cutoff_factor, $num_splits,
	   $bicscore/$num_points,
       $split_stat,
	   );

if($ANALYSE_CACHE) {
	dbg_print($D_EXTRA, "ch = %cache hit rate, cw=cache writes, co=cache hits+misses,"
			  ." cm=cache misses\n");
	if($cache_hit + $cache_miss > 0) {
		printf("%4.2f ch %d cw %d co %d cm\n", 100.0*$cache_hit/($cache_hit + $cache_miss),
			   $cache_write, $cache_hit + $cache_miss, $cache_miss);
	} else {
		printf("%4.2f ch %d cw %d co %d cm\n", 100.0*$cache_hit/1,
			   $cache_write, $cache_hit + $cache_miss, $cache_miss);
	}
	dbg_print($D_EXTRA, "hits\n");
	for my $nhits (sort {$a <=> $b} keys %hits) {
		print "h[$nhits] = $hits{$nhits}\n";
	}
	for my $nhits (sort {$a <=> $b} keys %unhits) {
		print "uh[$nhits] = $unhits{$nhits}\n";
	}
}

if($clusterfile) {
  # need to refine the clustering
  my $WRITER = gensym();
  my $READER = gensym();
  my $ERROR = gensym();
  open3($WRITER, $READER, $ERROR, "$bindir/$kmeans_binary membership -in $tmpclust > $clusterfile") or die("cant do process!");
  while(<$READER>) {
  }
  close($WRITER);
  close($READER);
  close($ERROR);
  unlink($tmpclust);
}
}

if($destroy_datafile) {
	unlink($datafile);
    unlink($unifile);
}
print "$datafile\n" if($create_datafile && !$run_kmeans);

#sub find_nearest_neighbor(@$point, @$candidates) {
sub find_nearest_neighbor {
	my ($point, $candidates) = @_;
	my $min_dist = d_sqd($point, $candidates->[0]);
	my $min_idx = 0;
 	for(my $i=0; $i<@$candidates; $i++) {
 		if(d_sqd($point, $candidates->[$i]) < $min_dist) {
 			$min_dist = &d_sqd($point, $candidates->[$i]);
 			$min_idx = $i;
 		}
 	}
	return $min_idx;
}

#sub d_sqd(@$v1, @$v2) {
sub d_sqd {
	my $sum = 0;
	die("Vectors not equal length!") if ($#{$_[0]} != $#{$_[1]});
	for(my $i=0; $i<=$#{$_[0]}; $i++) {
		my $diff = $_[0][$i] - $_[1][$i];
		$sum += $diff*$diff;
	}
	return $sum;
}

sub dbg_print
{
	local($level, $msg) = @_;
	if($level & $debug_level) {
		print STDERR $msg;
	}
}
