fitting.h

/***************************************************************************
                          fitting.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2001 by Thies Jochimsen & Michael von Mengershausen
    email                : jochimse@cns.mpg.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 <odindata/data.h>
#include <odindata/linalg.h>
#include <odindata/utils.h>


struct fitpar {

  fitpar() : val(0.0), err(0.0) {}

  float val;

  float 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;


  bool fit(const Array<float,1>& yvals,
           const Array<float,1>& ysigma=defaultArray,
           const Array<float,1>& xvals=defaultArray,
           unsigned int max_iterations=500, double tolerance=1.0e-4);
  
  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 ModelData; // forward declaration
class GslData4Fit; // forward declaration

class FunctionFit {

 public:

  FunctionFit(ModelFunction& model_func, unsigned int nvals, unsigned int max_iterations=500, double tolerance=1.0e-4);

  ~FunctionFit();

  bool fit(const Array<float,1>& yvals,
           const Array<float,1>& ysigma=defaultArray,
           const Array<float,1>& xvals=defaultArray);


  // dummy array used for default arguments
  static const Array<float,1> defaultArray;

 private:

  void print_state (size_t iter);

  ModelFunction& func;
  unsigned int max_it;
  double tol;

  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);
};


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(ysigma.size()==npts) sigma=ysigma;
  else sigma=1.0;

  Array<float,1> x(npts);
  if(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;
};





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

---