|
1 #====================================================================== |
|
2 # . L S F I T . B I L I N |
|
3 # doc: Wed Feb 24 09:40:06 1999 |
|
4 # dlm: Fri Jul 28 13:36:36 2006 |
|
5 # (c) 1999 A.M. Thurnherr |
|
6 # uE-Info: 32 41 NIL 0 0 72 2 2 4 NIL ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # What you need to provide if you wanna fit a different |
|
10 # linear model function to your data: |
|
11 # - a number of global variables to be set during loading |
|
12 # - a number of subs to perform admin tasks (usage, init, ...) |
|
13 # - a sub to evaluate the basis funs at a given x value; each |
|
14 # y value must be stored in @A (beginning with |
|
15 # A[1]!!!). |
|
16 # - the interface is documented between +++++++ lines |
|
17 |
|
18 # fit bi-linear function to data, i.e. y = A + B*x1 + C*x2 |
|
19 |
|
20 # HISTORY: |
|
21 # Aug 01, 1999: - adapted from [.lsfit.poly] |
|
22 # Aug 02, 1999: - added &antsDescription() |
|
23 # Sep 26, 1999: - cosmetics |
|
24 # - added vars & covars |
|
25 # Sep 27, 1999: - changed from covar to sigmas |
|
26 # Oct 01, 1999: - cosmetics |
|
27 # Oct 06, 1999: - added -l |
|
28 # Mar 17, 2001: - param->arg |
|
29 # May 2, 2001: - updated doc |
|
30 # Nov 17, 2005: - commented out antsDescription() |
|
31 # - updated stats on -p |
|
32 # Jul 28, 2006: - Version 3.3 [HISTORY] |
|
33 |
|
34 # NOTES: |
|
35 # - could be easily extended to multidimensional linear fit but |
|
36 # what's the use? |
|
37 # - -p zeroes param C for bilinear spice method (T = A + B*sig + C*neph) |
|
38 # so that residual field becomes spice anomaly (linearly correlated |
|
39 # with neph) |
|
40 |
|
41 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
42 # |
|
43 # THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading) |
|
44 # |
|
45 # $modelOpts string of allowed options |
|
46 # $modelOptsUsage usage information string for options |
|
47 # $modelMinArgs min # of arguments of model |
|
48 # $modelArgsUsage usage information string for arguments |
|
49 # |
|
50 # The following variables may be set later but not after &modelInit() |
|
51 # |
|
52 # $modelNFit number of params to fit in model |
|
53 # @nameA symbolic names of model parameters |
|
54 # |
|
55 # You should call &antsDescription() for the -ct options here |
|
56 # |
|
57 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
58 |
|
59 $modelOpts = "pl:"; |
|
60 $modelOptsUsage = "[-p)artial f/r] [-l)imit r_BC <min,max>]"; |
|
61 $modelMinArgs = 1; |
|
62 $modelArgsUsage = "<2nd x-field>"; |
|
63 $modelNFit = 3; |
|
64 $nameA[1] = "A"; $nameA[2] = "B"; $nameA[3] = "C"; |
|
65 $A[1] = nan; $A[2] = nan; $A[3] = nan; |
|
66 #&antsDescription("c","bilin_$nameA[1]", |
|
67 # "c","bilin_$nameA[2]", |
|
68 # "c","bilin_$nameA[3]"); |
|
69 #&antsDescription("t","bilin_sigma_$nameA[1]", |
|
70 # "t","bilin_sigma_$nameA[2]", |
|
71 # "t","bilin_sigma_$nameA[3]", |
|
72 # "t","bilin_ccc_$nameA[1]_$nameA[2]", |
|
73 # "t","bilin_ccc_$nameA[1]_$nameA[3]", |
|
74 # "t","bilin_ccc_$nameA[2]_$nameA[3]"); |
|
75 |
|
76 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
77 # |
|
78 # &modelUsage() mangle parameters; NB: there may be `infinite' # of |
|
79 # filenames after model arguments; this usually sets |
|
80 # @A (the model parameters) but these can later be |
|
81 # calculated heuristically during &modelInit() |
|
82 # |
|
83 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
84 |
|
85 sub modelUsage() |
|
86 { |
|
87 if (defined($opt_l)) { |
|
88 ($minBC,$maxBC) = split(',',$opt_l); |
|
89 &antsUsageError("\n>>> error with -l") |
|
90 unless (defined($maxBC) && $maxBC > $minBC); |
|
91 } |
|
92 $x2fnr = &antsFieldArg(); |
|
93 &antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]); |
|
94 } |
|
95 |
|
96 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
97 # |
|
98 # &modelInit() initializes model after reading of data |
|
99 # |
|
100 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
101 |
|
102 sub modelInit() {} |
|
103 |
|
104 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
105 # |
|
106 # &modelEvaluate(idx,xfnr,vals) evaluate basis funs at x (NB: x1, x2) |
|
107 # idx current index in @ants_ |
|
108 # xfnr field number of x field |
|
109 # vals reference to return values (1-relative!) |
|
110 # |
|
111 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
112 |
|
113 sub modelEvaluate($$$) |
|
114 { |
|
115 my($idx,$xfnr,$valsR) = @_; |
|
116 my($i); |
|
117 |
|
118 $valsR->[1] = 1; |
|
119 $valsR->[2] = $ants_[$idx][$xfnr]; |
|
120 $valsR->[3] = $ants_[$idx][$x2fnr]; |
|
121 } |
|
122 |
|
123 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
124 # &modelCleanup() cleans up after fitting but before output |
|
125 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
126 |
|
127 sub modelCleanup() |
|
128 { |
|
129 if (defined($opt_l)) { |
|
130 my($ccc) = $covar[2][3]/(sqrt($covar[2][2])*sqrt($covar[3][3])); |
|
131 if ($ccc < $minBC || $ccc > $maxBC) { |
|
132 &antsInfo("CCC B/C = %.3g out of range, fit discarded",$ccc); |
|
133 $suppressFit = 1; |
|
134 } |
|
135 } |
|
136 if ($opt_p) { |
|
137 &antsInfo("parameter $nameA[3] = %.3g discarded",$A[3]); |
|
138 $A[3] = 'discarded'; |
|
139 $covar[3][3] = $RMS = $sig = nan; |
|
140 } |
|
141 } |