FittingFunctionSigmoid Class Reference

#include <FittingFunctionSigmoid.h>

Inherits FittingFunction.

Inheritance diagram for FittingFunctionSigmoid:

Inheritance graph
[legend]
Collaboration diagram for FittingFunctionSigmoid:

Collaboration graph
[legend]
List of all members.

Public Member Functions

 FittingFunctionSigmoid ()
 FittingFunctionSigmoid (long lNbIteration, double dLimitSearch)
virtual ~FittingFunctionSigmoid ()
virtual PtrFunction getFunction () const
virtual PtrFunctionDerivate getFunctionDerivate () const
virtual PtrFunctionAndDerivate getFunctionAndDerivate () const
virtual int doFit (long points, gsl_vector *x, gsl_vector *y, gsl_vector *sigma, gsl_vector *guess, bool bEstimateGuess, int iSearchStoppingMethod=1, bool bUseScaled=true)
virtual int generateFunctionFit (double dValMinX, double dResolutionX, long lNbPoints, gsl_vector *funcX, gsl_vector *funcY)
virtual int generateFunctionFit (long lNbPoints, gsl_vector *x, gsl_vector *funcX, gsl_vector *funcY)
void printState (size_t iter, gsl_multifit_fdfsolver *s)
void printResults ()
double computeValue (double dX, double dPosition, double dWidth, double dHeight, double dBackground)
double computeDerivateValue (double dX, double dPosition, double dWidth, double dHeight)
void initializeInitialsParameters (gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY, gsl_vector *vInitialGuess, bool bEstimate)
void estimateInitialGuess (gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
double estimateInitialPosition (gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
double estimateInitialSigma (gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
double estimateInitialHeight (gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
double estimateInitialBackground (gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
void estimateSigmoidIntersections ()
double getXLow () const
double getXHigh () const
std::string getEquation () const

Static Public Member Functions

static int Function (const gsl_vector *peak_params, void *calc_data, gsl_vector *func)
 Implements the Sigmoid function for the model :

\[ Y_i = \frac{height}{1 + \exp \left[ - \frac{x_i - pos}{width} \right]} \]

The HWHM is equal to :

\[ HWHM = ??? \]

ModelId=43849A020287.

static int FunctionDerivate (const gsl_vector *peak_params, void *calc_data, gsl_matrix *Jacobian)
 Partial Derivates of the Sigmoid function according each of the variable parameters

\[ Y_i = \frac{height}{1 + \exp \left( - \frac{x_i - pos}{width} \right)} + background \]

\[ \frac{\partial Y_i}{\partial pos} = - \frac{height}{width}* \frac{\exp \left( - \frac{x_i - pos}{width}\right)}{ \left[ 1 + \exp \left( - \frac{x_i - pos}{width}\right) \right]^2} \]

\[ \frac{\partial Y_i}{\partial width} = - \frac{height}{width} * \left( \frac{x_i - pos}{width}\right) * \frac{\exp \left( - \frac{x_i - pos}{width}\right)}{ \left[ 1 + \exp \left( - \frac{x_i - pos}{width}\right) \right]^2} \]

\[ \frac{\partial Y_i}{\partial height} = \frac{1}{1 + \exp \left( - \frac{x_i - pos}{width} \right)} \]

The jacobian matrix is filled with the derivate of Yi with the parameters :

\[ \frac{\partial Y_i}{\partial pos} / \sigma_i \]

\[ \frac{\partial Y_i}{\partial width} / \sigma_i \]

\[ \frac{\partial Y_i}{\partial height} / \sigma_i \]

with $ \sigma_i $ the error on data number $ i $
ModelId=43849A02028C.

static int FunctionAndDerivate (const gsl_vector *peak_params, void *calc_data, gsl_vector *func, gsl_matrix *Jacobian)

Private Attributes

double _dXLow
double _dXHigh
long _lPositionIndex

Detailed Description

Definition at line 18 of file FittingFunctionSigmoid.h.


Constructor & Destructor Documentation

FittingFunctionSigmoid::FittingFunctionSigmoid  ) 
 

Definition at line 20 of file FittingFunctionSigmoid.cpp.

References FittingFunction::_dSearchLimit, FittingFunction::_iNbParameters, FittingFunction::_lNbMaxIteration, FittingFunction::_mFunctionMap, FittingFunction::_pParameters, FittingFunction::_pParametersErrors, FittingFunction::_vInitialGuess, NB_ITERATION_MAX, and SEARCH_LIMIT.

00020                                               : _dXLow(0),_dXHigh(0)
00021 {       
00022         _iNbParameters                          =       3;
00023         _lNbMaxIteration                        =       NB_ITERATION_MAX;
00024         _dSearchLimit                           =       SEARCH_LIMIT;
00025         _pParameters                            =       new double[_iNbParameters];
00026         _pParametersErrors                      =       new double[_iNbParameters];
00027         _mFunctionMap["Position"]       =       0;
00028         _mFunctionMap["Width"]          =       1;
00029         _mFunctionMap["Height"]         =       2;      
00030         _vInitialGuess = gsl_vector_alloc(3);
00031 }

FittingFunctionSigmoid::FittingFunctionSigmoid long  lNbIteration,
double  dLimitSearch
 

Definition at line 34 of file FittingFunctionSigmoid.cpp.

References FittingFunction::_dSearchLimit, FittingFunction::_iNbParameters, FittingFunction::_lNbMaxIteration, FittingFunction::_mFunctionMap, FittingFunction::_pParameters, FittingFunction::_pParametersErrors, and FittingFunction::_vInitialGuess.

00034                                                                                    : _dXLow(0),_dXHigh(0) 
00035 {
00036         _lNbMaxIteration                        =       lNbIteration;
00037         _dSearchLimit                           =       dLimitSearch;
00038         _iNbParameters                          =       3;
00039         _pParameters                            =       new double[_iNbParameters];
00040         _pParametersErrors                      =       new double[_iNbParameters];
00041         _mFunctionMap["Position"]       =       0;
00042         _mFunctionMap["Width"]          =       1;
00043         _mFunctionMap["Height"]         =       2;      
00044         _vInitialGuess = gsl_vector_alloc(3);
00045 }

FittingFunctionSigmoid::~FittingFunctionSigmoid  )  [virtual]
 

Definition at line 48 of file FittingFunctionSigmoid.cpp.

References FittingFunction::_pParameters, FittingFunction::_pParametersErrors, and FittingFunction::_vInitialGuess.

00049 {
00050         printf ("Destruction FittingFunctionSigmoid !!!\n");
00051         delete[] _pParameters;
00052         delete[] _pParametersErrors;
00053         _pParameters = NULL;
00054         _pParametersErrors = NULL;
00055         gsl_vector_free(_vInitialGuess);
00056 }


Member Function Documentation

double FittingFunctionSigmoid::computeDerivateValue double  dX,
double  dPosition,
double  dWidth,
double  dHeight
[virtual]
 

Implements FittingFunction.

Definition at line 434 of file FittingFunctionSigmoid.cpp.

Referenced by estimateSigmoidIntersections().

00435 {
00436         double dXPos,dYi;
00437         
00438         if (dWidth == 0) dWidth = 1E-20;        //avoid to divide by zero
00439 
00440         dXPos                           = dX - dPosition;               
00441         dYi                                     = dHeight * exp(-dXPos/dWidth) /  ((1 + exp(-dXPos/dWidth)) * (1 + exp(-dXPos/dWidth)) * dWidth);
00442                         
00443         return dYi;
00444 }

double FittingFunctionSigmoid::computeValue double  dX,
double  dPosition,
double  dWidth,
double  dHeight,
double  dBackground
[virtual]
 

Implements FittingFunction.

Definition at line 421 of file FittingFunctionSigmoid.cpp.

Referenced by estimateSigmoidIntersections().

00422 {
00423 
00424         double dXPos,dYi;
00425         
00426         if (dWidth == 0) dWidth = 1E-20;        //avoid to divide by zero
00427 
00428         dXPos                           = dX - dPosition;               
00429         dYi                                     = dHeight / (1 + exp(-dXPos/dWidth));
00430                         
00431         return dYi;
00432 }

int FittingFunctionSigmoid::doFit long  points,
gsl_vector *  x,
gsl_vector *  y,
gsl_vector *  sigma,
gsl_vector *  guess,
bool  bEstimateGuess,
int  iSearchStoppingMethod = 1,
bool  bUseScaled = true
[virtual]
 

Implements FittingFunction.

Definition at line 178 of file FittingFunctionSigmoid.cpp.

References FittingFunction::_dSearchLimit, FittingFunction::_iNbParameters, FittingFunction::_vInitialGuess, getFunction(), getFunctionAndDerivate(), getFunctionDerivate(), initializeInitialsParameters(), FittingData::n, printState(), FittingData::sigma, FittingData::x, and FittingData::y.

00179 {
00180 initializeInitialsParameters(x,y,guess,bEstimateGuess);
00181 
00182         long            i;
00183         long            iter;
00184         long            status;
00185         
00186         double* datax = new double[points];
00187         double* datay = new double[points];
00188         double* datas = new double[points];
00189                         
00190         // fill the data arrays 
00191         for (i=0; i<points; i++)
00192         {               
00193                 datax[i] = gsl_vector_get (x, i);               // x = index            
00194                 datay[i] = gsl_vector_get (y, i);               // y = data to be fit    
00195                 datas[i] = gsl_vector_get (sigma, i);   // s = sigma (weight) of data values
00196         }
00197         
00198         // gsl data structure   
00199         struct FittingData gsl_data;
00200         gsl_data.n      = points;
00201         gsl_data.x      = datax;
00202         gsl_data.y              = datay;
00203         gsl_data.sigma  = datas;
00204 
00205         cout << "\n---------------------------------------------" << endl;
00206         cout << "Initials Parameters For The Fit" << endl;
00207         // set-up the solver 
00208         const gsl_multifit_fdfsolver_type *T;
00209         if (bUseScaled)
00210         {
00211                 T       = gsl_multifit_fdfsolver_lmsder;
00212                 cout << "Jacobian Data Scaled" << endl;
00213         }
00214         else
00215         {
00216                 T       = gsl_multifit_fdfsolver_lmder;
00217                 cout << "Jacobian Data Not Scaled" << endl;
00218         }
00219         gsl_multifit_fdfsolver *s                               = gsl_multifit_fdfsolver_alloc (T, points, _iNbParameters);
00220         gsl_multifit_function_fdf f;
00221 
00222         f.f             = getFunction();            
00223         f.df            = getFunctionDerivate();    
00224         f.fdf           = getFunctionAndDerivate(); 
00225         f.n             = points;
00226         f.p             = _iNbParameters;
00227         f.params        = &gsl_data;
00228   
00229         gsl_multifit_fdfsolver_set ( s,  &f,  _vInitialGuess );
00230         switch(iSearchStoppingMethod)
00231         {
00232                 case 1 :
00233                         cout << "Test Delta Search Stopping Function Used" << endl;
00234                         break;
00235 
00236                 case 2 :
00237                         cout << "Test Gradient Search Stopping Function Used" << endl;
00238                         break;
00239         }
00240 
00241         cout << "\n---------------------------------------------" << endl;
00242         cout << "Initials Values For The Fit" << endl;
00243         cout << "Position = " << gsl_vector_get (_vInitialGuess, 0) << endl;
00244         cout << "Width    = " << gsl_vector_get (_vInitialGuess, 1) << endl;
00245         cout << "Height   = " << gsl_vector_get (_vInitialGuess, 2) << endl;
00246         
00247   
00248         iter = 0;
00249         cout << "\n---------------------------------------------" << endl;
00250         cout << "Iterations" << endl;
00251         #ifdef DEBUG_FIT
00252           printState (iter, s);
00253         #endif
00254 
00255         do
00256     {
00257                 iter++;
00258                 status = gsl_multifit_fdfsolver_iterate (s);
00259 
00260                 #ifdef DEBUG_FIT
00261                           printf ("status = %s\n", gsl_strerror (status));
00262                           printState (iter, s);
00263                 #endif
00264                 
00265                 if (status)
00266                 {
00267                 break;
00268                 }
00269 
00270 
00271                 switch(iSearchStoppingMethod)
00272                 {
00273                 case 1 :
00274                         status = gsl_multifit_test_delta (s->dx, s->x,  /*_dSearchLimit*/0.0, _dSearchLimit);
00275                         break;
00276 
00277                 case 2 :
00278                         gsl_vector* gradient = gsl_vector_alloc(_iNbParameters);
00279                         gsl_multifit_gradient(s->J,s->f,gradient);
00280                         status = gsl_multifit_test_gradient (gradient, _dSearchLimit);
00281                         gsl_vector_free(gradient);
00282                         break;
00283                 }
00284 
00285     }
00286         while (status == GSL_CONTINUE && iter < _lNbMaxIteration);
00287 
00288         _lNbIterations = iter;
00289         // Check whether the fitting has converged or failed  
00290         //if ( status == GSL_SUCCESS )
00291         //      {
00292         gsl_matrix *covar = gsl_matrix_alloc (_iNbParameters, _iNbParameters);
00293         gsl_multifit_covar (s->J, 0.0, covar);
00294         double chi                              = gsl_blas_dnrm2(s->f);
00295 
00296         _pParameters[0]                 =       FIT(0);
00297         _pParameters[1]                 =       FIT(1);
00298         _pParameters[2]                 =       FIT(2);
00299 
00300         _pParametersErrors[0]   =       ERR(0);
00301         _pParametersErrors[1]   =       ERR(1);
00302         _pParametersErrors[2]   =       ERR(2);
00303 
00304         double dHeight  = _pParameters[2];
00305         _dQualityFactor =       pow(chi, 2.0)/ (points - _iNbParameters);
00306         _dSigma                 =       1.0/ (sqrt(2*M_PI)*dHeight);
00307         _dFHWM                  =       _pParameters[1]; // width //2.355*_dSigma;
00308         _dHWHM                  =       0.5*_dFHWM;     
00309         #ifdef DEBUG_FIT
00310                         cout << "\n---------------------------------------------" << endl;
00311                         cout << "Results For Sigmoid Fit" << endl;
00312                         printResults(); 
00313         #endif  
00314                 
00315         gsl_matrix_free(covar);         
00316         //      }
00317 
00318         #ifdef DEBUG_FIT  
00319                 printf ("status = %s\n", gsl_strerror (status));
00320         #endif
00321         
00322         // Free data 
00323         gsl_multifit_fdfsolver_free (s);
00324         
00325         estimateSigmoidIntersections();
00326 
00327         delete[] datax;
00328         delete[] datay;
00329         delete[] datas;
00330         
00331         datax = NULL;
00332         datax = NULL;
00333         datax = NULL;
00334 
00335         return status;
00336 }

Here is the call graph for this function:

double FittingFunctionSigmoid::estimateInitialBackground gsl_vector *  vExperimentalDataX,
gsl_vector *  vExperimentalDataY
 

Definition at line 496 of file FittingFunctionSigmoid.cpp.

Referenced by estimateInitialHeight(), estimateInitialPosition(), and estimateInitialSigma().

00497 {
00498         double dBackground = 0.0;
00499         return dBackground;
00500 }

void FittingFunctionSigmoid::estimateInitialGuess gsl_vector *  vExperimentalDataX,
gsl_vector *  vExperimentalDataY
[virtual]
 

Implements FittingFunction.

Definition at line 632 of file FittingFunctionSigmoid.cpp.

Referenced by initializeInitialsParameters().

00633 {
00634         // fill the initials parameters values estimation vector for the initialisation
00635         // [0] = position 
00636         // [1] = width  
00637         // [2] = height
00638         // [3] = background
00639         // Thanks to the max index of all the points, we can get the
00640         // position_value       = x_value at index_max
00641         // width                        = standard deviation
00642         // height_value         = y_value - min_value 
00643         //min_index = gsl_vector_min_index (vExperimentalDataY);
00644         //max_index = gsl_vector_max_index (vExperimentalDataY);
00645 
00646         double dInitialPosition = estimateInitialPosition(vExperimentalDataX,vExperimentalDataY);
00647         double dInitialSigma    = estimateInitialSigma(vExperimentalDataX,vExperimentalDataY);
00648         double dInitialHeight   = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00649 
00650         gsl_vector_set(_vInitialGuess, 0, dInitialPosition);
00651         gsl_vector_set(_vInitialGuess, 1, dInitialSigma);
00652         gsl_vector_set(_vInitialGuess, 2, dInitialHeight);
00653 }

double FittingFunctionSigmoid::estimateInitialHeight gsl_vector *  vExperimentalDataX,
gsl_vector *  vExperimentalDataY
 

Definition at line 484 of file FittingFunctionSigmoid.cpp.

References estimateInitialBackground().

Referenced by estimateInitialPosition(), and estimateInitialSigma().

00485 {
00486         size_t max_index        = gsl_vector_max_index (vExperimentalDataY);
00487 
00488         double background       = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00489         double maxValue         = gsl_vector_get(vExperimentalDataY, max_index);
00490         double height           = maxValue - background;
00491 
00492         return height;
00493 }

Here is the call graph for this function:

double FittingFunctionSigmoid::estimateInitialPosition gsl_vector *  vExperimentalDataX,
gsl_vector *  vExperimentalDataY
 

Definition at line 447 of file FittingFunctionSigmoid.cpp.

References estimateInitialBackground(), and estimateInitialHeight().

00448 {
00449         double dHeight                          = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00450         double dBackground                      = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00451         double dHalfMaximum                     = 0.5 * dHeight + dBackground;
00452           
00453         size_t lVectorSize                      = vExperimentalDataX->size;//lSizeX;
00454         size_t iCurrentIndex            = 0;
00455         double dYCurrentValue           = gsl_vector_get(vExperimentalDataY, iCurrentIndex);
00456 
00457         double dPosition;
00458 
00459         while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue < dHalfMaximum)
00460         {
00461                 dYCurrentValue = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00462         }
00463           
00464     if (iCurrentIndex != lVectorSize)
00465     {
00466                 dPosition                       = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00467                 _lPositionIndex         = (long)iCurrentIndex;
00468     }
00469     else
00470     {
00471                 //- try an arbitrary value that should be not so far from reality
00472                 cout << "Arbitrary value !!!!" << endl;
00473                 dPosition = gsl_vector_get(vExperimentalDataX, lVectorSize / 2);
00474                 _lPositionIndex         = (long)(lVectorSize/2);
00475     }
00476 /*      cout    << "\n--------------------- Estimation Initial Position -------------------" << endl;
00477         cout    << "\tindex=" << _lPositionIndex 
00478                         << " --> (x,y)=(" << gsl_vector_get(vExperimentalDataX,_lPositionIndex) 
00479                         << "," << gsl_vector_get(vExperimentalDataY,_lPositionIndex) << ")" << endl;
00480 */
00481         return dPosition;
00482 }

Here is the call graph for this function:

double FittingFunctionSigmoid::estimateInitialSigma gsl_vector *  vExperimentalDataX,
gsl_vector *  vExperimentalDataY
 

Definition at line 503 of file FittingFunctionSigmoid.cpp.

References estimateInitialBackground(), and estimateInitialHeight().

00504 {
00505         double dSigma;
00506         double dCovariance;
00507         double dHeight                          = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00508         double dBackground                      = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00509         
00510         size_t iIndexMax                        = gsl_vector_max_index(vExperimentalDataY);
00511         double dYMaxValue                       = gsl_vector_get(vExperimentalDataY, iIndexMax);
00512   
00513         double  dYCurrentValue          = dYMaxValue;
00514 
00515         size_t iCurrentIndex            = iIndexMax;
00516         size_t lVectorSize                      = vExperimentalDataX->size;//lSizeX;
00517         double lHalfVectorSize          = 0.5 * lVectorSize;
00518         double dHalfMaximum                     = 0.5 * dHeight + dBackground;
00519 
00520         double d34Maximum                       = 3/4 * dHeight + dBackground;
00521         double d14Maximum                       = 1/4 * dHeight + dBackground;
00522 
00523         size_t iIndexLow                        = iCurrentIndex;
00524         size_t iIndexHigh                       = iCurrentIndex;
00525 
00526         size_t iQuarterIndexLow         = iCurrentIndex;
00527         size_t iQuarterIndexHigh        = iCurrentIndex;
00528 
00529 
00530     if (iIndexMax < lHalfVectorSize)
00531     {
00532                 //cout << "\n--> Look in increasing direction" << endl;
00533                 //- look in increasing direction
00534                 while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue > dHalfMaximum)
00535                 {
00536                         
00537                         dYCurrentValue  = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00538                         iIndexHigh              = iCurrentIndex;
00539                         iIndexLow               = iIndexMax - iCurrentIndex;
00540                 }
00541     }
00542     else
00543     {
00544                 //cout << "\n--> Look in decreasing direction" << endl;
00545                 //- look in decreasing direction
00546                 while( iCurrentIndex > 0 && dYCurrentValue > dHalfMaximum)
00547                 {
00548                         dYCurrentValue =  gsl_vector_get(vExperimentalDataY, --iCurrentIndex);
00549                         iIndexLow               = iCurrentIndex;
00550                         iIndexHigh              = iIndexMax + iCurrentIndex;
00551                 }
00552     }
00553 
00554     if (iCurrentIndex != iIndexMax)
00555     {
00556                 double dPosition        = gsl_vector_get(vExperimentalDataX, iIndexMax);
00557                 double dX2                      = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00558                 dSigma = ::fabs(dX2 - dPosition);
00559     }
00560     else
00561     {
00562                 //- try an arbitrary value that should be not so far from reality
00563                 dSigma = lVectorSize / 6;
00564     }
00565 
00566 /*      cout << "-------------------------------" << endl;
00567         cout << "Peak : (x,y) --> ("    << gsl_vector_get(vExperimentalDataX, iIndexMax) << "," 
00568                                                                         << gsl_vector_get(vExperimentalDataY, iIndexMax) << ")" << endl;
00569                 
00570         
00571 
00572         cout << "Low : index=" << iIndexLow << endl;
00573         cout << "\txlow=" << gsl_vector_get(vExperimentalDataX, iIndexLow) << " | ylow=" << gsl_vector_get(vExperimentalDataY, iIndexLow) << endl;
00574         cout << "\txhigh=" << gsl_vector_get(vExperimentalDataX, iIndexHigh) << " | yhigh=" << gsl_vector_get(vExperimentalDataY, iIndexHigh) << endl;
00575         
00576         cout << "\txhigh - xlow= " << gsl_vector_get(vExperimentalDataX, iIndexHigh) - gsl_vector_get(vExperimentalDataX, iIndexLow) << endl;
00577         cout << "\tsigma=" << dSigma << endl;
00578 */
00579 
00580         cout << "\n---------------------------------------------" << endl;
00581         cout << "Estimation Initial Width" << endl;
00582         double d34Value = 3.0/4*dHeight + dBackground;
00583         double d12Value = 1.0/2*dHeight + dBackground;
00584         double d14Value = 1.0/4*dHeight + dBackground;
00585 
00586         int i34Index = getIndexValue(d34Value,vExperimentalDataY);
00587         int i12Index = getIndexValue(d12Value,vExperimentalDataY);
00588         int i14Index = getIndexValue(d14Value,vExperimentalDataY);
00589 
00590         double d34XValue = gsl_vector_get(vExperimentalDataX,i34Index);
00591         double d34YValue = gsl_vector_get(vExperimentalDataY,i34Index);
00592         
00593         double d12XValue = gsl_vector_get(vExperimentalDataX,i12Index);
00594         double d12YValue = gsl_vector_get(vExperimentalDataY,i12Index);
00595         
00596         double d14XValue = gsl_vector_get(vExperimentalDataX,i14Index);
00597         double d14YValue = gsl_vector_get(vExperimentalDataY,i14Index);
00598 
00599         cout << "3/4A=" << "(" << d34Value << ")-->" << i34Index << " (x,y)=(" << d34XValue << "," << d34YValue << ")" << endl;
00600         cout << "1/2A=" << "(" << d12Value << ")-->" << i12Index << " (x,y)=(" << d12XValue << "," << d12YValue << ")" << endl;
00601         cout << "1/4A=" << "(" << d14Value << ")-->" << i14Index << " (x,y)=(" << d14XValue << "," << d14YValue << ")" << endl;
00602 
00603         double dWidthEstimated = ::fabs(d34XValue - d14XValue);
00604 
00605         cout << "Width estimated = " << dWidthEstimated << endl;
00606 
00607     dCovariance = dSigma * dSigma;
00608 
00609 /*      cout << "------------- ESTIMATION DU WIDTH DE LA SIGMOID --------------" << endl;
00610         
00611         double xpos = gsl_vector_get(vExperimentalDataX, _lPositionIndex);
00612         double ypos = gsl_vector_get(vExperimentalDataY, _lPositionIndex);
00613         double xpos1 = gsl_vector_get(vExperimentalDataX, _lPositionIndex+1);
00614         double ypos1 = gsl_vector_get(vExperimentalDataY, _lPositionIndex+1);
00615 
00616 
00617         cout << "Index pos=" << _lPositionIndex << " --> (xpos,ypos)=(" << xpos << "," << ypos << ")" << endl;
00618         cout << "Index pos+1=" << _lPositionIndex+1 << " --> (xpos+1,ypos+1)=(" << xpos1 << "," << ypos1 << ")" << endl;
00619 
00620         double dPente = (ypos1 - ypos)/(xpos1 - xpos);
00621         cout << "Pente=" << dPente << endl;
00622 
00623         dSigma = dPente;
00624 */
00625         dSigma = dWidthEstimated/2.0;           //better than pente
00626 
00627         return dSigma;
00628 }

Here is the call graph for this function:

void FittingFunctionSigmoid::estimateSigmoidIntersections  ) 
 

Definition at line 671 of file FittingFunctionSigmoid.cpp.

References _dXHigh, _dXLow, computeDerivateValue(), computeValue(), and FittingFunction::getParameter().

00672 {
00673 
00674         double xlow;
00675         double xhigh;
00676         double dBackground;
00677 //      if (fitFunction->getNbParameters() == 4)
00678 //      {
00679 //              dBackground     = fitFunction->getParameter(3);
00680 //      }
00681 //      else 
00682                 dBackground = 0.0;
00683         
00684         double ylow                     = dBackground;                                                  //ylow = background
00685         double yhigh            = ylow + getParameter(2);       //yhigh = ylow + height;
00686 
00687         double dPosition        = getParameter(0);
00688         double dWidth           = getParameter(1);
00689         double dHeight          = getParameter(2);
00690 
00691         //At this time the function is estimated position/width/height and background are known
00692 
00693         //the tangente equation : y = f'(position)(x - position)+ f(position)  -->f'
00694         //f is equal to sigmoid function  --> 
00695         //f' is the derivate of the sigmoid function
00696         // on reevalue yhigh avec la fonction sigmoid 
00697         // so to get x it is : (y - f(position))/(f'(position)) + position
00698         double fa               = computeValue(dPosition,dPosition,dWidth,dHeight,dBackground); 
00699         double fprimea  = computeDerivateValue(dPosition,dPosition,dWidth,dHeight);
00700 
00701         xlow                    = (ylow - fa)/fprimea + dPosition;
00702         xhigh                   = (yhigh - fa)/fprimea + dPosition;
00703 
00704         cout << "\n---------------------------------------------" << endl;
00705         cout << "Intersections" << endl;
00706         cout << "(xlow ,ylow )=(" << xlow << "," << ylow << ")" << endl;  
00707         cout << "(xhigh,yhigh)=(" << xhigh << "," << yhigh << ")" << endl;
00708 
00709         _dXLow  = xlow;
00710         _dXHigh = xhigh;
00711 }

Here is the call graph for this function:

int FittingFunctionSigmoid::Function const gsl_vector *  peak_params,
void *  calc_data,
gsl_vector *  func
[static]
 

Implements the Sigmoid function for the model :

\[ Y_i = \frac{height}{1 + \exp \left[ - \frac{x_i - pos}{width} \right]} \]

The HWHM is equal to :

\[ HWHM = ??? \]

ModelId=43849A020287.

Definition at line 79 of file FittingFunctionSigmoid.cpp.

References FittingData::sigma, FittingData::x, and FittingData::y.

Referenced by FunctionAndDerivate(), and getFunction().

00080 {
00081         size_t i;
00082         size_t iNbPoints        = ((struct FittingData *)calc_data)->n;
00083         double *x                       = ((struct FittingData *)calc_data)->x;
00084         double *y                       = ((struct FittingData *)calc_data)->y;
00085         double *sigma           = ((struct FittingData *) calc_data)->sigma;
00086 
00087         double pos                      = gsl_vector_get (peak_params, 0);
00088         double width            = gsl_vector_get (peak_params, 1);
00089         double height           = gsl_vector_get (peak_params, 2);
00090 
00091 //#ifdef DEBUG_FIT      
00092 //  printf("pos = %g, width = %g, height = %g\n", pos, width, height);
00093 //#endif        
00094         double dXPos,dYi;
00095         
00096         if (width == 0) width = 1E-20;  //avoid to divide by zero
00097 
00098         for (i=0; i < iNbPoints; i++)
00099         {               
00100                 dXPos                           = x[i] - pos;           
00101                 dYi                                     = height / (1 + exp(-dXPos/width));
00102 
00103                 gsl_vector_set (func, i,( dYi - y[i]) / sigma[i]);              
00104         }
00105         return GSL_SUCCESS;
00106 }

int FittingFunctionSigmoid::FunctionAndDerivate const gsl_vector *  peak_params,
void *  calc_data,
gsl_vector *  func,
gsl_matrix *  Jacobian
[static]
 

Definition at line 168 of file FittingFunctionSigmoid.cpp.

References Function(), and FunctionDerivate().

Referenced by getFunctionAndDerivate().

00169 {       
00170         Function  (peak_params, calc_data, func);
00171         FunctionDerivate (peak_params, calc_data, Jacobian);
00172         return GSL_SUCCESS;
00173 }

Here is the call graph for this function:

int FittingFunctionSigmoid::FunctionDerivate const gsl_vector *  peak_params,
void *  calc_data,
gsl_matrix *  Jacobian
[static]
 

Partial Derivates of the Sigmoid function according each of the variable parameters

\[ Y_i = \frac{height}{1 + \exp \left( - \frac{x_i - pos}{width} \right)} + background \]

\[ \frac{\partial Y_i}{\partial pos} = - \frac{height}{width}* \frac{\exp \left( - \frac{x_i - pos}{width}\right)}{ \left[ 1 + \exp \left( - \frac{x_i - pos}{width}\right) \right]^2} \]

\[ \frac{\partial Y_i}{\partial width} = - \frac{height}{width} * \left( \frac{x_i - pos}{width}\right) * \frac{\exp \left( - \frac{x_i - pos}{width}\right)}{ \left[ 1 + \exp \left( - \frac{x_i - pos}{width}\right) \right]^2} \]

\[ \frac{\partial Y_i}{\partial height} = \frac{1}{1 + \exp \left( - \frac{x_i - pos}{width} \right)} \]

The jacobian matrix is filled with the derivate of Yi with the parameters :

\[ \frac{\partial Y_i}{\partial pos} / \sigma_i \]

\[ \frac{\partial Y_i}{\partial width} / \sigma_i \]

\[ \frac{\partial Y_i}{\partial height} / \sigma_i \]

with $ \sigma_i $ the error on data number $ i $
ModelId=43849A02028C.

Definition at line 120 of file FittingFunctionSigmoid.cpp.

References FittingData::n, FittingData::sigma, and FittingData::x.

Referenced by FunctionAndDerivate(), and getFunctionDerivate().

00121 {
00122 
00123         size_t i;
00124         
00125         size_t  n               = ((struct FittingData *)calc_data)->n;
00126         double *x               = ((struct FittingData *)calc_data)->x;
00127         double *sigma   = ((struct FittingData *)calc_data)->sigma;
00128 
00129         double pos              = gsl_vector_get (peak_params, 0);
00130         double width    = gsl_vector_get (peak_params, 1);
00131         double height   = gsl_vector_get (peak_params, 2);
00132   
00133         double d_pos, d_width, d_height;
00134         double s, alpha, e,e2,expo,xpos;
00135         
00136         // Fill the Jacobian Matrix     
00137         // Jacobian matrix J(i,j) = dfi / dxj, 
00138         // where fi = (Yi - yi)/sigma[i],  
00139         // Yi = height / [ 1 + exp(-(x-pos)/width)] 
00140         // and the xj are the parameters (pos,width,height)
00141         for (i=0; i< n; i++)
00142         {
00143                 s                               = sigma[i]; 
00144                 xpos                    = x[i] - pos;
00145                 alpha                   = -xpos/width;
00146                 expo                    = exp(alpha);
00147                 e                               = 1 + expo;
00148                 e2                              = e * e;
00149 
00150                 // derivate dYi/dPos
00151                 d_pos                   = -((height*expo)/(width*e2))/s;                
00152                 // derivate dYi/dWidth
00153                 d_width                 = -((height*xpos*expo)/(e2*width*width)) / s;           
00154                 // derivate dYi/dHeight
00155                 d_height                = (1/e)/s;
00156                 
00157                 // fill the jacobion matrix J(points, parameters)
00158                 gsl_matrix_set (Jacobian, i, 0, d_pos); 
00159                 gsl_matrix_set (Jacobian, i, 1, d_width);
00160                 gsl_matrix_set (Jacobian, i, 2, d_height);
00161         }
00162          
00163         return GSL_SUCCESS;
00164 }

int FittingFunctionSigmoid::generateFunctionFit long  lNbPoints,
gsl_vector *  x,
gsl_vector *  funcX,
gsl_vector *  funcY
[virtual]
 

Implements FittingFunction.

Definition at line 366 of file FittingFunctionSigmoid.cpp.

References FittingFunction::_pParameters.

00367 {
00368         int i;
00369 
00370         double dPos             =       _pParameters[0];
00371         double dWidth   =       _pParameters[1];
00372         double dHeight  =       _pParameters[2];
00373 
00374         double dXi,dXPos,dExp,dYi;
00375 
00376         for (i=0; i < lNbPoints; i++)
00377     {           
00378                 dXi             = gsl_vector_get(x,i);
00379                 dXPos   = dXi  - dPos;
00380                 dExp    = exp(-dXPos/dWidth);
00381                                 
00382                 dYi             = dHeight / (1.0 + dExp);
00383 
00384                 gsl_vector_set (funcX, i, dXi);
00385                 gsl_vector_set (funcY, i, dYi);         
00386     }   
00387         return 1;
00388 }

int FittingFunctionSigmoid::generateFunctionFit double  dValMinX,
double  dResolutionX,
long  lNbPoints,
gsl_vector *  funcX,
gsl_vector *  funcY
[virtual]
 

Implements FittingFunction.

Definition at line 340 of file FittingFunctionSigmoid.cpp.

References FittingFunction::_pParameters.

00340                                                                                                                                      {
00341 
00342         int i;
00343 
00344         double dPos             =       _pParameters[0];
00345         double dWidth   =       _pParameters[1];
00346         double dHeight  =       _pParameters[2];
00347 
00348         double dXi,dXPos,dExp,dYi;
00349 
00350         for (i=0; i < lNbPoints; i++)
00351     {           
00352                 dXi             = dValMinX + i*dResolutionX;
00353                 dXPos   = dXi  - dPos;
00354                 dExp    = exp(-dXPos/dWidth);
00355                                 
00356                 dYi             = dHeight / (1.0 + dExp);
00357 
00358                 gsl_vector_set (funcX, i, dXi);
00359                 gsl_vector_set (funcY, i, dYi);         
00360     }   
00361         return 1;
00362 }

std::string FittingFunctionSigmoid::getEquation  )  const [virtual]
 

Implements FittingFunction.

Definition at line 13 of file FittingFunctionSigmoid.cpp.

00014 {
00015         return "f(x) = Height / [1 + exp[ - (x - Position) / Width]";
00016 }

PtrFunction FittingFunctionSigmoid::getFunction  )  const [virtual]
 

Implements FittingFunction.

Definition at line 59 of file FittingFunctionSigmoid.cpp.

References Function().

Referenced by doFit().

00059                                                      {
00060         return FittingFunctionSigmoid::Function;
00061 }

Here is the call graph for this function:

PtrFunctionAndDerivate FittingFunctionSigmoid::getFunctionAndDerivate  )  const [virtual]
 

Implements FittingFunction.

Definition at line 69 of file FittingFunctionSigmoid.cpp.

References FunctionAndDerivate().

Referenced by doFit().

00069                                                                             {
00070         return FittingFunctionSigmoid::FunctionAndDerivate;
00071 }

Here is the call graph for this function:

PtrFunctionDerivate FittingFunctionSigmoid::getFunctionDerivate  )  const [virtual]
 

Implements FittingFunction.

Definition at line 64 of file FittingFunctionSigmoid.cpp.

References FunctionDerivate().

Referenced by doFit().

00064                                                                      {
00065         return FittingFunctionSigmoid::FunctionDerivate;
00066 }

Here is the call graph for this function:

double FittingFunctionSigmoid::getXHigh  )  const [inline]
 

Definition at line 66 of file FittingFunctionSigmoid.h.

References _dXHigh.

00066 {return _dXHigh;};

double FittingFunctionSigmoid::getXLow  )  const [inline]
 

Definition at line 65 of file FittingFunctionSigmoid.h.

References _dXLow.

00065 {return _dXLow;};

void FittingFunctionSigmoid::initializeInitialsParameters gsl_vector *  vExperimentalDataX,
gsl_vector *  vExperimentalDataY,
gsl_vector *  vInitialGuess,
bool  bEstimate
[virtual]
 

Implements FittingFunction.

Definition at line 656 of file FittingFunctionSigmoid.cpp.

References FittingFunction::_vInitialGuess, and estimateInitialGuess().

Referenced by doFit().

00657 {
00658         if (bEstimate) 
00659         {
00660                 estimateInitialGuess(vExperimentalDataX,vExperimentalDataY);
00661         }
00662         else
00663         {
00664                 gsl_vector_set(_vInitialGuess,0,gsl_vector_get(vInitialGuess,0));
00665                 gsl_vector_set(_vInitialGuess,1,gsl_vector_get(vInitialGuess,1));
00666                 gsl_vector_set(_vInitialGuess,2,gsl_vector_get(vInitialGuess,2));
00667         }
00668 }

Here is the call graph for this function:

void FittingFunctionSigmoid::printResults  ) 
 

Definition at line 409 of file FittingFunctionSigmoid.cpp.

References FittingFunction::_dFHWM, FittingFunction::_dHWHM, FittingFunction::_dQualityFactor, FittingFunction::_dSigma, FittingFunction::_pParametersErrors, and FittingFunction::getParameter().

00410 {
00411                 printf ("pos            = %g \t+/- %g\n",       getParameter("Position"),_pParametersErrors[0]);
00412                 printf ("width          = %g \t+/- %g\n",       getParameter("Width"),  _pParametersErrors[1]);
00413                 printf ("height         = %g \t+/- %g\n",       getParameter("Height"), _pParametersErrors[2]);
00414         printf ("chisq/dof      = %g\n",                                _dQualityFactor);
00415                 printf ("Sigma          = %g\n",                                _dSigma);
00416                 printf ("HWHM           = %g\n",                                _dHWHM);
00417                 printf ("FHWM           = %g\n",                                _dFHWM);
00418 }

Here is the call graph for this function:

void FittingFunctionSigmoid::printState size_t  iter,
gsl_multifit_fdfsolver *  s
 

Definition at line 398 of file FittingFunctionSigmoid.cpp.

Referenced by doFit().

00399 {
00400   printf ("iter: %3u x : %g %g %g |f(x)|=%g\n",
00401           iter,
00402           gsl_vector_get (s->x, 0), 
00403           gsl_vector_get (s->x, 1),
00404           gsl_vector_get (s->x, 2), 
00405           gsl_blas_dnrm2 (s->f));
00406 }


Member Data Documentation

double FittingFunctionSigmoid::_dXHigh [private]
 

Definition at line 74 of file FittingFunctionSigmoid.h.

Referenced by estimateSigmoidIntersections(), and getXHigh().

double FittingFunctionSigmoid::_dXLow [private]
 

Definition at line 73 of file FittingFunctionSigmoid.h.

Referenced by estimateSigmoidIntersections(), and getXLow().

long FittingFunctionSigmoid::_lPositionIndex [private]
 

Definition at line 75 of file FittingFunctionSigmoid.h.


The documentation for this class was generated from the following files:
Generated on Tue Apr 14 09:35:14 2009 for Data Fitting Library by  doxygen 1.4.5