Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
software/scripts/gradient.sh
Go to the documentation of this file.
1#!/usr/bin/env zsh
2script=${0##*/}
3
4zmodload zsh/mathfunc
5
6if [[ ! $ZSH_EVAL_CONTEXT =~ :file$ ]]; then
7 echo "Auxiliary script for running the conjugate gradient fit in zsh."
8 echo "This script should be sourced by the steering script."
9 exit
10fi
11
12# This script contains auxiliary functions for running the conjugate gradient fit.
13# It should be sourced by the steering script.
14# The steering script should provide for the method getChi2 and set the fit parameters before calling method gradient.
15# On return, the method getChi2 should set the current value of chi2 to its first positional parameter.
16# The second positional parameter corresponds to an option that can be used within method getChi2.
17# The value of the option corresponds to the following cases.
18# 0 => step wise improvement of the chi2;
19# 1 => evaluation of the chi2 before the determination of the gradient of the chi2; and
20# 2 => evaluation of the derivative of the chi2 to each fit parameter.
21# The fit parameters are defined by the associative array PARAMETERS.
22# In this, each key corresponds to a valid zsh command that can be used for changing a specific fit parameter.
23# The command should contain the wild card '%' which is repeatedly replaced by a value, referred to as step size.
24# The associated value is used to evaluate the derivative of chi2 per unit step size of each fit parameter.
25# These step sizes should be tuned beforehand to yield similar absolute values of the derivatives.
26
27typeset -a CHI2 # internal chi2 values; index [1] corresponds to best value
28typeset -a GRADIENT # derivative of chi2 per unit step size of each fit parameter
29let "EPSILON = 5.0E-4" # default maximal distance to minimum
30let "NUMBER_OF_ITERATIONS = 1000" # default maximum number of iterations
31typeset -A PARAMETERS # fit parameters
32
33if [[ -z "$DEBUG" ]]; then
34 DEBUG=0
35fi
36
37
38#
39# Auxiliary function to evaluate chi2 and gradient thereof at current position.
40#
41# The external method getChi2 is used to evaluate the chi2.
42# The value corresponding to each key in the associative array PARAMETERS
43# is used as unit step size for the evaluation of the derivative of the chi2.
44# The results are stored in the internal variables CHI2[1] and GRADIENT, respectively.
45#
46function evaluate()
47{
48 GRADIENT=()
49
50 getChi2 CHI2\[1\] 1
51
52 let "WIDTH = 1"
53
54 for KEY in ${(k)PARAMETERS}; do
55 if (( ${#KEY} > $WIDTH )); then
56 let "WIDTH = ${#KEY}"
57 fi
58 done
59
60 let "V = 0.0"
61 let "i = 1"
62
63 for KEY in ${(ko)PARAMETERS}; do
64
65 BUFFER=(${(s/;/)${(Q)KEY}})
66
67 VALUE=$PARAMETERS[$KEY]
68
69 if (( $VALUE != 0.0 )); then
70
71 for K1 in ${(o)BUFFER}; do
72 eval ${(Q)${(qq)K1//\%/$((+0.5 * $VALUE))}}
73 done
74
75 getChi2 CHI2\[2\] 2
76
77 for K1 in ${(O)BUFFER}; do
78 eval ${(Q)${(qq)K1//\%/$((-0.5 * $VALUE))}}
79 done
80
81 for K1 in ${(o)BUFFER}; do
82 eval ${(Q)${(qq)K1//\%/$((-0.5 * $VALUE))}}
83 done
84
85 getChi2 CHI2\[3\] 2
86
87 for K1 in ${(O)BUFFER}; do
88 eval ${(Q)${(qq)K1//\%/$((+0.5 * $VALUE))}}
89 done
90
91 GRADIENT[$i]=$(($CHI2[2] - $CHI2[3]))
92 else
93
94 GRADIENT[$i]=0.0
95 fi
96
97 if (( $DEBUG >= 3 )); then
98 printf "%-${WIDTH}s %12.5f %12.5f\n" $KEY $PARAMETERS[$KEY] $GRADIENT[$i]
99 fi
100
101 let "V += $GRADIENT[$i]*$GRADIENT[$i]"
102 let "i += 1"
103 done
104
105 if (( $DEBUG >= 3 )); then
106 printf "%-${WIDTH}s %12s %12.5f\n" "|gradient|" "" $((sqrt($V)))
107 fi
108}
109
110
111#
112# Move current position along gradient by given scaling factor of the unit step sizes.
113#
114# \param 1 scaling factor
115#
116function move()
117{
118 if (( $1 > 0.0 )); then
119
120 let "i = 1"
121
122 for KEY in ${(ko)PARAMETERS}; do
123
124 BUFFER=(${(s/;/)${(Q)KEY}})
125
126 for K1 in ${(o)BUFFER}; do
127 eval ${(Q)${(qq)K1//\%/$(($GRADIENT[$i] * $1 * $PARAMETERS[$KEY]))}}
128 done
129
130 let "i += 1"
131 done
132
133 else
134
135 let "i = ${#GRADIENT}"
136
137 for KEY in ${(kO)PARAMETERS}; do
138
139 BUFFER=(${(s/;/)${(Q)KEY}})
140
141 for K1 in ${(O)BUFFER}; do
142 eval ${(Q)${(qq)K1//\%/$(($GRADIENT[$i] * $1 * $PARAMETERS[$KEY]))}}
143 done
144
145 let "i -= 1"
146 done
147 fi
148}
149
150
151#
152# Main fit function.
153#
154# The algorithm is based on the conjugate gradient method.
155# The position is moved to the optimal value and the final chi2 is returned.
156#
157# If the option is one, the step sizes are fine-tuned according to the gradient of the chi2.
158# The value corresponding to each key in the associative array PARAMETERS is then
159# multiplied by the corresponding element in GRADIENT and only the sign thereof is maintained.
160#
161# \param 1 maximum number of extra steps
162# \param 2 variable containing final chi2 value on return
163#
164function gradient()
165{
166 typeset -a G
167 typeset -a H
168
169 if (( ${#PARAMETERS} == 0 )); then
170 return
171 fi
172
173 evaluate
174
175 # initialise
176
177 for (( i = 1; i <= ${#GRADIENT}; ++i)); do
178 G[$i]=$((-1.0 * $GRADIENT[$i]))
179 H[$i]=$G[$i]
180 GRADIENT[$i]=$H[$i]
181 done
182
183 for (( N = 0; $N != $NUMBER_OF_ITERATIONS; N += 1 )); do
184
185 if (( $DEBUG >= 3 )); then
186 printf "chi2[1] %4d %12.5f\n" $N $CHI2[1]
187 fi
188
189 # minimise chi2 in direction of gradient
190
191 CHI2[2]=$CHI2[1]
192 CHI2[3]=$CHI2[2]
193
194 for (( DS = 1.0, M = 0 ; $DS > 1.0e-3; )); do
195
196 move $DS
197
198 getChi2 CHI2\[4\] 0
199
200 if (( $DEBUG >= 3 )); then
201 printf "chi2[4] %4d %12.5f %12.5e\n" $M $CHI2[4] $DS
202 fi
203
204 if (( $CHI2[4] < $CHI2[3] )); then
205
206 # continue
207
208 CHI2[2]=$CHI2[3]
209 CHI2[3]=$CHI2[4]
210
211 let "M = 0"
212
213 continue
214 fi
215
216 # if chi2 increases, try additional steps first to overcome possible hurdle, then try single step with reduced size
217
218 if (( $DS == 1.0 )); then
219
220 if (( $M == 0 )); then
221 CHI2[5]=$CHI2[4]
222 fi
223
224 if (( $M != $1 )); then
225
226 let "M += 1"
227
228 continue;
229
230 else
231
232 for (( ; $M != 0; --M )); do
233 move $((-1.0 * $DS))
234 done
235
236 CHI2[4]=$CHI2[5]
237 fi
238 fi
239
240 move $((-1.0 * $DS))
241
242 if (( $CHI2[3] < $CHI2[2] )); then
243
244 # final step based on parabolic interpolation through following points
245 #
246 # X1 = -1 * DS -> CHI2[2]
247 # X2 = 0 * DS -> CHI2[3]
248 # X3 = +1 * DS -> CHI2[4]
249
250 let "F21 = $CHI2[3] - $CHI2[2]" # F(X2) - F(X1)
251 let "F23 = $CHI2[3] - $CHI2[4]" # F(X2) - F(X3)
252
253 let "XS = 0.5 * ($F21 - $F23) / ($F23 + $F21)"
254
255 move $((+1.0 * $XS * $DS))
256
257 getChi2 CHI2\[4\] 0
258
259 if (( $CHI2[4] < $CHI2[3] )); then
260
261 CHI2[3]=$CHI2[4]
262
263 else
264
265 move $((-1.0 * $XS * $DS))
266
267 getChi2 CHI2\[3\] 0
268 fi
269
270 if (( $DEBUG >= 3 )); then
271 printf "chi2[3] %4d %12.5f %12.5e\n" $N $CHI2[3] $DS
272 fi
273
274 break
275
276 else
277
278 # reduce step size
279
280 let "DS *= 0.5"
281 fi
282 done
283
284 if (( fabs($CHI2[3] - $CHI2[1]) < $EPSILON * 0.5 * (fabs($CHI2[1]) + fabs($CHI2[3])) )); then
285
286 CHI2[1]=$CHI2[3]
287
288 break;
289 fi
290
291 evaluate
292
293 let "GG = 0.0"
294 let "DGG = 0.0"
295
296 for (( i = 1; $i <= ${#GRADIENT}; ++i )); do
297 let "GG += $G[$i]*$G[$i]"
298 let "DGG += ($GRADIENT[$i] + $G[$i]) * $GRADIENT[$i]"
299 done
300
301 if (( $GG == 0.0 )); then
302 break
303 fi
304
305 let "GAM = $DGG / $GG"
306
307 for (( i = 1; $i <= ${#GRADIENT}; ++i )); do
308 G[$i]=$((-1.0 * $GRADIENT[$i]))
309 H[$i]=$(($G[$i] + $GAM * $H[$i]))
310 GRADIENT[$i]=$H[$i]
311 done
312 done
313
314 if (( $DEBUG >= 3 )); then
315 printf "chi2[1] %4d %12.5f\n" $N $CHI2[1]
316 fi
317
318 eval $2=$CHI2[1]
319}
320
321
322#
323# Additional fit function to be called after gradient.
324#
325# The algorithm is based on sorting the derivatives of the chi2 and
326# moving the position one-by-one to the optimal value.
327# The final chi2 is returned.
328#
329# \param 1 variable containing final chi2 value on return
330#
331function simplex()
332{
333 typeset -A G
334
335
336 let "WIDTH = 1"
337
338 for KEY in ${(k)PARAMETERS}; do
339 if (( ${#KEY} > $WIDTH )); then
340 let "WIDTH = ${#KEY}"
341 fi
342 done
343
344
345 for (( N = 0; $N != $NUMBER_OF_ITERATIONS; N += 1 )); do
346
347 G=()
348
349 GRADIENT=()
350
351 getChi2 CHI2\[1\] 1
352
353 let "i = 1"
354
355 for KEY in ${(ko)PARAMETERS}; do
356
357 VALUE=$PARAMETERS[$KEY]
358
359 BUFFER=(${(s/;/)${(Q)KEY}})
360
361 for K1 in ${(o)BUFFER}; do
362 eval ${(Q)${(qq)K1//\%/$((+1.0 * $VALUE))}}
363 done
364
365 getChi2 CHI2\[2\] 2
366
367 for K1 in ${(O)BUFFER}; do
368 eval ${(Q)${(qq)K1//\%/$((-1.0 * $VALUE))}}
369 done
370
371 for K1 in ${(o)BUFFER}; do
372 eval ${(Q)${(qq)K1//\%/$((-1.0 * $VALUE))}}
373 done
374
375 getChi2 CHI2\[3\] 2
376
377 for K1 in ${(O)BUFFER}; do
378 eval ${(Q)${(qq)K1//\%/$((+1.0 * $VALUE))}}
379 done
380
381 if (( $DEBUG >= 3 )); then
382 printf "%-${WIDTH}s %12.5f %12.5f -> %+9.5f <- %+9.5f\n" $KEY $VALUE $CHI2[1] $(($CHI2[2] - $CHI2[1])) $(($CHI2[3] - $CHI2[1]))
383 fi
384
385 # map index to derivative of chi2 and maintain sign of step
386
387 if (($CHI2[2] < $CHI2[1] && $CHI2[2] < $CHI2[3])); then
388
389 G[$i]=$(($CHI2[2] - $CHI2[1]))
390
391 GRADIENT[$i]=+1.0
392
393 elif (($CHI2[3] < $CHI2[1] && $CHI2[3] < $CHI2[2])); then
394
395 G[$i]=$(($CHI2[3] - $CHI2[1]))
396
397 GRADIENT[$i]=-1.0
398#
399# elif (($CHI2[3] - $CHI2[2] > 0.5*($CHI2[2] + $CHI2[3]) - $CHI2[1])); then
400#
401# G[$i]=$(($CHI2[2] - $CHI2[1]))
402#
403# GRADIENT[$i]=+0.33
404#
405# elif (($CHI2[2] - $CHI2[3] > 0.5*($CHI2[2] + $CHI2[3]) - $CHI2[1])); then
406#
407# G[$i]=$(($CHI2[3] - $CHI2[1]))
408#
409# GRADIENT[$i]=-0.33
410 fi
411
412 let "i += 1"
413 done
414
415
416 if (( ${#G} == 0 )); then
417 break
418 fi
419
420
421 # sort derivatives of chi2 and optimise position one-by-one
422
423 CHI2[2]=$CHI2[1]
424
425 for i VALUE in `printf "%d %12.3e\n" ${(kv)G} | sort -g -k2`; do
426
427 KEY=${${(ko)PARAMETERS}[$i]}
428
429 BUFFER=(${(s/;/)${(Q)KEY}})
430
431 # initial step size
432
433 let "DS = $GRADIENT[$i] * $PARAMETERS[$KEY]"
434
435 if (($DEBUG >= 3)); then
436 printf "%-${WIDTH}s %12.5f %12.5f\n" $KEY $DS $G[$i]
437 fi
438
439 CHI2[3]=$CHI2[2]
440
441 for (( n = 1; ; ++n )); do
442
443 for K1 in ${(o)BUFFER}; do
444 eval ${(Q)${(qq)K1//\%/$((+1.0 * $DS))}}
445 done
446
447 getChi2 CHI2\[4\] 0
448
449 if (( $DEBUG >= 3 )); then
450 printf "chi2[4] %4d %12.5f\n" $n $CHI2[4]
451 fi
452
453 if (( $CHI2[4] < $CHI2[3] )); then
454
455 CHI2[3]=$CHI2[4]
456
457 continue
458 fi
459
460 for K1 in ${(O)BUFFER}; do
461 eval ${(Q)${(qq)K1//\%/$((-1.0 * $DS))}}
462 done
463
464 if (( $CHI2[3] < $CHI2[2] )); then
465
466 # final step based on parabolic interpolation through following points
467 #
468 # X1 = -1 * DS -> CHI2[2]
469 # X2 = 0 * DS -> CHI2[3]
470 # X3 = +1 * DS -> CHI2[4]
471
472 let "F21 = $CHI2[3] - $CHI2[2]" # F(X2) - F(X1)
473 let "F23 = $CHI2[3] - $CHI2[4]" # F(X2) - F(X3)
474
475 let "XS = 0.5 * ($F21 - $F23) / ($F23 + $F21)"
476
477 for K1 in ${(o)BUFFER}; do
478 eval ${(Q)${(qq)K1//\%/$((+1.0 * $XS * $DS))}}
479 done
480
481 getChi2 CHI2\[4\] 0
482
483 if (( $CHI2[4] < $CHI2[3] )); then
484
485 CHI2[3]=$CHI2[4]
486
487 else
488
489 for K1 in ${(O)BUFFER}; do
490 eval ${(Q)${(qq)K1//\%/$((-1.0 * $XS * $DS))}}
491 done
492 fi
493 fi
494
495 if (( $DEBUG >= 3 )); then
496 printf "chi2[3] %4d %12.5f\n" $n $CHI2[3]
497 fi
498
499 break
500 done
501
502 CHI2[2]=$CHI2[3]
503 done
504
505 if (( fabs($CHI2[2] - $CHI2[1]) < $EPSILON * 0.5 * (fabs($CHI2[1]) + fabs($CHI2[2])) )); then
506
507 CHI2[1]=$CHI2[2]
508
509 break;
510 fi
511 done
512
513 if (( $DEBUG >= 3 )); then
514 printf "chi2[1] %4d %12.5f\n" $N $CHI2[1]
515 fi
516
517 eval $1=$CHI2[1]
518}