fitting.h

/***************************************************************************
                          fitting.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2000-2021 by Thies Jochimsen & Michael von Mengershausen
    email                : thies@jochimsen.de  mengers@cns.mpg.de
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef FITTING_H
#define FITTING_H
 
#include <tjutils/tjnumeric.h> // for MinimizationFunction
#include <odindata/data.h>
#include <odindata/linalg.h>
#include <odindata/utils.h>
 
 
#define DEFAULT_MAX_ITER 1000
#define DEFAULT_TOLERANCE 1e-4
 
struct fitpar {
 
  fitpar() : val(0.0), err(0.0) {}
 
  float val;
 
  float err;
 
  friend STD_ostream& operator << (STD_ostream& s, const fitpar& fp) {
    return s << fp.val << " +/- " << fp.err;
  }
 
};
 
 
 
class ModelFunction {
 
 public:
 
  virtual float evaluate_f(float x) const = 0;
 
  virtual fvector evaluate_df(float x) const = 0;
 
  virtual unsigned int numof_fitpars() const = 0;
 
  virtual fitpar& get_fitpar(unsigned int i) = 0;
 
  Array<float,1> get_function(const Array<float,1>& xvals) const;
  
 
  // dummy array used for default arguments
  static const Array<float,1> defaultArray;
 
 protected:
   ModelFunction() {}
   virtual ~ModelFunction() {}
 
   fitpar dummy_fitpar;
   
};
 
 
class FunctionFitInterface {
 
 public:
 
  virtual ~FunctionFitInterface() {}
 
  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;
 
  // dummy array used for default arguments
  static const Array<float,1> defaultArray;
};
 
 
class ModelData; // forward declaration
class GslData4Fit; // forward declaration
 
class FunctionFitDerivative : public virtual FunctionFitInterface {
 
 public:
 
   FunctionFitDerivative() : gsldata(0), data4fit(0) {}
 
  ~FunctionFitDerivative();
 
  // overloading virtual functions of FunctionFitInterface
  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);
 
 
 private:
 
  void print_state (size_t iter);
 
  GslData4Fit* gsldata;
  ModelData* data4fit;
};
 
 
 
 
 
struct ExponentialFunction : public ModelFunction {
 
  fitpar A;
  fitpar lambda;
 
  // implementing virtual functions of ModelFunction
  float evaluate_f(float x) const;
  fvector evaluate_df(float x) const;
  unsigned int numof_fitpars() const;
  fitpar& get_fitpar(unsigned int i);
};
 
 
 
struct ExponentialFunctionWithOffset : public ModelFunction {
 
  fitpar A;
  fitpar lambda;
  fitpar C;
 
  // implementing virtual functions of ModelFunction
  float evaluate_f(float x) const;
  fvector evaluate_df(float x) const;
  unsigned int numof_fitpars() const;
  fitpar& get_fitpar(unsigned int i);
};
 
 
 
struct GaussianFunction : public ModelFunction {
 
  fitpar A;
  fitpar x0;
  fitpar fwhm;
 
  // implementing virtual functions of ModelFunction
  float evaluate_f(float x) const;
  fvector evaluate_df(float x) const;
  unsigned int numof_fitpars() const;
  fitpar& get_fitpar(unsigned int i);
};
 
 
 
 
struct SinusFunction : public ModelFunction {
 
  fitpar A;
  fitpar m;
  fitpar c;
 
  // implementing virtual functions of ModelFunction
  float evaluate_f(float x) const;
  fvector evaluate_df(float x) const;
  unsigned int numof_fitpars() const;
  fitpar& get_fitpar(unsigned int i);
};
 
 
 
struct GammaVariateFunction : public ModelFunction {
 
  void set_pars(float alphaval, float xmax, float ymax);
 
  fitpar A;
  fitpar alpha;
  fitpar beta;
 
  // implementing virtual functions of ModelFunction
  float evaluate_f(float x) const;
  fvector evaluate_df(float x) const;
  unsigned int numof_fitpars() const;
  fitpar& get_fitpar(unsigned int i);
};
 
 
template <int N_rank>
struct PolynomialFunction {
 
  fitpar a[N_rank+1];
 
  bool fit(const Array<float,1>& yvals,
           const Array<float,1>& ysigma,
           const Array<float,1>& xvals);
 
/*
  bool fit(const Array<float,1>& yvals,
           const Array<float,1>& ysigma){
    firstIndex fi;
    Array<float,1> xvals(yvals.size());
    xvals=fi;
    return fit(yvals,ysigma,xvals);
  };
*/
 
  bool fit(const Array<float,1>& yvals){
    Array<float,1> ysigma(yvals.size());
    ysigma=1.;
    return fit(yvals,ysigma);
  };
  
  Array<float,1> get_function(const Array<float,1>& xvals) const;
 
};
 
 
template <int N_rank>
bool PolynomialFunction<N_rank>::fit(const Array<float,1>& yvals, const Array<float,1>& ysigma, const Array<float,1>& xvals) {
 
  int npol=N_rank+1;
  for(int i=0; i<npol; i++) a[i]=fitpar(); // reset
 
  int npts=yvals.size();
 
  Array<float,1> sigma(npts);
  if(int(ysigma.size())==npts) sigma=ysigma;
  else sigma=1.0;
 
  Array<float,1> x(npts);
  if(int(xvals.size())==npts) x=xvals;
  else for(int ipt=0; ipt<npts; ipt++) x(ipt)=ipt;
 
 
  Array<float,2> A(npts,npol);
  Array<float,1> b(npts);
 
 
  for(int ipt=0; ipt<npts; ipt++) {
    float weight=secureInv( sigma(ipt));
 
    b(ipt)=weight*yvals(ipt);
 
    for(int ipol=0; ipol<npol; ipol++) {
      A(ipt,ipol)=weight*pow(x(ipt),ipol);
    }
  }
 
  Array<float,1> coeff(solve_linear(A,b));
 
  for(int ipol=0; ipol<npol; ipol++) a[ipol].val=coeff(ipol);
 
  return true;
}
 
 
template <int N_rank>
Array<float,1> PolynomialFunction<N_rank>::get_function(const Array<float,1>& xvals) const {
  int npts=xvals.size();
  Array<float,1> result(npts); result=0.0;
 
  for(int ipt=0; ipt<npts; ipt++) {
    for(int ipol=0; ipol<(N_rank+1); ipol++) {
      result(ipt)+=a[ipol].val*pow(xvals(ipt),ipol);
    }
  }
 
  return result;
}
 
 
 
 
struct LinearFunction {
 
  fitpar m;
  fitpar c;
 
  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;
 
 
  // dummy array used for default arguments
  static const Array<float,1> defaultArray;
};
 
 
 
 
class GslData4DownhillSimplex; // forward declaration
 
class DownhillSimplex {
 
 public:
 
  DownhillSimplex(MinimizationFunction& function);
 
  ~DownhillSimplex();
 
  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);
 
 
 private:
  unsigned int ndim;
  GslData4DownhillSimplex* gsldata;
 
};
 
 
 
class FunctionFitDownhillSimplex : public virtual FunctionFitInterface, public MinimizationFunction {
 
 public:
 
   FunctionFitDownhillSimplex();
 
  ~FunctionFitDownhillSimplex();
 
  // overloading virtual functions of FunctionFitInterface
  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);
 
 
  // overloading virtual functions of MinimizationFunction
  unsigned int numof_fitpars() const;
  float evaluate(const fvector& pars) const;
 
 
 private:
  ModelFunction* func;
  DownhillSimplex* ds;
  Array<float,1> yvals_cache;
  Array<float,1> ysigma_cache;
  Array<float,1> xvals_cache;
};
 
 
template <int N_rank>
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) {
  Log<OdinData> odinlog("","polyniomial_fit");
 
  Data<float,N_rank> result(value_map.shape());
  result=0.0;
 
  if(!same_shape(value_map,reliability_map)) {
    ODINLOG(odinlog,errorLog) << "size mismatch (value_map.shape()=" << value_map.shape() << ") != (reliability_map.shape()=" << reliability_map.shape() << ")" << STD_endl;
    return result;
  }
 
  if(min(reliability_map)<0.0) {
    ODINLOG(odinlog,errorLog) << "reliability_map must be non-negative" << STD_endl;
    return result;
  }
 
  int minsize=max(value_map.shape());
  for(int idim=0; idim<N_rank; idim++) {
    int dimsize=value_map.shape()(idim);
    if( (dimsize>1) && (dimsize<minsize) ) minsize=dimsize;
  }
  if(minsize<=0) {
    return result;
  }
 
  if((minsize-1)<int(polynom_order)) {
    polynom_order=minsize-1;
    ODINLOG(odinlog,warningLog) << "array size too small, restricting polynom_order to " << polynom_order << STD_endl;
  }
 
  TinyVector<int,N_rank> valshape(value_map.shape());
  int nvals=value_map.numElements();
 
  TinyVector<int,N_rank> polsize; polsize=polynom_order+1;
  Data<int,N_rank> polarr(polsize);
  int npol=polarr.numElements();
 
  if(pow(kernel_size,float(N_rank))<float(npol)) {
    kernel_size=pow(double(npol),double(1.0/float(N_rank)));
    ODINLOG(odinlog,warningLog) << "kernel_size too small for polynome, increasing to " << kernel_size << STD_endl;
  }
 
 
 
  int neighb_pixel=int(kernel_size);
  if(neighb_pixel<=0) neighb_pixel=1;
  TinyVector<int,N_rank> neighbsize; neighbsize=2*neighb_pixel+1;
  TinyVector<int,N_rank> neighboffset; neighboffset=-neighb_pixel;
  Data<int,N_rank> neighbarr(neighbsize); // neighbour grid around root pixel
  int nneighb=neighbarr.numElements();
 
  ODINLOG(odinlog,normalDebug) << "nvals/npol/nneighb=" << nvals << "/" << npol << "/" << nneighb << STD_endl;
 
  if(npol>nneighb) {
    ODINLOG(odinlog,warningLog) << "polynome order (" << npol << ") larger than number of neighbours (" << nneighb << ")" << STD_endl;
  }
 
 
  Array<float,2> A(npol,npol);
  Array<float,1> c(npol);
  Array<float,1> b(npol);
 
  TinyVector<int,N_rank> valindex;
  TinyVector<int,N_rank> neighbindex;
  TinyVector<int,N_rank> currindex;
  TinyVector<int,N_rank> diffindex;
  TinyVector<int,N_rank> polindex;
  TinyVector<int,N_rank> polindex_sum;
 
  float epsilon=0.01;
  float relevant_radius=0.5*kernel_size*sqrt(double(N_rank))+epsilon;
 
  // iterate through pixels of value_map
  for(int ival=0; ival<nvals; ival++) {
    valindex=result.create_index(ival);
 
    if( (!only_zero_reliability) || (reliability_map(valindex)<=0.0) ) { // fit only pixel with zero reliability
 
      A=0.0;
      b=0.0;
 
      int n_relevant_neighb_pixel=0;
 
      // iterate through neigbourhood of pixel and accumulate them in a single
      // set of equations, weighted by their reliability
      for(int ineighb=0; ineighb<nneighb; ineighb++) {
        neighbindex=neighbarr.create_index(ineighb);
        currindex=valindex+neighboffset+neighbindex;
 
        bool valid_pixel=true;
 
        // is the pixel within value_map ?
        for(int irank=0; irank<N_rank; irank++) {
          if(currindex(irank)<0 || currindex(irank)>=valshape(irank)) valid_pixel=false;
        }
 
        // does the pixel have non-vanishing reliability
        float reliability=0.0;
        if(valid_pixel) reliability=reliability_map(currindex);
        if(reliability<=0.0) valid_pixel=false;
 
 
        if(valid_pixel) {
 
          diffindex=currindex-valindex; // (xk-x0,yk-y0,...)
 
          float radiussqr=sum(diffindex*diffindex);
          float weight=reliability*exp(-2.0*radiussqr/(kernel_size*kernel_size));
 
          if(weight>0.0) {
            if(sqrt(radiussqr)<=relevant_radius) n_relevant_neighb_pixel++;
 
            // create b_i,j
            for(int ipol=0; ipol<npol; ipol++) {
              polindex=polarr.create_index(ipol); // (i,j,..)
              float polproduct=1.0;
              for(int irank=0; irank<N_rank; irank++) polproduct*=pow(float(diffindex(irank)),float(polindex(irank)));
              b(ipol)+=weight*value_map(currindex)*polproduct;
            }
 
            // create A_ii',jj'
            for(int ipol=0; ipol<npol; ipol++) {
              for(int ipol_prime=0; ipol_prime<npol; ipol_prime++) {
                polindex_sum=polarr.create_index(ipol)+polarr.create_index(ipol_prime); 
                float polproduct=1.0;
                for(int irank=0; irank<N_rank; irank++) polproduct*=pow(float(diffindex(irank)),float(polindex_sum(irank)));
                A(ipol,ipol_prime)+=weight*polproduct;
              }
            }
 
          }
 
 
        }
      }
 
      if(n_relevant_neighb_pixel>=npol) { // do we have enough pixel for the fit ?
        c=solve_linear(A,b);
        result(valindex)=c(0);
      }
 
    } else result(valindex)=value_map(valindex);
 
  }
 
  return result;
}
 
 
#endif