diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FitPeak.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FitPeak.h index 82b084da5617337eb94356a8aed3bedcf6bd7b55..bbc9576fd01cb385f07394ca997958a8bd092537 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FitPeak.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FitPeak.h @@ -47,6 +47,8 @@ namespace Algorithms void setFitPeakCriteria(bool usepeakpostol, double peakpostol); + std::string getDebugMessage(); + /// Fit peak and background together bool simpleFit(); @@ -208,6 +210,11 @@ namespace Algorithms /// Final goodness value (Rwp/Chi-square) double m_finalGoodnessValue; + /// + size_t m_numFitCalls; + + /// String stream + std::stringstream m_sstream; }; @@ -410,6 +417,9 @@ namespace Algorithms // std::map<std::string, double> m_fitErrorPeakFunc; // std::map<std::string, double> m_fitErrorBkgdFunc; + /// Option on output + bool m_lightWeightOutput; + }; diff --git a/Code/Mantid/Framework/Algorithms/src/FitPeak.cpp b/Code/Mantid/Framework/Algorithms/src/FitPeak.cpp index 9ba7ad82a60116c17bab206a6e0fe44675d2af93..3e912a6f3c8ca2dd282b9560171a952a51cdbdfb 100644 --- a/Code/Mantid/Framework/Algorithms/src/FitPeak.cpp +++ b/Code/Mantid/Framework/Algorithms/src/FitPeak.cpp @@ -67,7 +67,7 @@ namespace Algorithms */ FitOneSinglePeak::FitOneSinglePeak() : m_fitMethodSet(false), m_peakRangeSet(false), m_peakWidthSet(false), m_peakWindowSet(false), m_usePeakPositionTolerance(false), - m_wsIndex(0) + m_wsIndex(0), m_numFitCalls(0), m_sstream("") { } @@ -167,7 +167,24 @@ namespace Algorithms void FitOneSinglePeak::setFittingMethod(std::string minimizer, std::string costfunction) { m_minimizer = minimizer; - m_costFunction = costfunction; + if (costfunction == "Chi-Square") + { + m_costFunction = "Least squares"; + } + else if (costfunction == "Rwp") + { + m_costFunction = "Rwp"; + } + else if (costfunction == "Least squares") + { + m_costFunction = costfunction; + } + else + { + stringstream errss; + errss << "FitOneSinglePeak: cost function " << costfunction << " is not supported. "; + throw runtime_error(errss.str()); + } m_fitMethodSet = true; @@ -192,8 +209,10 @@ namespace Algorithms if (usrwidth <= 0) { // Set up default FWHM if user does not give reasonable peak width - g_log.warning() << "Client inputs user-defined peak width = " << usrwidth - << "; Automatically reset to 4 as default." << "\n"; + + m_sstream << "Client inputs user-defined peak width = " << usrwidth + << "; Automatically reset to 4 as default." << "\n"; + if (!fitwithsteppedfwhm) { fitwithsteppedfwhm = true; @@ -224,8 +243,8 @@ namespace Algorithms int i_centre = static_cast<int>(getVectorIndex(m_dataWS->readX(m_wsIndex), m_peakFunc->centre())); int i_maxindex = static_cast<int>(vecX.size())-1; - g_log.information() << "FWHM to guess. Range = " << minfwhm << ", " << maxfwhm - << "; Step = " << stepsize << "\n"; + m_sstream << "FWHM to guess. Range = " << minfwhm << ", " << maxfwhm + << "; Step = " << stepsize << "\n"; if (stepsize == 0 || maxfwhm < minfwhm) throw runtime_error("FWHM is not given right."); @@ -243,17 +262,17 @@ namespace Algorithms double in_fwhm = vecX[irightside] - vecX[ileftside]; if (in_fwhm < 1.0E-20) - { - g_log.warning() << "It is impossible to have zero peak width as iCentre = " - << i_centre << ", iWidth = " << iwidth << "\n" - << "More information: Spectrum = " << m_wsIndex << "; Range of X is " - << vecX.front() << ", " << vecX.back() - << "; Peak centre = " << vecX[i_centre]; + { + m_sstream << "It is impossible to have zero peak width as iCentre = " + << i_centre << ", iWidth = " << iwidth << "\n" + << "More information: Spectrum = " << m_wsIndex << "; Range of X is " + << vecX.front() << ", " << vecX.back() + << "; Peak centre = " << vecX[i_centre] << "\n"; } else { - g_log.debug() << "Fx330 i_width = " << iwidth << ", i_left = " << ileftside << ", i_right = " - << irightside << ", FWHM = " << in_fwhm << ".\n"; + m_sstream << "Fx330 i_width = " << iwidth << ", i_left = " << ileftside << ", i_right = " + << irightside << ", FWHM = " << in_fwhm << ".\n"; } m_vecFWHM.push_back(in_fwhm); @@ -313,11 +332,20 @@ namespace Algorithms return true; } + //---------------------------------------------------------------------------------------------- + /** Get debug message + */ + std::string FitOneSinglePeak::getDebugMessage() + { + return m_sstream.str(); + } + //---------------------------------------------------------------------------------------------- /** Fit peak with simple schemem */ bool FitOneSinglePeak::simpleFit() - { + { + m_numFitCalls = 0; string errmsg; if (!hasSetupToFitPeak(errmsg)) { @@ -333,7 +361,7 @@ namespace Algorithms compfunc->addFunction(m_peakFunc); compfunc->addFunction(m_bkgdFunc); - g_log.information() << "One-Step-Fit Function: " << compfunc->asString() << "\n"; + m_sstream << "One-Step-Fit Function: " << compfunc->asString() << "\n"; // Store starting setup push(m_peakFunc, m_bkupPeakFunc); @@ -344,7 +372,7 @@ namespace Algorithms for (size_t i = 0; i < numfits; ++i) { // set FWHM - g_log.debug() << "[SingleStepFit] FWHM = " << m_vecFWHM[i] << "\n"; + m_sstream << "[SingleStepFit] FWHM = " << m_vecFWHM[i] << "\n"; m_peakFunc->setFwhm(m_vecFWHM[i]); // fit and process result @@ -365,8 +393,9 @@ namespace Algorithms m_finalGoodnessValue = m_bestRwp; - g_log.information() << "One-Step-Fit Best (Chi^2 = " << m_bestRwp << ") Fitted Function: " - << compfunc->asString() << "\n"; + m_sstream << "One-Step-Fit Best (Chi^2 = " << m_bestRwp << ") Fitted Function: " + << compfunc->asString() << "\n" + << "Number of calls of Fit = " << m_numFitCalls << "\n"; return false; } @@ -409,8 +438,9 @@ namespace Algorithms double curpeakheight = svvalues[0]; // Get maximum peak value among - const MantidVec& vecX = dataws->readX(wsindex); + // const MantidVec& vecX = dataws->readX(wsindex); + const MantidVec& vecX = dataws->readX(wsindex); const MantidVec& vecY = dataws->readY(wsindex); double ymax = vecY[ixmin+1]; size_t iymax = ixmin+1; @@ -423,10 +453,11 @@ namespace Algorithms iymax = i; } } - g_log.debug() << "Estimate-Peak-Height: Current peak height = " << curpeakheight - << ". Estimate-Peak-Height: Maximum Y value between " << vecX[ixmin] << " and " - << vecX[ixmax] << " is " - << ymax << " at X = " << vecX[iymax] << ".\n"; + + m_sstream << "Estimate-Peak-Height: Current peak height = " << curpeakheight + << ". Estimate-Peak-Height: Maximum Y value between " << vecX[ixmin] << " and " + << vecX[ixmax] << " is " + << ymax << " at X = " << vecX[iymax] << ".\n"; // Compute peak height (not the maximum peak intensity) double estheight = ymax/curpeakheight*peakfunc->height(); @@ -475,12 +506,13 @@ namespace Algorithms // Check validity and debug output if (!peakfunc) throw std::runtime_error("fitPeakFunction's input peakfunc has not been initialized."); - else - g_log.debug() << "Function (to fit): " << peakfunc->asString() << " From " - << startx << " to " << endx << ".\n"; + + m_sstream << "Function (to fit): " << peakfunc->asString() << " From " + << startx << " to " << endx << ".\n"; double goodness = fitFunctionSD(peakfunc, dataws, wsindex, startx, endx, false); - g_log.debug() << "Peak parameter goodness-Fit = " << goodness << "\n"; + + m_sstream << "Peak parameter goodness-Fit = " << goodness << "\n"; return goodness; } @@ -494,6 +526,8 @@ namespace Algorithms */ bool FitOneSinglePeak::highBkgdFit() { + m_numFitCalls = 0; + // Check and initialization string errmsg; if (!hasSetupToFitPeak(errmsg)) @@ -501,6 +535,10 @@ namespace Algorithms g_log.error(errmsg); throw runtime_error("Object has not been set up completely to fit peak."); } + else + { + m_sstream << "F1158: Well-setup and good to go!\n"; + } m_bestRwp = DBL_MAX; @@ -529,8 +567,8 @@ namespace Algorithms // Set FWHM m_peakFunc->setFwhm(m_vecFWHM[i]); - g_log.debug() << "Round " << i << " of " << m_vecFWHM.size() << ". Using proposed FWHM = " - << m_vecFWHM[i] << "\n"; + m_sstream << "Round " << i << " of " << m_vecFWHM.size() << ". Using proposed FWHM = " + << m_vecFWHM[i] << "\n"; // Fit double rwp = fitPeakFunction(m_peakFunc, purePeakWS, 0, m_minFitX, m_maxFitX); @@ -546,8 +584,9 @@ namespace Algorithms m_finalFitGoodness = fitCompositeFunction(m_peakFunc, m_bkgdFunc, m_dataWS, m_wsIndex, m_minFitX, m_maxFitX); - g_log.information() << "MultStep-Fit: Best Fitted Peak: " << m_peakFunc->asString() - << "Final " << m_costFunction << " = " << m_finalFitGoodness << "\n"; + m_sstream << "MultStep-Fit: Best Fitted Peak: " << m_peakFunc->asString() + << "Final " << m_costFunction << " = " << m_finalFitGoodness << "\n" + << "Number of calls on Fit = " << m_numFitCalls << "\n"; return false; } @@ -646,7 +685,6 @@ namespace Algorithms IAlgorithm_sptr fit; try { - fit = createChildAlgorithm("Fit", -1, -1, true); } catch (Exception::NotFoundError &) @@ -657,6 +695,8 @@ namespace Algorithms throw std::runtime_error(errss.str()); } + m_sstream << "F1200 Cost function is " << m_costFunction << "\n"; + // Set the properties fit->setProperty("Function", fitfunc); fit->setProperty("InputWorkspace", dataws); @@ -669,8 +709,8 @@ namespace Algorithms fit->setProperty("CalcErrors", true); // Execute fit and get result of fitting background - g_log.debug() << "FitSingleDomain: Fit " << fit->asString() << "; StartX = " << xmin - << ", EndX = " << xmax << ".\n"; + m_sstream << "FitSingleDomain: Fit " << fit->asString() << "; StartX = " << xmin + << ", EndX = " << xmax << ".\n"; fit->executeAsChildAlg(); if (!fit->isExecuted()) @@ -678,6 +718,7 @@ namespace Algorithms g_log.error("Fit for background is not executed. "); throw std::runtime_error("Fit for background is not executed. "); } + ++ m_numFitCalls; // Retrieve result std::string fitStatus = fit->getProperty("OutputStatus"); @@ -696,9 +737,9 @@ namespace Algorithms } // Debug information - g_log.debug() << "FitSingleDomain Fitted-Function " << fitfunc->asString() - << ": Fit-status = " << fitStatus - << ", chi^2 = " << chi2 << ".\n"; + m_sstream << "[F1201] FitSingleDomain Fitted-Function " << fitfunc->asString() + << ": Fit-status = " << fitStatus + << ", chi^2 = " << chi2 << ".\n"; return chi2; } @@ -759,7 +800,7 @@ namespace Algorithms fit->setProperty("Minimizer", m_minimizer); fit->setProperty("CostFunction", "Least squares"); - g_log.debug() << "FitMultiDomain: Funcion " << funcmd->asString() << "\n"; + m_sstream << "FitMultiDomain: Funcion " << funcmd->asString() << "\n"; // Execute fit->execute(); @@ -767,17 +808,18 @@ namespace Algorithms { throw runtime_error("Fit is not executed on multi-domain function/data. "); } + ++ m_numFitCalls; // Retrieve result std::string fitStatus = fit->getProperty("OutputStatus"); - g_log.debug() << "[DB] Multi-domain fit status: " << fitStatus << ".\n"; + m_sstream << "[DB] Multi-domain fit status: " << fitStatus << ".\n"; double chi2 = EMPTY_DBL(); if (fitStatus == "success") { chi2 = fit->getProperty("OutputChi2overDoF"); - g_log.debug() << "FitMultidomain: Successfully-Fitted Function " <<fitfunc->asString() - << ", Chi^2 = "<< chi2 << "\n"; + m_sstream << "FitMultidomain: Successfully-Fitted Function " <<fitfunc->asString() + << ", Chi^2 = "<< chi2 << "\n"; } return chi2; @@ -805,7 +847,7 @@ namespace Algorithms // Do calculation for starting chi^2/Rwp bool modecal = true; double goodness_init = fitFunctionSD(compfunc, dataws, wsindex, startx, endx, modecal); - g_log.debug() << "Peak+Backgruond: Pre-fit Goodness = " << goodness_init << "\n"; + m_sstream << "Peak+Backgruond: Pre-fit Goodness = " << goodness_init << "\n"; map<string, double> bkuppeakmap, bkupbkgdmap; push(peakfunc, bkuppeakmap); @@ -820,8 +862,9 @@ namespace Algorithms // Check fit result goodness = checkFittedPeak(peakfunc, goodness, errorreason); + if (errorreason.size() > 0) - g_log.information() << "Error reason: " << errorreason << "\n"; + m_sstream << "Error reason: " << errorreason << "\n"; double goodness_final = DBL_MAX; if (goodness <= goodness_init) @@ -834,15 +877,17 @@ namespace Algorithms { // A worse result is got. Revert to original function parameters goodness_final = goodness_init; - g_log.information() << "Fit peak/background composite function FAILS to render a better solution. " - << "Input cost function value = " << goodness_init << ", output cost function value = " - << goodness << "\n"; + + m_sstream << "Fit peak/background composite function FAILS to render a better solution. " + << "Input cost function value = " << goodness_init << ", output cost function value = " + << goodness << "\n"; + pop(bkuppeakmap, peakfunc); pop(bkupbkgdmap, bkgdfunc); } else { - g_log.error("Fit peak-background function fails in all approaches! "); + m_sstream << "Fit peak-background function fails in all approaches! \n"; } return goodness_final; @@ -975,8 +1020,8 @@ namespace Algorithms fitsuccess = false; } - g_log.debug() << "Process and Store is called. " << "Rwp = " << rwp << ", best Rwp = " - << m_bestRwp << ", Fit success = " << fitsuccess << "\n"; + m_sstream << "Process and Store is called. " << "Rwp = " << rwp << ", best Rwp = " + << m_bestRwp << ", Fit success = " << fitsuccess << "\n"; // Store result if if (rwp < m_bestRwp && fitsuccess) @@ -992,7 +1037,7 @@ namespace Algorithms } else if (!fitsuccess) { - g_log.debug() << "Reason of fit's failure: " << failreason << "\n"; + m_sstream << "Reason of fit's failure: " << failreason << "\n"; } return; @@ -1205,6 +1250,9 @@ namespace Algorithms { fit1peakalg.simpleFit(); } + string dbmsg = fit1peakalg.getDebugMessage(); + g_log.information(dbmsg); + m_finalGoodnessValue = fit1peakalg.getFitCostFunctionValue(); map<string, double> peakfuncfiterrormap, bkgdfuncfiterrormap; @@ -1390,7 +1438,7 @@ namespace Algorithms // Generate background function m_bkgdFunc = boost::dynamic_pointer_cast<IBackgroundFunction>( FunctionFactory::Instance().createFunction(bkgdtype)); - g_log.debug() << "Created background function of type " << bkgdtype << "\n"; + // g_log.debug() << "Created background function of type " << bkgdtype << "\n"; // Set background function parameter values m_bkgdParameterNames = getProperty("BackgroundParameterNames"); @@ -1429,7 +1477,7 @@ namespace Algorithms string peaktype = parseFunctionTypeFull(peaktypeprev, defaultparorder); m_peakFunc = boost::dynamic_pointer_cast<IPeakFunction>( FunctionFactory::Instance().createFunction(peaktype)); - g_log.debug() << "Create peak function of type " << peaktype << "\n"; + // g_log.debug() << "Create peak function of type " << peaktype << "\n"; // Peak parameters' names m_peakParameterNames = getProperty("PeakParameterNames"); @@ -1464,7 +1512,7 @@ namespace Algorithms m_peakFunc->setParameter(m_peakParameterNames[i], vec_peakparvalues[i]); } - g_log.information() << "Created peak function " << m_peakFunc->asString() << "\n"; + // g_log.information() << "Created peak function " << m_peakFunc->asString() << "\n"; return; } @@ -1504,9 +1552,11 @@ namespace Algorithms // Check validity on peak centre double centre_guess = m_peakFunc->centre(); +#if 0 g_log.debug() << "Fit Peak with given window: Guessed center = " << centre_guess << " x-min = " << m_minFitX << ", x-max = " << m_maxFitX << "\n"; +#endif if (m_minFitX >= centre_guess || m_maxFitX <= centre_guess) { g_log.error("Peak centre is out side of fit window."); @@ -1530,7 +1580,7 @@ namespace Algorithms size_t i_minFitX = getVectorIndex(m_dataWS->readX(m_wsIndex), m_minFitX); size_t i_maxFitX = getVectorIndex(m_dataWS->readX(m_wsIndex), m_maxFitX); - g_log.debug() << "[DB] i_min/max Fit X = " << i_minFitX << ", " << i_maxFitX << "\n"; + // g_log.debug() << "[DB] i_min/max Fit X = " << i_minFitX << ", " << i_maxFitX << "\n"; // Data workspace size_t nspec = 3; diff --git a/Code/Mantid/Framework/Algorithms/test/FitPeakTest.h b/Code/Mantid/Framework/Algorithms/test/FitPeakTest.h index 3f694f8f649985645a5c969f18a32bc23cac4d4a..29a371a9e311a2b8cb322d84354db5b60d94a200 100644 --- a/Code/Mantid/Framework/Algorithms/test/FitPeakTest.h +++ b/Code/Mantid/Framework/Algorithms/test/FitPeakTest.h @@ -34,7 +34,7 @@ public: //---------------------------------------------------------------------------------------------- /** Test on init and setup */ - void test_Init() + void Ptest_Init() { // Generate input workspace MatrixWorkspace_sptr dataws = gen_4866P5Data(); @@ -79,7 +79,7 @@ public: //---------------------------------------------------------------------------------------------- /** Test on fit a peak with significantly high background */ - void test_FitPeakWithHighBkgd() + void Ptest_FitPeakWithHighBkgd() { // Generate input workspace MatrixWorkspace_sptr dataws = gen_4866P5Data(); @@ -114,6 +114,7 @@ public: fitpeak.setProperty("MinGuessedPeakWidth", 2); fitpeak.setProperty("MaxGuessedPeakWidth", 20); fitpeak.setProperty("GuessedPeakWidthStep", 2); + // fitpeak.setProperty("OutputFitFunctionOnly", false); // Execute fitpeak.execute();