# This library file contains Amoeba n-D Minimisation routine for Perl
# $Id: Amoeba.pm,v 1.2 1995/12/24 12:37:46 willijar Exp $
# $Id: Amoeba.pm,v 1.2 1995/12/24 12:37:46 willijar Exp $
package Math::Amoeba;
use strict;
use warnings;
our $VERSION = 0.05;
use Carp;
use constant TINY => 1e-16;
use Exporter;
our @ISA=qw(Exporter);
our @EXPORT_OK=qw(ConstructVertices EvaluateVertices Amoeba MinimiseND);
=head1 NAME
Math::Amoeba - Multidimensional Function Minimisation
=head1 SYNOPSIS
use Math::Amoeba qw(ConstructVertices EvaluateVertices Amoeba MinimiseND);
my ($vertice,$y)=MinimiseND(\@guess,\@scales,\&func,$tol,$itmax,$verbose);
my @vertices=ConstructVertices(\@vector,\@offsets);
my @y=EvaluateVertices(\@vertices,\&func);
my ($vertice,$y)=Amoeba(\@vertices,\@y,\&func,$tol,$itmax,$verbose);
=head1 DESCRIPTION
This is an implimenation of the Downhill Simpex Method in
Multidimensions (Nelder and Mead) for finding the (local) minimum of a
function. Doing this in Perl makes it easy for that function to
actually be the output of another program such as a simulator.
Arrays and the function are passed by reference to the routines.
=over
=item C<MinimiseND>
The simplest use is the B<MinimiseND> function. This takes a reference
to an array of guess values for the parameters at the function
minimum, a reference to an array of scales for these parameters
(sensible ranges around the guess in which to look), a reference to
the function, a convergence tolerence for the minimum, the maximum
number of iterations to be taken and the verbose flag (default ON).
It returns an array consisting of a reference to the function parameters
at the minimum and the value there.
=item C<Amoeba>
The B<Amoeba> function is the actual implimentation of the Downhill
Simpex Method in Multidimensions. It takes a reference to an array of
references to arrays which are the initial n+1 vertices (where n is
the number of function parameters), a reference to the function
valuation at these vertices, a reference to the function, a
convergence tolerence for the minimum, the maximum number of
iterations to be taken and the verbose flag (default ON).
It returns an array consisting of a reference to the function parameters
at the minimum and the value there.
=item C<ConstructVertices>
The B<ConstructVertices> is used by B<MinimiseND> to construct the
initial vertices for B<Amoeba> as the initial guess plus the parameter
scale parameters as vectors along the parameter axis.
=item C<EvaluateVertices>
The B<EvaluateVertices> takes these set of vertices, calling the
function for each one and returning the vector of results.
=back
=head1 EXAMPLE
use Math::Amoeba qw(MinimiseND);
sub afunc {
my ($a,$b)=@_;
print "$a\t$b\n";
return ($a-7)**2+($b+3)**2;
}
my @guess=(1,1);
my @scale=(1,1);
($p,$y)=MinimiseND(\@guess,\@scale,\&afunc,1e-7,100);
print "(",join(',',@{$p}),")=$y\n";
produces the output
(6.99978191653352,-2.99981241563247)=1.00000008274829
=head1 LICENSE
PERL
=cut
my ($ALPHA,$BETA,$GAMMA)=(1.0,0.5,2.0);
sub MinimiseND {
my ($guesses,$scales,$func,$tol,$itmax, $verbose)=@_;
my @p=ConstructVertices($guesses,$scales);
my @y=EvaluateVertices(\@p,$func);
return Amoeba(\@p,\@y,$func,$tol,$itmax, $verbose);
}
sub ConstructVertices {
# given 2 vector references constructs an amoeba
# returning the vertices
my ($vector,$ofs)=@_;
my $n=$#{$vector};
my @vector=@{$vector};
my (@p,@y,$i);
$p[0]=[]; @{$p[0]}=@{$vector};
for($i=0; $i<=$n; $i++) {
my $v=[]; @{$v}=@{$vector};
$v->[$i]+=$ofs->[$i];
$p[$i+1]=$v;
}
return @p;
}
sub EvaluateVertices {
# evaluates functions for all vertices of the amoeba
my ($p,$func)=@_;
my ($i,@y);
for($i=0; $i<=$#{$p}; $i++) {
$y[$i]=&$func(@{$p->[$i]});
}
return @y;
}
sub Amoeba {
my ($p,$y,$func,$ftol,$itmax, $verbose)=@_;
my $n=$#{$p}; # no points
# Default parameters
$verbose = (defined($verbose)) ? $verbose : 1;
if (!$itmax) { $itmax=200; }
if (!$ftol) { $ftol=1e-6; }
# Member variables
my ($i,$j);
my $iter=0;
my ($ilo, $inhi, $ihi);
my ($pbar, $pr, $pe, $pc, $ypr, $ype, $ypc);
# To control the recalculation of centroid
my $recalc = 1;
my $ihi_o;
# Loop until any of stopping conditions hit
while (1)
{
($ilo, $inhi, $ihi) = _FindMarkers($y);
# Stopping conditions
my $rtol = 2*abs($y->[$ihi]-$y->[$ilo])/(abs($y->[$ihi])+abs($y->[$ilo])+TINY);
if ($rtol<$ftol) { last; }
if ($iter++>$itmax) {
carp "Amoeba exceeded maximum iterations\n" if ($verbose);
last;
}
# Determine the Centroid
if ($recalc) {
$pbar = _CalcCentroid($p, $ihi);
} else {
_AdjustCentroid($pbar, $p, $ihi_o, $ihi);
}
# Reset the re-calculation flag, and remember the current highest
$recalc = 0;
# Determine the reflection point, evaluate its value
$pr = _CalcReflection($pbar, $p->[$ihi], $ALPHA);
$ypr = &$func(@$pr);
# if it gives a better value than best point, try an
# additional extrapolation by a factor gamma, accept best
if ($ypr < $y->[$ilo]) {
$pe = _CalcReflection($pbar, $pr, -$GAMMA);
$ype=&$func(@$pe);
if ($ype < $y->[$ilo]) {
$p->[$ihi] = $pe; $y->[$ihi] = $ype;
}
else {
$p->[$ihi] = $pr; $y->[$ihi] = $ypr;
}
}
# if reflected point worse than 2nd highest
elsif ($ypr >= $y->[$inhi]) {
# if it is better than highest, replace it
if ($ypr < $y->[$ihi] ) {
$p->[$ihi] = $pr; $y->[$ihi] = $ypr;
}
# look for an intermediate lower point by performing a
# contraction of the simplex along one dimension
$pc = _CalcReflection($pbar, $p->[$ihi], -$BETA);
$ypc = &$func(@$pc);
# if contraction gives an improvement, accept it
if ($ypc < $y->[$ihi]) {
$p->[$ihi] = $pc; $y->[$ihi] = $ypc;
}
# otherwise cant seem to remove high point
# so contract around lo (best) point
else {
for($i=0; $i<=$n; $i++) {
if ($i!=$ilo) {
$p->[$i] = _CalcReflection($p->[$ilo], $p->[$i], -$BETA);
$y->[$i] = &$func(@{$p->[$i]});
}
}
$recalc = 1;
}
}
# if reflected point is in-between lowest and 2nd highest
else {
$p->[$ihi] = $pr; $y->[$ihi] = $ypr;
}
# Remember the replacing position and its value
$ihi_o = $ihi;
}
return ($p->[$ilo],$y->[$ilo]);
}
# Helper function - find the lowest, 2nd highest and highest position
sub _FindMarkers
{
my $y = shift;
my ($ilo, $inhi, $ihi);
my ($i, $n);
$n = @$y - 1;
$ilo=0;
if ($y->[0]>$y->[1]) {
$ihi=0; $inhi=1;
}
else {
$ihi=1; $inhi=0;
}
for($i = 0; $i <= $n; $i++) {
if ($y->[$i] < $y->[$ilo]) { $ilo = $i; }
if ($y->[$i] > $y->[$ihi]) { $inhi = $ihi; $ihi = $i; }
elsif ($y->[$i] > $y->[$inhi] && $ihi != $i) { $inhi = $i; }
}
return ($ilo, $inhi, $ihi);
}
# Determine the centroid (except the highest point)
sub _CalcCentroid
{
my ($p, $ihi) = @_;
my ($i, $j, $n);
$n = @$p - 1;
my $pbar = [];
for($j=0; $j<$n; $j++) {
for($i=0; $i<=$n; $i++) {
if ($i!=$ihi) {
$pbar->[$j] += $p->[$i][$j];
}
}
$pbar->[$j] /= $n;
}
return $pbar;
}
# Adjust the centroid only
sub _AdjustCentroid
{
my ($pbar, $p, $ihi_o, $ihi) = @_;
my ($j, $n);
$n = @$pbar;
if ($ihi_o != $ihi) {
for($j=0; $j<$n; $j++) {
$pbar->[$j] += ($p->[$ihi_o][$j] - $p->[$ihi][$j]) / $n;
}
}
}
# Determine the reflection point
sub _CalcReflection
{
my ($p1, $p2, $scale) = @_;
my $j;
my $n = @$p1;
my $pr = [];
for($j=0; $j<$n; $j++) {
$pr->[$j] = $p1->[$j] + $scale*($p1->[$j]-$p2->[$j]);
}
return $pr;
}
return 1;
__END__
=head1 HISTORY
See "REAME".
=head1 BUGS
Let me know.
=head1 AUTHOR
John A.R. Williams <J.A.R.Williams@aston.ac.uk>
Tom Chau <tom@cpan.org>
=head1 SEE ALSO
"Numerical Recipies: The Art of Scientific Computing"
W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling.
Cambridge University Press. ISBN 0 521 30811 9.
=head1 COPYRIGHT
Copyright (c) 1995 John A.R. Williams. All rights reserved.
This program is free software; you can redistribute it and/or
modify it under the same terms as Perl itself.
Since 2005, this module was co-developed with Tom.