Skip to content

Commit

Permalink
Refs #10623. Changed PoldiFitPeaks1D
Browse files Browse the repository at this point in the history 
The algorithm now tries to find out which Chebyshev-polynomial (as background) gives the best fit.
  • Loading branch information
Michael Wedel committed Nov 26, 2014
1 parent eea77dc commit 9e04b33
Show file tree
Hide file tree
Showing 3 changed files with 144 additions and 104 deletions.
18 changes: 11 additions & 7 deletions Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks1D.h
View file
Expand Up @@ -84,32 +84,36 @@ namespace Poldi
virtual const std::string category() const;

protected:
PoldiPeakCollection_sptr fitPeaks(const PoldiPeakCollection_sptr &peaks);

int getBestChebyshevPolynomialDegree(const DataObjects::Workspace2D_sptr &dataWorkspace, const RefinedRange_sptr &range);

PoldiPeakCollection_sptr getReducedPeakCollection(const PoldiPeakCollection_sptr &peaks) const;
bool peakIsAcceptable(const PoldiPeak_sptr &peak) const;

void setPeakFunction(const std::string &peakFunction);
PoldiPeakCollection_sptr getInitializedPeakCollection(const DataObjects::TableWorkspace_sptr &peakTable) const;

std::vector<RefinedRange_sptr> getRefinedRanges(const PoldiPeakCollection_sptr &peaks) const;
std::vector<RefinedRange_sptr> getReducedRanges(const std::vector<RefinedRange_sptr> &ranges) const;

API::IFunction_sptr getRangeProfile(const RefinedRange_sptr &range, int n) const;
API::IFunction_sptr getPeakProfile(const PoldiPeak_sptr &poldiPeak) const;
void setValuesFromProfileFunction(PoldiPeak_sptr poldiPeak, const API::IFunction_sptr &fittedFunction) const;
double getFwhmWidthRelation(API::IPeakFunction_sptr peakFunction) const;

API::IAlgorithm_sptr getFitAlgorithm(const DataObjects::Workspace2D_sptr &dataWorkspace, const PoldiPeak_sptr &peak, const API::IFunction_sptr &profile);

void addPeakFitCharacteristics(const API::ITableWorkspace_sptr &fitResult);
void initializeFitResultWorkspace(const API::ITableWorkspace_sptr &fitResult);
API::IAlgorithm_sptr getFitAlgorithm(const DataObjects::Workspace2D_sptr &dataWorkspace, const RefinedRange_sptr &range, int n);

void initializePeakResultWorkspace(const DataObjects::TableWorkspace_sptr &peakResultWorkspace) const;
void storePeakResult(API::TableRow tableRow, const PoldiPeak_sptr &peak) const;
DataObjects::TableWorkspace_sptr generateResultTable(const PoldiPeakCollection_sptr &peaks) const;

PoldiPeakCollection_sptr m_peaks;
std::string m_profileTemplate;
API::IFunction_sptr m_backgroundTemplate;
std::string m_profileTies;

DataObjects::TableWorkspace_sptr m_fitCharacteristics;
DataObjects::TableWorkspace_sptr m_peakResultOutput;
API::WorkspaceGroup_sptr m_fitplots;


double m_fwhmMultiples;

Expand Down
216 changes: 129 additions & 87 deletions Code/Mantid/Framework/SINQ/src/PoldiFitPeaks1D.cpp
View file
Expand Up @@ -16,22 +16,25 @@
#include "MantidAPI/CompositeFunction.h"
#include <boost/make_shared.hpp>

#include <functional>
#include <boost/math/distributions/normal.hpp>


namespace Mantid
{
namespace Poldi
{

using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::CurveFitting;
using namespace Kernel;
using namespace API;
using namespace DataObjects;
using namespace CurveFitting;

RefinedRange::RefinedRange(const PoldiPeak_sptr &peak, double fwhmMultiples) :
m_peaks(1, peak)
{
double width = peak->fwhm();
double extent = std::min(0.05, std::max(0.002, width)) * fwhmMultiples;
double extent = std::max(0.002, width) * fwhmMultiples;

m_xStart = peak->q() - extent;
m_xEnd = peak->q() + extent;
Expand Down Expand Up @@ -86,10 +89,8 @@ DECLARE_ALGORITHM(PoldiFitPeaks1D)
PoldiFitPeaks1D::PoldiFitPeaks1D() :
m_peaks(),
m_profileTemplate(),
m_backgroundTemplate(),
m_profileTies(),
m_fitCharacteristics(),
m_peakResultOutput(),
m_fitplots(new WorkspaceGroup),
m_fwhmMultiples(1.0)
{

Expand Down Expand Up @@ -124,11 +125,7 @@ void PoldiFitPeaks1D::init()

declareProperty(new WorkspaceProperty<TableWorkspace>("OutputWorkspace","RefinedPeakTable",Direction::Output), "Output workspace with refined peak data.");
declareProperty(new WorkspaceProperty<TableWorkspace>("ResultTableWorkspace","ResultTable",Direction::Output), "Fit results.");
declareProperty(new WorkspaceProperty<TableWorkspace>("FitCharacteristicsWorkspace","FitCharacteristics",Direction::Output), "Fit characteristics for each peak.");
declareProperty(new WorkspaceProperty<Workspace>("FitPlotsWorkspace","FitPlots",Direction::Output), "Plots of all peak fits.");

m_backgroundTemplate = FunctionFactory::Instance().createInitialized("name=UserFunction, Formula=A0 + A1*(x - x0)^2");
m_profileTies = "f1.x0 = f0.PeakCentre";
}

void PoldiFitPeaks1D::setPeakFunction(const std::string &peakFunction)
Expand Down Expand Up @@ -176,38 +173,40 @@ std::vector<RefinedRange_sptr> PoldiFitPeaks1D::getReducedRanges(const std::vect
return reducedRanges;
}

IFunction_sptr PoldiFitPeaks1D::getPeakProfile(const PoldiPeak_sptr &poldiPeak) const {
IPeakFunction_sptr clonedProfile = boost::dynamic_pointer_cast<IPeakFunction>(FunctionFactory::Instance().createFunction(m_profileTemplate));
clonedProfile->setCentre(poldiPeak->q());
clonedProfile->setFwhm(poldiPeak->fwhm(PoldiPeak::AbsoluteQ));
clonedProfile->setHeight(poldiPeak->intensity());

IFunction_sptr clonedBackground = m_backgroundTemplate->clone();

API::IFunction_sptr PoldiFitPeaks1D::getRangeProfile(const RefinedRange_sptr &range, int n) const
{
CompositeFunction_sptr totalProfile(new CompositeFunction);
totalProfile->initialize();
totalProfile->addFunction(clonedProfile);
totalProfile->addFunction(clonedBackground);

if(!m_profileTies.empty()) {
totalProfile->addTies(m_profileTies);
std::vector<PoldiPeak_sptr> peaks = range->getPeaks();
for(auto it = peaks.begin(); it != peaks.end(); ++it) {
totalProfile->addFunction(getPeakProfile(*it));
}

totalProfile->addFunction(FunctionFactory::Instance().createInitialized("name=Chebyshev,n=" + boost::lexical_cast<std::string>(n)
+ ",StartX=" + boost::lexical_cast<std::string>(range->getXStart())
+ ",EndX=" + boost::lexical_cast<std::string>(range->getXEnd())));

return totalProfile;
}

IFunction_sptr PoldiFitPeaks1D::getPeakProfile(const PoldiPeak_sptr &poldiPeak) const {
IPeakFunction_sptr clonedProfile = boost::dynamic_pointer_cast<IPeakFunction>(FunctionFactory::Instance().createFunction(m_profileTemplate));
clonedProfile->setCentre(poldiPeak->q());
clonedProfile->setFwhm(poldiPeak->fwhm(PoldiPeak::AbsoluteQ));
clonedProfile->setHeight(poldiPeak->intensity());

return clonedProfile;
}

void PoldiFitPeaks1D::setValuesFromProfileFunction(PoldiPeak_sptr poldiPeak, const IFunction_sptr &fittedFunction) const
{
CompositeFunction_sptr totalFunction = boost::dynamic_pointer_cast<CompositeFunction>(fittedFunction);

if(totalFunction) {
IPeakFunction_sptr peakFunction = boost::dynamic_pointer_cast<IPeakFunction>(totalFunction->getFunction(0));
IPeakFunction_sptr peakFunction = boost::dynamic_pointer_cast<IPeakFunction>(fittedFunction);

if(peakFunction) {
poldiPeak->setIntensity(UncertainValue(peakFunction->height(), peakFunction->getError(0)));
poldiPeak->setQ(UncertainValue(peakFunction->centre(), peakFunction->getError(1)));
poldiPeak->setFwhm(UncertainValue(peakFunction->fwhm(), getFwhmWidthRelation(peakFunction) * peakFunction->getError(2)));
}
if(peakFunction) {
poldiPeak->setIntensity(UncertainValue(peakFunction->height(), peakFunction->getError(0)));
poldiPeak->setQ(UncertainValue(peakFunction->centre(), peakFunction->getError(1)));
poldiPeak->setFwhm(UncertainValue(peakFunction->fwhm(), getFwhmWidthRelation(peakFunction) * peakFunction->getError(2)));
}
}

Expand All @@ -216,95 +215,138 @@ double PoldiFitPeaks1D::getFwhmWidthRelation(IPeakFunction_sptr peakFunction) co
return peakFunction->fwhm() / peakFunction->getParameter(2);
}

void PoldiFitPeaks1D::exec()
PoldiPeakCollection_sptr PoldiFitPeaks1D::fitPeaks(const PoldiPeakCollection_sptr &peaks)
{
setPeakFunction(getProperty("PeakFunction"));

// Number of points around the peak center to use for the fit
m_fwhmMultiples = getProperty("FwhmMultiples");

// try to construct PoldiPeakCollection from provided TableWorkspace
TableWorkspace_sptr poldiPeakTable = getProperty("PoldiPeakTable");
m_peaks = getInitializedPeakCollection(poldiPeakTable);
g_log.information() << "Peaks to fit: " << peaks->peakCount() << std::endl;

g_log.information() << "Peaks to fit: " << m_peaks->peakCount() << std::endl;

std::vector<RefinedRange_sptr> rawRanges = getRefinedRanges(m_peaks);
std::vector<RefinedRange_sptr> rawRanges = getRefinedRanges(peaks);
std::vector<RefinedRange_sptr> reducedRanges = getReducedRanges(rawRanges);

g_log.information() << "Ranges used for fitting: " << reducedRanges.size() << std::endl;

Workspace2D_sptr dataWorkspace = getProperty("InputWorkspace");
m_fitplots->removeAll();

for(size_t i = 0; i < reducedRanges.size(); ++i) {
RefinedRange_sptr currentRange = reducedRanges[i];
int nMin = getBestChebyshevPolynomialDegree(dataWorkspace, currentRange);

if(nMin > -1) {
IAlgorithm_sptr fit = getFitAlgorithm(dataWorkspace, currentRange, nMin);
fit->execute();

IFunction_sptr fitFunction = fit->getProperty("Function");
CompositeFunction_sptr composite = boost::dynamic_pointer_cast<CompositeFunction>(fitFunction);

if(!composite) {
throw std::runtime_error("Not a composite function!");
}

m_fitCharacteristics = boost::dynamic_pointer_cast<TableWorkspace>(WorkspaceFactory::Instance().createTable());
WorkspaceGroup_sptr fitPlotGroup(new WorkspaceGroup);
std::vector<PoldiPeak_sptr> peaks = currentRange->getPeaks();
for(size_t i = 0; i < peaks.size(); ++i) {
setValuesFromProfileFunction(peaks[i], composite->getFunction(i));
MatrixWorkspace_sptr fpg = fit->getProperty("OutputWorkspace");
m_fitplots->addWorkspace(fpg);
}
}
}

for(size_t i = 0; i < m_peaks->peakCount(); ++i) {
PoldiPeak_sptr currentPeak = m_peaks->peak(i);
IFunction_sptr currentProfile = getPeakProfile(currentPeak);
return getReducedPeakCollection(peaks);
}

IAlgorithm_sptr fit = getFitAlgorithm(dataWorkspace, currentPeak, currentProfile);
int PoldiFitPeaks1D::getBestChebyshevPolynomialDegree(const Workspace2D_sptr &dataWorkspace, const RefinedRange_sptr &range)
{
int n = 0;
double chiSquareMin = 1e10;
int nMin = -1;

while((n < 3)) {
IAlgorithm_sptr fit = getFitAlgorithm(dataWorkspace, range, n);
bool fitSuccess = fit->execute();

if(fitSuccess) {
setValuesFromProfileFunction(currentPeak, fit->getProperty("Function"));
addPeakFitCharacteristics(fit->getProperty("OutputParameters"));
ITableWorkspace_sptr fitCharacteristics = fit->getProperty("OutputParameters");
TableRow row = fitCharacteristics->getRow(fitCharacteristics->rowCount() - 1);

double chiSquare = row.Double(1);

MatrixWorkspace_sptr fpg = fit->getProperty("OutputWorkspace");
fitPlotGroup->addWorkspace(fpg);

if(fabs(chiSquare - 1) < fabs(chiSquareMin - 1)) {
chiSquareMin = chiSquare;
nMin = n;
}
}

++n;
}

m_peakResultOutput = generateResultTable(m_peaks);
g_log.information() << "Chi^2 for range [" << range->getXStart() << " - " << range->getXEnd() << "] is minimal at n = " << nMin << " with Chi^2 = " << chiSquareMin << std::endl;

setProperty("OutputWorkspace", m_peaks->asTableWorkspace());
setProperty("FitCharacteristicsWorkspace", m_fitCharacteristics);
setProperty("ResultTableWorkspace", m_peakResultOutput);
setProperty("FitPlotsWorkspace", fitPlotGroup);
return nMin;
}

IAlgorithm_sptr PoldiFitPeaks1D::getFitAlgorithm(const Workspace2D_sptr &dataWorkspace, const PoldiPeak_sptr &peak, const IFunction_sptr &profile)
PoldiPeakCollection_sptr PoldiFitPeaks1D::getReducedPeakCollection(const PoldiPeakCollection_sptr &peaks) const
{
double width = peak->fwhm();
double extent = std::min(0.05, std::max(0.002, width)) * m_fwhmMultiples;
PoldiPeakCollection_sptr reducedPeaks = boost::make_shared<PoldiPeakCollection>();
reducedPeaks->setProfileFunctionName(peaks->getProfileFunctionName());

std::pair<double, double> xBorders(peak->q() - extent, peak->q() + extent);
for(size_t i = 0; i < peaks->peakCount(); ++i) {
PoldiPeak_sptr currentPeak = peaks->peak(i);

IAlgorithm_sptr fitAlgorithm = createChildAlgorithm("Fit", -1, -1, false);
fitAlgorithm->setProperty("CreateOutput", true);
fitAlgorithm->setProperty("Output", "FitPeaks1D");
fitAlgorithm->setProperty("CalcErrors", true);
fitAlgorithm->setProperty("Function", profile);
fitAlgorithm->setProperty("InputWorkspace", dataWorkspace);
fitAlgorithm->setProperty("WorkspaceIndex", 0);
fitAlgorithm->setProperty("StartX", xBorders.first);
fitAlgorithm->setProperty("EndX", xBorders.second);
if(peakIsAcceptable(currentPeak)) {
reducedPeaks->addPeak(currentPeak);
}
}

return fitAlgorithm;
return reducedPeaks;
}

void PoldiFitPeaks1D::addPeakFitCharacteristics(const ITableWorkspace_sptr &fitResult)
bool PoldiFitPeaks1D::peakIsAcceptable(const PoldiPeak_sptr &peak) const
{
if(m_fitCharacteristics->columnCount() == 0) {
initializeFitResultWorkspace(fitResult);
}
return peak->intensity() > 0 && peak->fwhm(PoldiPeak::Relative) < 0.02;
}

void PoldiFitPeaks1D::exec()
{
setPeakFunction(getProperty("PeakFunction"));

TableRow newRow = m_fitCharacteristics->appendRow();
// Number of points around the peak center to use for the fit
m_fwhmMultiples = getProperty("FwhmMultiples");

for(size_t i = 0; i < fitResult->rowCount(); ++i) {
TableRow currentRow = fitResult->getRow(i);
// try to construct PoldiPeakCollection from provided TableWorkspace
TableWorkspace_sptr poldiPeakTable = getProperty("PoldiPeakTable");
m_peaks = getInitializedPeakCollection(poldiPeakTable);

newRow << UncertainValueIO::toString(UncertainValue(currentRow.Double(1), currentRow.Double(2)));
PoldiPeakCollection_sptr fittedPeaksNew = fitPeaks(m_peaks);
PoldiPeakCollection_sptr fittedPeaksOld = m_peaks;
while(fittedPeaksNew->peakCount() < fittedPeaksOld->peakCount()) {
fittedPeaksOld = fittedPeaksNew;
fittedPeaksNew = fitPeaks(fittedPeaksOld);
}

m_peakResultOutput = generateResultTable(fittedPeaksNew);

setProperty("OutputWorkspace", m_peaks->asTableWorkspace());
setProperty("ResultTableWorkspace", m_peakResultOutput);
setProperty("FitPlotsWorkspace", m_fitplots);
}

void PoldiFitPeaks1D::initializeFitResultWorkspace(const API::ITableWorkspace_sptr &fitResult)
IAlgorithm_sptr PoldiFitPeaks1D::getFitAlgorithm(const Workspace2D_sptr &dataWorkspace, const RefinedRange_sptr &range, int n)
{
for(size_t i = 0; i < fitResult->rowCount(); ++i) {
TableRow currentRow = fitResult->getRow(i);
m_fitCharacteristics->addColumn("str", currentRow.cell<std::string>(0));
}
IFunction_sptr rangeProfile = getRangeProfile(range, n);

IAlgorithm_sptr fitAlgorithm = createChildAlgorithm("Fit", -1, -1, false);
fitAlgorithm->setProperty("CreateOutput", true);
fitAlgorithm->setProperty("Output", "FitPeaks1D");
fitAlgorithm->setProperty("CalcErrors", true);
fitAlgorithm->setProperty("OutputCompositeMembers", true);
fitAlgorithm->setProperty("Function", rangeProfile);
fitAlgorithm->setProperty("InputWorkspace", dataWorkspace);
fitAlgorithm->setProperty("WorkspaceIndex", 0);
fitAlgorithm->setProperty("StartX", range->getXStart());
fitAlgorithm->setProperty("EndX", range->getXEnd());

return fitAlgorithm;
}

void PoldiFitPeaks1D::initializePeakResultWorkspace(const DataObjects::TableWorkspace_sptr &peakResultWorkspace) const
Expand Down
14 changes: 4 additions & 10 deletions Code/Mantid/Framework/SINQ/test/PoldiFitPeaks1DTest.h
View file
Expand Up @@ -59,18 +59,13 @@ class PoldiFitPeaks1DTest : public CxxTest::TestSuite
void testGetPeakProfile()
{
TestablePoldiFitPeaks1D poldiFitPeaks;
poldiFitPeaks.m_backgroundTemplate = m_backgroundTestFunction;
poldiFitPeaks.initialize();
poldiFitPeaks.setPeakFunction(m_profileTestFunction);

IFunction_sptr totalProfile = poldiFitPeaks.getPeakProfile(m_testPeak);

// make sure that we get back a composite of peak and background
CompositeFunction_sptr composite = boost::dynamic_pointer_cast<CompositeFunction>(totalProfile);
TS_ASSERT(composite);
IFunction_sptr peakFunction = poldiFitPeaks.getPeakProfile(m_testPeak);

// make sure that the profile is the first function in the composite
IPeakFunction_sptr profile = boost::dynamic_pointer_cast<IPeakFunction>(composite->getFunction(0));
// make sure that the profile is correct
IPeakFunction_sptr profile = boost::dynamic_pointer_cast<IPeakFunction>(peakFunction);
TS_ASSERT(profile);

TS_ASSERT_EQUALS(profile->centre(), m_testPeak->q());
Expand Down Expand Up @@ -100,7 +95,7 @@ class PoldiFitPeaks1DTest : public CxxTest::TestSuite
Mantid::Poldi::PoldiFitPeaks1D fitPeaks1D;
fitPeaks1D.initialize();

TS_ASSERT_EQUALS(fitPeaks1D.propertyCount(), 8);
TS_ASSERT_EQUALS(fitPeaks1D.propertyCount(), 7);

std::vector<Property *> properties = fitPeaks1D.getProperties();
std::set<std::string> names;
Expand All @@ -115,7 +110,6 @@ class PoldiFitPeaks1DTest : public CxxTest::TestSuite
TS_ASSERT_EQUALS(names.count("PoldiPeakTable"), 1);
TS_ASSERT_EQUALS(names.count("OutputWorkspace"), 1);
TS_ASSERT_EQUALS(names.count("ResultTableWorkspace"), 1);
TS_ASSERT_EQUALS(names.count("FitCharacteristicsWorkspace"), 1);
TS_ASSERT_EQUALS(names.count("FitPlotsWorkspace"), 1);
}

Expand Down

0 comments on commit 9e04b33

Please sign in to comment.