#!/afs/athena/contrib/perl5/perl

use Plotter;

require Gaussian;
import Gaussian;

$N = 100;
$alpha = $ARGV[0];
$alpha = .1 unless $alpha;
$beta = $alpha;
$npops = 50;
srand(17);
#srand;

$pl = new Plotter(500,200);
$pl->addPlot("p", 0);
$pl->addPlot("a", 0);

for $i (0..($N-1)) {
    $j[$i] = &gauss_rand;
}

for $pop (1..$npops){
    my @spin = [];
    for $i (0..($N-1)) {
	if(rand() > .5){
	    $spin[$i] = 1;
	}
	else{
	    $spin[$i] = -1;
	}
    }
    push(@pops, \@spin);
}

    
$emin = &emin;
print("Emin: $emin\n");

while(1){
    $prob = 0;
    $sume = $ect = 0;
    foreach $pn (0..$#pops){
#	print("Population $pn with @{$pops[$pn]}\n");
	$e = &energy(@{$pops[$pn]});
#	print("Population $pn has energy $e\n");
	$p = exp(-1*$beta*($e-$emin));
	$prob+= $p;
	$prob[$pn] = $prob;
	$beste = ($e < $beste) ? $e:$beste;
	$sume+=$e;
	$ect++;
    }
    $avge = $sume/$ect;
    print("$beste:$avge ($emin)\n");
    $pl->plot("p", -1*(($beste-$emin)/$emin));
    $pl->plot("a", -1*(($avge-$emin)/$emin));
    $pl->update();

    for $pop (1..$npops){
	$mom = &pickparent;
	$dad = &pickparent;
	$cross = $N*rand();
	my @ch = [];
	@ch = (@{$mom}[0..$cross], @{$dad}[$cross+1..$N-1]);
	$mutate = $N*rand();
	$ch[$mutate] = ($ch[$mutate]==1) ? -1:1;
	push(@newpops, \@ch);
    }
    undef @pops;
    @pops = @newpops;
    undef @newpops;
}

sub pickparent {
    $pv = $prob*rand(1);
    for $pn (0..$#pops){
	if($pv <= $prob[$pn]){
	    return $pops[$pn];
	}
    }
    return $pops[$#pops];
}

sub emin {
    my $emin = 0;
    for $i (0..($N-1)){
	$emin+= abs($j[$i]);
    }
    return -1*$emin;
}

sub energy {
    my @spin = @_;
    my $energy=0;
    for $i (0..($N-2)){
	$energy += $j[$i]*$spin[$i]*$spin[$i+1];
    }
    $energy += $j[$N-1]*$spin[0]*$spin[$N-1];
    return -1*$energy;
}