fitting.h
1 /***************************************************************************
2  fitting.h - description
3  -------------------
4  begin : Fri Apr 6 2001
5  copyright : (C) 2000-2015 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  bool fit(const Array<float,1>& yvals){
346  Array<float,1> ysigma(yvals.size());
347  ysigma=1.;
348  return fit(yvals,ysigma);
349  };
350 
355  Array<float,1> get_function(const Array<float,1>& xvals) const;
356 
357 };
358 
359 
360 template <int N_rank>
361 bool PolynomialFunction<N_rank>::fit(const Array<float,1>& yvals, const Array<float,1>& ysigma, const Array<float,1>& xvals) {
362 
363  int npol=N_rank+1;
364  for(int i=0; i<npol; i++) a[i]=fitpar(); // reset
365 
366  int npts=yvals.size();
367 
368  Array<float,1> sigma(npts);
369  if(int(ysigma.size())==npts) sigma=ysigma;
370  else sigma=1.0;
371 
372  Array<float,1> x(npts);
373  if(int(xvals.size())==npts) x=xvals;
374  else for(int ipt=0; ipt<npts; ipt++) x(ipt)=ipt;
375 
376 
377  Array<float,2> A(npts,npol);
378  Array<float,1> b(npts);
379 
380 
381  for(int ipt=0; ipt<npts; ipt++) {
382  float weight=secureInv( sigma(ipt));
383 
384  b(ipt)=weight*yvals(ipt);
385 
386  for(int ipol=0; ipol<npol; ipol++) {
387  A(ipt,ipol)=weight*pow(x(ipt),ipol);
388  }
389  }
390 
391  Array<float,1> coeff(solve_linear(A,b));
392 
393  for(int ipol=0; ipol<npol; ipol++) a[ipol].val=coeff(ipol);
394 
395  return true;
396 }
397 
398 
399 template <int N_rank>
400 Array<float,1> PolynomialFunction<N_rank>::get_function(const Array<float,1>& xvals) const {
401  int npts=xvals.size();
402  Array<float,1> result(npts); result=0.0;
403 
404  for(int ipt=0; ipt<npts; ipt++) {
405  for(int ipol=0; ipol<(N_rank+1); ipol++) {
406  result(ipt)+=a[ipol].val*pow(xvals(ipt),ipol);
407  }
408  }
409 
410  return result;
411 }
412 
413 
415 
416 
426 
427  fitpar m;
428  fitpar c;
429 
440  bool fit(const Array<float,1>& yvals,
441  const Array<float,1>& ysigma=defaultArray,
442  const Array<float,1>& xvals=defaultArray);
443 
448  Array<float,1> get_function(const Array<float,1>& xvals) const;
449 
450 
451  // dummy array used for default arguments
452  static const Array<float,1> defaultArray;
453 };
454 
455 
456 
458 
459 class GslData4DownhillSimplex; // forward declaration
460 
465 
466  public:
467 
473 
477  ~DownhillSimplex();
478 
486  fvector get_minimum_parameters(const fvector& starting_point, const fvector& step_size, unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE);
487 
488 
489  private:
490  unsigned int ndim;
491  GslData4DownhillSimplex* gsldata;
492 
493 };
494 
495 
497 
502 
503  public:
504 
509 
514 
515  // overloading virtual functions of FunctionFitInterface
516  bool init(ModelFunction& model_func, unsigned int nvals);
517  bool fit(const Array<float,1>& yvals,
518  const Array<float,1>& ysigma=defaultArray,
519  const Array<float,1>& xvals=defaultArray,
520  unsigned int max_iterations=DEFAULT_MAX_ITER, double tolerance=DEFAULT_TOLERANCE);
521 
522 
523  // overloading virtual functions of MinimizationFunction
524  unsigned int numof_fitpars() const;
525  float evaluate(const fvector& pars) const;
526 
527 
528  private:
529  ModelFunction* func;
530  DownhillSimplex* ds;
531  Array<float,1> yvals_cache;
532  Array<float,1> ysigma_cache;
533  Array<float,1> xvals_cache;
534 };
535 
537 
551 template <int N_rank>
552 Array<float,N_rank> polyniomial_fit(const Array<float,N_rank>& value_map, const Array<float,N_rank>& reliability_map,
553  unsigned int polynom_order, float kernel_size, bool only_zero_reliability=false) {
554  Log<OdinData> odinlog("","polyniomial_fit");
555 
556  Data<float,N_rank> result(value_map.shape());
557  result=0.0;
558 
559  if(!same_shape(value_map,reliability_map)) {
560  ODINLOG(odinlog,errorLog) << "size mismatch (value_map.shape()=" << value_map.shape() << ") != (reliability_map.shape()=" << reliability_map.shape() << ")" << STD_endl;
561  return result;
562  }
563 
564  if(min(reliability_map)<0.0) {
565  ODINLOG(odinlog,errorLog) << "reliability_map must be non-negative" << STD_endl;
566  return result;
567  }
568 
569  int minsize=max(value_map.shape());
570  for(int idim=0; idim<N_rank; idim++) {
571  int dimsize=value_map.shape()(idim);
572  if( (dimsize>1) && (dimsize<minsize) ) minsize=dimsize;
573  }
574  if(minsize<=0) {
575  return result;
576  }
577 
578  if((minsize-1)<int(polynom_order)) {
579  polynom_order=minsize-1;
580  ODINLOG(odinlog,warningLog) << "array size too small, restricting polynom_order to " << polynom_order << STD_endl;
581  }
582 
583  TinyVector<int,N_rank> valshape(value_map.shape());
584  int nvals=value_map.numElements();
585 
586  TinyVector<int,N_rank> polsize; polsize=polynom_order+1;
587  Data<int,N_rank> polarr(polsize);
588  int npol=polarr.numElements();
589 
590  if(pow(kernel_size,float(N_rank))<float(npol)) {
591  kernel_size=pow(double(npol),double(1.0/float(N_rank)));
592  ODINLOG(odinlog,warningLog) << "kernel_size too small for polynome, increasing to " << kernel_size << STD_endl;
593  }
594 
595 
596 
597  int neighb_pixel=int(kernel_size);
598  if(neighb_pixel<=0) neighb_pixel=1;
599  TinyVector<int,N_rank> neighbsize; neighbsize=2*neighb_pixel+1;
600  TinyVector<int,N_rank> neighboffset; neighboffset=-neighb_pixel;
601  Data<int,N_rank> neighbarr(neighbsize); // neighbour grid around root pixel
602  int nneighb=neighbarr.numElements();
603 
604  ODINLOG(odinlog,normalDebug) << "nvals/npol/nneighb=" << nvals << "/" << npol << "/" << nneighb << STD_endl;
605 
606  if(npol>nneighb) {
607  ODINLOG(odinlog,warningLog) << "polynome order (" << npol << ") larger than number of neighbours (" << nneighb << ")" << STD_endl;
608  }
609 
610 
611  Array<float,2> A(npol,npol);
612  Array<float,1> c(npol);
613  Array<float,1> b(npol);
614 
615  TinyVector<int,N_rank> valindex;
616  TinyVector<int,N_rank> neighbindex;
617  TinyVector<int,N_rank> currindex;
618  TinyVector<int,N_rank> diffindex;
619  TinyVector<int,N_rank> polindex;
620  TinyVector<int,N_rank> polindex_sum;
621 
622  float epsilon=0.01;
623  float relevant_radius=0.5*kernel_size*sqrt(double(N_rank))+epsilon;
624 
625  // iterate through pixels of value_map
626  for(int ival=0; ival<nvals; ival++) {
627  valindex=result.create_index(ival);
628 
629  if( (!only_zero_reliability) || (reliability_map(valindex)<=0.0) ) { // fit only pixel with zero reliability
630 
631  A=0.0;
632  b=0.0;
633 
634  int n_relevant_neighb_pixel=0;
635 
636  // iterate through neigbourhood of pixel and accumulate them in a single
637  // set of equations, weighted by their reliability
638  for(int ineighb=0; ineighb<nneighb; ineighb++) {
639  neighbindex=neighbarr.create_index(ineighb);
640  currindex=valindex+neighboffset+neighbindex;
641 
642  bool valid_pixel=true;
643 
644  // is the pixel within value_map ?
645  for(int irank=0; irank<N_rank; irank++) {
646  if(currindex(irank)<0 || currindex(irank)>=valshape(irank)) valid_pixel=false;
647  }
648 
649  // does the pixel have non-vanishing reliability
650  float reliability=0.0;
651  if(valid_pixel) reliability=reliability_map(currindex);
652  if(reliability<=0.0) valid_pixel=false;
653 
654 
655  if(valid_pixel) {
656 
657  diffindex=currindex-valindex; // (xk-x0,yk-y0,...)
658 
659  float radiussqr=sum(diffindex*diffindex);
660  float weight=reliability*exp(-2.0*radiussqr/(kernel_size*kernel_size));
661 
662  if(weight>0.0) {
663  if(sqrt(radiussqr)<=relevant_radius) n_relevant_neighb_pixel++;
664 
665  // create b_i,j
666  for(int ipol=0; ipol<npol; ipol++) {
667  polindex=polarr.create_index(ipol); // (i,j,..)
668  float polproduct=1.0;
669  for(int irank=0; irank<N_rank; irank++) polproduct*=pow(float(diffindex(irank)),float(polindex(irank)));
670  b(ipol)+=weight*value_map(currindex)*polproduct;
671  }
672 
673  // create A_ii',jj'
674  for(int ipol=0; ipol<npol; ipol++) {
675  for(int ipol_prime=0; ipol_prime<npol; ipol_prime++) {
676  polindex_sum=polarr.create_index(ipol)+polarr.create_index(ipol_prime);
677  float polproduct=1.0;
678  for(int irank=0; irank<N_rank; irank++) polproduct*=pow(float(diffindex(irank)),float(polindex_sum(irank)));
679  A(ipol,ipol_prime)+=weight*polproduct;
680  }
681  }
682 
683  }
684 
685 
686  }
687  }
688 
689  if(n_relevant_neighb_pixel>=npol) { // do we have enough pixel for the fit ?
690  c=solve_linear(A,b);
691  result(valindex)=c(0);
692  }
693 
694  } else result(valindex)=value_map(valindex);
695 
696  }
697 
698  return result;
699 }
700 
701 
706 #endif
707 
TinyVector< int, N_rank > create_index(unsigned long index) const
Definition: tjlog.h:218
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
double secureInv(double denominator)
Definition: tjtools.h:271
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:552
bool fit(const Array< float, 1 > &yvals, const Array< float, 1 > &ysigma, const Array< float, 1 > &xvals)
Definition: fitting.h:361
friend STD_ostream & operator<<(STD_ostream &s, const fitpar &fp)
Definition: fitting.h:55
Definition: fitting.h:38
Data< float, 1 > solve_linear(const Data< float, 2 > &A, const Data< float, 1 > &b, float sv_truncation=0.0)
Definition: tjarray.h:41
float val
Definition: fitting.h:45
Array< float, 1 > get_function(const Array< float, 1 > &xvals) const
Definition: fitting.h:400
float err
Definition: fitting.h:50