.lsfit.poly
author A.M. Thurnherr <athurnherr@yahoo.com>
Sat, 18 Nov 2017 18:55:17 -0500
changeset 32 4c6434bb5124
parent 30 1a1a12d5edc1
child 36 04e8cb4f8073
permissions -rw-r--r--
before Hamburg
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
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     4
#                    dlm: Thu Jan 10 16:47:06 2013
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     5
#                    (c) 1999 A.M. Thurnherr
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
     6
#                    uE-Info: 118 16 NIL 0 0 72 2 2 4 NIL ofnI
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)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    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
# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    34
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    35
#	$modelOpts			string of allowed options
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    36
#	$modelOptsUsage		usage information string for options
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    37
#	$modelMinArgs		min # of arguments of model
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    38
#	$modelArgsUsage		usage information string for arguments
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    39
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    40
# The following variables may be set later but not after &modelInit()
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    41
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    42
#	$modelNFit			number of params to fit in model
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    43
#	@nameA				symbolic names of model parameters
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    44
#
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    45
# You should call &antsDescription() for the -c option here
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    46
#
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
$modelOpts = "o:";
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    50
$modelOptsUsage = "-o)rder <n>";
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    51
$modelMinArgs = 0;
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    52
$modelArgsUsage = "";
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    53
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
# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    57
#					filenames after model arguments; this usually sets
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    58
#					@A (the model parameters) but these can later be
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    59
#					calculated heuristically during &modelInit()
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    60
#
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
sub modelUsage()
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    64
{
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    65
	my($c);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    66
	
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    67
	$modelNFit = &antsCardOpt($opt_o) + 1;						# order of polynomial
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    68
	die("$0 (.lsfit.poly): ERROR! -o required\n")
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    69
		unless (defined($opt_o) && $opt_o >= 0);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    70
	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    71
}
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
# &modelInit()		initializes model after reading of data
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    76
#
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
sub modelInit()
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    80
{
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    81
	for ($c=0; $c<$modelNFit; $c++) {				# init coefficients
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    82
		$A[$c+1] = nan;
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    83
		$nameA[$c+1] = "c$c";						# and names
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    84
	}
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
# &modelEvaluate(idx,xfnr,vals)	evaluate polynomial basis funs at x
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    90
#		idx				       	current index in @ants_
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    91
#		xfnr					field number of x field
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    92
#		vals					reference to return values (1-relative!)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    93
#
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
sub modelEvaluate($$$)
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    97
{
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    98
	my($idx,$xfnr,$valsR) = @_;
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
    99
	my($i);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   100
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   101
	$valsR->[1] = 1;
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   102
	for ($i=2; $i<=$#{$valsR}; $i++) {
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   103
		$valsR->[$i] = $valsR->[$i-1] * $ants_[$idx][$xfnr];
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   104
    }
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
# &modelCleanup()	cleans up after fitting but before output
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   109
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   110
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   111
sub modelCleanup() 
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   112
{
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   113
	return unless ($opt_o == 2);		# calculate loc of extremum on parabolic fits
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   114
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   115
	my($extX) = -$A[2] / (2 * $A[3]);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   116
	if ($A[3] > 0) {
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   117
		&antsInfo(".lsfit.poly: minimum at %.1f",$extX);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   118
	} elsif ($A[3] < 0) {
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   119
		&antsInfo(".lsfit.poly: maximum at %.1f",$extX);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   120
	} else {
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   121
		&antsInfo(".lsfit.poly: saddle point at %.1f",$extX);
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
diff changeset
   122
	}
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
1;