FittingFunctionSigmoid.cpp

Go to the documentation of this file.
00001 // FittingFunctionSigmoid.cpp: implementation of the FittingFunctionSigmoid class.
00002 //
00004 #ifdef WIN32 
00005 #pragma warning(disable:4786)
00006 #endif
00007 
00008 #include "FittingFunctionSigmoid.h"
00010 // Construction/Destruction
00012 
00013 std::string     FittingFunctionSigmoid::getEquation() const
00014 {
00015         return "f(x) = Height / [1 + exp[ - (x - Position) / Width]";
00016 }
00017 
00018 
00019 //##ModelId=43849A020277
00020 FittingFunctionSigmoid::FittingFunctionSigmoid(): _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 }
00032 
00033 //##ModelId=43849A020278
00034 FittingFunctionSigmoid::FittingFunctionSigmoid(long lNbIteration,double dLimitSearch): _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 }
00046 
00047 //##ModelId=43849A02027B
00048 FittingFunctionSigmoid::~FittingFunctionSigmoid()
00049 {
00050         printf ("Destruction FittingFunctionSigmoid !!!\n");
00051         delete[] _pParameters;
00052         delete[] _pParametersErrors;
00053         _pParameters = NULL;
00054         _pParametersErrors = NULL;
00055         gsl_vector_free(_vInitialGuess);
00056 }
00057 
00058 //##ModelId=43849A02027D
00059 PtrFunction FittingFunctionSigmoid::getFunction() const{
00060         return FittingFunctionSigmoid::Function;
00061 }
00062 
00063 //##ModelId=43849A02027F
00064 PtrFunctionDerivate FittingFunctionSigmoid::getFunctionDerivate() const{
00065         return FittingFunctionSigmoid::FunctionDerivate;
00066 }
00067 
00068 //##ModelId=43849A020281
00069 PtrFunctionAndDerivate FittingFunctionSigmoid::getFunctionAndDerivate() const {
00070         return FittingFunctionSigmoid::FunctionAndDerivate;
00071 }
00072 
00073 
00078 //##ModelId=43849A020287
00079 int FittingFunctionSigmoid::Function(const gsl_vector *peak_params, void *calc_data, gsl_vector *func)
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 }
00107 
00108 
00119 //##ModelId=43849A02028C
00120 int FittingFunctionSigmoid::FunctionDerivate(const gsl_vector *peak_params, void *calc_data, gsl_matrix *Jacobian)
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 }
00165 
00166 // Interface function for the gsl_multifit_fdfsolver
00167 //##ModelId=43849A020296
00168 int FittingFunctionSigmoid::FunctionAndDerivate(const gsl_vector *peak_params, void *calc_data,gsl_vector *func, gsl_matrix *Jacobian) 
00169 {       
00170         Function  (peak_params, calc_data, func);
00171         FunctionDerivate (peak_params, calc_data, Jacobian);
00172         return GSL_SUCCESS;
00173 }
00174 
00175 
00176 
00177 //##ModelId=43849A02029C
00178 int FittingFunctionSigmoid::doFit(long points, gsl_vector *x, gsl_vector *y,gsl_vector *sigma ,gsl_vector *guess, bool bEstimateGuess, int iSearchStoppingMethod,bool bUseScaled)
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 }
00337 
00338 
00339 //##ModelId=43849A0202A8
00340 int FittingFunctionSigmoid::generateFunctionFit(double dValMinX,double dResolutionX,long lNbPoints,gsl_vector* funcX,gsl_vector* funcY){
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 }
00363 
00364 
00365 
00366 int FittingFunctionSigmoid::generateFunctionFit(long lNbPoints, gsl_vector *x, gsl_vector* funcX,gsl_vector* funcY)
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 }
00389 
00390 
00391 
00392 
00393 
00394 
00395 
00396 // print the state of the fitting during iterations 
00397 //##ModelId=43849A0202AF
00398 void FittingFunctionSigmoid::printState (size_t iter, gsl_multifit_fdfsolver * s)
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 }
00407 
00408 //##ModelId=43849A0202B7
00409 void FittingFunctionSigmoid::printResults()
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 }
00419 
00420 
00421 double  FittingFunctionSigmoid::computeValue(double dX,double dPosition, double dWidth, double dHeight, double dBackground=0.0)
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 }
00433 
00434 double  FittingFunctionSigmoid::computeDerivateValue(double dX,double dPosition, double dWidth, double dHeight)
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 }
00445 
00446 
00447 double FittingFunctionSigmoid::estimateInitialPosition(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
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 }
00483 
00484 double FittingFunctionSigmoid::estimateInitialHeight(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
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 }
00494 
00495 
00496 double FittingFunctionSigmoid::estimateInitialBackground(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00497 {
00498         double dBackground = 0.0;
00499         return dBackground;
00500 }
00501 
00502 
00503 double FittingFunctionSigmoid::estimateInitialSigma(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
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 }
00629 
00630 
00631 
00632 void    FittingFunctionSigmoid::estimateInitialGuess(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY)
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 }
00654 
00655 
00656 void    FittingFunctionSigmoid::initializeInitialsParameters(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY,gsl_vector* vInitialGuess,bool bEstimate)
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 }
00669 
00670 
00671 void FittingFunctionSigmoid::estimateSigmoidIntersections()
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 }
00712 
00713 

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