ODIN
fitting.h
1 /***************************************************************************
2  fitting.h - description
3  -------------------
4  begin : Fri Apr 6 2001
5  copyright : (C) 2000-2021 by Thies Jochimsen & Michael von Mengershausen
6  email : thies@jochimsen.de mengers@cns.mpg.de
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #ifndef FITTING_H
19 #define FITTING_H
20 
21 #include <tjutils/tjnumeric.h> // for MinimizationFunction
22 #include <odindata/data.h>
23 #include <odindata/linalg.h>
24 #include <odindata/utils.h>
25 
26 
27 #define DEFAULT_MAX_ITER 1000
28 #define DEFAULT_TOLERANCE 1e-4
29 
38 struct fitpar {
39 
40  fitpar() : val(0.0), err(0.0) {}
41 
45  float val;
46 
50  float err;
51 
55  friend STD_ostream& operator << (STD_ostream& s, const fitpar& fp) {
56  return s << fp.val << " +/- " << fp.err;
57  }
58 
59 };
60 
61 
63 
77 
78  public:
79 
83  virtual float evaluate_f(float x) const = 0;
84 
88  virtual fvector evaluate_df(float x) const = 0;
89 
93  virtual unsigned int numof_fitpars() const = 0;
94 
98  virtual fitpar& get_fitpar(unsigned int i) = 0;
99 
103  Array<float,1> get_function(const Array<float,1>& xvals) const;
104 
105 
106  // dummy array used for default arguments
107  static const Array<float,1> defaultArray;
108 
109  protected:
110  ModelFunction() {}
111  virtual ~ModelFunction() {}
112 
113  fitpar dummy_fitpar;
114 
115 };
116 
118 
123 
124  public:
125 
126  virtual ~FunctionFitInterface() {}
127 
131  virtual bool init(ModelFunction& model_func, unsigned int nvals) = 0;
132 
142  virtual bool fit(const Array<float,1>& yvals,
143  const Array<float,1>& ysigma=defaultArray,
144  const Array<float,1>& xvals=defaultArray,
145  unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE) = 0;
146 
147  // dummy array used for default arguments
148  static const Array<float,1> defaultArray;
149 };
150 
152 
153 class ModelData; // forward declaration
154 class GslData4Fit; // forward declaration
155 
160 
161  public:
162 
166  FunctionFitDerivative() : gsldata(0), data4fit(0) {}
167 
172 
173  // overloading virtual functions of FunctionFitInterface
174  bool init(ModelFunction& model_func, unsigned int nvals);
175  bool fit(const Array<float,1>& yvals,
176  const Array<float,1>& ysigma=defaultArray,
177  const Array<float,1>& xvals=defaultArray,
178  unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE);
179 
180 
181  private:
182 
183  void print_state (size_t iter);
184 
185  GslData4Fit* gsldata;
186  ModelData* data4fit;
187 };
188 
189 
190 
192 
193 
201 
202  fitpar A;
203  fitpar lambda;
204 
205  // implementing virtual functions of ModelFunction
206  float evaluate_f(float x) const;
207  fvector evaluate_df(float x) const;
208  unsigned int numof_fitpars() const;
209  fitpar& get_fitpar(unsigned int i);
210 };
211 
212 
214 
222 
223  fitpar A;
224  fitpar lambda;
225  fitpar C;
226 
227  // implementing virtual functions of ModelFunction
228  float evaluate_f(float x) const;
229  fvector evaluate_df(float x) const;
230  unsigned int numof_fitpars() const;
231  fitpar& get_fitpar(unsigned int i);
232 };
233 
235 
236 
244 
245  fitpar A;
246  fitpar x0;
247  fitpar fwhm;
248 
249  // implementing virtual functions of ModelFunction
250  float evaluate_f(float x) const;
251  fvector evaluate_df(float x) const;
252  unsigned int numof_fitpars() const;
253  fitpar& get_fitpar(unsigned int i);
254 };
255 
256 
258 
259 
266 struct SinusFunction : public ModelFunction {
267 
268  fitpar A;
269  fitpar m;
270  fitpar c;
271 
272  // implementing virtual functions of ModelFunction
273  float evaluate_f(float x) const;
274  fvector evaluate_df(float x) const;
275  unsigned int numof_fitpars() const;
276  fitpar& get_fitpar(unsigned int i);
277 };
278 
280 
281 
289 
294  void set_pars(float alphaval, float xmax, float ymax);
295 
296  fitpar A;
297  fitpar alpha;
298  fitpar beta;
299 
300  // implementing virtual functions of ModelFunction
301  float evaluate_f(float x) const;
302  fvector evaluate_df(float x) const;
303  unsigned int numof_fitpars() const;
304  fitpar& get_fitpar(unsigned int i);
305 };
306 
308 
317 template <int N_rank>
319 
320  fitpar a[N_rank+1];
321 
332  bool fit(const Array<float,1>& yvals,
333  const Array<float,1>& ysigma,
334  const Array<float,1>& xvals);
335 
336 /*
337  bool fit(const Array<float,1>& yvals,
338  const Array<float,1>& ysigma){
339  firstIndex fi;
340  Array<float,1> xvals(yvals.size());
341  xvals=fi;
342  return fit(yvals,ysigma,xvals);
343  };
344 */
345 
346  bool fit(const Array<float,1>& yvals){
347  Array<float,1> ysigma(yvals.size());
348  ysigma=1.;
349  return fit(yvals,ysigma);
350  };
351 
356  Array<float,1> get_function(const Array<float,1>& xvals) const;
357 
358 };
359 
360 
361 template <int N_rank>
362 bool PolynomialFunction<N_rank>::fit(const Array<float,1>& yvals, const Array<float,1>& ysigma, const Array<float,1>& xvals) {
363 
364  int npol=N_rank+1;
365  for(int i=0; i<npol; i++) a[i]=fitpar(); // reset
366 
367  int npts=yvals.size();
368 
369  Array<float,1> sigma(npts);
370  if(int(ysigma.size())==npts) sigma=ysigma;
371  else sigma=1.0;
372 
373  Array<float,1> x(npts);
374  if(int(xvals.size())==npts) x=xvals;
375  else for(int ipt=0; ipt<npts; ipt++) x(ipt)=ipt;
376 
377 
378  Array<float,2> A(npts,npol);
379  Array<float,1> b(npts);
380 
381 
382  for(int ipt=0; ipt<npts; ipt++) {
383  float weight=secureInv( sigma(ipt));
384 
385  b(ipt)=weight*yvals(ipt);
386 
387  for(int ipol=0; ipol<npol; ipol++) {
388  A(ipt,ipol)=weight*pow(x(ipt),ipol);
389  }
390  }
391 
392  Array<float,1> coeff(solve_linear(A,b));
393 
394  for(int ipol=0; ipol<npol; ipol++) a[ipol].val=coeff(ipol);
395 
396  return true;
397 }
398 
399 
400 template <int N_rank>
401 Array<float,1> PolynomialFunction<N_rank>::get_function(const Array<float,1>& xvals) const {
402  int npts=xvals.size();
403  Array<float,1> result(npts); result=0.0;
404 
405  for(int ipt=0; ipt<npts; ipt++) {
406  for(int ipol=0; ipol<(N_rank+1); ipol++) {
407  result(ipt)+=a[ipol].val*pow(xvals(ipt),ipol);
408  }
409  }
410 
411  return result;
412 }
413 
414 
416 
417 
427 
428  fitpar m;
429  fitpar c;
430 
441  bool fit(const Array<float,1>& yvals,
442  const Array<float,1>& ysigma=defaultArray,
443  const Array<float,1>& xvals=defaultArray);
444 
449  Array<float,1> get_function(const Array<float,1>& xvals) const;
450 
451 
452  // dummy array used for default arguments
453  static const Array<float,1> defaultArray;
454 };
455 
456 
457 
459 
460 class GslData4DownhillSimplex; // forward declaration
461 
466 
467  public:
468 
474 
479 
488  bool get_minimum_parameters(fvector& result, const fvector& starting_point, const fvector& step_size, unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE);
489 
490 
491  private:
492  unsigned int ndim;
493  GslData4DownhillSimplex* gsldata;
494 
495 };
496 
497 
499 
504 
505  public:
506 
511 
516 
517  // overloading virtual functions of FunctionFitInterface
518  bool init(ModelFunction& model_func, unsigned int nvals);
519  bool fit(const Array<float,1>& yvals,
520  const Array<float,1>& ysigma=defaultArray,
521  const Array<float,1>& xvals=defaultArray,
522  unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE);
523 
524 
525  // overloading virtual functions of MinimizationFunction
526  unsigned int numof_fitpars() const;
527  float evaluate(const fvector& pars) const;
528 
529 
530  private:
531  ModelFunction* func;
532  DownhillSimplex* ds;
533  Array<float,1> yvals_cache;
534  Array<float,1> ysigma_cache;
535  Array<float,1> xvals_cache;
536 };
537 
539 
553 template <int N_rank>
554 Array<float,N_rank> polyniomial_fit(const Array<float,N_rank>& value_map, const Array<float,N_rank>& reliability_map,
555  unsigned int polynom_order, float kernel_size, bool only_zero_reliability=false) {
556  Log<OdinData> odinlog("","polyniomial_fit");
557 
558  Data<float,N_rank> result(value_map.shape());
559  result=0.0;
560 
561  if(!same_shape(value_map,reliability_map)) {
562  ODINLOG(odinlog,errorLog) << "size mismatch (value_map.shape()=" << value_map.shape() << ") != (reliability_map.shape()=" << reliability_map.shape() << ")" << STD_endl;
563  return result;
564  }
565 
566  if(min(reliability_map)<0.0) {
567  ODINLOG(odinlog,errorLog) << "reliability_map must be non-negative" << STD_endl;
568  return result;
569  }
570 
571  int minsize=max(value_map.shape());
572  for(int idim=0; idim<N_rank; idim++) {
573  int dimsize=value_map.shape()(idim);
574  if( (dimsize>1) && (dimsize<minsize) ) minsize=dimsize;
575  }
576  if(minsize<=0) {
577  return result;
578  }
579 
580  if((minsize-1)<int(polynom_order)) {
581  polynom_order=minsize-1;
582  ODINLOG(odinlog,warningLog) << "array size too small, restricting polynom_order to " << polynom_order << STD_endl;
583  }
584 
585  TinyVector<int,N_rank> valshape(value_map.shape());
586  int nvals=value_map.numElements();
587 
588  TinyVector<int,N_rank> polsize; polsize=polynom_order+1;
589  Data<int,N_rank> polarr(polsize);
590  int npol=polarr.numElements();
591 
592  if(pow(kernel_size,float(N_rank))<float(npol)) {
593  kernel_size=pow(double(npol),double(1.0/float(N_rank)));
594  ODINLOG(odinlog,warningLog) << "kernel_size too small for polynome, increasing to " << kernel_size << STD_endl;
595  }
596 
597 
598 
599  int neighb_pixel=int(kernel_size);
600  if(neighb_pixel<=0) neighb_pixel=1;
601  TinyVector<int,N_rank> neighbsize; neighbsize=2*neighb_pixel+1;
602  TinyVector<int,N_rank> neighboffset; neighboffset=-neighb_pixel;
603  Data<int,N_rank> neighbarr(neighbsize); // neighbour grid around root pixel
604  int nneighb=neighbarr.numElements();
605 
606  ODINLOG(odinlog,normalDebug) << "nvals/npol/nneighb=" << nvals << "/" << npol << "/" << nneighb << STD_endl;
607 
608  if(npol>nneighb) {
609  ODINLOG(odinlog,warningLog) << "polynome order (" << npol << ") larger than number of neighbours (" << nneighb << ")" << STD_endl;
610  }
611 
612 
613  Array<float,2> A(npol,npol);
614  Array<float,1> c(npol);
615  Array<float,1> b(npol);
616 
617  TinyVector<int,N_rank> valindex;
618  TinyVector<int,N_rank> neighbindex;
619  TinyVector<int,N_rank> currindex;
620  TinyVector<int,N_rank> diffindex;
621  TinyVector<int,N_rank> polindex;
622  TinyVector<int,N_rank> polindex_sum;
623 
624  float epsilon=0.01;
625  float relevant_radius=0.5*kernel_size*sqrt(double(N_rank))+epsilon;
626 
627  // iterate through pixels of value_map
628  for(int ival=0; ival<nvals; ival++) {
629  valindex=result.create_index(ival);
630 
631  if( (!only_zero_reliability) || (reliability_map(valindex)<=0.0) ) { // fit only pixel with zero reliability
632 
633  A=0.0;
634  b=0.0;
635 
636  int n_relevant_neighb_pixel=0;
637 
638  // iterate through neigbourhood of pixel and accumulate them in a single
639  // set of equations, weighted by their reliability
640  for(int ineighb=0; ineighb<nneighb; ineighb++) {
641  neighbindex=neighbarr.create_index(ineighb);
642  currindex=valindex+neighboffset+neighbindex;
643 
644  bool valid_pixel=true;
645 
646  // is the pixel within value_map ?
647  for(int irank=0; irank<N_rank; irank++) {
648  if(currindex(irank)<0 || currindex(irank)>=valshape(irank)) valid_pixel=false;
649  }
650 
651  // does the pixel have non-vanishing reliability
652  float reliability=0.0;
653  if(valid_pixel) reliability=reliability_map(currindex);
654  if(reliability<=0.0) valid_pixel=false;
655 
656 
657  if(valid_pixel) {
658 
659  diffindex=currindex-valindex; // (xk-x0,yk-y0,...)
660 
661  float radiussqr=sum(diffindex*diffindex);
662  float weight=reliability*exp(-2.0*radiussqr/(kernel_size*kernel_size));
663 
664  if(weight>0.0) {
665  if(sqrt(radiussqr)<=relevant_radius) n_relevant_neighb_pixel++;
666 
667  // create b_i,j
668  for(int ipol=0; ipol<npol; ipol++) {
669  polindex=polarr.create_index(ipol); // (i,j,..)
670  float polproduct=1.0;
671  for(int irank=0; irank<N_rank; irank++) polproduct*=pow(float(diffindex(irank)),float(polindex(irank)));
672  b(ipol)+=weight*value_map(currindex)*polproduct;
673  }
674 
675  // create A_ii',jj'
676  for(int ipol=0; ipol<npol; ipol++) {
677  for(int ipol_prime=0; ipol_prime<npol; ipol_prime++) {
678  polindex_sum=polarr.create_index(ipol)+polarr.create_index(ipol_prime);
679  float polproduct=1.0;
680  for(int irank=0; irank<N_rank; irank++) polproduct*=pow(float(diffindex(irank)),float(polindex_sum(irank)));
681  A(ipol,ipol_prime)+=weight*polproduct;
682  }
683  }
684 
685  }
686 
687 
688  }
689  }
690 
691  if(n_relevant_neighb_pixel>=npol) { // do we have enough pixel for the fit ?
692  c=solve_linear(A,b);
693  result(valindex)=c(0);
694  }
695 
696  } else result(valindex)=value_map(valindex);
697 
698  }
699 
700  return result;
701 }
702 
703 
708 #endif
709 
TinyVector< int, N_rank > create_index(unsigned long index) const
DownhillSimplex(MinimizationFunction &function)
bool get_minimum_parameters(fvector &result, const fvector &starting_point, const fvector &step_size, unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE)
bool fit(const Array< float, 1 > &yvals, const Array< float, 1 > &ysigma=defaultArray, const Array< float, 1 > &xvals=defaultArray, unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE)
bool init(ModelFunction &model_func, unsigned int nvals)
bool fit(const Array< float, 1 > &yvals, const Array< float, 1 > &ysigma=defaultArray, const Array< float, 1 > &xvals=defaultArray, unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE)
float evaluate(const fvector &pars) const
bool init(ModelFunction &model_func, unsigned int nvals)
unsigned int numof_fitpars() const
virtual bool init(ModelFunction &model_func, unsigned int nvals)=0
virtual bool fit(const Array< float, 1 > &yvals, const Array< float, 1 > &ysigma=defaultArray, const Array< float, 1 > &xvals=defaultArray, unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE)=0
Definition: tjlog.h:218
virtual unsigned int numof_fitpars() const =0
virtual fvector evaluate_df(float x) const =0
virtual fitpar & get_fitpar(unsigned int i)=0
virtual float evaluate_f(float x) const =0
Array< float, 1 > get_function(const Array< float, 1 > &xvals) const
Definition: tjarray.h:41
Data< float, 1 > solve_linear(const Data< float, 2 > &A, const Data< float, 1 > &b, float sv_truncation=0.0)
Array< float, N_rank > polyniomial_fit(const Array< float, N_rank > &value_map, const Array< float, N_rank > &reliability_map, unsigned int polynom_order, float kernel_size, bool only_zero_reliability=false)
Definition: fitting.h:554
bool same_shape(const Array< T, N_rank > &a1, const Array< T, N_rank > &a2, const TinyVector< int, N_rank > &dimmask=1)
Definition: utils.h:168
Array< float, 1 > get_function(const Array< float, 1 > &xvals) const
Definition: fitting.h:401
bool fit(const Array< float, 1 > &yvals, const Array< float, 1 > &ysigma, const Array< float, 1 > &xvals)
Definition: fitting.h:362
double secureInv(double denominator)
Definition: tjtools.h:271
unsigned int numof_fitpars() const
fitpar & get_fitpar(unsigned int i)
float evaluate_f(float x) const
fvector evaluate_df(float x) const
float evaluate_f(float x) const
fitpar & get_fitpar(unsigned int i)
fvector evaluate_df(float x) const
unsigned int numof_fitpars() const
fitpar & get_fitpar(unsigned int i)
float evaluate_f(float x) const
void set_pars(float alphaval, float xmax, float ymax)
unsigned int numof_fitpars() const
fvector evaluate_df(float x) const
float evaluate_f(float x) const
fvector evaluate_df(float x) const
fitpar & get_fitpar(unsigned int i)
unsigned int numof_fitpars() const
bool fit(const Array< float, 1 > &yvals, const Array< float, 1 > &ysigma=defaultArray, const Array< float, 1 > &xvals=defaultArray)
Array< float, 1 > get_function(const Array< float, 1 > &xvals) const
unsigned int numof_fitpars() const
fvector evaluate_df(float x) const
fitpar & get_fitpar(unsigned int i)
float evaluate_f(float x) const
Definition: fitting.h:38
float err
Definition: fitting.h:50
float val
Definition: fitting.h:45
friend STD_ostream & operator<<(STD_ostream &s, const fitpar &fp)
Definition: fitting.h:55