diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/LeBailFit.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/LeBailFit.h index 498c8a8d06a16b74931a472f22ef64b1c6262575..cc4b2144a8adbc3e75a3a33921c0b855e56b933e 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/LeBailFit.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/LeBailFit.h @@ -222,7 +222,7 @@ namespace CurveFitting void addParameterToMCMinimize(vector<string>& parnamesforMC, string parname); /// Calculate diffraction pattern in Le Bail algorithm for MC Random walk - bool calculateDiffractionPatternMC(const MantidVec &vecX, const MantidVec &vecY, + bool calculateDiffractionPattern(const MantidVec &vecX, const MantidVec &vecY, bool inputraw, bool outputwithbkgd, MantidVec& vecBkgd, MantidVec& values, Rfactor& rfactor); @@ -266,7 +266,7 @@ namespace CurveFitting void storeBackgroundParameters(vector<double> &bkgdparamvec); /// Restore/recover the buffered background parameters to m_background function - void recoverBackgroundParameters(vector<double> bkgdparamvec); + void recoverBackgroundParameters(const vector<double> &bkgdparamvec); /// Propose new background parameters void proposeNewBackgroundValues(); @@ -408,7 +408,7 @@ namespace CurveFitting vector<string> m_bkgdParameterNames; size_t m_numberBkgdParameters; vector<double> m_bkgdParameterBuffer; - vector<double> m_bkgdParameterBest; + vector<double> m_bestBkgdParams; int m_roundBkgd; vector<double> m_bkgdParameterStepVec; diff --git a/Code/Mantid/Framework/CurveFitting/src/LeBailFit.cpp b/Code/Mantid/Framework/CurveFitting/src/LeBailFit.cpp index 0ded7fe2dd74d6d3387aa3bcaf69f50e9f6decd1..11ee2cf0f64d1d457baaca5ba877d5e10389dc7d 100644 --- a/Code/Mantid/Framework/CurveFitting/src/LeBailFit.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/LeBailFit.cpp @@ -28,6 +28,12 @@ as peak heights are calculated but not refined but highly correlated to backgrou See [[Le Bail Fit]]. *WIKI*/ +/* COMMIT NOTES * + 1. Rename calculateDiffractionPatternMC to calculateDiffractionPattern + 2. + + * COMMIT NOTES */ + #include "MantidCurveFitting/LeBailFit.h" #include "MantidKernel/ListValidator.h" #include "MantidAPI/TableRow.h" @@ -416,8 +422,8 @@ namespace CurveFitting // Calculate background vector<double> emptyvec; - bool resultphysical = calculateDiffractionPatternMC(m_dataWS->readX(m_wsIndex), m_dataWS->readY(m_wsIndex), - true, true, emptyvec, vecY, rfactor); + bool resultphysical = calculateDiffractionPattern(m_dataWS->readX(m_wsIndex), m_dataWS->readY(m_wsIndex), + true, true, emptyvec, vecY, rfactor); if (!resultphysical) { @@ -465,11 +471,11 @@ namespace CurveFitting */ void LeBailFit::execRefineBackground() { - // 0. Set up + // Set up class variables for background m_bkgdParameterNames = m_backgroundFunction->getParameterNames(); m_numberBkgdParameters = m_bkgdParameterNames.size(); m_bkgdParameterBuffer.resize(m_numberBkgdParameters); - m_bkgdParameterBest.resize(m_numberBkgdParameters); + m_bestBkgdParams.resize(m_numberBkgdParameters); m_roundBkgd = 0; m_bkgdParameterStepVec.resize(m_numberBkgdParameters, 0.01); for (size_t i = 1; i < m_numberBkgdParameters; ++i) @@ -478,8 +484,8 @@ namespace CurveFitting } // 1. Generate domain and value - const vector<double> vecX = m_dataWS->readX(m_wsIndex); - const vector<double> vecY = m_dataWS->readY(m_wsIndex); + const MantidVec& vecX = m_dataWS->readX(m_wsIndex); + const MantidVec& vecY = m_dataWS->readY(m_wsIndex); vector<double> valueVec(vecX.size(), 0); size_t numpts = vecX.size(); @@ -498,14 +504,14 @@ namespace CurveFitting } map<string, double> parammap = convertToDoubleMap(m_funcParameters); m_lebailFunction->setProfileParameterValues(parammap); - calculateDiffractionPatternMC(m_outputWS->readX(INPUTPUREPEAKINDEX), m_outputWS->readY(INPUTPUREPEAKINDEX), + calculateDiffractionPattern(m_outputWS->readX(INPUTPUREPEAKINDEX), m_outputWS->readY(INPUTPUREPEAKINDEX), false, true, backgroundvalues, valueVec, currR); Rfactor bestR = currR; - storeBackgroundParameters(m_bkgdParameterBest); + storeBackgroundParameters(m_bestBkgdParams); stringstream bufss; bufss << "Starting background parameter "; - for (size_t i = 0; i < m_bkgdParameterBest.size(); ++i) - bufss << "[" << i << "] = " << m_bkgdParameterBest[i] << ", "; + for (size_t i = 0; i < m_bestBkgdParams.size(); ++i) + bufss << "[" << i << "] = " << m_bestBkgdParams[i] << ", "; bufss << ". Starting Rwp = " << currR.Rwp; g_log.notice(bufss.str()); @@ -525,7 +531,7 @@ namespace CurveFitting } map<string, double> parammap = convertToDoubleMap(m_funcParameters); m_lebailFunction->setProfileParameterValues(parammap); - calculateDiffractionPatternMC(m_outputWS->readX(INPUTPUREPEAKINDEX), m_outputWS->readY(INPUTPUREPEAKINDEX), + calculateDiffractionPattern(m_outputWS->readX(INPUTPUREPEAKINDEX), m_outputWS->readY(INPUTPUREPEAKINDEX), false, true, backgroundvalues, valueVec, newR); g_log.information() << "[DBx800] New Rwp = " << newR.Rwp << ", Rp = " << newR.Rp << ".\n"; @@ -546,12 +552,12 @@ namespace CurveFitting { // Is it the best? bestR = newR; - storeBackgroundParameters(m_bkgdParameterBest); + storeBackgroundParameters(m_bestBkgdParams); stringstream bufss; bufss << "Temp best background parameter "; - for (size_t i = 0; i < m_bkgdParameterBest.size(); ++i) - bufss << "[" << i << "] = " << m_bkgdParameterBest[i] << ", "; + for (size_t i = 0; i < m_bestBkgdParams.size(); ++i) + bufss << "[" << i << "] = " << m_bestBkgdParams[i] << ", "; g_log.information(bufss.str()); } } @@ -561,7 +567,7 @@ namespace CurveFitting } // 3. Recover the best - recoverBackgroundParameters(m_bkgdParameterBest); + recoverBackgroundParameters(m_bestBkgdParams); stringstream bufss1; bufss1 << "Best background parameter "; @@ -578,7 +584,7 @@ namespace CurveFitting } parammap = convertToDoubleMap(m_funcParameters); m_lebailFunction->setProfileParameterValues(parammap); - calculateDiffractionPatternMC(m_outputWS->readX(INPUTPUREPEAKINDEX), m_outputWS->readY(INPUTPUREPEAKINDEX), + calculateDiffractionPattern(m_outputWS->readX(INPUTPUREPEAKINDEX), m_outputWS->readY(INPUTPUREPEAKINDEX), false, true, backgroundvalues, valueVec, outputR); g_log.notice() << "[DBx604] Best Rwp = " << bestR.Rwp << ", vs. recovered best Rwp = " << outputR.Rwp << ".\n"; @@ -635,7 +641,7 @@ namespace CurveFitting /** Restore/recover the buffered background parameters to m_background function * @param bkgdparamvec :: vector holding the background parameters whose order is same in background function */ - void LeBailFit::recoverBackgroundParameters(vector<double> bkgdparamvec) + void LeBailFit::recoverBackgroundParameters(const vector<double>& bkgdparamvec) { for (size_t i = 0; i < m_numberBkgdParameters; ++i) { @@ -751,7 +757,8 @@ namespace CurveFitting else { g_log.debug() << "DBx307: Cropped Workspace... Range From " << cropws->readX(wsindex)[0] << " To " - << cropws->readX(wsindex).back() << "\n"; + << cropws->readX(wsindex).back() << " of size " << cropws->readX(wsindex).size() + << "\n"; } return cropws; @@ -1489,7 +1496,7 @@ namespace CurveFitting map<string, double> pardblmap = convertToDoubleMap(parammap); m_lebailFunction->setProfileParameterValues(pardblmap); - bool startvaluevalid = calculateDiffractionPatternMC(vecX, vecPurePeak, false, false, vecBkgd, purepeakvalues, startR); + bool startvaluevalid = calculateDiffractionPattern(vecX, vecPurePeak, false, false, vecBkgd, purepeakvalues, startR); if (!startvaluevalid) { // Throw exception if starting values are not valid for all @@ -1542,7 +1549,7 @@ namespace CurveFitting // ii. Evaluate map<string, double> newpardblmap = convertToDoubleMap(newparammap); m_lebailFunction->setProfileParameterValues(newpardblmap); - bool validparams = calculateDiffractionPatternMC(vecX, vecPurePeak, false, false, vecBkgd, + bool validparams = calculateDiffractionPattern(vecX, vecPurePeak, false, false, vecBkgd, purepeakvalues, newR); g_log.information() << "[Calculation] Rwp = " << newR.Rwp << ", Rp = " << newR.Rp << ".\n"; @@ -1680,7 +1687,7 @@ namespace CurveFitting // c) Calculate again map<string, double> bestparams = convertToDoubleMap(m_bestParameters); m_lebailFunction->setProfileParameterValues(bestparams); - calculateDiffractionPatternMC(vecX, vecPurePeak, false, false, vecBkgd, purepeakvalues, currR); + calculateDiffractionPattern(vecX, vecPurePeak, false, false, vecBkgd, purepeakvalues, currR); MantidVec& vecCalY = m_outputWS->dataY(CALDATAINDEX); MantidVec& vecDiff = m_outputWS->dataY(DATADIFFINDEX); @@ -1961,10 +1968,10 @@ namespace CurveFitting * * @return :: boolean value. whether all the peaks' parameters are physical. */ - bool LeBailFit::calculateDiffractionPatternMC(const MantidVec& vecX, const MantidVec &vecY, - bool inputraw, bool outputwithbkgd, - MantidVec& vecBkgd, MantidVec& values, - Rfactor& rfactor) + bool LeBailFit::calculateDiffractionPattern(const MantidVec& vecX, const MantidVec &vecY, + bool inputraw, bool outputwithbkgd, + MantidVec& vecBkgd, MantidVec& values, + Rfactor& rfactor) { vector<double> veccalbkgd; @@ -2091,7 +2098,14 @@ namespace CurveFitting { // Find out the i-th parameter to be refined or not string paramname = mcgroup[i]; - Parameter param = curparammap[paramname]; +#if 0 + Parameter& param = curparammap[paramname]; +#else + map<string, Parameter>::iterator mapiter = curparammap.find(paramname); + if (mapiter == curparammap.end()) + throw runtime_error("Parameter to update is not in the pool of parameters to get updated."); + Parameter& param = mapiter->second; +#endif if (param.fit) anyparamtorefine = true; else @@ -2168,12 +2182,23 @@ namespace CurveFitting } // Apply to new parameter map +#if 0 newparammap[paramname].curvalue = newvalue; +#else + map<string, Parameter>::iterator newmiter = newparammap.find(paramname); + if (newmiter == newparammap.end()) + throw runtime_error("New parameter map does not contain parameter that is updated."); + newmiter->second.curvalue = newvalue; +#endif g_log.information() << "[ProposeNewValue] " << paramname << " --> " << newvalue << "; random number = " << randomnumber << "\n"; // g) record some trace +#if 0 Parameter& p = curparammap[paramname]; +#else + Parameter& p = param; +#endif if (stepsize > 0) { p.movedirection = 1; diff --git a/Code/Mantid/Framework/CurveFitting/src/LeBailFunction.cpp b/Code/Mantid/Framework/CurveFitting/src/LeBailFunction.cpp index 2669c2dc152d91735f179318eac5762b751122f2..ccc04b71d7974b7dea544e7890665d32e1442bd8 100644 --- a/Code/Mantid/Framework/CurveFitting/src/LeBailFunction.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/LeBailFunction.cpp @@ -386,7 +386,7 @@ namespace CurveFitting } else { - g_log.debug() << "[Fx155] Peaks group size = " << peakgroup.size() << "\n"; + g_log.information() << "[Fx155] Peaks group size = " << peakgroup.size() << "\n"; } if (peakgroup.size() > 1) sort(peakgroup.begin(), peakgroup.end()); @@ -486,8 +486,8 @@ namespace CurveFitting g_log.error(errmsg.str()); throw runtime_error(errmsg.str()); } - g_log.debug() << "[DBx356] Number of data points = " << ndata << " index from " << ileft - << " to " << iright << "; Size(datax, datay) = " << datax.size() << "\n"; + g_log.information() << "[DBx356] Number of data points = " << ndata << " index from " << ileft + << " to " << iright << "; Size(datax, datay) = " << datax.size() << "\n"; // Prepare to integrate dataY to calculate peak intensity vector<double> sumYs(ndata, 0.0); @@ -846,15 +846,20 @@ namespace CurveFitting } ++ ipk; - } + } // still in bound else { // Peak is get out of boundary + inbound = false; + g_log.information() << "[Fx301] Group peak: peak @ " << thispeak->centre() << " causes grouping " + << "peak over at maximum TOF = " << xmax << ".\n"; - // Right most peak in the boundary. Add current peak group vector to - vector<pair<double, IPowderDiffPeakFunction_sptr> > peakgroupcopy = peakgroup; - peakgroupvec.push_back(peakgroupcopy); - } + if (peakgroup.size() > 0) + { + vector<pair<double, IPowderDiffPeakFunction_sptr> > peakgroupcopy = peakgroup; + peakgroupvec.push_back(peakgroupcopy); + } + } // FIRST out of boundary } // ENDWHILE while (ipk < m_numPeaks) @@ -865,8 +870,8 @@ namespace CurveFitting ipk += 1; } - g_log.debug() << "[Calculate Peak Intensity]: Number of Peak Groups = " << peakgroupvec.size() - << "\n"; + g_log.information() << "[Calculate Peak Intensity]: Number of Peak Groups = " << peakgroupvec.size() + << "\n"; return; }