Skip to content

Commit

Permalink
refs #10384 added couple of unit tests and perf test for new Algo
Browse files Browse the repository at this point in the history 
and couple of minor bugs discovered when running these tests.
  • Loading branch information
@abuts
abuts committed Oct 31, 2014
1 parent d9f646a commit f8e1aeb
Show file tree
Hide file tree
Showing 3 changed files with 212 additions and 70 deletions.
2 changes: 1 addition & 1 deletion Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RemoveBackground.h
View file
Expand Up @@ -45,7 +45,7 @@ namespace Mantid

void initialize(const API::MatrixWorkspace_const_sptr &bkgWS,const API::MatrixWorkspace_sptr &sourceWS,
int emode,Kernel::Logger *pLog=NULL,int nTreads=1,bool inPlace=true);
void removeBackground(int hist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data,
void removeBackground(int hist,MantidVec &XValues,MantidVec &y_data,MantidVec &e_data,
int tread_num=0)const;

private:
Expand Down
17 changes: 11 additions & 6 deletions Code/Mantid/Framework/Algorithms/src/RemoveBackground.cpp
View file
Expand Up @@ -113,8 +113,9 @@ namespace Mantid
for (int hist=0; hist < histnumber;++hist)
{
PARALLEL_START_INTERUPT_REGION
// get const references to input Workspace arrays (no copying)
const MantidVec& XValues = outputWS->readX(hist);
// get references to output Workspace X-arrays.
MantidVec& XValues = outputWS->dataX(hist);
// get references to output workspace data and error. If it is new workspace, data will be copied there, if old, modified in-place
MantidVec& YValues = outputWS->dataY(hist);
MantidVec& YErrors = outputWS->dataE(hist);

Expand Down Expand Up @@ -228,7 +229,7 @@ namespace Mantid
* @param e_data -- the spectra errors
* @param threadNum -- number of thread doing conversion (by default 0, single thread)
*/
void BackgroundHelper::removeBackground(int nHist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data,int threadNum)const
void BackgroundHelper::removeBackground(int nHist,MantidVec &x_data,MantidVec &y_data,MantidVec &e_data,int threadNum)const
{

double dtBg,IBg;
Expand Down Expand Up @@ -256,6 +257,7 @@ namespace Mantid
double L2 = detector->getDistance(*m_Sample);
double delta(std::numeric_limits<double>::quiet_NaN());
// get access to source workspace in case if target is different from source
const MantidVec& XValues = m_wkWS->readX(nHist);
const MantidVec& YValues = m_wkWS->readY(nHist);
const MantidVec& YErrors = m_wkWS->readE(nHist);

Expand All @@ -264,10 +266,12 @@ namespace Mantid
unitConv->initialize(m_L1, L2,twoTheta, m_Emode, m_Efix,delta);


double tof1 = unitConv->singleToTOF(XValues[0]);
x_data[0] = XValues[0];
double tof1 = unitConv->singleToTOF(x_data[0]);
for(size_t i=0;i<y_data.size();i++)
{
double tof2=unitConv->singleToTOF(XValues[i+1]);
double X=XValues[i+1];
double tof2=unitConv->singleToTOF(X);
double Jack = std::fabs((tof2-tof1)/dtBg);
double normBkgrnd = IBg*Jack;
tof1=tof2;
Expand All @@ -280,8 +284,9 @@ namespace Mantid
}
else
{
x_data[i+1] =X;
y_data[i] =YValues[i]-normBkgrnd;
e_data[i] =std::sqrt((normBkgrnd+YErrors[i]*YErrors[i])/2.); // needs further clarification -- Gaussian error summation?
e_data[i] =std::sqrt((normBkgrnd+YErrors[i]*YErrors[i])/2.);
}

}
Expand Down
263 changes: 200 additions & 63 deletions Code/Mantid/Framework/Algorithms/test/RemoveBackgroundTest.h
View file
Expand Up @@ -11,70 +11,69 @@
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
using namespace Mantid;

class RemoveBackgroundTest : public CxxTest::TestSuite
{
public:

// This pair of boilerplate methods prevent the suite being created statically
// This means the constructor isn't called when running other tests
static RemoveBackgroundTest *createSuite() { return new RemoveBackgroundTest(); }
static void destroySuite( RemoveBackgroundTest *suite ) { delete suite; }
void init_workspaces(int nHist,int nBins, API::MatrixWorkspace_sptr &BgWS,API::MatrixWorkspace_sptr &SourceWS)
{
DataObjects::Workspace2D_sptr theWS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(nHist, nBins);
// Add incident energy necessary for unit conversion
theWS->mutableRun().addProperty("Ei",13.,"meV",true);

RemoveBackgroundTest()
{
API::AnalysisDataService::Instance().addOrReplace("sourceWS",theWS);

DataObjects::Workspace2D_sptr theWS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(1, 15000);
Algorithms::Rebin rebinner;
rebinner.initialize();
rebinner.setPropertyValue("InputWorkspace",theWS->getName());
rebinner.setPropertyValue("OutputWorkspace","Background");
rebinner.setPropertyValue("Params","10000,5000,15000");

std::string wsName = theWS->getName();
if(wsName.empty())
{
wsName = "sourceWS";
}
// Add incident energy necessary for unit conversion
theWS->mutableRun().addProperty("Ei",13.,"meV",true);
rebinner.execute();

API::AnalysisDataService::Instance().addOrReplace(wsName,theWS);
BgWS = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("Background");

Algorithms::Rebin rebinner;
rebinner.initialize();
rebinner.setPropertyValue("InputWorkspace",theWS->getName());
rebinner.setPropertyValue("OutputWorkspace","Background");
rebinner.setPropertyValue("Params","10000,5000,15000");
Algorithms::ConvertUnits unitsConv;
unitsConv.initialize();
unitsConv.setPropertyValue("InputWorkspace",theWS->getName());
unitsConv.setPropertyValue("OutputWorkspace","sourceWSdE");
unitsConv.setPropertyValue("Target","DeltaE");
unitsConv.setPropertyValue("EMode","Direct");

rebinner.execute();
unitsConv.execute();

BgWS = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("Background");
Algorithms::CalculateFlatBackground bgRemoval;

Algorithms::ConvertUnits unitsConv;
unitsConv.initialize();
unitsConv.setPropertyValue("InputWorkspace",theWS->getName());
unitsConv.setPropertyValue("OutputWorkspace","sourceWSdE");
unitsConv.setPropertyValue("Target","DeltaE");
unitsConv.setPropertyValue("EMode","Direct");
bgRemoval.initialize();
bgRemoval.setPropertyValue("InputWorkspace",theWS->getName());
bgRemoval.setPropertyValue("OutputWorkspace",theWS->getName());
bgRemoval.setPropertyValue("StartX","10000");
bgRemoval.setPropertyValue("EndX","15000");
bgRemoval.setPropertyValue("Mode","Mean");

unitsConv.execute();
bgRemoval.execute();

Algorithms::CalculateFlatBackground bgRemoval;
unitsConv.setPropertyValue("InputWorkspace",theWS->getName());
unitsConv.setPropertyValue("OutputWorkspace","sampleWSdE");
unitsConv.setPropertyValue("Target","DeltaE");
unitsConv.setPropertyValue("EMode","Direct");

bgRemoval.initialize();
bgRemoval.setPropertyValue("InputWorkspace",theWS->getName());
bgRemoval.setPropertyValue("OutputWorkspace",theWS->getName());
bgRemoval.setPropertyValue("StartX","10000");
bgRemoval.setPropertyValue("EndX","15000");
bgRemoval.setPropertyValue("Mode","Mean");
unitsConv.execute();

bgRemoval.execute();
SourceWS = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("sourceWSdE");

unitsConv.setPropertyValue("InputWorkspace",theWS->getName());
unitsConv.setPropertyValue("OutputWorkspace","sampleWSdE");
unitsConv.setPropertyValue("Target","DeltaE");
unitsConv.setPropertyValue("EMode","Direct");

unitsConv.execute();
}

SourceWS = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("sourceWSdE");
class RemoveBackgroundTest : public CxxTest::TestSuite
{
public:

// This pair of boilerplate methods prevent the suite being created statically
// This means the constructor isn't called when running other tests
static RemoveBackgroundTest *createSuite() { return new RemoveBackgroundTest(); }
static void destroySuite( RemoveBackgroundTest *suite ) { delete suite; }

RemoveBackgroundTest()
{
init_workspaces(1,15000,BgWS,SourceWS);
}

~RemoveBackgroundTest()
Expand All @@ -84,33 +83,37 @@ class RemoveBackgroundTest : public CxxTest::TestSuite
}
void testWrongInit()
{
Algorithms::BackgroundHelper bgRemoval;
// create workspace with units of energy transfer
auto bkgWS = WorkspaceCreationHelper::createProcessedInelasticWS(std::vector<double>(1,1.), std::vector<double>(1,20.), std::vector<double>(1,10.));
TSM_ASSERT_THROWS("Should throw if background workspace is not in TOF units",bgRemoval.initialize(bkgWS,SourceWS,0),std::invalid_argument);
Algorithms::BackgroundHelper bgRemoval;
// create workspace with units of energy transfer
auto bkgWS = WorkspaceCreationHelper::createProcessedInelasticWS(std::vector<double>(1,1.), std::vector<double>(1,20.), std::vector<double>(1,10.));
TSM_ASSERT_THROWS("Should throw if background workspace is not in TOF units",bgRemoval.initialize(bkgWS,SourceWS,0),std::invalid_argument);


bkgWS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(2, 15);
TSM_ASSERT_THROWS("Should throw if background is not 1 or equal to source",bgRemoval.initialize(bkgWS,SourceWS,0),std::invalid_argument);
bkgWS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(2, 15);
TSM_ASSERT_THROWS("Should throw if background is not 1 or equal to source",bgRemoval.initialize(bkgWS,SourceWS,0),std::invalid_argument);

auto sourceWS = WorkspaceCreationHelper::Create2DWorkspace(5,10);
TSM_ASSERT_THROWS("Should throw if source workspace does not have units",bgRemoval.initialize(BgWS,sourceWS,0),std::invalid_argument);
auto sourceWS = WorkspaceCreationHelper::Create2DWorkspace(5,10);
TSM_ASSERT_THROWS("Should throw if source workspace does not have units",bgRemoval.initialize(BgWS,sourceWS,0),std::invalid_argument);

sourceWS ->getAxis(0)->setUnit("TOF");
TSM_ASSERT_THROWS("Should throw if source workspace does not have proper instrument",bgRemoval.initialize(BgWS,sourceWS,0),std::invalid_argument);
sourceWS ->getAxis(0)->setUnit("TOF");
TSM_ASSERT_THROWS("Should throw if source workspace does not have proper instrument",bgRemoval.initialize(BgWS,sourceWS,0),std::invalid_argument);
}


void testBackgroundHelper()
{
Algorithms::BackgroundHelper bgRemoval;
auto clone = cloneSourceWS();

API::AnalysisDataService::Instance().addOrReplace("TestWS",clone);


int emode = static_cast<int>(Kernel::DeltaEMode().fromString("Direct"));
bgRemoval.initialize(BgWS,SourceWS,emode);

MantidVec& dataX = SourceWS->dataX(0);
MantidVec& dataY = SourceWS->dataY(0);
MantidVec& dataE = SourceWS->dataE(0);
MantidVec& dataX = clone->dataX(0);
MantidVec& dataY = clone->dataY(0);
MantidVec& dataE = clone->dataE(0);

bgRemoval.removeBackground(0,dataX,dataY,dataE);

Expand All @@ -126,13 +129,147 @@ class RemoveBackgroundTest : public CxxTest::TestSuite
}

}
void testRemoveBkgInPlace()
{
auto clone = cloneSourceWS();
API::AnalysisDataService::Instance().addOrReplace("TestWS",clone);

Algorithms::RemoveBackground bkgRem;
bkgRem.initialize();
bkgRem.setPropertyValue("InputWorkspace","TestWS");
bkgRem.setPropertyValue("OutputWorkspace","TestWS");
bkgRem.setPropertyValue("FlatBkgWorkspace",BgWS->getName());
bkgRem.setPropertyValue("EMode","Direct");

TS_ASSERT_THROWS_NOTHING(bkgRem.execute());


auto SampleWS = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("sampleWSdE");
auto result = clone;

const MantidVec & sampleX = SampleWS->readX(0);
const MantidVec & sampleY = SampleWS->readY(0);

const MantidVec & resultX = result->readX(0);
const MantidVec & resultY = result->readY(0);

//const MantidVec & sampleE = SampleWS->readE(0);
for(size_t i=0;i<sampleY.size();i++)
{
TS_ASSERT_DELTA(resultX[i],sampleX[i],1.e-7);
TS_ASSERT_DELTA(resultY[i],sampleY[i],1.e-7);
}
}

void testRemoveBkgNewRez()
{
auto clone = cloneSourceWS();
API::AnalysisDataService::Instance().addOrReplace("TestWS",clone);

Algorithms::RemoveBackground bkgRem;
bkgRem.initialize();
bkgRem.setPropertyValue("InputWorkspace","TestWS");
bkgRem.setPropertyValue("OutputWorkspace","TestWS2");
bkgRem.setPropertyValue("FlatBkgWorkspace",BgWS->getName());
bkgRem.setPropertyValue("EMode","Direct");

TS_ASSERT_THROWS_NOTHING(bkgRem.execute());


auto SampleWS = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("sampleWSdE");
auto result = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("TestWS2");

const MantidVec & sampleX = SampleWS->readX(0);
const MantidVec & sampleY = SampleWS->readY(0);

const MantidVec & resultX = result->readX(0);
const MantidVec & resultY = result->readY(0);

//const MantidVec & sampleE = SampleWS->readE(0);
for(size_t i=0;i<sampleY.size();i++)
{
TS_ASSERT_DELTA(resultX[i],sampleX[i],1.e-7);
TS_ASSERT_DELTA(resultY[i],sampleY[i],1.e-7);
}
}

private:
API::MatrixWorkspace_sptr cloneSourceWS()
{
auto cloneWS = API::WorkspaceFactory::Instance().create(SourceWS);

const MantidVec &X = SourceWS->readX(0);
const MantidVec &Y = SourceWS->readY(0);
const MantidVec &E = SourceWS->readE(0);
cloneWS->setX(0,X);
cloneWS->getSpectrum(0)->setData(Y,E);

return cloneWS;
}

API::MatrixWorkspace_sptr BgWS;
API::MatrixWorkspace_sptr SourceWS;

};


class RemoveBackgroundTestPerformance : public CxxTest::TestSuite
{
public:

// This pair of boilerplate methods prevent the suite being created statically
// This means the constructor isn't called when running other xests
static RemoveBackgroundTestPerformance *createSuite() { return new RemoveBackgroundTestPerformance (); }
static void destroySuite( RemoveBackgroundTestPerformance *suite ) { delete suite; }

RemoveBackgroundTestPerformance()
{

init_workspaces(1000,15000,BgWS,SourceWS);
}

void testRemoveBkgInPlace()
{

Algorithms::RemoveBackground bkgRem;
bkgRem.initialize();
bkgRem.setPropertyValue("InputWorkspace","sourceWSdE");
bkgRem.setPropertyValue("OutputWorkspace","sourceWSdE");
bkgRem.setPropertyValue("FlatBkgWorkspace",BgWS->getName());
bkgRem.setPropertyValue("EMode","Direct");

TS_ASSERT_THROWS_NOTHING(bkgRem.execute());


auto SampleWS = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("sampleWSdE");
auto result = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("sourceWSdE");

size_t spectra[]={0,10,100,999};
std::vector<size_t> list_to_check(spectra,spectra+4);


for(size_t i=0;i<list_to_check.size();i++)
{
const MantidVec & sampleX = SampleWS->readX(list_to_check[i]);
const MantidVec & sampleY = SampleWS->readY(list_to_check[i]);

const MantidVec & resultX = result->readX(list_to_check[i]);
const MantidVec & resultY = result->readY(list_to_check[i]);


//const MantidVec & sampleE = SampleWS->readE(0);
for(size_t i=0;i<sampleY.size();i++)
{
TS_ASSERT_DELTA(resultX[i],sampleX[i],1.e-7);
TS_ASSERT_DELTA(resultY[i],sampleY[i],1.e-7);
}
}
}

private:
API::MatrixWorkspace_sptr BgWS;
API::MatrixWorkspace_sptr SourceWS;
};



#endif /*ALGORITHMTEST_H_*/
#endif /*REMOVE_BACKGROUD_TEST_H_*/

0 comments on commit f8e1aeb

Please sign in to comment.