#======================================================================
# L F I T . P L
# doc: Sat Jul 31 11:24:47 1999
# dlm: Mon May 18 10:07:47 2015
# (c) 1999 A.M. Thurnherr
# uE-Info: 20 61 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# LFIT routine from Numerical Recipes adapted to ANTS
# HISTORY:
# Jul 31, 1999: - manually converted from c-source
# Aug 01, 1999: - changed funcs() interface
# Sep 26, 1999: - made sure right version of covsrt is used
# Jun 28, 2001: - re-added commented out code 'cause it's required
# on some perl versions
# Jan 5, 2012: - BUG: non-numeric x/y were not handled correctly;
# this was only easily apparent when the last
# record contained non-numeric values
# May 18, 2015: - added firstRec, lastRec parameters (V6.1)
# Notes:
# - x,y,sig are field numbers for data in $ants_
# - funcs is passed the current index & xfnr instead of the x-value
# - if sig is a negative number, -sig is used as constant input stddev
# - @a, @ia, @covar, &funcs passed as refs
# - chi square is returned
require "$ANTS/nrutil.pl";
require "$ANTS/covsrt.pl";
require "$ANTS/gaussj.pl";
sub lfit(@)
{
my($xfnr,$yfnr,$sig,$aR,$iaR,$covarR,$funcsR,$firstRec,$lastRec) = @_;
$firstRec = 0 unless defined($firstRec);
$lastRec = $#ants_ unless defined($lastRec);
my($i,$j,$k,$l,$m,$mfit); # int
my($ym,$wt,$sum,$sig2i,$chisq); # float
my(@beta,@afunc); # float[]
&matrix(\@beta,1,$#{$aR},1,1);
&vector(\@afunc,1,$#{$aR});
for ($j=1; $j<=$#{$aR}; $j++) {
$mfit++ if ($iaR->[$j]);
}
croak("lfit: no parameters to be fitted") if ($mfit == 0);
for ($j=1; $j<=$mfit; $j++) { # REQUIRED FOR SOME PERL VERSIONS!!!
for ($k=1;$ k<=$mfit; $k++) {
$covarR->[$j][$k] = 0;
}
$beta[$j][1] = 0;
}
for ($i=$firstRec; $i<=$lastRec; $i++) {
next if ($antsFlagged[$i]);
next unless numberp($ants_[$i][$xfnr]) && numberp($ants_[$i][$yfnr]);
&$funcsR($i,$xfnr,\@afunc);
$ym = $ants_[$i][$yfnr];
if ($mfit < $#{$aR}) {
for ($j=1; $j<=$#{$aR}; $j++) {
$ym -= $aR->[$j]*$afunc[$j] if (!$iaR->[$j]);
}
}
if ($sig > 0) { # field number
$sig2i = 1.0/($ants_[$i][$sig]*$ants_[$i][$sig]);
} else { # const value
$sig2i = 1.0/($sig*$sig);
}
for ($j=0,$l=1; $l<=$#{$aR}; $l++) {
if ($iaR->[$l]) {
$wt = $afunc[$l]*$sig2i;
for ($j++,$k=0,$m=1; $m<=$l; $m++) {
$covarR->[$j][++$k] += $wt*$afunc[$m] if ($iaR->[$m]);
}
$beta[$j][1] += $ym*$wt;
}
}
}
for ($j=2; $j<=$mfit; $j++) {
for ($k=1;$k<$j;$k++) {
$covarR->[$k][$j] = $covarR->[$j][$k];
# print(STDERR "covarR[$k][$j] = $covarR->[$k][$j]\n");
}
}
&gaussj($covarR,\@beta);
for ($j=0,$l=1;$l<=$#{$aR};$l++) {
$aR->[$l]=$beta[++$j][1] if ($iaR->[$l]);
}
for ($i=$firstRec; $i<=$lastRec; $i++) {
next if ($antsFlagged[$i]);
next unless numberp($ants_[$i][$xfnr]) && numberp($ants_[$i][$yfnr]);
&$funcsR($i,$xfnr,\@afunc);
for ($sum=0,$j=1; $j<=$#{$aR}; $j++) {
$sum += $aR->[$j]*$afunc[$j];
}
my($tmpval) = ($sig > 0) ?
($ants_[$i][$yfnr] - $sum) / $ants_[$i][$sig] :
($ants_[$i][$yfnr] - $sum) / -$sig;
$chisq += $tmpval * $tmpval;
}
&covsrt($covarR,$iaR);
return $chisq;
}