00001
00002
00004 #ifdef WIN32
00005 #pragma warning(disable:4786)
00006 #endif
00007
00008 #include "FittingFunctionSigmoid.h"
00010
00012
00013 std::string FittingFunctionSigmoid::getEquation() const
00014 {
00015 return "f(x) = Height / [1 + exp[ - (x - Position) / Width]";
00016 }
00017
00018
00019
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
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
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
00059 PtrFunction FittingFunctionSigmoid::getFunction() const{
00060 return FittingFunctionSigmoid::Function;
00061 }
00062
00063
00064 PtrFunctionDerivate FittingFunctionSigmoid::getFunctionDerivate() const{
00065 return FittingFunctionSigmoid::FunctionDerivate;
00066 }
00067
00068
00069 PtrFunctionAndDerivate FittingFunctionSigmoid::getFunctionAndDerivate() const {
00070 return FittingFunctionSigmoid::FunctionAndDerivate;
00071 }
00072
00073
00078
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
00092
00093
00094 double dXPos,dYi;
00095
00096 if (width == 0) width = 1E-20;
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
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
00137
00138
00139
00140
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
00151 d_pos = -((height*expo)/(width*e2))/s;
00152
00153 d_width = -((height*xpos*expo)/(e2*width*width)) / s;
00154
00155 d_height = (1/e)/s;
00156
00157
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
00167
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
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
00191 for (i=0; i<points; i++)
00192 {
00193 datax[i] = gsl_vector_get (x, i);
00194 datay[i] = gsl_vector_get (y, i);
00195 datas[i] = gsl_vector_get (sigma, i);
00196 }
00197
00198
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
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, 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
00290
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];
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
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
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
00397
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
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;
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;
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;
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
00472 cout << "Arbitrary value !!!!" << endl;
00473 dPosition = gsl_vector_get(vExperimentalDataX, lVectorSize / 2);
00474 _lPositionIndex = (long)(lVectorSize/2);
00475 }
00476
00477
00478
00479
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;
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
00533
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
00545
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
00563 dSigma = lVectorSize / 6;
00564 }
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
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
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625 dSigma = dWidthEstimated/2.0;
00626
00627 return dSigma;
00628 }
00629
00630
00631
00632 void FittingFunctionSigmoid::estimateInitialGuess(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY)
00633 {
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
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
00678
00679
00680
00681
00682 dBackground = 0.0;
00683
00684 double ylow = dBackground;
00685 double yhigh = ylow + getParameter(2);
00686
00687 double dPosition = getParameter(0);
00688 double dWidth = getParameter(1);
00689 double dHeight = getParameter(2);
00690
00691
00692
00693
00694
00695
00696
00697
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