FittingFunction.cpp

Go to the documentation of this file.
00001 // FittingFunction.cpp: implementation of the FittingFunction class.
00002 //
00004 #ifdef WIN32 
00005 #pragma warning(disable:4786)
00006 #endif
00007 
00008 #include "FittingFunction.h"
00009 #include "FileWriter.h"
00010 #include <iostream>
00012 // Construction/Destruction
00014 
00015 
00016 //##ModelId=43849A03012E
00017 FittingFunction::FittingFunction()
00018 {
00019 }
00020 
00021 //##ModelId=43849A03012F
00022 FittingFunction::~FittingFunction()
00023 {
00024 
00025 }
00026 
00027 
00028 //##ModelId=43849A030150
00029 double FittingFunction::getParameter(std::string sParameterName) const{
00030 
00031         STRING2INT::const_iterator MyIterator;
00032         if((MyIterator = _mFunctionMap.find(sParameterName)) != _mFunctionMap.end()) {
00033                 return _pParameters[(*MyIterator).second];
00034         }
00035         else return 0;
00036 }
00037 
00038 //##ModelId=43849A03015D
00039 double FittingFunction::getParameterError(std::string sParameterName) const {
00040 
00041         STRING2INT::const_iterator MyIterator;
00042         if((MyIterator = _mFunctionMap.find(sParameterName)) != _mFunctionMap.end()) {
00043                 return _pParametersErrors[(*MyIterator).second];
00044         }
00045         else return 0;
00046 }
00047 
00048 //##ModelId=43849A030153
00049 double FittingFunction::getParameter(int iIndiceParameter) const{
00050         
00051         if ((0 <= iIndiceParameter) && (iIndiceParameter < _iNbParameters)) {
00052                 return _pParameters[iIndiceParameter];
00053         }
00054         return 0;
00055 }
00056 
00057 //##ModelId=43849A030156
00058 double* FittingFunction::getAllParameters() const{
00059         
00060         return _pParameters;
00061 }
00062 
00063 //##ModelId=43849A030160
00064 double  FittingFunction::getParameterError(int iIndiceParameter) const {
00065         
00066         if ((0 <= iIndiceParameter) && (iIndiceParameter < _iNbParameters)) {
00067                 return _pParametersErrors[iIndiceParameter];
00068         }
00069         return 0;
00070 }
00071 
00072 //##ModelId=43849A030163
00073 double* FittingFunction::getAllParametersErrors() const {
00074         return _pParametersErrors;
00075 }
00076 
00077 
00078 double FittingFunction::getStandardDeviation()
00079 {
00080         return _dStandardDeviation;
00081 }
00082 
00083 
00084 double          FittingFunction::computeStandardDeviation()
00085 {
00086 
00087         size_t lNbdata  = _cData->n;
00088 
00089         double dSum = 0.0;
00090         double dSumSquare = 0.0;
00091         double dXValue;
00092         for (unsigned int i=0;i<lNbdata;i++)
00093         {
00094                 dXValue         = *(_cData->x+i);
00095                 dSum            += dXValue;
00096                 dSumSquare      += dXValue*dXValue;
00097         }
00098 
00099         double dInvNbData = 1/lNbdata;
00100         
00101         double dSqrt = dInvNbData*dSumSquare - dInvNbData*dInvNbData*dSum*dSum;
00102         
00103         _dStandardDeviation = sqrt(dSqrt);
00104 
00105         return getStandardDeviation();
00106 
00107 }
00108 
00109 double          FittingFunction::computeMean(gsl_vector* vVector)
00110 {
00111         double dMean    = 0.0;
00112         double dSum             = 0.0;
00113         size_t dNbData  = vVector->size;
00114 
00115         for (unsigned int i=0;i<dNbData;i++)
00116         {
00117                 dSum+= gsl_vector_get(vVector,i);
00118         }
00119         dMean = dSum / dNbData;
00120         return dMean;
00121 }
00122 
00123 double          FittingFunction::computeVariance(gsl_vector* vVector)
00124 {
00125         double dMean            = computeMean(vVector);
00126         double dSum                     = 0.0;
00127         double dDiff            = 0.0;
00128         double dDeviation       = 0.0;
00129         size_t dNbData  = vVector->size;
00130 
00131         for (unsigned int i=0;i<dNbData;i++)
00132         {
00133                 dDiff=gsl_vector_get(vVector,i) - dMean;
00134                 dSum+= dDiff*dDiff;
00135         }
00136         dDeviation = dSum;// / dNbData;
00137         return dDeviation;
00138 }
00139 
00140 double          FittingFunction::computeResidual(gsl_vector* vVectorExp,gsl_vector* vVectorPred)
00141 {
00142         double dSum                     = 0.0;
00143         double dDiff            = 0.0;
00144         double dResidual        = 0.0;
00145         size_t dNbData  = vVectorExp->size;
00146 
00147         for (unsigned int i=0;i<dNbData;i++)
00148         {
00149                 dDiff=gsl_vector_get(vVectorExp,i) - gsl_vector_get(vVectorPred,i);
00150                 dSum+= dDiff*dDiff;
00151         }
00152         dResidual = dSum;// / dNbData;
00153         return dResidual;
00154 } 
00155 
00156 double          FittingFunction::computeDeterminationCoefficient(gsl_vector* vExperimentalDataY,gsl_vector* vFittedDataY)
00157 {
00158         double dDeviation       = computeVariance(vExperimentalDataY);
00159         double dResidual        = computeResidual(vExperimentalDataY,vFittedDataY);
00160 
00161         double dRSQ                     = 1 - dResidual/dDeviation;
00162         return dRSQ;
00163 }
00164 
00165 
00166 double          FittingFunction::computeObservedFStatistic(gsl_vector* vExperimentalDataY,gsl_vector* vFittedDataY)
00167 {
00168         double dDeviation               = computeVariance(vExperimentalDataY);
00169         double dResidual                = computeResidual(vExperimentalDataY,vFittedDataY);
00170         size_t lNbData                  = vExperimentalDataY->size;
00171 
00172         double dSSR                             = dDeviation - dResidual;
00173         double dNbParameters    = getNbParameters();
00174         double dDegreeOfFreedom = lNbData - 1 - dNbParameters;
00175         
00176         double dMSE                             = dResidual/dDegreeOfFreedom;
00177         double dMSR                             = dSSR/dNbParameters;
00178         double dFStatistic              = dMSR/dMSE;
00179 
00180         return dFStatistic;
00181 }
00182 
00183 
00184 int FittingFunction::getIndexValue(double dWantedValue,gsl_vector* vVector)
00185 {
00186         size_t  lVectorSize                     = vVector->size;
00187         double  lHalfVectorSize         = 0.5*lVectorSize;
00188 
00189         size_t iIndexMax                        = gsl_vector_max_index(vVector);
00190         double dMaxValue                        = gsl_vector_get(vVector, iIndexMax);
00191   
00192         double dCurrentValue            = dMaxValue;
00193         double dLowValue                        = dMaxValue;
00194         double dHighValue                       = dMaxValue;
00195 
00196         size_t  iIndexHigh                      = iIndexMax;
00197         size_t  iIndexLow                       = iIndexMax;
00198 
00199         size_t iCurrentIndex            = iIndexMax;
00200         size_t iGoodIndex                       = iIndexMax;
00201 
00202         double dLowStep                         = dLowValue;
00203         double dHighStep                        = dHighValue;
00204 
00205 
00206 
00207     if (iIndexMax < lHalfVectorSize)
00208     {
00209                 //- look in increasing direction
00210                 while(iCurrentIndex < lVectorSize - 1 && dCurrentValue > dWantedValue)
00211                 {
00212                         dCurrentValue   = gsl_vector_get(vVector, ++iCurrentIndex);
00213                         iIndexHigh              = iCurrentIndex;
00214                         iIndexLow               = iCurrentIndex - 1;
00215                         dLowValue               = gsl_vector_get(vVector, iIndexLow);
00216                         dHighValue              = gsl_vector_get(vVector, iIndexHigh);
00217         
00218 //                      cout << "icurrentIndex = " << iCurrentIndex << " | dlowvalue = " << dLowValue << " | dhighvalue = " << dHighValue << endl;
00219 
00220                 }
00221                 dLowStep                        = ::fabs(dLowValue - dWantedValue); 
00222                 dHighStep                       = ::fabs(dHighValue - dWantedValue); 
00223                 if (dLowStep < dHighStep) iGoodIndex = iIndexLow;
00224                 else iGoodIndex = iIndexHigh;
00225     }
00226     else
00227     {
00228                 //- look in decreasing direction
00229                 while( iCurrentIndex > 0 && dCurrentValue > dWantedValue)
00230                 {
00231                         dCurrentValue =  gsl_vector_get(vVector, --iCurrentIndex);
00232                         iIndexLow               = iCurrentIndex;
00233                         iIndexHigh              = iCurrentIndex + 1;
00234                         dLowValue               = gsl_vector_get(vVector, iIndexLow);
00235                         dHighValue              = gsl_vector_get(vVector, iIndexHigh);
00236 //                      cout << "icurrentIndex = " << iCurrentIndex << " | dlowvalue = " << dLowValue << " | dhighvalue = " << dHighValue << endl;
00237                 }
00238                 dLowStep                        = ::fabs(dLowValue - dWantedValue); 
00239                 dHighStep                       = ::fabs(dHighValue - dWantedValue); 
00240                 if (dLowStep < dHighStep) iGoodIndex = iIndexLow;
00241                 else iGoodIndex = iIndexHigh;
00242 
00243     }
00244 
00245     if (iCurrentIndex != iIndexMax)
00246     {
00247                 //cout << "\nFind nearest index for " << dWantedValue << " --> " << iGoodIndex << endl;
00248                 return (int)iGoodIndex;
00249                 
00250     }
00251     else
00252     {
00253                 //iGoodIndex = -1;
00254                 cout << "\nFind nearest index for " << dWantedValue << " --> " << -1 << endl;
00255                 return -1;
00256     }
00257 
00258         
00259 
00260         return -1;
00261 }
00262 
00263 
00264 
00265 /*void                  FittingFunction::WriteFittedDataInFile(const std::string sFilePath,const std::string sFileName,FittingData* pFittedData)
00266 {
00267         Interpolator::FileWriter  mFile(sFilePath,sFileName,_cFittedData.n,_cFittedData.x,cFittedData.y);
00268 }
00269 
00270 
00271 
00272 void                    FittingFunction::WriteExperimentalDataInFile(const std::string sFilePath,const std::string sFileName//,FittingData* pFittedData)
00273 {
00274         Interpolator::FileWriter  mFile(sFilePath,sFileName,_cEData.n,_cFittedData.x,cFittedData.y);
00275 }
00276 
00277 */
00278 
00279 
00280 /*
00281 
00282 double FittingFunction::estimateInitialPosition()
00283 {
00284         size_t max_index        = gsl_vector_max_index (vExperimentalDataY);
00285         double dPosition        = gsl_vector_get(vExperimentalDataX,max_index);
00286 
00287         return dPosition;
00288 }
00289 
00290 
00291 
00292 
00293 double FittingFunction::estimateInitialHeight()
00294 {
00295         size_t max_index        = gsl_vector_max_index (vExperimentalDataY);
00296 
00297         double background       = estimateInitialBackground();
00298         double maxValue         = gsl_vector_get(vExperimentalDataY, max_index);
00299         double height           = maxValue - background;
00300         return height;
00301 }
00302 
00303 
00304 double FittingFunction::estimateInitialBackground()
00305 {
00306         double dBackground = 0.0;
00307         if (vInitialsParameters->size == 4) 
00308         {       
00309                 size_t min_index = gsl_vector_min_index (vExperimentalDataY);
00310                 //TODO do a mean with more points to better estimate background
00311                 dBackground =  gsl_vector_get(vExperimentalDataY, min_index);
00312         }
00313         return dBackground;
00314 }
00315 
00316 
00317 double FittingFunction::estimateInitialSigma()
00318 {
00319         double dSigma;
00320         double dCovariance;
00321         double dHeight                          = estimateInitialHeight();
00322         double dBackground                      = estimateInitialBackground();
00323         
00324         size_t iIndexMax                        = gsl_vector_max_index(vExperimentalDataY);
00325         double dYMaxValue                       = gsl_vector_get(vExperimentalDataY, iIndexMax);
00326   
00327         double  dYCurrentValue          = dYMaxValue;
00328 
00329         size_t iCurrentIndex            = iIndexMax;
00330         size_t lVectorSize                      = lSizeX;
00331         double lHalfVectorSize          = 0.5 * lVectorSize;
00332         double dHalfMaximum                     = 0.5 * dHeight + dBackground;
00333   
00334     if (iIndexMax < lHalfVectorSize)
00335     {
00336                 //- look in increasing direction
00337                 while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue > dHalfMaximum)
00338                 {
00339                         dYCurrentValue = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00340                 }
00341     }
00342     else
00343     {
00344                 //- look in decreasing direction
00345                 while( iCurrentIndex > 0 && dYCurrentValue > dHalfMaximum)
00346                 {
00347                         dYCurrentValue =  gsl_vector_get(vExperimentalDataY, --iCurrentIndex);
00348                 }
00349     }
00350 
00351     if (iCurrentIndex != iIndexMax)
00352     {
00353                 double dPosition        = gsl_vector_get(vExperimentalDataX, iIndexMax);
00354                 double dX2                      = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00355                 dSigma = ::fabs(dX2 - dPosition);
00356     }
00357     else
00358     {
00359                 //- try an arbitrary value that should be not so far from reality
00360                 dSigma = lVectorSize / 6;
00361     }
00362 
00363     dCovariance = dSigma * dSigma;
00364         return dSigma;
00365 }
00366 
00367 
00368 
00369 
00370 double FittingFunction::estimateInitialPositionSigmoid()
00371 {
00372         double dHeight                          = estimateInitialHeight();
00373         double dBackground                      = estimateInitialBackground();
00374         double dHalfMaximum                     = 0.5 * dHeight + dBackground;
00375           
00376         size_t lVectorSize                      = lSizeX;
00377         size_t iCurrentIndex            = 0;
00378         double dYCurrentValue           = gsl_vector_get(vExperimentalDataY, iCurrentIndex);
00379 
00380         double dPosition;
00381 
00382         while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue < dHalfMaximum)
00383         {
00384                 dYCurrentValue = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00385         }
00386           
00387     if (iCurrentIndex != lVectorSize)
00388     {
00389                 dPosition                       = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00390     }
00391     else
00392     {
00393                 //- try an arbitrary value that should be not so far from reality
00394                 dPosition = gsl_vector_get(vExperimentalDataX, lVectorSize / 2);
00395     }
00396 
00397 
00398         return dPosition;
00399 }
00400 
00401 double FittingFunction::estimateInitialSigmaSigmoid()
00402 {
00403         double dSigma;
00404         double dCovariance;
00405         double dHeight                          = estimateInitialHeight();
00406         double dBackground                      = estimateInitialBackground();
00407         
00408         size_t iIndexMax                        = gsl_vector_max_index(vExperimentalDataY);
00409         double dYMaxValue                       = gsl_vector_get(vExperimentalDataY, iIndexMax);
00410   
00411         double  dYCurrentValue          = dYMaxValue;
00412 
00413         size_t iCurrentIndex            = iIndexMax;
00414         size_t lVectorSize                      = lSizeX;
00415         double lHalfVectorSize          = 0.5 * lVectorSize;
00416         double dHalfMaximum                     = 0.5 * dHeight + dBackground;
00417   
00418     if (iIndexMax < lHalfVectorSize)
00419     {
00420                 //- look in increasing direction
00421                 while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue > dHalfMaximum)
00422                 {
00423                         dYCurrentValue = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00424                 }
00425     }
00426     else
00427     {
00428                 //- look in decreasing direction
00429                 while( iCurrentIndex > 0 && dYCurrentValue > dHalfMaximum)
00430                 {
00431                         dYCurrentValue =  gsl_vector_get(vExperimentalDataY, --iCurrentIndex);
00432                 }
00433     }
00434 
00435     if (iCurrentIndex != iIndexMax)
00436     {
00437                 double dPosition        = gsl_vector_get(vExperimentalDataX, iIndexMax);
00438                 double dX2                      = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00439                 dSigma = ::fabs(dX2 - dPosition);
00440     }
00441     else
00442     {
00443                 //- try an arbitrary value that should be not so far from reality
00444                 dSigma = lVectorSize / 6;
00445     }
00446 
00447     dCovariance = dSigma * dSigma;
00448         return dSigma;
00449 }
00450 
00451 */
00452 
00453 

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