|
1 #====================================================================== |
|
2 # . L M F I T . G A U S S |
|
3 # doc: Wed Feb 24 09:40:06 1999 |
|
4 # dlm: Fri Jul 28 13:32:35 2006 |
|
5 # (c) 1999 A.M. Thurnherr |
|
6 # uE-Info: 35 51 NIL 0 0 72 2 2 4 NIL ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # What you need to provide if you wanna fit a different |
|
10 # 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 model function which is to be fitted using |
|
14 # a number of pararams which are all stored in @A (beginning at |
|
15 # A[1]!!!). You also need to return the partial derivatives of |
|
16 # the model function wrt all params. |
|
17 # - the interface is documented between +++++++ lines |
|
18 |
|
19 # Gauss data model (i.e. fit Gaussian curve) |
|
20 # NB: - fitting is rather sensitive to the input parameters, thus |
|
21 # a heuristic has been added to guess them (by setting them |
|
22 # to NaN) |
|
23 # - another fickle parameter is the y-offset (zero line); thus |
|
24 # a heuristics has been added for this one as well |
|
25 # - the parameters are peak, mean, standard deviation |
|
26 |
|
27 # HISTORY: |
|
28 # Feb 24, 1999: - created together with [./cfit] |
|
29 # Feb 25, 1999: - cosmetic changes |
|
30 # Jul 31, 1999: - parameter typecheck |
|
31 # Oct 04, 1999: - changed param names |
|
32 # Oct 05, 1999: - improved heuristics |
|
33 # - changed e-scale to sigma |
|
34 # Mar 17, 2001: - param->arg |
|
35 # Jul 28, 2006: - Version 3.3 [HISTORY]; &isnan() |
|
36 |
|
37 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
38 # |
|
39 # THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading) |
|
40 # |
|
41 # $modelOpts string of allowed options |
|
42 # $modelOptsUsage usage information string for options |
|
43 # $modelMinArgs min # of arguments of model |
|
44 # $modelArgsUsage usage information string for arguments |
|
45 # |
|
46 # The following variables may be set later but not after &modelInit() |
|
47 # |
|
48 # $modelNFit number of params to fit in model |
|
49 # @nameA symbolic names of model parameters |
|
50 # |
|
51 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
52 |
|
53 $modelOpts = "y"; |
|
54 $modelOptsUsage = "[-y)shift]"; |
|
55 $modelMinArgs = 0; |
|
56 $modelArgsUsage = "[peak guess [mean guess [sigma guess]]]"; |
|
57 $modelNFit = 3; |
|
58 $nameA[1] = "peak"; |
|
59 $nameA[2] = "mean"; |
|
60 $nameA[3] = "sigma"; |
|
61 |
|
62 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
63 # |
|
64 # &modelUsage() mangle parameters; NB: there may be `infinite' # of |
|
65 # filenames after model arguments; this usually sets |
|
66 # @A (the model parameters) but these can later be |
|
67 # calculated heuristically during &modelInit() |
|
68 # |
|
69 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
70 |
|
71 sub modelUsage() |
|
72 { |
|
73 $A[1] = nan; $A[2] = nan; $A[3] = nan; # usage |
|
74 $A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]); |
|
75 $A[2] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]); |
|
76 $A[3] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]); |
|
77 &antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]); |
|
78 } |
|
79 |
|
80 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
81 # |
|
82 # &modelInit() initializes model after reading of data |
|
83 # |
|
84 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
85 |
|
86 sub modelInit() |
|
87 { |
|
88 my($i,$j,$ymin,$ymax,$xatymax); |
|
89 |
|
90 # -------------------------------------------------- |
|
91 # heuristics for initial model param values |
|
92 # -------------------------------------------------- |
|
93 |
|
94 $ymin = 1e33, $ymax = -1e33, $xatymax = 0; |
|
95 for ($i=0; $i<=$#ants_; $i++) { |
|
96 next if ($antsFlagged[$i]); |
|
97 $ymin = $ants_[$i][$yfnr] |
|
98 if ($ants_[$i][$yfnr] < $ymin); |
|
99 $ymax = $ants_[$i][$yfnr], $xatymax = $ants_[$i][$xfnr] |
|
100 if ($ants_[$i][$yfnr] > $ymax); |
|
101 } |
|
102 $A[1] = $ymax - $ymin if isnan($A[1]); # peak guess |
|
103 $A[2] = $xatymax if isnan($A[2]); # mean guess |
|
104 if (isnan($A[3])) { # sigma guess |
|
105 for ($i=1; |
|
106 $i<=$#ants_ && !$antsFlagged[$i] |
|
107 && $ants_[$i][$yfnr]-$ymin<0.36*$A[1]; |
|
108 $i++) {} |
|
109 for ($j=$#ants_; |
|
110 $j>=1 && !$antsFlagged[$i] |
|
111 && $ants_[$j][$yfnr]-$ymin < 0.36*$A[1]; |
|
112 $j--) {} |
|
113 $A[3] = abs($ants_[$i][$xfnr]-$ants_[$j][$xfnr]) / 2.0; |
|
114 $A[3] *= 0.71; # scale by 1/sqrt(2) |
|
115 if ($A[3] == 0) { |
|
116 &antsInfo("$model: sigma heuristic failed (set to 1)!"); |
|
117 $A[3] = 1; |
|
118 } |
|
119 } |
|
120 |
|
121 # -------------------------------------------------- |
|
122 # y shift (-y option) |
|
123 # -------------------------------------------------- |
|
124 |
|
125 if ($opt_y) { |
|
126 for ($i=1; $i<=$#ants_; $i++) { |
|
127 $ants_[$i][$yfnr] -= $ymin; |
|
128 } |
|
129 } |
|
130 |
|
131 } |
|
132 |
|
133 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
134 # |
|
135 # &modelEvaluate(x,A,dyda) evaluate sum of Gaussians (p.528) at x |
|
136 # x x value (NOT xfnr) |
|
137 # A reference to @A |
|
138 # dyda reference to array for partial derivatives |
|
139 # (wrt individaul params in @A) |
|
140 # <ret val> y value |
|
141 # |
|
142 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
143 |
|
144 sub modelEvaluate($$$) |
|
145 { |
|
146 my($x,$AR,$dydaR) = @_; |
|
147 my($i,$fac,$ex,$arg,$sqrt2sig); |
|
148 my($y) = 0; |
|
149 |
|
150 for ($i=1; $i < $#{$AR}; $i+=3) { |
|
151 $sqrt2sig = (1.4142135623731*$AR->[$i+2]); |
|
152 $arg = ($x - $AR->[$i+1]) / $sqrt2sig; |
|
153 $ex = exp(-$arg*$arg); |
|
154 $fac = $AR->[$i] * $ex * 2*$arg; |
|
155 $y += $AR->[$i] * $ex; |
|
156 |
|
157 $dydaR->[$i] = $ex; |
|
158 $dydaR->[$i+1] = $fac / $sqrt2sig; |
|
159 $dydaR->[$i+2] = $fac * $arg / $sqrt2sig; |
|
160 } |
|
161 return $y; |
|
162 } |
|
163 |
|
164 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
165 # &modelCleanup() cleans up after fitting but before output |
|
166 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
167 |
|
168 sub modelCleanup() |
|
169 { |
|
170 if ($opt_y) { |
|
171 $A[1] += $ymin; |
|
172 for ($i=1; $i<=$#ants_; $i++) { |
|
173 $ants_[$i][$yfnr] += $ymin; |
|
174 } |
|
175 } |
|
176 } |