#======================================================================
#                    . L S F I T . P O L Y 
#                    doc: Wed Feb 24 09:40:06 1999
#                    dlm: Mon May 11 11:54:03 2020
#                    (c) 1999 A.M. Thurnherr
#                    uE-Info: 31 55 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

# What you need to provide if you wanna fit a different
# linear 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 basis funs at a given x value; each
#	  y value must be stored in @A (beginning with
#	  A[1]!!!).
#	- the interface is documented between +++++++ lines

# fit polynomial (sum of A_i x^i) to data
# NB: - preferable to [./.lmfit.poly]

# HISTORY:
#	Jul 31, 1999: - adapted from [./.lmfit.poly]
#	Aug 01, 1999: - changed &modelEvaluate() interface
#	Aug 02, 1999: - added &antsDescription()
#   Mar 17, 2001: - param->arg
#	Jul 12, 2004: - made poly-order argument into -o option
#	Jul 28, 2006: - Version 3.3 [HISTORY]
#	Sep 19, 2011: - moved part of the usage code into init() to allow use in [pgram]
#	Jan 10, 2013: - added extremum output when fitting parabola (-o 2)
#	May 13, 2018: - BUG: replaced opt_o with modelNFit in &modelCleanup()
#	May 11, 2020: - increased extremum output precision

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# 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
#
# You should call &antsDescription() for the -c option here
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

$modelOpts = "o:";
$modelOptsUsage = "-o)rder <n>";
$modelMinArgs = 0;
$modelArgsUsage = "";

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &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()
{
	my($c);
	
	$modelNFit = &antsCardOpt($opt_o) + 1;						# order of polynomial
	die("$0 (.lsfit.poly): ERROR! -o required\n")
		unless (defined($opt_o) && $opt_o >= 0);
	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
}

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

sub modelInit()
{
	for ($c=0; $c<$modelNFit; $c++) {				# init coefficients
		$A[$c+1] = nan;
		$nameA[$c+1] = "c$c";						# and names
	}
}

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelEvaluate(idx,xfnr,vals)	evaluate polynomial basis funs at x
#		idx				       	current index in @ants_
#		xfnr					field number of x field
#		vals					reference to return values (1-relative!)
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub modelEvaluate($$$)
{
	my($idx,$xfnr,$valsR) = @_;
	my($i);

	$valsR->[1] = 1;
	for ($i=2; $i<=$#{$valsR}; $i++) {
		$valsR->[$i] = $valsR->[$i-1] * $ants_[$idx][$xfnr];
    }
}

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

sub modelCleanup() 
{
	return unless ($0 eq 'lsfit' && $modelNFit == 3);		# calculate loc of extremum on parabolic fits with lsfit only

	my($extX) = -$A[2] / (2 * $A[3]);
	if ($A[3] > 0) {
		&antsInfo(".lsfit.poly: minimum at %g",$extX);
	} elsif ($A[3] < 0) {
		&antsInfo(".lsfit.poly: maximum at %g",$extX);
	} else {
		&antsInfo(".lsfit.poly: saddle point at %g",$extX);
	}
}

1;