TBCI Numerical high perf. C++ Library 2.8.0
my_nr.h
Go to the documentation of this file.
1
5/*
6 * last update 29 dec 98
7 * all matrix and vector index ranges are from 0..n-1 !
8 */
9
10/* Changed 98/12/29 by AMB */
11/* 98/01/16 by KG */
12/* $Id: my_nr.h,v 1.10.2.14 2022/11/03 17:28:11 garloff Exp $ */
13
14#ifndef TBCI_MY_NR_H
15#define TBCI_MY_NR_H
16
17#include "tbci/vector.h"
18#include "tbci/matrix.h"
19#include "tbci/constants.h"
20#include "tbci/solver/gauss_jordan.h"
21
22#ifdef HAVE_SIGNALS
23# include <signal.h>
24#endif
25/*
26#include <cstdlib>
27#include <cmath>
28#include <iostream>
29#include <iomanip>
30*/
31
32
33NAMESPACE_TBCI
34
35#ifdef PRAGMA_I
36# pragma interface
37#endif
38
39INST (template <typename T> class dummy friend void do_nothing(T);)
40template <typename T>
41inline void do_nothing(T dummy) {}
42
43// the kronecker delta
44#define DELTA(i,j) ((i==j)?1:0)
45
46/*
47inline int delta (const int i, const int j)
48{
49 return ((i==j) ? 1 : 0);
50}
51 */
52
53#if defined(SIGHANDLE) && defined(HAVE_SIGNALS)
54volatile int sigint = 0;
55
56void breakhandler (int signum)
57{
58 STD__ cerr << "\n SIGINT !\n";
59 sigint++;
60 if (sigint > 2)
61 {
62 STD__ cerr << "\n ... aborted due to too many ignored SIGINTs !\n";
63 exit (1);
64 }
65 else signal (signum, breakhandler);
66}
67#endif
68
69#define DELTAX 1e-30
70
71INST(template <typename T> class Vector friend int basis_trafo (long int, long int, Vector<T>&);)
72template <typename T>
73int basis_trafo (long int value, long int basis, Vector<T>& index)
74/* transforms the integer 'value' into the basis 'basis';
75 index contain the resulting digits with lowest significant first */
76{
77 int n=index.size();
78 int i=0;
79
80 index=(T)0;
81 if (value==0) return 1; //special zero treatment
82 if (n<(int)ceil(MATH__ sqrt(2*log((double)value)/log((double)basis))) || basis<=1) return 0;
83 while( value>0 )
84 {
85 index(i)=value%basis;
86 value/=basis;
87 i++;
88 }
89 return 1;
90}
91
92
93#define CON2 2.0
94#define CON SQRT2
95#define BIG 1.0e30
96#define NTAB 10
97#define SAFE 2.0
98
99INST(template <typename T> class Vector friend T partial_del (T (*func)(const Vector<T>&, const Vector<T>&), const Vector<T>&, const Vector<T>&, int, T, T&);)
100template <typename T>
101T partial_del (T (*func)(const Vector<T>&, const Vector<T>&),
102 const Vector<T>& x, const Vector<T>& p,
103 int mu, T h, T& err)
104 // calculate partial derivative of func with respect to mu
105 // Acc. to Num. recipes, ch. 5.7
106{
107 //static T maxarg1 ALIGN(MIN_ALIGN), maxarg2 ALIGN(MIN_ALIGN);
108 T errt ALIGN(MIN_ALIGN), fac, hh, ans(0);
109 int xlen = x.size();
110 Matrix<T> a(NTAB,NTAB);
111 Vector<T> buf_plus(xlen);
112 Vector<T> buf_minus(xlen);
113#ifdef FOURTHORDER
114 Vector<T> buf_plus2(xlen);
115 Vector<T> buf_minus2(xlen);
116#endif
117
118 if (mu >= xlen)
119 STD__ cerr << "\n error: partial_del: derivate variable index out of range \n";
120 hh = h ;
121 buf_plus = x; buf_minus = x; buf_plus(mu) += hh; buf_minus(mu) -= hh;
122#ifdef FOURTHORDER
123 buf_plus2 = x; buf_minus2 = x; buf_plus2(mu) += hh*2; buf_minus2 -= hh*2;
124 a(0,0) = (-(*func)(buf_plus2, p) + 8.0*(*func)(buf_plus, p)
125 - 8.0*(*func)(buf_minus, p) + (*func)(buf_minus2, p)) / (hh*12.0);
126#else
127 a(0,0) = ((*func)(buf_plus, p) - (*func)(buf_minus, p)) / (2.0*hh);
128#endif
129 err = BIG;
130 for (int i = 1; i < NTAB; i++)
131 {
132 hh /= CON;
133 buf_plus = x; buf_minus = x; buf_plus(mu) += hh; buf_minus(mu) -= hh;
134#ifdef FOURTHORDER
135 buf_plus2 = x; buf_minus2 = x; buf_plus2(mu) += hh*2; buf_minus2 -= hh*2;
136 a(0,i) = (-(*func)(buf_plus2, p) + 8.0*(*func)(buf_plus, p)
137 - 8.0*(*func)(buf_minus, p) + (*func)(buf_minus2, p)) / (hh*12.0);
138#else
139 a(0,i) = ((*func)(buf_plus, p) - (*func)(buf_minus, p)) / (2.0*hh);
140#endif
141 fac = CON2;
142 for (REGISTER int j = 1; j <= i; j++)
143 {
144 a(j,i) = (a(j-1,i)*fac - a(j-1,i-1))/(fac-1.0);
145 fac *= CON2;
146 errt = MAX (MATH__ fabs(a(j,i)-a(j-1,i)), MATH__ fabs(a(j,i)-a(j-1,i-1)));
147 if (errt <= err) {
148 err = errt;
149 ans = a(j,i);
150 }
151 }
152 if (MATH__ fabs (a(i,i)-a(i-1,i-1)) >= MATH__ fabs (SAFE*err)) break;
153 }
154 return ans;
155}
156
157INST(template <typename T> class Vector friend T chisq (Vector<T>&, Vector<T>&, Vector<T>&);)
158template <typename T>
159T chisq (Vector<T>& data, Vector<T>& func, Vector<T>& sig)
160 // calculates the chi-squared sum
161
162{
163 T chi ALIGN(MIN_ALIGN) = 0;
164 int len=data.size();
165
166 for(REGISTER int i=0;i<len;i++) chi+=sqr((data(i)-func(i))/sig(i));
167 return(chi);
168}
169
170INST(template <typename T> class Vector friend T chisquare (T (*func)(const Vector<T>&, const Vector<T>&), Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&);)
171template <typename T>
172T chisquare (T (*func)(const Vector<T>&, const Vector<T>&),
173 Vector<T>& param, Vector<T>& x,
174 Vector<T>& y, Vector<T>& sigma)
175 // calculates function values of func at position x with params param
176 // and return chisq with given y values weighted with sigmas
177{
178 Vector<T> arg(1);
179
180 T chi ALIGN(MIN_ALIGN) = 0, wfvaldiff;
181 int anzpkt = x.size();
182
183 for (REGISTER int i = 0; i < anzpkt; i++)
184 {
185 arg(0) = x(i);
186 wfvaldiff = (y(i) - (*func)(param, arg)) / sigma(i);
187 chi += sqr (wfvaldiff);
188 }
189 return chi;
190}
191
192#if 0
193//class nslist<minimum<T>,<T>> minima;
194
195template <typename T>
196class minimum
197{
198 public:
199 Vector<T> koord;
200 T value ALIGN(MIN_ALIGN);
201 T getval () { return value; }
202 minimum () { value = 1e38; }
203};
204
205INST(template <typename T> class Vector friend TVector<T> grid_mins(T (*func)(const Vector<T>&, const Vector<T>&), Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, int, char, int, nsList<minimum<T>,T>& );)
206template <typename T>
207TVector<T> grid_mins(T (*func)(const Vector<T>&, const Vector<T>&),
208 Vector<T>& x, Vector<T>& y,
209 Vector<T>& sigma, Vector<T>& pmin, Vector<T>& pmax,
210 int res, char v, int anzmin,
211 nsList<minimum<T>,T>& minima)
212 //return absmin;
213 // find minima of func on grid with resolution res
214{
215 int anzpar = pmin.size();
216 minimum<T>* min;
217 minimum<T> min2; min2.koord.resize (anzpar);
218 Vector<T> actpar(anzpar);
219 TVector<T> absmin(anzpar);
220 //absmin.resize(anzpar);
221 Vector<T> factor(anzpar);
222 Vector<T> plim(anzpar);
223 long int len = (long int) MATH__ pow ((double)res, anzpar);
224 bool overflow; int k;
225
226 actpar = pmin;
227 T actval ALIGN(MIN_ALIGN) = chisquare (func, actpar, x, y, sigma);
228 min = minima.setlast();
229 STD__ cerr << "Parameter space discretisation in process ...\n";
230 //factors
231 for (REGISTER int j = 0; j < anzpar; j++)
232 {
233 factor(j) = (pmax(j) - pmin(j)) / (res - 1);
234 plim(j) = pmax(j) + 0.5 * factor(j);
235 }
236 //discretisation
237 for (long int i = 0; i < len; i++)
238 {
239 actval = chisquare (func, actpar, x, y, sigma);
240 if (actval < min->value)
241 {
242 if (v=='y' || v=='Y')
243 STD__ cerr << "Found new min " << actval << " at " << actpar << '\n';
244#ifdef USEQSORT
245 min->koord = actpar;
246 min->value = actval;
247 if (anzmin > 1) minima.qsort ();
248 min = minima.setlast ();
249#else
250 min2.koord = actpar; min2.value = actval;
251 min = minima.sortin_r ( min2 );
252#endif
253 }
254 if (!(i % 1000)) STD__ cerr << i << '\r';
255 //avoid basis_trafo
256 overflow = true; k = 0;
257 while (overflow && (k < anzpar))
258 {
259 actpar(k) += factor(k);
260 if (actpar(k) > plim(k)) actpar(k) = pmin(k);
261 else overflow = false;
262 k++;
263 }
264#ifdef SIGHANDLE
265 if (sigint) { sigint--; return 0; }
266#endif
267 }
268 //STD__ cerr << "\nLast par: " << actpar << "\n";
269 STD__ cerr << "\n ... done.\n\n";
270 absmin = minima.getfirst()->koord;
271 return absmin;
272}
273#endif /* 0 */
274
275INST(template <typename T> class Vector friend TVector<T> grid_min(T (*func)(const Vector<T>&, const Vector<T>&), Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, long int, int, char);)
276template <typename T>
277TVector<T> grid_min(T (*func)(const Vector<T>&, const Vector<T>&),
278 Vector<T>& x, Vector<T>& y,
279 Vector<T>& sig, Vector<T>& pmin,
280 Vector<T>& pmax, Vector<T>& lambda,
281 long int res, int max_it, char v)
282 // find minimun of func on grid with resolution res
283{
284
285 int plen=pmin.size();
286 int xlen=x.size();
287 long int len=0;
288 Vector<T> p(pmin);
289 TVector<T> min_pos(pmin);
290 Vector<T> factor(plen);
291 Vector<T> arg(1);
292 Vector<T> fval(x.size());
293 Vector<T> index(plen);
294 Vector<T> b(plen,0);
295 Vector<T> a(plen,0);
296 T min ALIGN(MIN_ALIGN) = 0, chi=0, err=0.0, fac=1e-8;
297
298 len=(long int) pow((double)res, plen);
299 for (int s=0; s<xlen; s++) { arg(0)=x(s); fval(s)=(*func)(p, arg); }
300 min=chisq(y, fval, sig);
301 STD__ cerr << "\n parameter space discretisation in progress ... \n ";
302 for (REGISTER int j = 0; j < plen; j++)
303 factor(j) = (pmax(j) - pmin(j)) / ((T)res - (T)1);
304 for (long int i=0;i<len;i++)
305 {
306 // linearize recursive for structure by basis trafo
307 if ( !basis_trafo(i, res, index) )
308 STD__ cerr << "\n error: basis trafo failed\n";
309 for (REGISTER int j=0;j<plen;j++)
310 p(j) = factor(j) * index(j) + pmin(j); //p-space discr.
311 //do max_it steepest gradient steps
312 for (int l=0;l<max_it;l++)
313 {
314 // setup Vector
315 err=0.0;
316 for (int k=0;k<plen;k++)
317 {
318 for (int j=0;j<xlen;j++)
319 {
320 arg(0)=x(j);
321 b(k)+=(y(j)-(*func)(p, arg))/(sig(j)*sig(j))*
322 partial_del(func, p, arg, k, x(j)*fac, err);
323 }
324 }
325 // setup diagonal Matrix for steepest gradient const
326 err = 0.0;
327 for (int t=0; t<plen; t++)
328 {
329 for (REGISTER int u=0; u<xlen; u++)
330 {
331 arg(0) = x(u);
332 a(t) += sqr(partial_del(func, p, arg, t, x(u)*fac, err))/
333 (sig(u)*sig(u));
334 }
335 }
336 for (int q=0; q<plen; q++)
337 p(q) += b(q)*lambda(q)/(a(q)+DELTAX);
338 }
339 // calc func values with new parameter
340 for (int q=0; q<xlen; q++) { arg(0)=x(q); fval(q)=(*func)(p, arg); }
341 // calc chi2 sqared
342 chi=chisq(y, fval, sig);
343 if ( chi<min ) { min=chi; min_pos=p; } //choose new min
344 if ( i%1000==0 ) STD__ cerr << i << "\r";
345#ifdef SIGHANDLE
346 if (sigint) { sigint--; return 0; }
347#endif
348 }
349 STD__ cerr << " \n done! \n\n\n";
350 return min_pos;
351}
352
353INST(template <typename T> class Vector friend T lev_mar(Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, T (*func)(const Vector<T>&, const Vector<T>&), Vector<T>, Vector<T>&, T, int, char);)
354template <typename T>
355T lev_mar(Vector<T>& x, Vector<T>& y, Vector<T>& sig, Vector<T>& param,
356 Vector<T>& plow, Vector<T>& phigh,
357 T (*func)(const Vector<T>&, const Vector<T>&),
358 Vector<T> lambda, Vector<T>& lscale,
359 T tol, int iter, char ver)
360/* fit method: levenberg-marquardt
361 leads to an linear equ
362 derivatives are calculated numerically by num_derivatives
363 parameter space is restricted by plow and phigh */
364
365{
366 T chi ALIGN(MIN_ALIGN), chi_new, chi_min, chi_dif, err=0.0, fac=1.0e-8;
367 int len=x.size();
368 int plen=param.size();
369 Vector<T> param2(plen); // Vector<T> for new parameter set
370 Vector<T> param_best(plen); // Vector<T> for best parameter set
371 Vector<T> fval(len); // Vector for fit-func values
372 Vector<T> arg(1); // Vector for arguments of fit-func (now only 1d)
373
374 // comp chi-squared
375 chi=chisquare(func,param,x,y,sig);
376 STD__ cerr << " \n chi-squared with start parameters is: " << chi;
377 chi_min=chi;
378 // ini matrices
379 Matrix<T> a(plen,plen);
380 Matrix<T> b(plen,1);
381 STD__ cerr << "\n initializing matrices and starting loop... \n";
382 int j=0; // iteration counter
383 param_best=param;
384 do
385 {
386 // setup Matrix
387 for (int k=0;k<plen;k++)
388 for (int l=0;l<plen;l++)
389 {
390 for (int i=0;i<len;i++)
391 {
392 arg(0)=x(i);
393 a(k,l)+=partial_del(func, param, (Vector<T>&)arg, k, x(i)*fac, err)*
394 partial_del(func, param, (Vector<T>&)arg, l, x(i)*fac, err) /(sig(i)*sig(i));
395 }
396 a(k,l)*=((T)1+lambda(k)*(T)DELTA(k,l));
397 }
398 // setup Vector<T>
399 for (int q=0; q<plen; q++)
400 {
401 for (int i=0;i<len;i++)
402 {
403 arg(0)=x(i);
404 b(q,0)+=(y(i)-(*func)(param, (Vector<T>&)arg))/(sig(i)*sig(i))*
405 partial_del(func, param, (Vector<T>&)arg, q, x(i)*fac, err);
406 }
407 }
408 // solve equ a.x=b
409 char exceed = gaussj (a, b);
410 // calc new parameter
411 for (REGISTER int t=0; t<plen; t++) param2(t)=param(t)+b(t,0);
412 // check parameter range
413 exceed = 0;
414 for (REGISTER int i=0; i<plen; i++)
415 {
416 if ((param2(i) > phigh(i)) && (phigh(i)!=plow(i)) ) exceed = 1;
417 if ((param2(i) < plow(i)) && (phigh(i)!=plow(i)) ) exceed = 1;
418 }
419 // calc chi2 sqared
420 chi_new = chisquare(func, param2, x, y, sig);
421 chi_dif = MATH__ fabs (chi_new - chi);
422 if ((chi_new >= chi) || (exceed == 1)) lambda *= lscale(0);
423 //increase lambda cause downhill is exhausted or limits are exceeded
424 else {
425 // decrease lambda and update trial solution
426 lambda*=(T)1/lscale(1);
427 param=param2;
428 // update best est. param if new min is reached
429 if (chi_new<chi_min) {
430 chi_min=chi_new;
431 param_best=param;
432 }
433 }
434 chi=chi_new; // update chi
435 if (ver=='y' || ver=='Y')
436 STD__ cerr << "\n iteration " << j << " " << param <<
437 " " << param2 << " " << chi_new;
438#ifdef SIGHANDLE
439 if (sigint) { sigint--; param = param_best; return 1e38; }
440#endif
441 j++;
442 if (j%10 == 0) STD__ cerr << ".";
443 // end loop if tolerance is reached or number of iterations exceeded
444 } while( (chi_dif >= tol) && (j < iter) );
445 param=param_best; // return best parameter set found
446 STD__ cerr << "\n loop finished after " << j << " iterations with last delta chi: "
447 << chi_dif << "\n and chi-squared: " << chi_new << "\n";
448 return (chi_new);
449}
450
451INST(template <typename T> class Matrix friend void hpsort (Matrix<T>&, int);)
452template <typename T>
453void hpsort (Matrix<T>& x, int pos)
454/* sorts rows of input Matrix x into ascending numerical order of values
455 in colomn pos;
456 uses heap sort algorithm, rf. nr p 336 */
457
458{
459 unsigned long i, ir, j, l, n;
460 Vector<T> rra(x.columns());
461
462 n=x.rows();
463 if ( n<2 ) return;
464 l=(n >> 1)+1;
465 ir=n;
466 for (;;)
467 {
468 if ( l>1 ) {
469 rra=x(--l-1);
470 }
471 else {
472 rra=x(ir-1);
473 x(ir-1)=x(0);
474 if (--ir==1) {
475 x(0)=rra;
476 break;
477 }
478 }
479 i=l;
480 j=l+l;
481 while ( j<=ir ) {
482 if ( j<ir && x(j-1,pos) < x(j,pos) ) j++;
483 if ( rra(pos)<x(j-1,pos) ) {
484 x(i-1)=x(j-1);
485 i=j;
486 j <<= 1;
487 }
488 else j=ir+1;
489 }
490 x(i-1)=rra;
491 }
492}
493
494NAMESPACE_END
495
496#endif /* TBCI_MY_NR_H */
497
NTAB
#define NTAB
Definition LM_fit.h:89
mu
const Vector< T > const Vector< T > const Vector< T > int mu
Definition LM_fit.h:99
SAFE
#define SAFE
Definition LM_fit.h:90
BIG
#define BIG
Definition LM_fit.h:88
basis
long int basis
Definition LM_fit.h:66
p
const Vector< T > const Vector< T > const Vector< T > & p
Definition LM_fit.h:98
index
long int Vector< T > & index
Definition LM_fit.h:69
DELTAX
#define DELTAX
Definition LM_fit.h:49
xlen
int xlen
Definition LM_fit.h:105
x
const Vector< T > const Vector< T > & x
Definition LM_fit.h:97
CON
#define CON
Definition LM_fit.h:87
buf_minus
Vector< T > buf_minus(xlen)
h
const Vector< T > const Vector< T > const Vector< T > int T h
Definition LM_fit.h:99
err
const Vector< T > const Vector< T > const Vector< T > int T T & err
Definition LM_fit.h:102
buf_plus
Vector< T > buf_plus(xlen)
i
int i
Definition LM_fit.h:71
CON2
#define CON2
Definition LM_fit.h:86
y
doublereal y
Definition TOMS_707.C:27
STD__
#define STD__
Definition basics.h:338
NAMESPACE_END
#define NAMESPACE_END
Definition basics.h:323
INST
#define INST(x)
Definition basics.h:238
NAMESPACE_TBCI
#define NAMESPACE_TBCI
Definition basics.h:317
REGISTER
#define REGISTER
Definition basics.h:108
MATH__
#define MATH__
Definition basics.h:339
ALIGN
#define ALIGN(x)
Definition basics.h:444
MIN_ALIGN
#define MIN_ALIGN
Definition basics.h:421
MAX
#define MAX(a, b)
Definition basics.h:656
T
#define T
Definition bdmatlib.cc:20
TVector
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition vector.h:73
TVector::size
unsigned long size() const
Definition vector.h:104
nsList
A numerically sortable List.
Definition list.h:771
fabs
double fabs(const TBCI__ cplx< T > &c)
Definition cplx.h:746
arg
T arg(const TBCI__ cplx< T > &c)
Definition cplx.h:690
pow
TBCI__ cplx< T > pow(const TBCI__ cplx< T > &z, const double n)
Definition cplx.h:761
sqr
cplx< T > sqr(const cplx< T > &c)
Definition cplx.h:449
log
TBCI__ cplx< T > log(const TBCI__ cplx< T > &z)
Definition cplx.h:771
sqrt
NAMESPACE_END NAMESPACE_CSTD TBCI__ cplx< T > sqrt(const TBCI__ cplx< T > &z)
Definition cplx.h:751
min
#define min(a, b)
Definition f2c.h:180
fac
tm fac
Definition f_matrix.h:1052
b
F_TMatrix< T > b
Definition f_matrix.h:736
gaussj
NAMESPACE_TBCI char gaussj(Matrix< T > &a, Matrix< T > &b)
Definition gauss_jordan.h:46
a
const unsigned TMatrix< T > const Matrix< T > * a
Definition matrix_kernels.h:26
res
const unsigned TMatrix< T > * res
Definition matrix_kernels.h:26
do_nothing
NAMESPACE_TBCI void do_nothing(T dummy)
Definition my_nr.h:41
p
const Vector< T > const Vector< T > const Vector< T > & p
Definition my_nr.h:102
x
const Vector< T > const Vector< T > & x
Definition my_nr.h:101
buf_minus
Vector< T > buf_minus(xlen)
buf_plus
Vector< T > buf_plus(xlen)
DELTA
#define DELTA(i, j)
Definition my_nr.h:44
basis_trafo
int basis_trafo(long int value, long int basis, Vector< T > &index)
Definition my_nr.h:73
Generated by doxygen 1.16.1