#======================================================================
#                    . L M F I T . G A U S S 
#                    doc: Wed Feb 24 09:40:06 1999
#                    dlm: Fri Jul 28 13:32:35 2006
#                    (c) 1999 A.M. Thurnherr
#                    uE-Info: 35 51 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

# What you need to provide if you wanna fit a different
# model function to your data:
#	- a number of global variables to be set during loading
#	- a number of subs to perform admin tasks (usage, init, ...)
#	- a sub to evaluate the model function which is to be fitted using
#	  a number of pararams which are all stored in @A (beginning at
#	  A[1]!!!). You also need to return the partial derivatives of
#	  the model function wrt all params.
#	- the interface is documented between +++++++ lines

# Gauss data model (i.e. fit Gaussian curve)
# NB: - fitting is rather sensitive to the input parameters, thus
# 	    a heuristic has been added to guess them (by setting them
# 	    to NaN)
# 	  - another fickle parameter is the y-offset (zero line); thus
# 	    a heuristics has been added for this one as well
#	  - the parameters are peak, mean, standard deviation

# HISTORY:
#	Feb 24, 1999: - created together with [./cfit]
#	Feb 25, 1999: - cosmetic changes
#	Jul 31, 1999: - parameter typecheck
#	Oct 04, 1999: - changed param names
#	Oct 05, 1999: - improved heuristics
#				  - changed e-scale to sigma
#   Mar 17, 2001: - param->arg
#	Jul 28, 2006: - Version 3.3 [HISTORY]; &isnan()

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
#
#	$modelOpts			string of allowed options
#	$modelOptsUsage		usage information string for options
#	$modelMinArgs		min # of arguments of model
#	$modelArgsUsage		usage information string for arguments
#
# The following variables may be set later but not after &modelInit()
#
#	$modelNFit			number of params to fit in model
#	@nameA				symbolic names of model parameters
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

$modelOpts = "y";
$modelOptsUsage = "[-y)shift]";
$modelMinArgs = 0;
$modelArgsUsage = "[peak guess [mean guess [sigma guess]]]";
$modelNFit = 3;			
$nameA[1] = "peak";
$nameA[2] = "mean";		
$nameA[3] = "sigma";	

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
#					filenames after model arguments; this usually sets
#					@A (the model parameters) but these can later be
#					calculated heuristically during &modelInit()
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub modelUsage()
{
	$A[1] = nan; $A[2] = nan; $A[3] = nan;				# usage
	$A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
	$A[2] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
	$A[3] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
}

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelInit()		initializes model after reading of data
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub modelInit()
{
	my($i,$j,$ymin,$ymax,$xatymax);

#	--------------------------------------------------
#	heuristics for initial model param values
#	--------------------------------------------------

	$ymin = 1e33, $ymax = -1e33, $xatymax = 0;
	for ($i=0; $i<=$#ants_; $i++) {
		next if ($antsFlagged[$i]);
		$ymin = $ants_[$i][$yfnr]
			if ($ants_[$i][$yfnr] < $ymin);
		$ymax = $ants_[$i][$yfnr], $xatymax = $ants_[$i][$xfnr]
			if ($ants_[$i][$yfnr] > $ymax);
	}
	$A[1] = $ymax - $ymin if isnan($A[1]);			# peak guess
	$A[2] = $xatymax if isnan($A[2]);				# mean guess
	if (isnan($A[3])) {								# sigma guess
		for ($i=1;
			 $i<=$#ants_ && !$antsFlagged[$i]
						 && $ants_[$i][$yfnr]-$ymin<0.36*$A[1];
			 $i++) {}				   
		for ($j=$#ants_;
			 $j>=1 && !$antsFlagged[$i]
				   && $ants_[$j][$yfnr]-$ymin < 0.36*$A[1];
			 $j--) {}
		$A[3] = abs($ants_[$i][$xfnr]-$ants_[$j][$xfnr]) / 2.0;
		$A[3] *= 0.71;								# scale by 1/sqrt(2)
		if ($A[3] == 0) {
			&antsInfo("$model: sigma heuristic failed (set to 1)!");
			$A[3] = 1;
	    }
	}

#	--------------------------------------------------
#	y shift (-y option)
#	--------------------------------------------------

	if ($opt_y) {
		for ($i=1; $i<=$#ants_; $i++) {
			$ants_[$i][$yfnr] -= $ymin;
		}
	}

}

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelEvaluate(x,A,dyda)	evaluate sum of Gaussians (p.528) at x
#		x					x value (NOT xfnr)
#		A					reference to @A
#		dyda				reference to array for partial derivatives
#							(wrt individaul params in @A)
#		<ret val>			y value
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub modelEvaluate($$$)
{
	my($x,$AR,$dydaR) = @_;
	my($i,$fac,$ex,$arg,$sqrt2sig);
	my($y) = 0;

	for ($i=1; $i < $#{$AR}; $i+=3) {
		$sqrt2sig = (1.4142135623731*$AR->[$i+2]);
		$arg = ($x - $AR->[$i+1]) / $sqrt2sig;
		$ex  = exp(-$arg*$arg);
		$fac = $AR->[$i] * $ex * 2*$arg;
		$y += $AR->[$i] * $ex;
		
		$dydaR->[$i]   = $ex;
		$dydaR->[$i+1] = $fac / $sqrt2sig;
		$dydaR->[$i+2] = $fac * $arg / $sqrt2sig;
	}
	return $y;
}

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# &modelCleanup()	cleans up after fitting but before output
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub modelCleanup()
{
	if ($opt_y) {
		$A[1] += $ymin;
		for ($i=1; $i<=$#ants_; $i++) {
			$ants_[$i][$yfnr] += $ymin;
		}
	}
}