#======================================================================
#                    A M O E B A . P L 
#                    doc: Wed Aug 23 05:11:48 2006
#                    dlm: Wed Aug 23 23:52:12 2006
#                    (c) 2006 A.M. Thurnherr
#                    uE-Info: 88 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

# perlified amoeba implementation of NR code

# NOTES:
#	- 0-based arrays
#	- amoeba returns undef if NMAX is exceeded and # of evals otherwise

use strict;

sub amotry($$$$$$)
{
	my($pR,$yR,$psumR,$funR,$ihi,$fac) = @_;
	my(@ptry);

	my($ndim) = scalar(@{$pR->[0]});
	my($fac1) = (1-$fac) / $ndim;
	my($fac2) = $fac1 - $fac;
	
	for (my($j)=0; $j<$ndim; $j++) {
		$ptry[$j] = $psumR->[$j]*$fac1 - $pR->[$ihi][$j]*$fac2;
	}
	my($ytry) = &$funR(@ptry);
	if ($ytry < $yR->[$ihi]) {
		$yR->[$ihi] = $ytry;
		for (my($j)=0; $j<$ndim; $j++) {
			$psumR->[$j] += $ptry[$j] - $pR->[$ihi][$j];
			$pR->[$ihi][$j] = $ptry[$j];
		}
	}
	return $ytry;
}

sub get_psum($$)
{
	my($pR,$psumR) = @_;

	for (my($j)=0; $j<@{$pR->[0]}; $j++) {
		my($sum);
		for ($sum=my($i)=0; $i<=@{$pR->[0]}; $i++) {
			$sum += $pR->[$i][$j];
		}
		$psumR->[$j] = $sum;
	}
}

sub amoeba($$$$)
{
	my($pR,$yR,$ftol,$funR,$NMAX) = @_;
	my($nfunk) = 0;
	my($ndim) = scalar(@{$pR->[0]});
	my(@psum);
	
	&get_psum($pR,\@psum);

	while (1) {
		my($i,$ihi,$inhi,$j);
		my($sum);
		
		my($ilo) = 0;
		if ($yR->[0] > $yR->[1]) {
			$ihi = 0; $inhi = 1;
		} else {
			$ihi = 1; $inhi = 0;
		}

		for ($i=0; $i<$ndim+1; $i++) {
			if ($yR->[$i] <= $yR->[$ilo]) {
				$ilo = $i;
			}
			if ($yR->[$i] > $yR->[$ihi]) {
				$inhi = $ihi;
				$ihi  = $i;
			} elsif ($yR->[$i] > $yR->[$inhi] && $i != $ihi) {
				$inhi = $i;
			}
		}
		print(STDERR "best = $yR->[$ilo]\n");

		my($rtol) = 2 * abs($yR->[$ihi] - $yR->[$ilo]) /
						(abs($yR->[$ihi]) + abs($yR->[$ilo]));
		if ($rtol < $ftol) {
			my($tmp) = $yR->[0]; $yR->[0] = $yR->[$ilo]; $yR->[$ilo] = $tmp;
			for ($i=0; $i<$ndim; $i++) {
				my($tmp) = $pR->[1][$i]; $pR->[1][$i] = $pR->[$ilo][$i];
				$pR->[$ilo][$i] = $tmp;
			}
			return $nfunk;
		}
		
		return undef if ($nfunk >= $NMAX);
		$nfunk += 2;
		
		my($ytry) = amotry($pR,$yR,\@psum,$funR,$ihi,-1);
		if ($ytry <= $yR->[$ilo]) {
			$ytry = amotry($pR,$yR,\@psum,$funR,$ihi,2);
		} elsif ($ytry >= $yR->[$inhi]) {
			my($ysave) = $yR->[$ihi];
			$ytry = amotry($pR,$yR,\@psum,$funR,$ihi,0.5);
			if ($ytry >= $ysave) {
				for ($i=0; $i<$ndim+1; $i++) {
					if ($i != $ilo) {
						for ($j=0; $j<$ndim; $j++) {
							$pR->[$i][$j] = $psum[$j] =
								0.5 * ($pR->[$i][$j] + $pR->[$ilo][$j]);
						}
						$yR->[$i] = &$funR(@psum);
					}
				}
				$nfunk += $ndim;
				&get_psum($pR,\@psum);
			}
		} else {
			--$nfunk;
		}
	}
}

1;