#======================================================================

#                    M R Q M I N . P L 

#                    doc: Wed Feb 24 15:10:22 1999

#                    dlm: Tue Aug 22 22:05:43 2006

#                    (c) 1999 A.M. Thurnherr

#                    uE-Info: 15 67 NIL 0 0 72 2 2 4 NIL ofnI

#======================================================================

# MRQMIN routine from Numerical Recipes adapted to ANTS

# NB: based on 1st edtion of NR!!!!

# HISTORY:

#	Mar 11, 1999: - created

#	Sep 27, 1999: - adapted to allow for new version of covsrt.pl as well

#	Aug 22, 2006: - changed require from covsrt_old.pl to covsrt.pl

# Notes:

#	- x,y,sig are field numbers for data in $ants_

#	- if sig is a negative number, -sig is used as constant input stddev

#	- @A, @listA, @alpha, @covar, $chiSq, &funcs, $alamda passed as refs

require "$ANTS/mrqcof.pl";

require "$ANTS/gaussj.pl";

require "$ANTS/covsrt.pl";

require "$ANTS/nrutil.pl";

{													# static scope

	my(@da,@atry,@oneda,@beta,$oChiSq);

	sub mrqmin($$$$$$$$$$)

	{

		my($xfnr,$yfnr,$sig,$AR,$listAR,$covarR,

		   $alphaR,$chiSqR,$funcsR,$alamdaR) = @_;

		my($k,$kk,$j,$ihit);

		if ($$alamdaR < 0.0) {

			&matrix(\@oneda,1,$#{$AR},1,1);

			&vector(\@atry,1,$#{$AR});

			&vector(\@da,1,$#{$AR});

			&vector(\@beta,1,$#{$AR});

			$kk = $#{$listAR}+1;

			for ($j=1; $j<=$#{$AR}; $j++) {

				$ihit = 0;

				for ($k=1; $k<=$#{$listAR}; $k++) {

					if ($listAR->[$k] == $j) { $ihit++; }

				}

				if ($ihit == 0) {

					$listAR->[$kk++] = $j;

				} elsif ($ihit > 1) {

					croak("Bad listA permutation in MRQMIN-1");

				}

			}

			if ($kk != $#{$AR}+1) {

				for ($ei=1; $ei<=$#{$listAR}; $ei++) {

					print(STDERR "listA[$ei] = $listAR->[$ei]\n");

				}

				croak("Bad listA permutation in MRQMIN-2 " .

					"($kk != $#{$AR}+1)");

			}

			$$alamdaR = 0.001;

			&mrqcof($xfnr,$yfnr,$sig,$AR,$listAR,$alphaR,

					\@beta,$chiSqR,$funcsR);

			$oChiSq = $$chiSqR;

		}

		for ($j=1; $j<=$#{$listAR}; $j++) {

			for ($k=1; $k<=$#{$listAR}; $k++) {

				$covarR->[$j][$k] = $alphaR->[$j][$k];

#				print(STDERR "covar[$j][$k] = $covarR->[$j][$k]\n");

			}

			$covarR->[$j][$j] = $alphaR->[$j][$j]*(1.0+$$alamdaR);

			$oneda[$j][1] = $beta[$j];

		}

		&gaussj($covarR,\@oneda);

		for ($j=1; $j<=$#{$listAR}; $j++) {

			$da[$j] = $oneda[$j][1];

		}

		if ($$alamdaR == 0.0) {

			&covsrt($covarR,$listAR);

			return;

		}

		for ($j=1; $j<=$#{$AR}; $j++) { $atry[$j] = $AR->[$j]; }

		for ($j=1; $j<=$#{$listAR}; $j++) {

			$atry[$listAR->[$j]] = $AR->[$listAR->[$j]]+$da[$j];

		}

		&mrqcof($xfnr,$yfnr,$sig,\@atry,$listAR,$covarR,\@da,$chiSqR,$funcsR);

		if ($$chiSqR < $oChiSq) {

			$$alamdaR *= 0.1;

			$oChiSq = $$chiSqR;

			for ($j=1; $j<=$#{$listAR}; $j++) {

				for ($k=1; $k<=$#{$listAR}; $k++) {

					$alphaR->[$j][$k] = $covarR->[$j][$k];

				}

				$beta[$j] = $da[$j];

				$AR->[$listAR->[$j]] = $atry[$listAR->[$j]];

			}

		} else {

			$$alamdaR *= 10.0;

			$$chiSqR = $oChiSq;

		}

		return;

	}

} # end of static scope

1;