.lsfit.poly
author A.M. Thurnherr <athurnherr@yahoo.com>
Tue, 27 Nov 2018 12:58:21 -0500
changeset 36 04e8cb4f8073
parent 30 1a1a12d5edc1
child 42 22f5d5d35236
permissions -rw-r--r--
V7.1
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
30
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     1
#======================================================================
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     2
#                    . L S F I T . P O L Y 
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     3
#                    doc: Wed Feb 24 09:40:06 1999
36
A.M. Thurnherr <athurnherr@yahoo.com>
parents: 30
diff changeset
     4
#                    dlm: Sun May 13 08:25:37 2018
30
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     5
#                    (c) 1999 A.M. Thurnherr
36
A.M. Thurnherr <athurnherr@yahoo.com>
parents: 30
diff changeset
     6
#                    uE-Info: 114 36 NIL 0 0 72 2 2 4 NIL ofnI
30
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     7
#======================================================================
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     8
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     9
# What you need to provide if you wanna fit a different
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    10
# linear model function to your data:
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    11
#	- a number of global variables to be set during loading
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    12
#	- a number of subs to perform admin tasks (usage, init, ...)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    13
#	- a sub to evaluate the basis funs at a given x value; each
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    14
#	  y value must be stored in @A (beginning with
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    15
#	  A[1]!!!).
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    16
#	- the interface is documented between +++++++ lines
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    17
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    18
# fit polynomial (sum of A_i x^i) to data
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    19
# NB: - preferable to [./.lmfit.poly]
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    20
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    21
# HISTORY:
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    22
#	Jul 31, 1999: - adapted from [./.lmfit.poly]
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    23
#	Aug 01, 1999: - changed &modelEvaluate() interface
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    24
#	Aug 02, 1999: - added &antsDescription()
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    25
#   Mar 17, 2001: - param->arg
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    26
#	Jul 12, 2004: - made poly-order argument into -o option
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    27
#	Jul 28, 2006: - Version 3.3 [HISTORY]
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    28
#	Sep 19, 2011: - moved part of the usage code into init() to allow use in [pgram]
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    29
#	Jan 10, 2013: - added extremum output when fitting parabola (-o 2)
36
A.M. Thurnherr <athurnherr@yahoo.com>
parents: 30
diff changeset
    30
#	May 13, 2018: - BUG: replaced opt_o with modelNFit in &modelCleanup()
30
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    31
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    32
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    33
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    34
# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    35
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    36
#	$modelOpts			string of allowed options
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    37
#	$modelOptsUsage		usage information string for options
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    38
#	$modelMinArgs		min # of arguments of model
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    39
#	$modelArgsUsage		usage information string for arguments
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    40
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    41
# The following variables may be set later but not after &modelInit()
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    42
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    43
#	$modelNFit			number of params to fit in model
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    44
#	@nameA				symbolic names of model parameters
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    45
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    46
# You should call &antsDescription() for the -c option here
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    47
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    48
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    49
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    50
$modelOpts = "o:";
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    51
$modelOptsUsage = "-o)rder <n>";
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    52
$modelMinArgs = 0;
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    53
$modelArgsUsage = "";
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    54
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    55
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    56
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    57
# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    58
#					filenames after model arguments; this usually sets
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    59
#					@A (the model parameters) but these can later be
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    60
#					calculated heuristically during &modelInit()
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    61
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    62
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    63
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    64
sub modelUsage()
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    65
{
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    66
	my($c);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    67
	
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    68
	$modelNFit = &antsCardOpt($opt_o) + 1;						# order of polynomial
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    69
	die("$0 (.lsfit.poly): ERROR! -o required\n")
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    70
		unless (defined($opt_o) && $opt_o >= 0);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    71
	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    72
}
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    73
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    74
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    75
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    76
# &modelInit()		initializes model after reading of data
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    77
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    78
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    79
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    80
sub modelInit()
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    81
{
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    82
	for ($c=0; $c<$modelNFit; $c++) {				# init coefficients
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    83
		$A[$c+1] = nan;
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    84
		$nameA[$c+1] = "c$c";						# and names
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    85
	}
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    86
}
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    87
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    88
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    89
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    90
# &modelEvaluate(idx,xfnr,vals)	evaluate polynomial basis funs at x
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    91
#		idx				       	current index in @ants_
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    92
#		xfnr					field number of x field
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    93
#		vals					reference to return values (1-relative!)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    94
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    95
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    96
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    97
sub modelEvaluate($$$)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    98
{
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    99
	my($idx,$xfnr,$valsR) = @_;
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   100
	my($i);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   101
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   102
	$valsR->[1] = 1;
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   103
	for ($i=2; $i<=$#{$valsR}; $i++) {
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   104
		$valsR->[$i] = $valsR->[$i-1] * $ants_[$idx][$xfnr];
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   105
    }
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   106
}
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   107
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   108
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   109
# &modelCleanup()	cleans up after fitting but before output
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   110
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   111
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   112
sub modelCleanup() 
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   113
{
36
A.M. Thurnherr <athurnherr@yahoo.com>
parents: 30
diff changeset
   114
	return unless ($0 eq 'lsfit' && $modelNFit == 3);		# calculate loc of extremum on parabolic fits with lsfit only
30
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   115
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   116
	my($extX) = -$A[2] / (2 * $A[3]);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   117
	if ($A[3] > 0) {
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   118
		&antsInfo(".lsfit.poly: minimum at %.1f",$extX);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   119
	} elsif ($A[3] < 0) {
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   120
		&antsInfo(".lsfit.poly: maximum at %.1f",$extX);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   121
	} else {
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   122
		&antsInfo(".lsfit.poly: saddle point at %.1f",$extX);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   123
	}
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   124
}
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   125
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   126
1;