FittingFunctionGaussianWithBackground.cpp

Go to the documentation of this file.

// FittingFunctionGaussianWithBackground.cpp: implementation of the FittingFunctionGaussianWithBackground class.
//
#ifdef WIN32 
#pragma warning(disable:4786)
#endif

#include "FittingFunctionGaussianWithBackground.h"


// Construction/Destruction

std::string     FittingFunctionGaussianWithBackground::getEquation() const
{
        return "f(x) = Height * exp[ - (x - Position)² / (2 * Width²)] + Background";
}


//##ModelId=43849A020392
FittingFunctionGaussianWithBackground::FittingFunctionGaussianWithBackground()
{
        _iNbParameters                          =       4;
        _lNbMaxIteration                        =       NB_ITERATION_MAX;
        _dSearchLimit                           =       SEARCH_LIMIT;
        _pParameters                            =       new double[_iNbParameters];
        _pParametersErrors                      =       new double[_iNbParameters];
        _mFunctionMap["Position"]       =       0;
        _mFunctionMap["Width"]          =       1;
        _mFunctionMap["Height"]         =       2;
        _mFunctionMap["Background"]     =       3;
        _vInitialGuess = gsl_vector_alloc(4);
}

//##ModelId=43849A02039F
FittingFunctionGaussianWithBackground::FittingFunctionGaussianWithBackground(long lNbIteration,double dLimitSearch) 
{
        _iNbParameters                          =       4;
        _lNbMaxIteration                        =       lNbIteration;
        _dSearchLimit                           =       dLimitSearch;
        _pParameters                            =       new double[_iNbParameters];
        _pParametersErrors                      =       new double[_iNbParameters];
        _mFunctionMap["Position"]       =       0;
        _mFunctionMap["Width"]          =       1;
        _mFunctionMap["Height"]         =       2;
        _mFunctionMap["Background"]     =       3;
        _vInitialGuess = gsl_vector_alloc(4);
}

//##ModelId=43849A0203A2
FittingFunctionGaussianWithBackground::~FittingFunctionGaussianWithBackground()
{
        printf ("Destruction FittingFunctionGaussianWithBackground !!!\n");
        delete[] _pParameters;
        delete[] _pParametersErrors;
        _pParameters = NULL;
        _pParametersErrors = NULL;
        gsl_vector_free(_vInitialGuess);
}

//##ModelId=43849A0203A4
PtrFunction FittingFunctionGaussianWithBackground::getFunction() const{
        return FittingFunctionGaussianWithBackground::Function;
}

//##ModelId=43849A0203A6
PtrFunctionDerivate FittingFunctionGaussianWithBackground::getFunctionDerivate() const{
        return FittingFunctionGaussianWithBackground::FunctionDerivate;
}

//##ModelId=43849A0203A8
PtrFunctionAndDerivate FittingFunctionGaussianWithBackground::getFunctionAndDerivate() const {
        return FittingFunctionGaussianWithBackground::FunctionAndDerivate;
}


//##ModelId=43849A0203AF
int FittingFunctionGaussianWithBackground::Function(const gsl_vector *peak_params, void *calc_data, gsl_vector *func)
{
        size_t i;
        size_t iNbPoints        = ((struct FittingData *)calc_data)->n;
        double *x                       = ((struct FittingData *)calc_data)->x;
        double *y                       = ((struct FittingData *)calc_data)->y;
        double *sigma           = ((struct FittingData *)calc_data)->sigma;

        double pos                      = gsl_vector_get (peak_params, 0);
        double width            = gsl_vector_get (peak_params, 1);
        double height           = gsl_vector_get (peak_params, 2);
        double background       = gsl_vector_get (peak_params, 3);

//#ifdef DEBUG_FIT      
//  printf("pos = %g, width = %g, height = %g, background = %g\n", pos, width, height,background);
//#endif        

double dXPos,dXPos2,dWidth2,dYi;

        for (i=0; i < iNbPoints; i++)
        {       
                dXPos   = x[i] - pos;
                dXPos2  = dXPos * dXPos;
                dWidth2 = width * width;
        
                dYi = background + height * exp(-0.5 * dXPos2 / dWidth2);
        
                //gsl_vector_set (func, i,( dYi - y[i]) / sigma[i]);
                gsl_vector_set (func, i,( dYi - y[i]) / sigma[i]);                              
    }

  return GSL_SUCCESS;
}

//##ModelId=43849A0203B4
int FittingFunctionGaussianWithBackground::FunctionDerivate(const gsl_vector *peak_params, void *calc_data, gsl_matrix *Jacobian)
{

        size_t i;
        
        size_t  n                       = ((struct FittingData *)calc_data)->n;
        double *x                       = ((struct FittingData *)calc_data)->x;
        double *sigma           = ((struct FittingData *)calc_data)->sigma;

        double pos                      = gsl_vector_get (peak_params, 0);
        double width            = gsl_vector_get (peak_params, 1);
        double height           = gsl_vector_get (peak_params, 2);
        //double background     = gsl_vector_get (peak_params, 3);

        double d_pos, d_width, d_height, d_background;
        double s,xpos2,width2,e;
        // Fill the Jacobian Matrix
        for (i=0; i< n; i++)
    {
                s                               = sigma[i];
                xpos2                   = (x[i]-pos)*(x[i]-pos);
                width2                  = width * width;
                e                               = exp( -0.5 * xpos2 / width2);          
                 
                d_pos                   = (height/width2 * (x[i] - pos) * e) / s;       // derivate dYi/dpos    
                d_width                 = (height/width2/width * xpos2 * e) / s;        // derivate dYi/dwidth 
                d_height                = e/s;                                                                          // derivate dYi/dheight                  
                d_background    = 1/s;                                                                          // derivate dYi/dbackground

                // fill the jacobian matrix J(points, parameters) 
        gsl_matrix_set (Jacobian, i, 0, d_pos); 
                gsl_matrix_set (Jacobian, i, 1, d_width);
                gsl_matrix_set (Jacobian, i, 2, d_height);
                gsl_matrix_set (Jacobian, i, 3, d_background);
    }
         
        return GSL_SUCCESS;
}


// Interface function for the gsl_multifit_fdfsolver
//##ModelId=43849A0203B9
int FittingFunctionGaussianWithBackground::FunctionAndDerivate(const gsl_vector *peak_params, void *calc_data,gsl_vector *func, gsl_matrix *Jacobian) 
{       
        Function  (peak_params, calc_data, func);
        FunctionDerivate (peak_params, calc_data, Jacobian);
        return GSL_SUCCESS;
}

//The size of the 3 gsl_vector must be the same ! Not checked here !!!
//##ModelId=43849A0203C4
int FittingFunctionGaussianWithBackground::doFit(long points, gsl_vector *x, gsl_vector *y,gsl_vector *sigma,gsl_vector *guess,  bool bEstimateGuess,int iSearchStoppingMethod,bool bUseScaled)
{

initializeInitialsParameters(x,y,guess,bEstimateGuess);

        long            i;
        long            iter;
        long            status;
        
        double* datax = new double[points];
        double* datay = new double[points];
        double* datas = new double[points];
                
        // fill the data arrays 
        for (i=0; i<points; i++)
        {               
                datax[i] = gsl_vector_get (x, i);               // x = index            
                datay[i] = gsl_vector_get (y, i);               // y = data to be fit           
                datas[i] = gsl_vector_get (sigma, i);   // s = sigma (weight) of data values 
        }
        
        //gsl data structure    
        struct FittingData gsl_data;
        gsl_data.n      = points;
        gsl_data.x      = datax;
        gsl_data.y              = datay;
        gsl_data.sigma  = datas;

        cout << "\n---------------------------------------------" << endl;
        cout << "Initials Parameters For The Fit" << endl;
        // set-up the solver 
        const gsl_multifit_fdfsolver_type *T;
        if (bUseScaled)
        {
                T       = gsl_multifit_fdfsolver_lmsder;
                cout << "Jacobian Data Scaled" << endl;
        }
        else
        {
                T       = gsl_multifit_fdfsolver_lmder;
                cout << "Jacobian Data Not Scaled" << endl;
        }
        gsl_multifit_fdfsolver *s                               = gsl_multifit_fdfsolver_alloc (T, points, _iNbParameters);
        gsl_multifit_function_fdf f;

        f.f             = getFunction();              // gauss_f;
        f.df            = getFunctionDerivate();      // gauss_df;
        f.fdf           = getFunctionAndDerivate();   // gauss_fdf;
        f.n             = points;
        f.p             = _iNbParameters;
        f.params        = &gsl_data;
  
        gsl_multifit_fdfsolver_set ( s,  &f,  _vInitialGuess );

        switch(iSearchStoppingMethod)
        {
                case 1 :
                        cout << "Test Delta Search Stopping Function Used" << endl;
                        break;

                case 2 :
                        cout << "Test Gradient Search Stopping Function Used" << endl;
                        break;
        }

        cout << "\n----------------------------------------------" << endl;
        cout << "Initials Values For The Fit " << endl;
        cout << "Position   = " << gsl_vector_get (_vInitialGuess, 0) << endl;
        cout << "Width      = " << gsl_vector_get (_vInitialGuess ,1) << endl;
        cout << "Height     = " << gsl_vector_get (_vInitialGuess, 2) << endl;
        cout << "Background = " << gsl_vector_get (_vInitialGuess, 3) << endl;


                        cout << "\n---------------------------------------------" << endl;
                        cout << "Iterations" << endl;
        iter = 0;
        #ifdef DEBUG_FIT
          printState (iter, s);
        #endif

        do
    {
                iter++;
                status = gsl_multifit_fdfsolver_iterate (s);

                #ifdef DEBUG_FIT
                          printf ("status = %s\n", gsl_strerror (status));
                          printState (iter, s);
                #endif
                
                if (status)
                {
                break;
                }


                //
                switch(iSearchStoppingMethod)
                {
                case 1 :
                        status = gsl_multifit_test_delta (s->dx, s->x,  /*_dSearchLimit*/0.0, _dSearchLimit);
                        break;

                case 2 :
                        gsl_vector* gradient = gsl_vector_alloc(_iNbParameters);
                        gsl_multifit_gradient(s->J,s->f,gradient);
                        status = gsl_multifit_test_gradient (gradient, _dSearchLimit);
                        gsl_vector_free(gradient);
                        break;
                }

    }
        while (status == GSL_CONTINUE && iter < _lNbMaxIteration);
        _lNbIterations = iter;
        // Check whether the fitting has converged or failed  
        //if ( status == GSL_SUCCESS )
        //      {
        gsl_matrix *covar = gsl_matrix_alloc (_iNbParameters, _iNbParameters);
        gsl_multifit_covar (s->J, 0.0, covar);
        double chi                      = gsl_blas_dnrm2(s->f);

        _pParameters[0]                 =       FIT(0);
        _pParameters[1]                 =       FIT(1);
        _pParameters[2]                 =       FIT(2);
        _pParameters[3]                 =       FIT(3);

        _pParametersErrors[0]   =       ERR(0);
        _pParametersErrors[1]   =       ERR(1);
        _pParametersErrors[2]   =       ERR(2);
        _pParametersErrors[3]   =       ERR(3);

        //double dHeight = _pParameters[2];

        double dEstimatedWidth  =       _pParameters[1];
        _dQualityFactor                 =       pow(chi, 2.0)/ (points - _iNbParameters);
//      _dSigma                                 =       1.0/ (sqrt(2*M_PI)*dHeight);
//      _dFHWM                                  =       2.355*_dSigma;
        _dHWHM                                  =       sqrt(2*log(2.0))*dEstimatedWidth;               
        _dFHWM                                  =       2*_dHWHM;               
        cout << "2*sqrt(2*log(2.0))=" << 2*sqrt(2*log(2.0)) <<endl;

        #ifdef DEBUG_FIT        
                        cout << "\n---------------------------------------------" << endl;
                        cout << "Results For Gaussian With Background" << endl;
                        printResults();
        #endif          

        gsl_matrix_free(covar);         
        //      }

        #ifdef DEBUG_FIT  
                printf ("status = %s\n", gsl_strerror (status));
        #endif
        
        // Free data
        gsl_multifit_fdfsolver_free (s);
        delete[] datax;
        delete[] datay;
        delete[] datas;
        
        datax = NULL;
        datax = NULL;
        datax = NULL;

        return status;
}



void    FittingFunctionGaussianWithBackground::estimateInitialGuess(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY)
{
        // fill the initials parameters values estimation vector for the initialisation
        // [0] = position 
        // [1] = width  
        // [2] = height
        // [3] = background
        // Thanks to the max index of all the points, we can get the
        // position_value       = x_value at index_max
        // width                        = standard deviation
        // height_value         = y_value - min_value (=background)     at index_max
        // background_value = minimal value
        //min_index = gsl_vector_min_index (vExperimentalDataY);
        //max_index = gsl_vector_max_index (vExperimentalDataY);
        double dInitialPosition = estimateInitialPosition(vExperimentalDataX,vExperimentalDataY);
        double dInitialSigma    = estimateInitialSigma(vExperimentalDataX,vExperimentalDataY);
        double dInitialHeight   = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
        double dInitialBackground = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);

        gsl_vector_set(_vInitialGuess, 0, dInitialPosition);
        gsl_vector_set(_vInitialGuess, 1, dInitialSigma);
        gsl_vector_set(_vInitialGuess, 2, dInitialHeight);
        gsl_vector_set(_vInitialGuess, 3, dInitialBackground);  
}


void    FittingFunctionGaussianWithBackground::initializeInitialsParameters(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY,gsl_vector* vInitialGuess,bool bEstimate)
{
        if (bEstimate) 
        {
                estimateInitialGuess(vExperimentalDataX,vExperimentalDataY);
        }
        else
        {
                gsl_vector_set(_vInitialGuess,0,gsl_vector_get(vInitialGuess,0));
                gsl_vector_set(_vInitialGuess,1,gsl_vector_get(vInitialGuess,1));
                gsl_vector_set(_vInitialGuess,2,gsl_vector_get(vInitialGuess,2));
                gsl_vector_set(_vInitialGuess,3,gsl_vector_get(vInitialGuess,3));
        }
}





//##ModelId=43849A0203D0
int FittingFunctionGaussianWithBackground::generateFunctionFit(double dValMinX,double dResolutionX,long lNbPoints,gsl_vector* funcX,gsl_vector* funcY) {

        int i;
        double dPos                     =       _pParameters[0];
        double dWidth           =       _pParameters[1];
        double dHeight          =       _pParameters[2];
        double dBackground      =       _pParameters[3];
        double dXi                      =       0;
        double dXPos            =       0;
        double dXPos2           =       0;
        double dWidth2          =       1;
        double dYi                      =       0;

        for (i=0;i<lNbPoints;i++) 
        {
                dXi             = dValMinX+i*dResolutionX;
                dXPos   = dXi - dPos;
                dXPos2  = dXPos * dXPos;
                dWidth2 = dWidth * dWidth;
                dYi             = dBackground + dHeight * exp(-0.5 * dXPos2 / dWidth2);
        
                gsl_vector_set(funcX,i,dXi);
                gsl_vector_set(funcY,i,dYi);
        }
        return 1;
}


int FittingFunctionGaussianWithBackground::generateFunctionFit(long lNbPoints, gsl_vector *x, gsl_vector* funcX,gsl_vector* funcY)
{

        int i;
        double dPos                     =       _pParameters[0];
        double dWidth           =       _pParameters[1];
        double dHeight          =       _pParameters[2];
        double dBackground      =       _pParameters[3];
        double dXi                      =       0;
        double dXPos            =       0;
        double dXPos2           =       0;
        double dWidth2          =       1;
        double dYi                      =       0;

        for (i=0;i<lNbPoints;i++) 
        {
                dXi             = gsl_vector_get(x,i);
                dXPos   = dXi - dPos;
                dXPos2  = dXPos * dXPos;
                dWidth2 = dWidth * dWidth;
                dYi             = dBackground + dHeight * exp(-0.5 * dXPos2 / dWidth2);
        
                gsl_vector_set(funcX,i,dXi);
                gsl_vector_set(funcY,i,dYi);
        }
        return 1;
}



// print the state of the fitting during iterations 
//##ModelId=43849A0203D7
void FittingFunctionGaussianWithBackground::printState (size_t iter, gsl_multifit_fdfsolver * s)
{
  printf ("iter: %3u x : %g %g %g %g|f(x)|=%g\n",
          iter,
          gsl_vector_get (s->x, 0), 
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->x, 2),
                  gsl_vector_get (s->x, 3),
          gsl_blas_dnrm2 (s->f));
}

//##ModelId=43849A0203E0
void FittingFunctionGaussianWithBackground::printResults()
{
                printf ("pos            = %g \t+/- %g\n",       getParameter("Position")        ,_pParametersErrors[0]);
                printf ("width          = %g \t+/- %g\n",       getParameter("Width")           ,_pParametersErrors[1]);
                printf ("height         = %g \t+/- %g\n",       getParameter("Height")          ,_pParametersErrors[2]);
                printf ("background     = %g \t+/- %g\n",       getParameter("Background")      ,_pParametersErrors[3]);
        printf ("chisq/dof      = %g\n",                                _dQualityFactor);
                printf ("Sigma          = %g\n",                                _dSigma);
                printf ("HWHM           = %g\n",                                _dHWHM);
                printf ("FHWM           = %g\n",                                _dFHWM);
}


double  FittingFunctionGaussianWithBackground::computeValue(double dX,double dPosition, double dWidth, double dHeight, double dBackground)
{
        double dXPos,dXPos2,dWidth2,dYi;

        dXPos   = dX - dPosition;
        dXPos2  = dXPos * dXPos;
        dWidth2 = dWidth * dWidth;
        
        dYi = dBackground + dHeight * exp(-0.5 * dXPos2 / dWidth2);
        return dYi;
}

double  FittingFunctionGaussianWithBackground::computeDerivateValue(double dX,double dPosition, double dWidth, double dHeight)
{
        return 0.0;
}

double FittingFunctionGaussianWithBackground::estimateInitialPosition(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
{
        size_t max_index        = gsl_vector_max_index (vExperimentalDataY);
        double dPosition        = gsl_vector_get(vExperimentalDataX,max_index);

        return dPosition;
}




double FittingFunctionGaussianWithBackground::estimateInitialHeight(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
{
        size_t max_index        = gsl_vector_max_index (vExperimentalDataY);

        double background       = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
        double maxValue         = gsl_vector_get(vExperimentalDataY, max_index);
        double height           = maxValue - background;
        return height;
}


double FittingFunctionGaussianWithBackground::estimateInitialBackground(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
{
        double dBackground = 0.0;
        size_t min_index = gsl_vector_min_index (vExperimentalDataY);
        //TODO do a mean with more points to better estimate background
        dBackground =  gsl_vector_get(vExperimentalDataY, min_index);

        return dBackground;
}


double FittingFunctionGaussianWithBackground::estimateInitialSigma(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
{
        double dSigma;
        double dCovariance;
        double dHeight                          = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
        double dBackground                      = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
        
        size_t iIndexMax                        = gsl_vector_max_index(vExperimentalDataY);
        double dYMaxValue                       = gsl_vector_get(vExperimentalDataY, iIndexMax);
  
        double  dYCurrentValue          = dYMaxValue;

        size_t iCurrentIndex            = iIndexMax;
        size_t lVectorSize                      = vExperimentalDataX->size;//lSizeX;
        double lHalfVectorSize          = 0.5 * lVectorSize;
        double dHalfMaximum                     = 0.5 * dHeight + dBackground;
  
        size_t iIndexLow                        = iCurrentIndex;
        size_t iIndexHigh                       = iCurrentIndex;

    if (iIndexMax < lHalfVectorSize)
    {
//              cout << "\n\t--> Look in increasing direction" << endl;
                //- look in increasing direction
                while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue > dHalfMaximum)
                {
                        dYCurrentValue  = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
                        iIndexHigh              = iCurrentIndex;
                        iIndexLow               = iIndexMax - iCurrentIndex;
                }
    }
    else
    {
//              cout << "\n\t--> Look in decreasing direction" << endl;
                //- look in decreasing direction
                while( iCurrentIndex > 0 && dYCurrentValue > dHalfMaximum)
                {
                        dYCurrentValue =  gsl_vector_get(vExperimentalDataY, --iCurrentIndex);
                        iIndexLow               = iCurrentIndex;
                        iIndexHigh              = iIndexMax + iCurrentIndex;
                }
    }

    if (iCurrentIndex != iIndexMax)
    {
                double dPosition        = gsl_vector_get(vExperimentalDataX, iIndexMax);
                double dX2                      = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
                dSigma = ::fabs(dX2 - dPosition);
    }
    else
    {
                //- try an arbitrary value that should be not so far from reality
                dSigma = lVectorSize / 6;
    }

/*      cout << "----------------------------" << endl;
        cout << "Peak : (x,y) --> ("    << gsl_vector_get(vExperimentalDataX, iIndexMax) << "," 
                                                                        << gsl_vector_get(vExperimentalDataY, iIndexMax) << ")" << endl;
                
        
        cout << "Low : index=" << iIndexLow << endl;
        cout << "\txlow=" << gsl_vector_get(vExperimentalDataX, iIndexLow) << " | ylow=" << gsl_vector_get(vExperimentalDataY, iIndexLow) << endl;
        cout << "\txhigh=" << gsl_vector_get(vExperimentalDataX, iIndexHigh) << " | yhigh=" << gsl_vector_get(vExperimentalDataY, iIndexHigh) << endl;
        
        cout << "\txhigh - xlow= " << gsl_vector_get(vExperimentalDataX, iIndexHigh) - gsl_vector_get(vExperimentalDataX, iIndexLow) << endl;
        cout << "\tsigma=" << dSigma << endl;
*/
    dCovariance = dSigma * dSigma;
        return dSigma;
}

---