NewtonFitterGSL.h

Go to the documentation of this file.

#ifndef __NEWTONFITTERGSL_H
#define __NEWTONFITTERGSL_H

#include "BaseFitter.h"

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_eigen.h>

// Class NewtonFitterGSL

class NewtonFitterGSL : public BaseFitter {
  public:
    NewtonFitterGSL();
    virtual ~NewtonFitterGSL();
    
    virtual double fit();
    
    virtual int getError() const;
    
    virtual double getProbability() const;
    virtual double getChi2() const;
    virtual int    getDoF() const;
    virtual int  getIterations() const;

    virtual int    getNcon() const;
    
    virtual int    getNsoft() const;
    
    virtual int    getNpar() const;
    
    virtual int    getNunm() const;

    virtual bool initialize();
  
    virtual void setDebug (int debuglevel);
  
  protected:
    virtual double calcChi2();
    
    int calcDx ();
    
    int calcDxSVD ();
    
    void printMy (double M[], double y[], int idim);
    
    bool updateParams (gsl_vector *xnew);
    void fillxold();
    void fillperr();
    int calcM (bool errorpropagation = false);
    int calcy();
    int optimizeScale();
    int invertM();
  
    void calcCovMatrix();
    
    double meritFunction(double mu);
    double meritFunctionDeriv(double mu);
    
    enum {NPARMAX=50, NCONMAX=10, NUNMMAX=10};
    
    int npar;      
    int ncon;      
    int nsoft;     
    int nunm;      
    int ierr;      
    int nit;       

    double fitprob;   
    double chi2;      
    
    static void ini_gsl_permutation (gsl_permutation *&p, unsigned int size);
    static void ini_gsl_vector (gsl_vector *&v, int unsigned size);
    static void ini_gsl_matrix (gsl_matrix *&m, int unsigned size1, unsigned int size2);
    
    static void debug_print (gsl_matrix *m, const char *name);
    static void debug_print (gsl_vector *v, const char *name);

  private:
    unsigned int idim;
    gsl_vector *x;
    gsl_vector *xold;
    gsl_vector *xbest;
    gsl_vector *dx;
    gsl_vector *dxscal;
    gsl_vector *grad;
    gsl_vector *y;
    gsl_vector *yscal;
    gsl_vector *perr;
    gsl_vector *v1;
    gsl_vector *v2;
    gsl_vector *Meval;
  
    gsl_matrix *M;
    gsl_matrix *Mscal;
    gsl_matrix *M1;
    gsl_matrix *M2;
    gsl_matrix *M3;
    gsl_matrix *M4;
    gsl_matrix *M5;
    gsl_matrix *Mevec;
    gsl_matrix *CC;
    gsl_matrix *CC1;
    gsl_matrix *CCinv;

    gsl_permutation *permM;
    gsl_eigen_symmv_workspace *ws; 
    unsigned int wsdim;
    
    double chi2best;
    double chi2new;
    double chi2old;
    double fvalbest;
    double scale;
    double scalebest;
    double stepsize;
    double stepbest;
    enum {NITMAX = 100};
    double scalevals[NITMAX];
    double fvals[NITMAX];
    
    int imerit;
    
    int debug;
};

#endif // __NEWTONFITTERGSL_H