Skip to content

Commit

Permalink
Checkpointing work. Refs #7506.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdzhou
wdzhou committed Jul 26, 2013
1 parent dc27fc8 commit aaeec83
Show file tree
Hide file tree
Showing 2 changed files with 149 additions and 76 deletions.
10 changes: 8 additions & 2 deletions Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FindPeaks.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,8 @@ class DLLExport FindPeaks : public API::Algorithm
/// Add the fit record (failure) to output workspace
void addNonFitRecord(const size_t spectrum);

void updateFitResults(API::IAlgorithm_sptr fitAlg, std::vector<double> &bestEffparams, std::vector<double> &bestRawparams, double &mincost, const double expPeakPos, const double expPeakHeight);
void updateFitResults(API::IFunction_sptr fitfunc, std::string fitstatus, double chi2, std::vector<double> &bestEffparams,
std::vector<double> &bestRawparams, double &mincost, const double expPeakPos, const double expPeakHeight);

API::IFunction_sptr createFunction(const double height, const double centre, const double sigma, const double a0, const double a1, const double a2, const bool withPeak = true);
int getBackgroundOrder();
Expand Down Expand Up @@ -187,6 +188,11 @@ class DLLExport FindPeaks : public API::Algorithm
void processFitResult(PeakFittingRecord& r1, PeakFittingRecord& r2, API::IPeakFunction_sptr peak, API::IFunction_sptr bkgdfunc, size_t spectrum,
size_t imin, size_t imax, double windowsize);

bool fitFunction(const API::MatrixWorkspace_sptr ws, size_t wsindex, API::IFunction_sptr func,
int maxiteration, double startx, double endx, std::string minimizer,
std::string& fitstatus, double& chi2);


/// Get best result from a set of fitting result
int getBestResult(std::vector<double> vecRwp);

Expand All @@ -206,7 +212,7 @@ class DLLExport FindPeaks : public API::Algorithm
//Properties saved in the algo.
API::MatrixWorkspace_sptr m_dataWS; ///<workspace to check for peaks
int m_inputPeakFWHM; ///<holder for the requested peak FWHM
int index; ///<list of workspace indicies to check
int m_wsIndex; ///<list of workspace indicies to check
bool singleSpectrum; ///<flag for if only a single spectrum is present
bool m_highBackground; ///<flag for find relatively weak peak in high background
bool m_rawPeaksTable; ///<flag for whether the output is the raw peak parameters or effective (centre, width, height)
Expand Down
215 changes: 141 additions & 74 deletions Code/Mantid/Framework/Algorithms/src/FindPeaks.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -333,14 +333,14 @@ namespace Algorithms
m_dataWS = getProperty("InputWorkspace");

// WorkspaceIndex
index = getProperty("WorkspaceIndex");
singleSpectrum = !isEmpty(index);
if (singleSpectrum && index >= static_cast<int>(m_dataWS->getNumberHistograms()))
m_wsIndex = getProperty("WorkspaceIndex");
singleSpectrum = !isEmpty(m_wsIndex);
if (singleSpectrum && m_wsIndex >= static_cast<int>(m_dataWS->getNumberHistograms()))
{
g_log.error() << "The value of WorkspaceIndex provided (" << index
g_log.error() << "The value of WorkspaceIndex provided (" << m_wsIndex
<< ") is larger than the size of this workspace (" << m_dataWS->getNumberHistograms()
<< ")\n";
throw Kernel::Exception::IndexError(index, m_dataWS->getNumberHistograms() - 1,
throw Kernel::Exception::IndexError(m_wsIndex, m_dataWS->getNumberHistograms() - 1,
"FindPeaks WorkspaceIndex property");
}

Expand Down Expand Up @@ -389,6 +389,13 @@ namespace Algorithms
// Minimum peak height
m_minHeight = getProperty("MinimumPeakHeight");

// Create a workspace for fit
size_t sizex = m_dataWS->readX(0).size();
size_t sizey = m_dataWS->readY(0).size();

API::MatrixWorkspace_sptr m_tempWS = boost::dynamic_pointer_cast<MatrixWorkspace>(
WorkspaceFactory::Instance().create("Workspace2D", 1, sizex, sizey));

return;
}

Expand All @@ -404,8 +411,8 @@ namespace Algorithms
std::size_t numPeaks = peakcentres.size();

// Loop over the spectra searching for peaks
const int start = singleSpectrum ? index : 0;
const int end = singleSpectrum ? index + 1 : static_cast<int>(m_dataWS->getNumberHistograms());
const int start = singleSpectrum ? m_wsIndex : 0;
const int end = singleSpectrum ? m_wsIndex + 1 : static_cast<int>(m_dataWS->getNumberHistograms());
m_progress = new Progress(this, 0.0, 1.0, end - start);

for (int spec = start; spec < end; ++spec)
Expand Down Expand Up @@ -476,8 +483,8 @@ namespace Algorithms
// setProperty("SmoothedData",smoothedData);

// Loop over the spectra searching for peaks
const int start = singleSpectrum ? index : 0;
const int end = singleSpectrum ? index + 1 : static_cast<int>(smoothedData->getNumberHistograms());
const int start = singleSpectrum ? m_wsIndex : 0;
const int end = singleSpectrum ? m_wsIndex + 1 : static_cast<int>(smoothedData->getNumberHistograms());
m_progress = new Progress(this, 0.0, 1.0, end - start);
const int blocksize = static_cast<int>(smoothedData->blocksize());

Expand Down Expand Up @@ -1070,7 +1077,20 @@ namespace Algorithms
for (int width = m_minGuessedPeakWidth; width <= m_maxGuessedPeakWidth; width +=
m_stepGuessedPeakWidth)
{
// Set up function to fit
int vecsize = static_cast<int>(vecX.size());
const int i_right = (i_centre + width < vecsize) ? i_centre + width : vecsize-1;
const int i_left = (i_centre - width >= 0) ? i_centre - width : 0;
const double in_sigma = vecX[i_right] - vecX[i_left];

// Create function to fit
IFunction_sptr peakfunc = createFunction(in_height, in_centre, in_sigma, in_bg0, in_bg1,
in_bg2, true);
g_log.debug() << " Function: " << peakfunc->asString() << "; Background Type = "
<< m_backgroundType << std::endl;

// Set up Child Algorithm Fit
#if 0
IAlgorithm_sptr fit;
try
{
Expand All @@ -1084,19 +1104,6 @@ namespace Algorithms
throw std::runtime_error(errorstr);
}

// Calculate the guessed sigma
// const double in_sigma = (i0 + width < vecX.size()) ? vecX[i0 + width] - vecX[i0] : 0.;
int vecsize = static_cast<int>(vecX.size());
const int i_right = (i_centre + width < vecsize) ? i_centre + width : vecsize-1;
const int i_left = (i_centre - width >= 0) ? i_centre - width : 0;
const double in_sigma = vecX[i_right] - vecX[i_left];

// Create function to fit
IFunction_sptr fitFunction = createFunction(in_height, in_centre, in_sigma, in_bg0, in_bg1,
in_bg2, true);
g_log.debug() << " Function: " << fitFunction->asString() << "; Background Type = "
<< m_backgroundType << std::endl;

// Set up Fit algorithm
fit->setProperty("Function", fitFunction);
fit->setProperty("InputWorkspace", input);
Expand All @@ -1109,8 +1116,19 @@ namespace Algorithms

// e) Fit and get result
fit->executeAsChildAlg();

updateFitResults(fit, bestparams, bestRawParams, mincost, in_centre, in_height);
updateFitResults(peakfit, bestparams, bestRawParams, mincost, in_centre, in_height);
#else
double startx = vecX[i_min];
double endx = vecX[i_max];
std::string fitstatus("");
double chi2;
bool isexecuted = fitFunction(input, spectrum, peakfunc, 50, startx, endx, m_minimizer, fitstatus, chi2);
#endif

if (isexecuted)
updateFitResults(peakfunc, fitstatus, chi2, bestparams, bestRawParams, mincost, in_centre, in_height);
else
g_log.warning("Fit is not executed. ");

} // ENDFOR: Loop over "width"

Expand All @@ -1122,10 +1140,10 @@ namespace Algorithms
addInfoRow(spectrum, bestparams, bestRawParams, mincost, true);

// Update collection of peaks
IFunction_sptr fitFunction = createFunction(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, true);
IFunction_sptr fitfunc = createFunction(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, true);
for (size_t i = 0; i < bestRawParams.size(); ++i)
fitFunction->setParameter(i, bestRawParams[i]);
addFittedFunction(fitFunction, i_min, i_max);
fitfunc->setParameter(i, bestRawParams[i]);
addFittedFunction(fitfunc, i_min, i_max);

return;
}
Expand Down Expand Up @@ -2172,27 +2190,22 @@ namespace Algorithms
* @param expPeakPos :: double as expected peak position
* @param expPeakHeight :: double as expected peak height
*/
void FindPeaks::updateFitResults(API::IAlgorithm_sptr fitAlg, std::vector<double> &bestEffparams,
void FindPeaks::updateFitResults(IFunction_sptr fitfunc, std::string fitstatus, double chi2, std::vector<double> &bestEffparams,
std::vector<double> &bestRawparams, double &mincost, const double expPeakPos,
const double expPeakHeight)
{

// check the results of the fit status
std::string fitStatus = fitAlg->getProperty("OutputStatus");
bool isfitgood = isFitSuccessful(fitStatus);
if (!isfitgood)
if (!isFitSuccessful(fitstatus))
return;

// check that chi2 got better
const double chi2 = fitAlg->getProperty("OutputChi2overDoF");
g_log.debug() << "Fit Status = " << fitStatus << ", chi2 = " << chi2 << std::endl;
if (chi2 > mincost)
{
return;
}

// get out the parameter names
std::vector<double> tempEffectiveParams, tempRawParams;
getComponentFunctions(fitAlg->getProperty("Function"), tempEffectiveParams, tempRawParams);
getComponentFunctions(fitfunc, tempEffectiveParams, tempRawParams);

// check the height
double height = tempEffectiveParams[2];
Expand Down Expand Up @@ -2494,6 +2507,7 @@ namespace Algorithms
<< startx << " EndX = " << endx << ".\n";

// Set up the background fitting
#if 0
IAlgorithm_sptr fit;
try
{
Expand All @@ -2518,42 +2532,37 @@ namespace Algorithms

// Execute fit and get result of fitting background
fit->executeAsChildAlg();
if (!fit->isExecuted())
#else
std::string fitstatus;
bool fitexecuted = fitFunction(bkgdWS, 0, m_backgroundFunction, 50, startx, endx, m_minimizer, fitstatus, chi2);
#endif

if (!fitexecuted)
{
g_log.error("Fit for background is not executed. ");
throw std::runtime_error("Fit for background is not executed. ");
}

std::string fitStatus = fit->getProperty("OutputStatus");
m_backgroundFunction = fit->getProperty("Function");

g_log.debug() << "(HighBackground) Fit Background Function. Fit Status = " << fitStatus
<< std::endl;
// m_backgroundFunction = bkgdfit->getProperty("Function");

// Check fiting status
bool allowedfailure = fitStatus.find("cannot") > 0 && fitStatus.find("tolerance") > 0;
bool allowedfailure = fitstatus.find("cannot") > 0 && fitstatus.find("tolerance") > 0;
g_log.debug() << "(HighBackground) Fit Background Function. Fit Status = " << fitstatus
<< ", Chi2 = " << chi2 << ".\n";

double bkgdchi2;
if (fitStatus.compare("success") == 0 || allowedfailure)
{
// good fit assumed
bkgdchi2 = fit->getProperty("OutputChi2overDoF");
}
else
if (fitstatus.compare("success") != 0 && !allowedfailure)
{
// fit status is completely failed (neither success nor failed in an allowed manner):
// set background to zero background
// Create background function
m_backgroundFunction->setParameter("A0", in_bg0);
if (numparams >= 2)
m_backgroundFunction->setParameter("A1", in_bg1);
if (numparams >= 3)
m_backgroundFunction->setParameter("A2", in_bg2);

bkgdchi2 = DBL_MAX;
// Chi2 to maximum
chi2 = DBL_MAX;
}

chi2 = bkgdchi2;

return true;
}

Expand Down Expand Up @@ -2593,6 +2602,7 @@ namespace Algorithms
init_rwp = calculateFunctionRwp(peakbkgdfunc, dataws, wsindex, startx, endx);

// Create Child Algorithm Fit
#if 0
IAlgorithm_sptr gfit;
try
{
Expand All @@ -2618,40 +2628,97 @@ namespace Algorithms

// Fit
gfit->executeAsChildAlg();
if (!gfit->isExecuted())
#else
std::string fitstatus("");
double chi2;
bool fitexecuted = fitFunction(dataws, wsindex, peakbkgdfunc, 50, startx, endx, m_minimizer, fitstatus, chi2);
#endif

double final_rwp(-DBL_MAX);
if (fitexecuted)
{
g_log.error("Fit is not executed correctly.");
return DBL_MAX;
// Fit is executed
g_log.debug("Fit is executed. ");

// Analyze result
bool isfitgood = isFitSuccessful(fitstatus);
if (isfitgood)
{
// calculate Rwp
final_rwp = calculateFunctionRwp(peakbkgdfunc, dataws, wsindex, startx, endx);

std::stringstream dbss;
std::vector<std::string> parnames = peakbkgdfunc->getParameterNames();
dbss << "[Fx357] Fit Peak (+background) Status = " << fitstatus << ". Starting Rwp = "
<< init_rwp << ". Final Rwp = " << final_rwp << ".\n";
for (size_t i = 0; i < parnames.size(); ++i)
dbss << parnames[i] << "\t = " << peakbkgdfunc->getParameter(parnames[i]) << "\n";
g_log.debug(dbss.str());
}
else
{
final_rwp = DBL_MAX;
}
}
else
{
g_log.debug("Fit is executed. ");
final_rwp = DBL_MAX;
}

// Analyze result
std::string fitpeakstatus = gfit->getProperty("OutputStatus");
bool isfitgood = isFitSuccessful(fitpeakstatus);
return final_rwp;
}

double final_rwp;
if (isfitgood)
bool FindPeaks::fitFunction(const MatrixWorkspace_sptr ws, size_t wsindex, IFunction_sptr func,
int maxiteration, double startx, double endx, std::string minimizer,
std::string& fitstatus, double& chi2)
{
#if 1
IAlgorithm_sptr fit;
try
{
final_rwp = calculateFunctionRwp(peakbkgdfunc, dataws, wsindex, startx, endx);
// Fitting the candidate peaks to a Gaussian
fit = createChildAlgorithm("Fit", -1, -1, true);
}
catch (Exception::NotFoundError &)
{
std::string errorstr("FindPeaks algorithm requires the CurveFitting library");
g_log.error(errorstr);
throw std::runtime_error(errorstr);
}

std::stringstream dbss;
// Set up Fit algorithm
fit->setProperty("Function", func);
fit->setProperty("InputWorkspace", ws);
fit->setProperty("WorkspaceIndex", static_cast<int>(wsindex));
fit->setProperty("MaxIterations", maxiteration);
fit->setProperty("StartX", startx);
fit->setProperty("EndX", endx);
fit->setProperty("Minimizer", minimizer);
fit->setProperty("CostFunction", "Least squares");

std::vector<std::string> parnames = peakbkgdfunc->getParameterNames();
dbss << "[Fx357] Fit Peak (+background) Status = " << fitpeakstatus << ". Starting Rwp = "
<< init_rwp << ". Final Rwp = " << final_rwp << ".\n";
for (size_t i = 0; i < parnames.size(); ++i)
dbss << parnames[i] << "\t = " << peakbkgdfunc->getParameter(parnames[i]) << "\n";
g_log.debug(dbss.str());
// e) Fit and get result
fit->executeAsChildAlg();

bool isexecuted = fit->isExecuted();
fitstatus = fit->getPropertyValue("OutputStatus");;
if (isexecuted)
{
chi2 = fit->getProperty("OutputChi2overDoF");
}
else
{
final_rwp = DBL_MAX;
chi2 = DBL_MAX;
}
g_log.debug() << "Fit Status = " << fitstatus << ", chi2 = " << chi2 << std::endl;
#else
// FIXME - This is a fake!!!
fitstatus = "success";
chi2 = 1.0;
bool isexecuted = true;
return final_rwp;
#endif

return isexecuted;
}


Expand Down

0 comments on commit aaeec83

Please sign in to comment.