FittingFunctionLorentzian.cpp

Go to the documentation of this file.
00001 // FittingFunctionLorentzian.cpp: implementation of the FittingFunctionLorentzian class.
00002 //
00004 #ifdef WIN32 
00005 #pragma warning(disable:4786)
00006 #endif
00007 
00008 #include "FittingFunctionLorentzian.h"
00009 
00010 
00012 // Construction/Destruction
00014 
00015 std::string     FittingFunctionLorentzian::getEquation() const
00016 {
00017         return "f(x) = Height / [1 + (x - Position)˛ / Width˛]";
00018 }
00019 
00020 
00021 //##ModelId=43849A020334
00022 FittingFunctionLorentzian::FittingFunctionLorentzian()
00023 {       
00024         _iNbParameters                          =       3;
00025         _lNbMaxIteration                        =       NB_ITERATION_MAX;
00026         _dSearchLimit                           =       SEARCH_LIMIT;
00027         _pParameters                            =       new double[_iNbParameters];
00028         _pParametersErrors                      =       new double[_iNbParameters];
00029         _mFunctionMap["Position"]       =       0;
00030         _mFunctionMap["Width"]          =       1;
00031         _mFunctionMap["Height"]         =       2;
00032         _vInitialGuess = gsl_vector_alloc(3);
00033 }
00034 
00035 //##ModelId=43849A020335
00036 FittingFunctionLorentzian::FittingFunctionLorentzian(long lNbIteration,double dLimitSearch) 
00037 {
00038         _lNbMaxIteration                        =       lNbIteration;
00039         _dSearchLimit                           =       dLimitSearch;
00040         _iNbParameters                          =       3;
00041         _pParameters                            =       new double[_iNbParameters];
00042         _pParametersErrors                      =       new double[_iNbParameters];
00043         _mFunctionMap["Position"]       =       0;
00044         _mFunctionMap["Width"]          =       1;
00045         _mFunctionMap["Height"]         =       2;
00046         _vInitialGuess = gsl_vector_alloc(3);
00047 }
00048 
00049 //##ModelId=43849A020342
00050 FittingFunctionLorentzian::~FittingFunctionLorentzian()
00051 {
00052         printf ("Destruction FittingFunctionLorentzian !!!\n");
00053         delete[] _pParameters;
00054         delete[] _pParametersErrors;
00055         _pParameters = NULL;
00056         _pParametersErrors = NULL;
00057         gsl_vector_free(_vInitialGuess);
00058 }
00059 
00060 //##ModelId=43849A020344
00061 PtrFunction FittingFunctionLorentzian::getFunction() const{
00062         return FittingFunctionLorentzian::Function;
00063 }
00064 
00065 //##ModelId=43849A020346
00066 PtrFunctionDerivate FittingFunctionLorentzian::getFunctionDerivate() const{
00067         return FittingFunctionLorentzian::FunctionDerivate;
00068 }
00069 
00070 //##ModelId=43849A020348
00071 PtrFunctionAndDerivate FittingFunctionLorentzian::getFunctionAndDerivate() const {
00072         return FittingFunctionLorentzian::FunctionAndDerivate;
00073 }
00074 
00075 
00080 //##ModelId=43849A02034A
00081 int FittingFunctionLorentzian::Function(const gsl_vector *peak_params, void *calc_data, gsl_vector *func)
00082 {
00083         size_t i;
00084         size_t iNbPoints = ((struct FittingData *)calc_data)->n;
00085         double *x = ((struct FittingData *)calc_data)->x;
00086         double *y = ((struct FittingData *)calc_data)->y;
00087         double *sigma = ((struct FittingData *) calc_data)->sigma;
00088 
00089         double pos              = gsl_vector_get (peak_params, 0);
00090         double width    = gsl_vector_get (peak_params, 1);
00091         double height   = gsl_vector_get (peak_params, 2);
00092 
00093 //#ifdef DEBUG_FIT      
00094 //  printf("pos = %g, width = %g, height = %g\n", pos, width, height);
00095 //#endif        
00096 
00097         double dXPos,dXPos2,dWidth2,dYi;
00098 
00099   for (i=0; i < iNbPoints; i++)
00100   {             
00101         dXPos   = x[i] - pos;
00102         dXPos2  = dXPos * dXPos;
00103         dWidth2 = width * width;
00104                         
00105         dYi = height / (1 + dXPos2 / dWidth2);
00106 
00107         gsl_vector_set (func, i,( dYi - y[i]) / sigma[i]);              
00108   }
00109 
00110   return GSL_SUCCESS;
00111 }
00112 
00113 
00114 
00125 //##ModelId=43849A020352
00126 int FittingFunctionLorentzian::FunctionDerivate(const gsl_vector *peak_params, void *calc_data, gsl_matrix *Jacobian)
00127 {
00128 
00129         size_t i;
00130         
00131         size_t  n = ((struct FittingData *)calc_data)->n;
00132         double *x = ((struct FittingData *)calc_data)->x;
00133         double *sigma = ((struct FittingData *)calc_data)->sigma;
00134 
00135         double pos              = gsl_vector_get (peak_params, 0);
00136         double width    = gsl_vector_get (peak_params, 1);
00137         double height   = gsl_vector_get (peak_params, 2);
00138   
00139         double d_pos, d_width, d_height;
00140 
00141         // Fill the Jacobian Matrix     
00142         // Jacobian matrix J(i,j) = dfi / dxj, 
00143         // where fi = (Yi - yi)/sigma[i],  
00144         // Yi = Yi = c / [ 1 + ((x-a)/b)^2] 
00145         // and the xj are the parameters (a,b,c) 
00146 
00147         double s,alpha,beta,e,e2;
00148                 for (i=0; i< n; i++)
00149         {
00150 
00151                 s                               = sigma[i]; 
00152                 alpha                   = (x[i] - pos)/width;
00153                 beta                    = 2 * height / width;
00154                 e                               = (1+alpha*alpha);
00155                 e2                              = e * e;
00156 
00157                 d_pos                   = beta * alpha / e2 / s;                        // derivate dYi/dpos
00158                 d_width                 = beta * alpha * alpha / e2 / s;        // derivate dYi/dwidth
00159                 d_height                = 1/e/s;                                                        // derivate dYi/dheight
00160                 
00161                 // fill the jacobian matrix J(points, parameters) 
00162                 gsl_matrix_set (Jacobian, i, 0, d_pos); 
00163                 gsl_matrix_set (Jacobian, i, 1, d_width);
00164                 gsl_matrix_set (Jacobian, i, 2, d_height);
00165         }
00166          
00167         return GSL_SUCCESS;
00168 }
00169 
00170 // Interface function for the gsl_multifit_fdfsolver
00171 //##ModelId=43849A020357
00172 int FittingFunctionLorentzian::FunctionAndDerivate(const gsl_vector *peak_params, void *calc_data,gsl_vector *func, gsl_matrix *Jacobian) 
00173 {       
00174         Function  (peak_params, calc_data, func);
00175         FunctionDerivate (peak_params, calc_data, Jacobian);
00176         return GSL_SUCCESS;
00177 }
00178 
00179 
00180 
00181 //##ModelId=43849A020362
00182 int FittingFunctionLorentzian::doFit(long points, gsl_vector *x, gsl_vector *y,gsl_vector *sigma ,gsl_vector* guess, bool bEstimateGuess, int iSearchStoppingMethod,bool bUseScaled)
00183 {
00184 initializeInitialsParameters(x,y,guess,bEstimateGuess);
00185         long            i;
00186         long            iter;
00187         long            status;
00188         
00189         double* datax = new double[points];
00190         double* datay = new double[points];
00191         double* datas = new double[points];
00192                         
00193         // fill the data arrays 
00194         for (i=0; i<points; i++)
00195         {               
00196                 datax[i] = gsl_vector_get (x, i);               // x = index            
00197                 datay[i] = gsl_vector_get (y, i);               // y = data to be fit    
00198                 datas[i] = gsl_vector_get (sigma, i);   // s = sigma (weight) of data values
00199         }
00200         
00201         // gsl data structure   
00202         struct FittingData gsl_data;
00203         gsl_data.n      = points;
00204         gsl_data.x      = datax;
00205         gsl_data.y              = datay;
00206         gsl_data.sigma = datas;
00207 
00208 //      #ifdef DEBUG_FIT                
00209 //              printf ("Init = %g \t %g\t %g\t \n", gsl_vector_get(_vInitialGuess,0),gsl_vector_get(_vInitialGuess,1),gsl_vector_get(_vInitialGuess,2));
00210 //      #endif
00211 
00212                         cout << "\n---------------------------------------------" << endl;
00213                         cout << "Initials Parameters For The Fit" << endl;
00214         // set-up the solver 
00215         const gsl_multifit_fdfsolver_type *T;
00216         if (bUseScaled)
00217         {
00218                 T       = gsl_multifit_fdfsolver_lmsder;
00219                 cout << "Jacobian Data Scaled" << endl;
00220         }
00221         else
00222         {
00223                 T       = gsl_multifit_fdfsolver_lmder;
00224                         cout << "Jacobian Data Not Scaled" << endl;
00225         }
00226         gsl_multifit_fdfsolver *s                               = gsl_multifit_fdfsolver_alloc (T, points, _iNbParameters);
00227         gsl_multifit_function_fdf f;
00228 
00229         f.f             = getFunction();            // lorentz_f;
00230         f.df            = getFunctionDerivate();    // lorentz_df;
00231         f.fdf           = getFunctionAndDerivate(); // lorentz_fdf;
00232         f.n             = points;
00233         f.p             = _iNbParameters;
00234         f.params        = &gsl_data;
00235   
00236         switch(iSearchStoppingMethod)
00237         {
00238                 case 1 :
00239                         cout << "Test Delta Search Stopping Function Used" << endl;
00240                         break;
00241 
00242                 case 2 :
00243                         cout << "Test Gradient Search Stopping Function Used" << endl;
00244                         break;
00245         }
00246 
00247         cout << "\n----------------------------------------------" << endl;
00248         cout << "Initials Values For The Fit " << endl;
00249         cout << "Position = " << gsl_vector_get (_vInitialGuess, 0) << endl;
00250         cout << "Width    = " << gsl_vector_get (_vInitialGuess, 1) << endl;
00251         cout << "Height   = " << gsl_vector_get (_vInitialGuess, 2) << endl;
00252         
00253 
00254         gsl_multifit_fdfsolver_set ( s,  &f,  _vInitialGuess );
00255 
00256                         cout << "\n---------------------------------------------" << endl;
00257                         cout << "Iterations" << endl;
00258         iter = 0;
00259   
00260         #ifdef DEBUG_FIT
00261           printState (iter, s);
00262         #endif
00263 
00264         do
00265     {
00266                 iter++;
00267                 status = gsl_multifit_fdfsolver_iterate (s);
00268 
00269                 #ifdef DEBUG_FIT
00270                           printf ("status = %s\n", gsl_strerror (status));
00271                           printState (iter, s);
00272                 #endif
00273                 
00274                 if (status)
00275                 {
00276                 break;
00277                 }
00278 
00279 
00280                 switch(iSearchStoppingMethod)
00281                 {
00282                 case 1 :
00283                         status = gsl_multifit_test_delta (s->dx, s->x,  /*_dSearchLimit*/0.0, _dSearchLimit);
00284                         break;
00285 
00286                 case 2 :
00287                         gsl_vector* gradient = gsl_vector_alloc(_iNbParameters);
00288                         gsl_multifit_gradient(s->J,s->f,gradient);
00289                         status = gsl_multifit_test_gradient (gradient, _dSearchLimit);
00290                         gsl_vector_free(gradient);
00291                         break;
00292                 }
00293     }
00294         while (status == GSL_CONTINUE && iter < _lNbMaxIteration);
00295 
00296         _lNbIterations = iter;
00297         // Check whether the fitting has converged or failed  
00298         //if ( status == GSL_SUCCESS )
00299         //      {
00300         gsl_matrix *covar = gsl_matrix_alloc (_iNbParameters, _iNbParameters);
00301         gsl_multifit_covar (s->J, 0.0, covar);
00302         double chi                      = gsl_blas_dnrm2(s->f);
00303 
00304         _pParameters[0] =       FIT(0);
00305         _pParameters[1] =       FIT(1);
00306         _pParameters[2] =       FIT(2);
00307 
00308         _pParametersErrors[0]   =       ERR(0);
00309         _pParametersErrors[1]   =       ERR(1);
00310         _pParametersErrors[2]   =       ERR(2);
00311 
00312         //double dHeight = _pParameters[2];
00313         double dEstimatedWidth = _pParameters[1];
00314         _dQualityFactor =       pow(chi, 2.0)/ (points - _iNbParameters);
00315 //      _dSigma                 =       1.0/ (sqrt(2*M_PI)*dHeight);
00316 //      _dFHWM                  =       2.355*_dSigma;
00317                 
00318         _dHWHM                  = dEstimatedWidth;
00319         _dFHWM                  = 2*_dHWHM;
00320                 
00321         #ifdef DEBUG_FIT        
00322                         cout << "\n---------------------------------------------" << endl;
00323                         cout << "Results For Lorentzian" << endl;
00324                         printResults(); 
00325         #endif  
00326         
00327         gsl_matrix_free(covar);         
00328         //      }
00329 
00330         #ifdef DEBUG_FIT  
00331                 printf ("status = %s\n", gsl_strerror (status));
00332         #endif
00333         
00334         // Free data 
00335         gsl_multifit_fdfsolver_free (s);
00336 
00337         delete[] datax;
00338         delete[] datay;
00339         delete[] datas;
00340         
00341         datax = NULL;
00342         datax = NULL;
00343         datax = NULL;
00344 
00345         return status;
00346 }
00347 
00348 
00349 //##ModelId=43849A02036A
00350 int FittingFunctionLorentzian::generateFunctionFit(double dValMinX,double dResolutionX,long lNbPoints,gsl_vector* funcX,gsl_vector* funcY){
00351 
00352         int i;
00353 
00354         double dPos             =       _pParameters[0];
00355         double dWidth   =       _pParameters[1];
00356         double dHeight  =       _pParameters[2];
00357 
00358         double dXi,dXPos,dXPos2,dWidth2,dYi;
00359 
00360         for (i=0; i < lNbPoints; i++)
00361     {           
00362                 dXi             = dValMinX + i*dResolutionX;
00363                 dXPos   = dXi  - dPos;
00364                 dXPos2  = dXPos * dXPos;
00365                 dWidth2 = dWidth * dWidth;
00366                                 
00367                 dYi             = dHeight / (1.0 + dXPos2 / dWidth2);
00368 
00369                 gsl_vector_set (funcX, i, dXi);
00370                 gsl_vector_set (funcY, i, dYi);         
00371     }   
00372         return 1;
00373 }
00374 
00375 
00376 
00377 int FittingFunctionLorentzian::generateFunctionFit(long lNbPoints, gsl_vector *x, gsl_vector* funcX,gsl_vector* funcY)
00378 {
00379 
00380         int i;
00381 
00382         double dPos             =       _pParameters[0];
00383         double dWidth   =       _pParameters[1];
00384         double dHeight  =       _pParameters[2];
00385 
00386         double dXi,dXPos,dXPos2,dWidth2,dYi;
00387 
00388         for (i=0; i < lNbPoints; i++)
00389     {           
00390                 dXi             = gsl_vector_get(x,i);
00391                 dXPos   = dXi  - dPos;
00392                 dXPos2  = dXPos * dXPos;
00393                 dWidth2 = dWidth * dWidth;
00394                                 
00395                 dYi             = dHeight / (1.0 + dXPos2 / dWidth2);
00396 
00397                 gsl_vector_set (funcX, i, dXi);
00398                 gsl_vector_set (funcY, i, dYi);         
00399     }   
00400         return 1;
00401 }
00402 
00403 
00404 
00405 
00406 
00407 // print the state of the fitting during iterations 
00408 //##ModelId=43849A020376
00409 void FittingFunctionLorentzian::printState (size_t iter, gsl_multifit_fdfsolver * s)
00410 {
00411   printf ("iter: %3u x : %g %g %g |f(x)|=%g\n",
00412           iter,
00413           gsl_vector_get (s->x, 0), 
00414           gsl_vector_get (s->x, 1),
00415           gsl_vector_get (s->x, 2), 
00416           gsl_blas_dnrm2 (s->f));
00417 }
00418 
00419 //##ModelId=43849A020379
00420 void FittingFunctionLorentzian::printResults()
00421 {
00422                 printf ("pos            = %g \t+/- %g\n",       getParameter("Position"),_pParametersErrors[0]);
00423                 printf ("width          = %g \t+/- %g\n",       getParameter("Width"),  _pParametersErrors[1]);
00424                 printf ("height         = %g \t+/- %g\n",       getParameter("Height"), _pParametersErrors[2]);
00425         printf ("chisq/dof      = %g\n",                                _dQualityFactor);
00426                 printf ("Sigma          = %g\n",                                _dSigma);
00427                 printf ("HWHM           = %g\n",                                _dHWHM);
00428                 printf ("FHWM           = %g\n",                                _dFHWM);
00429 }
00430 
00431 
00432 double  FittingFunctionLorentzian::computeValue(double dX,double dPosition, double dWidth, double dHeight, double dBackground)
00433 {
00434 
00435         double dXPos,dXPos2,dWidth2,dYi;
00436 
00437         dXPos   = dX - dPosition;
00438         dXPos2  = dXPos * dXPos;
00439         dWidth2 = dWidth * dWidth;
00440                         
00441         dYi = dHeight / (1 + dXPos2 / dWidth2);
00442 
00443         return dYi;
00444 }
00445 
00447 double  FittingFunctionLorentzian::computeDerivateValue(double dX,double dPosition, double dWidth, double dHeight)
00448 {
00449         return 0.0;
00450 }
00451 
00452 
00453 double FittingFunctionLorentzian::estimateInitialPosition(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00454 {
00455         size_t max_index        = gsl_vector_max_index (vExperimentalDataY);
00456         double dPosition        = gsl_vector_get(vExperimentalDataX,max_index);
00457 
00458         return dPosition;
00459 }
00460 
00461 
00462 
00463 
00464 double FittingFunctionLorentzian::estimateInitialHeight(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00465 {
00466         size_t max_index        = gsl_vector_max_index (vExperimentalDataY);
00467 
00468         double background       = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00469         double maxValue         = gsl_vector_get(vExperimentalDataY, max_index);
00470         double height           = maxValue - background;
00471         return height;
00472 }
00473 
00474 
00475 double FittingFunctionLorentzian::estimateInitialBackground(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00476 {
00477         double dBackground = 0.0;
00478         return dBackground;
00479 }
00480 
00481 
00482 double FittingFunctionLorentzian::estimateInitialSigma(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00483 {
00484         double dSigma;
00485         double dCovariance;
00486         double dHeight                          = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00487         double dBackground                      = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00488         
00489         size_t iIndexMax                        = gsl_vector_max_index(vExperimentalDataY);
00490         double dYMaxValue                       = gsl_vector_get(vExperimentalDataY, iIndexMax);
00491   
00492         double  dYCurrentValue          = dYMaxValue;
00493 
00494         size_t iCurrentIndex            = iIndexMax;
00495         size_t lVectorSize                      = vExperimentalDataX->size;//lSizeX;
00496         double lHalfVectorSize          = 0.5 * lVectorSize;
00497         double dHalfMaximum                     = 0.5 * dHeight + dBackground;
00498   
00499         size_t iIndexLow                        = iCurrentIndex;
00500         size_t iIndexHigh                       = iCurrentIndex;
00501 
00502     if (iIndexMax < lHalfVectorSize)
00503     {
00504 //              cout << "\n\t--> Look in increasing direction" << endl;
00505                 //- look in increasing direction
00506                 while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue > dHalfMaximum)
00507                 {
00508                         dYCurrentValue  = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00509                         iIndexHigh              = iCurrentIndex;
00510                         iIndexLow               = iIndexMax - iCurrentIndex;
00511                 }
00512     }
00513     else
00514     {
00515 //              cout << "\n\t--> Look in decreasing direction" << endl;
00516                 //- look in decreasing direction
00517                 while( iCurrentIndex > 0 && dYCurrentValue > dHalfMaximum)
00518                 {
00519                         dYCurrentValue =  gsl_vector_get(vExperimentalDataY, --iCurrentIndex);
00520                         iIndexLow               = iCurrentIndex;
00521                         iIndexHigh              = iIndexMax + iCurrentIndex;
00522                 }
00523     }
00524 
00525     if (iCurrentIndex != iIndexMax)
00526     {
00527                 double dPosition        = gsl_vector_get(vExperimentalDataX, iIndexMax);
00528                 double dX2                      = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00529                 dSigma = ::fabs(dX2 - dPosition);
00530     }
00531     else
00532     {
00533                 //- try an arbitrary value that should be not so far from reality
00534                 dSigma = lVectorSize / 6;
00535     }
00536 
00537 /*      cout << "--------------------------------" << endl;
00538         cout << "Peak : (x,y) --> ("    << gsl_vector_get(vExperimentalDataX, iIndexMax) << "," 
00539                                                                         << gsl_vector_get(vExperimentalDataY, iIndexMax) << ")" << endl;
00540                 
00541         
00542         cout << "Low : index=" << iIndexLow << endl;
00543         cout << "\txlow=" << gsl_vector_get(vExperimentalDataX, iIndexLow) << " | ylow=" << gsl_vector_get(vExperimentalDataY, iIndexLow) << endl;
00544         cout << "\txhigh=" << gsl_vector_get(vExperimentalDataX, iIndexHigh) << " | yhigh=" << gsl_vector_get(vExperimentalDataY, iIndexHigh) << endl;
00545         
00546         cout << "\txhigh - xlow= " << gsl_vector_get(vExperimentalDataX, iIndexHigh) - gsl_vector_get(vExperimentalDataX, iIndexLow) << endl;
00547         cout << "\tsigma=" << dSigma << endl;
00548 */
00549     dCovariance = dSigma * dSigma;
00550         return dSigma;
00551 }
00552 
00553 
00554 void    FittingFunctionLorentzian::estimateInitialGuess(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY)
00555 {
00556         // fill the initials parameters values estimation vector for the initialisation
00557         // [0] = position 
00558         // [1] = width  
00559         // [2] = height
00560         // Thanks to the max index of all the points, we can get the
00561         // position_value = x_value at index_max
00562         // width                  = standard deviation
00563         // height_value   = y_value     at index_max
00564         //max_index = gsl_vector_max_index (vExperimentalDataY);
00565 
00566         double dInitialPosition = estimateInitialPosition(vExperimentalDataX,vExperimentalDataY);
00567         double dInitialSigma    = estimateInitialSigma(vExperimentalDataX,vExperimentalDataY);
00568         double dInitialHeight   = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00569 
00570         gsl_vector_set(_vInitialGuess, 0, dInitialPosition);
00571         gsl_vector_set(_vInitialGuess, 1, dInitialSigma);
00572         gsl_vector_set(_vInitialGuess, 2, dInitialHeight);
00573 }
00574 
00575 
00576 void    FittingFunctionLorentzian::initializeInitialsParameters(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY,gsl_vector* vInitialGuess,bool bEstimate)
00577 {
00578         if      (bEstimate)
00579         {
00580                 estimateInitialGuess(vExperimentalDataX,vExperimentalDataY);
00581         }
00582         else
00583         {
00584                 gsl_vector_set(_vInitialGuess,0,gsl_vector_get(vInitialGuess,0));
00585                 gsl_vector_set(_vInitialGuess,1,gsl_vector_get(vInitialGuess,1));
00586                 gsl_vector_set(_vInitialGuess,2,gsl_vector_get(vInitialGuess,2));
00587         }
00588 }

Generated on Tue Apr 14 09:34:28 2009 for Data Fitting Library by  doxygen 1.4.5