Commit 0fd1c3dc authored by Gigg, Martyn Anthony's avatar Gigg, Martyn Anthony
Browse files

Add a class that is able to produce faked MD data.

Refs #11209
parent 93b83c3b
......@@ -12,6 +12,7 @@ set ( SRC_FILES
src/EventWorkspaceHelpers.cpp
src/EventWorkspaceMRU.cpp
src/Events.cpp
src/FakeMD.cpp
src/GroupingWorkspace.cpp
src/Histogram1D.cpp
src/MDBin.cpp
......@@ -70,6 +71,7 @@ set ( INC_FILES
inc/MantidDataObjects/EventWorkspaceHelpers.h
inc/MantidDataObjects/EventWorkspaceMRU.h
inc/MantidDataObjects/Events.h
inc/MantidDataObjects/FakeMD.h
inc/MantidDataObjects/GroupingWorkspace.h
inc/MantidDataObjects/Histogram1D.h
inc/MantidDataObjects/MDBin.h
......@@ -123,6 +125,7 @@ set ( TEST_FILES
EventWorkspaceMRUTest.h
EventWorkspaceTest.h
EventsTest.h
FakeMDTest.h
GroupingWorkspaceTest.h
Histogram1DTest.h
LibraryManagerTest.h
......@@ -133,9 +136,9 @@ set ( TEST_FILES
MDBoxSaveableTest.h
MDBoxTest.h
MDDimensionStatsTest.h
MDEventTest.h
MDEventFactoryTest.h
MDEventInserterTest.h
MDEventTest.h
MDEventWorkspaceTest.h
MDGridBoxTest.h
MDHistoWorkspaceIteratorTest.h
......
#ifndef MANTID_DATAOBJECTS_FAKEMD_H_
#define MANTID_DATAOBJECTS_FAKEMD_H_
#include <vector>
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidDataObjects/DllConfig.h"
#include "MantidDataObjects/MDEventWorkspace.h"
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/variate_generator.hpp>
namespace Mantid {
namespace DataObjects {
/**
Provides a helper class to add fake data to an MD workspace
Copyright &copy; 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
National Laboratory & European Spallation Source
This file is part of Mantid.
Mantid is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
Mantid is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
File change history is stored at: <https://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class MANTID_DATAOBJECTS_DLL FakeMD {
public:
FakeMD(const std::vector<double> &uniformParams,
const std::vector<double> &peakParams, const int randomSeed,
const bool randomizeSignal);
void fill(API::IMDEventWorkspace_sptr workspace);
private:
void setupDetectorCache(const API::IMDEventWorkspace &workspace);
template <typename MDE, size_t nd>
void addFakePeak(typename MDEventWorkspace<MDE, nd>::sptr ws);
template <typename MDE, size_t nd>
void addFakeUniformData(typename MDEventWorkspace<MDE, nd>::sptr ws);
template <typename MDE, size_t nd>
void addFakeRandomData(const std::vector<double> &params,
typename MDEventWorkspace<MDE, nd>::sptr ws);
template <typename MDE, size_t nd>
void addFakeRegularData(const std::vector<double> &params,
typename MDEventWorkspace<MDE, nd>::sptr ws);
detid_t pickDetectorID();
//------------------ Member variables ------------------------------------
std::vector<double> m_uniformParams;
std::vector<double> m_peakParams;
const int m_randomSeed;
const bool m_randomizeSignal;
mutable std::vector<detid_t> m_detIDs;
boost::mt19937 m_randGen;
boost::uniform_int<size_t> m_uniformDist;
};
} // namespace DataObjects
} // namespace Mantid
#endif /* MANTID_DATAOBJECTS_FAKEMD_H_ */
//--------------------------------------------------------------------------------------------------
// Includes
//--------------------------------------------------------------------------------------------------
#include "MantidDataObjects/FakeMD.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidDataObjects/MDEventInserter.h"
#include "MantidKernel/ThreadPool.h"
#include "MantidKernel/ThreadScheduler.h"
#include "MantidKernel/Utils.h"
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
namespace Mantid {
namespace DataObjects {
using Kernel::ThreadPool;
using Kernel::ThreadSchedulerFIFO;
/**
* Constructor
* @param uniformParams Add a uniform, randomized distribution of events
* @param peakParams Add a peak with a normal distribution around a central
point
* @param randomSeed Seed int for the random number generator
* @param randomizeSignal If true, the events' signal and error values will be "
randomized around 1.0+-0.5
*/
FakeMD::FakeMD(const std::vector<double> &uniformParams,
const std::vector<double> &peakParams, const int randomSeed,
const bool randomizeSignal)
: m_uniformParams(uniformParams), m_peakParams(peakParams),
m_randomSeed(randomSeed), m_randomizeSignal(randomizeSignal), m_detIDs(),
m_randGen(1), m_uniformDist() {
if (uniformParams.empty() && peakParams.empty()) {
throw std::invalid_argument(
"You must specify at least one of peakParams or uniformParams");
}
}
/**
* Add the fake data to the given workspace
* @param workspace A pointer to MD event workspace to fill using the object
* parameters
*/
void FakeMD::fill(API::IMDEventWorkspace_sptr workspace) {
setupDetectorCache(*workspace);
CALL_MDEVENT_FUNCTION(this->addFakePeak, workspace)
CALL_MDEVENT_FUNCTION(this->addFakeUniformData, workspace)
// Mark that events were added, so the file back end (if any) needs updating
workspace->setFileNeedsUpdating(true);
}
/**
* Setup a detector cache for randomly picking IDs from the first
* instrument in the ExperimentInfo list.
* @param workspace The input workspace
*/
void FakeMD::setupDetectorCache(const API::IMDEventWorkspace &workspace) {
try {
auto inst = workspace.getExperimentInfo(0)->getInstrument();
m_detIDs = inst->getDetectorIDs(true); // true=skip monitors
size_t max = m_detIDs.size() - 1;
m_uniformDist = boost::uniform_int<size_t>(0, max); // Includes max
} catch (std::invalid_argument &) {
}
}
/** Function makes up a fake single-crystal peak and adds it to the workspace.
*
* @param ws A pointer to the workspace that receives the events
*/
template <typename MDE, size_t nd>
void FakeMD::addFakePeak(typename MDEventWorkspace<MDE, nd>::sptr ws) {
if (m_peakParams.empty())
return;
if (m_peakParams.size() != nd + 2)
throw std::invalid_argument("PeakParams needs to have ndims+2 arguments.");
if (m_peakParams[0] <= 0)
throw std::invalid_argument("PeakParams: number_of_events needs to be > 0");
size_t num = size_t(m_peakParams[0]);
// Width of the peak
double desiredRadius = m_peakParams.back();
boost::mt19937 rng;
boost::uniform_real<coord_t> u2(0, 1.0); // Random from 0 to 1.0
boost::variate_generator<boost::mt19937 &, boost::uniform_real<coord_t>>
genUnit(rng, u2);
rng.seed((unsigned int)(m_randomSeed));
// Inserter to help choose the correct event type
auto eventHelper =
MDEventInserter<typename MDEventWorkspace<MDE, nd>::sptr>(ws);
for (size_t i = 0; i < num; ++i) {
// Algorithm to generate points along a random n-sphere (sphere with not
// necessarily 3 dimensions)
// from http://en.wikipedia.org/wiki/N-sphere as of May 6, 2011.
// First, points in a hyper-cube of size 1.0, centered at 0.
coord_t centers[nd];
coord_t radiusSquared = 0;
for (size_t d = 0; d < nd; d++) {
centers[d] = genUnit() - 0.5f; // Distribute around +- the center
radiusSquared += centers[d] * centers[d];
}
// Make a unit vector pointing in this direction
coord_t radius = static_cast<coord_t>(sqrt(radiusSquared));
for (size_t d = 0; d < nd; d++)
centers[d] /= radius;
// Now place the point along this radius, scaled with ^1/n for uniformity.
coord_t radPos = genUnit();
radPos = static_cast<coord_t>(
pow(radPos, static_cast<coord_t>(1.0 / static_cast<coord_t>(nd))));
for (size_t d = 0; d < nd; d++) {
// Multiply by the scaling and the desired peak radius
centers[d] *= (radPos * static_cast<coord_t>(desiredRadius));
// Also offset by the center of the peak, as taken in Params
centers[d] += static_cast<coord_t>(m_peakParams[d + 1]);
}
// Default or randomized error/signal
float signal = 1.0;
float errorSquared = 1.0;
if (m_randomizeSignal) {
signal = float(0.5 + genUnit());
errorSquared = float(0.5 + genUnit());
}
// Create and add the event.
eventHelper.insertMDEvent(signal, errorSquared, 1, pickDetectorID(),
centers); // 1 = run number
}
ws->splitBox();
auto *ts = new ThreadSchedulerFIFO();
ThreadPool tp(ts);
ws->splitAllIfNeeded(ts);
tp.joinAll();
ws->refreshCache();
}
/**
* Function makes up a fake uniform event data and adds it to the workspace.
* @param ws
*/
template <typename MDE, size_t nd>
void FakeMD::addFakeUniformData(typename MDEventWorkspace<MDE, nd>::sptr ws) {
if (m_uniformParams.empty())
return;
bool randomEvents = true;
if (m_uniformParams[0] < 0) {
randomEvents = false;
m_uniformParams[0] = -m_uniformParams[0];
}
if (m_uniformParams.size() == 1) {
if (randomEvents) {
for (size_t d = 0; d < nd; ++d) {
m_uniformParams.push_back(ws->getDimension(d)->getMinimum());
m_uniformParams.push_back(ws->getDimension(d)->getMaximum());
}
} else // regular events
{
size_t nPoints = size_t(m_uniformParams[0]);
double Vol = 1;
for (size_t d = 0; d < nd; ++d)
Vol *= (ws->getDimension(d)->getMaximum() -
ws->getDimension(d)->getMinimum());
if (Vol == 0 || Vol > std::numeric_limits<float>::max())
throw std::invalid_argument(
" Domain ranges are not defined properly for workspace: " +
ws->getName());
double dV = Vol / double(nPoints);
double delta0 = std::pow(dV, 1. / double(nd));
for (size_t d = 0; d < nd; ++d) {
double min = ws->getDimension(d)->getMinimum();
m_uniformParams.push_back(min * (1 + FLT_EPSILON) - min + FLT_EPSILON);
double extent = ws->getDimension(d)->getMaximum() - min;
size_t nStrides = size_t(extent / delta0);
if (nStrides < 1)
nStrides = 1;
m_uniformParams.push_back(extent / static_cast<double>(nStrides));
}
}
}
if ((m_uniformParams.size() != 1 + nd * 2))
throw std::invalid_argument(
"UniformParams: needs to have ndims*2+1 arguments ");
if (randomEvents)
addFakeRandomData<MDE, nd>(m_uniformParams, ws);
else
addFakeRegularData<MDE, nd>(m_uniformParams, ws);
ws->splitBox();
auto *ts = new ThreadSchedulerFIFO();
ThreadPool tp(ts);
ws->splitAllIfNeeded(ts);
tp.joinAll();
ws->refreshCache();
}
/**
* Add fake randomized data to the workspace
* @param params A reference to the parameter vector
* @param ws The workspace to hold the data
*/
template <typename MDE, size_t nd>
void FakeMD::addFakeRandomData(const std::vector<double> &params,
typename MDEventWorkspace<MDE, nd>::sptr ws) {
size_t num = size_t(params[0]);
if (num == 0)
throw std::invalid_argument(
" number of distributed events can not be equal to 0");
boost::mt19937 rng;
rng.seed((unsigned int)(m_randomSeed));
// Unit-size randomizer
boost::uniform_real<double> u2(0, 1.0); // Random from 0 to 1.0
boost::variate_generator<boost::mt19937 &, boost::uniform_real<double>>
genUnit(rng, u2);
// Make a random generator for each dimensions
typedef boost::variate_generator<boost::mt19937 &,
boost::uniform_real<double>> gen_t;
// Inserter to help choose the correct event type
auto eventHelper =
MDEventInserter<typename MDEventWorkspace<MDE, nd>::sptr>(ws);
gen_t *gens[nd];
for (size_t d = 0; d < nd; ++d) {
double min = params[d * 2 + 1];
double max = params[d * 2 + 2];
if (max <= min)
throw std::invalid_argument(
"UniformParams: min must be < max for all dimensions.");
boost::uniform_real<double> u(min, max); // Range
gen_t *gen = new gen_t(rng, u);
gens[d] = gen;
}
// Create all the requested events
for (size_t i = 0; i < num; ++i) {
coord_t centers[nd];
for (size_t d = 0; d < nd; d++) {
centers[d] = static_cast<coord_t>(
(*gens[d])()); // use a different generator for each dimension
}
// Default or randomized error/signal
float signal = 1.0;
float errorSquared = 1.0;
if (m_randomizeSignal) {
signal = float(0.5 + genUnit());
errorSquared = float(0.5 + genUnit());
}
// Create and add the event.
eventHelper.insertMDEvent(signal, errorSquared, 1, pickDetectorID(),
centers); // 1 = run number
}
/// Clean up the generators
for (size_t d = 0; d < nd; ++d)
delete gens[d];
}
template <typename MDE, size_t nd>
void FakeMD::addFakeRegularData(const std::vector<double> &params,
typename MDEventWorkspace<MDE, nd>::sptr ws) {
// the parameters for regular distribution of events over the box
std::vector<double> startPoint(nd), delta(nd);
std::vector<size_t> indexMax(nd);
size_t gridSize(0);
size_t num = size_t(params[0]);
if (num == 0)
throw std::invalid_argument(
" number of distributed events can not be equal to 0");
// Inserter to help choose the correct event type
auto eventHelper =
MDEventInserter<typename MDEventWorkspace<MDE, nd>::sptr>(ws);
gridSize = 1;
for (size_t d = 0; d < nd; ++d) {
double min = ws->getDimension(d)->getMinimum();
double max = ws->getDimension(d)->getMaximum();
double shift = params[d * 2 + 1];
double step = params[d * 2 + 2];
if (shift < 0)
shift = 0;
if (shift >= step)
shift = step * (1 - FLT_EPSILON);
startPoint[d] = min + shift;
if ((startPoint[d] < min) || (startPoint[d] >= max))
throw std::invalid_argument("RegularData: starting point must be within "
"the box for all dimensions.");
if (step <= 0)
throw(std::invalid_argument(
"Step of the regular grid is less or equal to 0"));
indexMax[d] = size_t((max - min) / step);
if (indexMax[d] == 0)
indexMax[d] = 1;
// deal with round-off errors
while ((startPoint[d] + double(indexMax[d] - 1) * step) >= max)
step *= (1 - FLT_EPSILON);
delta[d] = step;
gridSize *= indexMax[d];
}
// Create all the requested events
std::vector<size_t> indexes;
size_t cellCount(0);
for (size_t i = 0; i < num; ++i) {
coord_t centers[nd];
Kernel::Utils::getIndicesFromLinearIndex(cellCount, indexMax, indexes);
++cellCount;
if (cellCount >= gridSize)
cellCount = 0;
for (size_t d = 0; d < nd; d++) {
centers[d] = coord_t(startPoint[d] + delta[d] * double(indexes[d]));
}
// Default or randomized error/signal
float signal = 1.0;
float errorSquared = 1.0;
// Create and add the event.
eventHelper.insertMDEvent(signal, errorSquared, 1, pickDetectorID(),
centers); // 1 = run number
}
}
/**
* Pick a detector ID for a particular event
* @returns A detector ID randomly selected from the instrument
*/
detid_t FakeMD::pickDetectorID() {
if (m_detIDs.empty()) {
return -1;
} else {
/// A variate generator to combine a random number generator with a
/// distribution
typedef boost::variate_generator<
boost::mt19937 &, boost::uniform_int<size_t>> uniform_generator;
uniform_generator uniformRand(m_randGen, m_uniformDist);
const size_t randIndex = uniformRand();
return m_detIDs[randIndex];
}
}
} // namespace DataObjects
} // namespace Mantid
#ifndef MANTID_DATAOBJECTS_FAKEMDTEST_H_
#define MANTID_DATAOBJECTS_FAKEMDTEST_H_
#include <vector>
#include "MantidDataObjects/FakeMD.h"
#include "MantidTestHelpers/ComponentCreationHelper.h"
#include "MantidTestHelpers/MDEventsTestHelper.h"
#include <boost/assign/list_of.hpp>
#include <cxxtest/TestSuite.h>
class FakeMDTest : 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 FakeMDTest *createSuite() { return new FakeMDTest(); }
static void destroySuite(FakeMDTest *suite) { delete suite; }
//---------------------------------------------------------------------------
// Failure cases
//---------------------------------------------------------------------------
void test_empty_peak_and_uniform_params_throws_exception() {
using Mantid::DataObjects::FakeMD;
const std::vector<double> peakParams;
const std::vector<double> uniformParams;
const int randomSeed(0);
const bool randomizeSignal(false);
TS_ASSERT_THROWS(
FakeMD(uniformParams, peakParams, randomSeed, randomizeSignal),
std::invalid_argument);
}
//---------------------------------------------------------------------------
// Success cases
//---------------------------------------------------------------------------
void test_no_randomize() {
using Mantid::DataObjects::FakeMD;
using Mantid::DataObjects::MDEventsTestHelper::makeMDEW;
using namespace boost::assign;
// Destination workspace
auto fakeData = makeMDEW<3>(10, 0.0, 10.0, 1);
// 1000 boxes with 1 event each
TS_ASSERT_EQUALS(fakeData->getNPoints(), 1000);
const std::vector<double> peakParams = list_of(1000.0)(5.0)(5.0)(5.0)(1.0);
const std::vector<double> uniformParams = list_of(10000.0);
const int randomSeed(0);
const bool randomizeSignal(false);
FakeMD faker(uniformParams, peakParams, randomSeed, randomizeSignal);
faker.fill(fakeData);
// Now there are 11000 more points.
TS_ASSERT_EQUALS(fakeData->getNPoints(), 12000);
}
void test_exec_randomizeSignal() {
using Mantid::DataObjects::FakeMD;
using Mantid::DataObjects::MDEventsTestHelper::makeMDEW;
using namespace boost::assign;
auto fakeData = makeMDEW<3>(10, 0.0, 10.0, 0);
TS_ASSERT_EQUALS(fakeData->getNPoints(), 0);
TS_ASSERT_DELTA(fakeData->getBox()->getSignal(), 0.0, 1e-5);
const std::vector<double> peakParams = list_of(100.0)(5.0)(5.0)(5.0)(1.0);
const std::vector<double> uniformParams = list_of(100.0);
const int randomSeed(0);
const bool randomizeSignal(true);
FakeMD faker(uniformParams, peakParams, randomSeed, randomizeSignal);
faker.fill(fakeData);
// Now there are 200 more points.
TS_ASSERT_EQUALS(fakeData->getNPoints(), 200);
// 200 +- 100 signal
TS_ASSERT_DELTA(fakeData->getBox()->getSignal(), 200.0, 100);
TS_ASSERT_DELTA(fakeData->getBox()->getErrorSquared(), 200.0, 100);
// But not exactly 200
TS_ASSERT_DIFFERS(fakeData->getBox()->getSignal(), 200.0);
TS_ASSERT_DIFFERS(fakeData->getBox()->getErrorSquared(), 200.0);
TSM_ASSERT("If the workspace is file-backed, then it needs updating.",
fakeData->fileNeedsUpdating());
}
void testExecRegularSignal() {
using Mantid::DataObjects::FakeMD;
using Mantid::DataObjects::MDEventsTestHelper::makeMDEW;
using namespace boost::assign;
auto fakeData = makeMDEW<3>(10, 0.0, 10.0, 0);
// No events
TS_ASSERT_EQUALS(fakeData->getNPoints(), 0);
TS_ASSERT_DELTA(fakeData->getBox()->getSignal(), 0.0, 1e-5);
const std::vector<double> peakParams;
const std::vector<double> uniformParams = list_of(-1000.0);
const int randomSeed(0);
const bool randomizeSignal(false);