12#include "MantidHistogramData/HistogramY.h"
15#include <boost/numeric/conversion/cast.hpp>
26constexpr int DEFAULT_ITERATIONS = 100;
27constexpr int DEFAULT_SEED = 89631139;
29std::string createRebinString(
double minimum,
double maximum,
double width) {
30 std::stringstream rebinStream;
31 rebinStream.precision(14);
32 rebinStream << minimum <<
", " << width <<
", " << maximum;
33 return rebinStream.str();
36template <
typename Generator>
37void randomizeHistogramWithinError(HistogramY &row,
const HistogramE &errors, Generator &generator) {
38 for (
auto i = 0u; i < row.size(); ++i)
39 row[i] += generator(errors[i]);
44 for (
auto i = 0u; i <
workspace->getNumberHistograms(); ++i)
45 randomizeHistogramWithinError(
workspace->mutableY(i),
workspace->e(i), randomNumberGenerator);
49double standardDeviation(
const std::vector<double> &inputValues) {
50 const auto inputSize = boost::numeric_cast<double>(inputValues.size());
51 const auto mean = std::accumulate(inputValues.begin(), inputValues.end(), 0.0) / inputSize;
52 const double sumOfXMinusMeanSquared =
53 std::accumulate(inputValues.cbegin(), inputValues.cend(), 0.,
54 [mean](
const auto sum,
const auto x) { return sum + std::pow(x - mean, 2); });
55 return sqrt(sumOfXMinusMeanSquared / (inputSize - 1));
58std::vector<double> standardDeviationArray(
const std::vector<std::vector<double>> &yValues) {
59 std::vector<double> standardDeviations;
60 standardDeviations.reserve(yValues.size());
61 std::transform(yValues.begin(), yValues.end(), std::back_inserter(standardDeviations), standardDeviation);
62 return standardDeviations;
70std::vector<std::vector<double>> allYValuesAtIndex(
const std::vector<MatrixWorkspace_sptr> &workspaces,
71 const std::size_t
index) {
72 std::vector<std::vector<double>> yValues(workspaces[0]->getDimension(0)->getNBins());
75 for (
auto j = 0u; j < values.size(); ++j)
76 yValues[j].emplace_back(values[j]);
82 return boost::numeric_cast<int>(
workspace->getNumberHistograms());
100 return "Calculates I(Q,t) from S(Q,w) and computes the errors using a "
101 "monte-carlo routine.";
107 "The name of the sample workspace.");
109 "The name of the resolution workspace.");
111 declareProperty(
"EnergyMin", -0.5,
"Minimum energy for fit. Default = -0.5.");
112 declareProperty(
"EnergyMax", 0.5,
"Maximum energy for fit. Default = 0.5.");
115 auto positiveInt = std::make_shared<Kernel::BoundedValidator<int>>();
116 positiveInt->setLower(1);
119 "Number of randomised simulations within error to run.");
121 "Seed the random number generator for monte-carlo error calculation.");
124 "The name to use for the output workspace.");
126 declareProperty(
"CalculateErrors",
true,
"Calculate monte-carlo errors.");
133 bool calculateErrors =
getProperty(
"CalculateErrors");
134 const int nIterations =
getProperty(
"NumberOfIterations");
138 auto outputWorkspace =
149 return createRebinString(e_min, e_max, e_width);
154 const std::string &rebinParams,
const int seed,
155 const bool calculateErrors,
const int nIterations) {
156 auto outputWorkspace =
calculateIqt(sample, resolution, rebinParams);
157 std::vector<MatrixWorkspace_sptr> simulatedWorkspaces;
158 simulatedWorkspaces.reserve(nIterations);
159 simulatedWorkspaces.emplace_back(outputWorkspace);
162 if (calculateErrors) {
163 Progress errorCalculationProg(
this, 0.0, 1.0, nIterations);
165 for (
auto i = 0; i < nIterations - 1; ++i) {
166 errorCalculationProg.
report(
"Calculating Monte Carlo errors...");
168 auto simulated =
doSimulation(sample->clone(), resolution, rebinParams, mTwister);
170 simulatedWorkspaces.emplace_back(simulated);
180 std::map<std::string, std::string> issues;
184 auto energy_swapped =
"EnergyMin is greater than EnergyMax";
185 issues[
"EnergyMin"] = energy_swapped;
186 issues[
"EnergyMax"] = energy_swapped;
193 rebinAlgorithm->initialize();
194 rebinAlgorithm->setProperty(
"InputWorkspace",
workspace);
195 rebinAlgorithm->setProperty(
"OutputWorkspace",
"_");
196 rebinAlgorithm->setProperty(
"Params", params);
197 rebinAlgorithm->execute();
198 return rebinAlgorithm->getProperty(
"OutputWorkspace");
203 integrationAlgorithm->initialize();
204 integrationAlgorithm->setProperty(
"InputWorkspace",
workspace);
205 integrationAlgorithm->setProperty(
"OutputWorkspace",
"_");
206 integrationAlgorithm->execute();
207 return integrationAlgorithm->getProperty(
"OutputWorkspace");
212 pointDataAlgorithm->initialize();
213 pointDataAlgorithm->setProperty(
"InputWorkspace",
workspace);
214 pointDataAlgorithm->setProperty(
"OutputWorkspace",
"_");
215 pointDataAlgorithm->execute();
216 return pointDataAlgorithm->getProperty(
"OutputWorkspace");
221 FFTAlgorithm->initialize();
222 FFTAlgorithm->setProperty(
"InputWorkspace",
workspace);
223 FFTAlgorithm->setProperty(
"OutputWorkspace",
"_");
224 FFTAlgorithm->setProperty(
"FFTPart", 2);
225 FFTAlgorithm->execute();
226 return FFTAlgorithm->getProperty(
"OutputWorkspace");
232 divideAlgorithm->initialize();
233 divideAlgorithm->setProperty(
"LHSWorkspace", lhsWorkspace);
234 divideAlgorithm->setProperty(
"RHSWorkspace", rhsWorkspace);
235 divideAlgorithm->setProperty(
"OutputWorkspace",
"_");
236 divideAlgorithm->execute();
237 return divideAlgorithm->getProperty(
"OutputWorkspace");
242 cropAlgorithm->initialize();
243 cropAlgorithm->setProperty(
"InputWorkspace",
workspace);
244 cropAlgorithm->setProperty(
"OutputWorkspace",
"_");
245 cropAlgorithm->setProperty(
"XMax", xMax);
246 cropAlgorithm->execute();
247 return cropAlgorithm->getProperty(
"OutputWorkspace");
252 specialValuesAlgorithm->initialize();
253 specialValuesAlgorithm->setProperty(
"InputWorkspace",
workspace);
254 specialValuesAlgorithm->setProperty(
"OutputWorkspace",
"_");
255 specialValuesAlgorithm->setProperty(
"InfinityValue", 0.0);
256 specialValuesAlgorithm->setProperty(
"BigNumberThreshold", 1.0001);
257 specialValuesAlgorithm->setProperty(
"NaNValue", 0.0);
258 specialValuesAlgorithm->execute();
259 return specialValuesAlgorithm->getProperty(
"OutputWorkspace");
263 const std::string &rebinParams) {
273 const std::string &rebinParams) {
280 auto simulatedWorkspace = randomizeWorkspaceWithinError(std::move(sample), mTwister);
281 return calculateIqt(simulatedWorkspace, resolution, rebinParams);
286 auto outputWorkspace = simulatedWorkspaces.front();
288 for (
auto i = 0; i < getWorkspaceNumberOfHistograms(outputWorkspace); ++i) {
290 outputWorkspace->mutableE(i) = standardDeviationArray(allYValuesAtIndex(simulatedWorkspaces, i));
294 return outputWorkspace;
298 auto outputWorkspace = simulatedWorkspaces.front();
300 for (
auto i = 0; i < getWorkspaceNumberOfHistograms(outputWorkspace); ++i) {
302 outputWorkspace->mutableE(i) = 0;
306 return outputWorkspace;
#define DECLARE_ALGORITHM(classname)
IPeaksWorkspace_sptr workspace
std::map< DeltaEMode::Type, std::string > index
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
#define PARALLEL_CRITICAL(name)
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Helper class for reporting progress from algorithms.
A property class for workspaces.
int version() const override
function to return a version of the algorithm, must be overridden in all algorithms
const std::string summary() const override
function returns a summary message that will be displayed in the default GUI, and in the help.
API::MatrixWorkspace_sptr calculateIqt(API::MatrixWorkspace_sptr workspace, const API::MatrixWorkspace_sptr &resolutionWorkspace, const std::string &rebinParams)
API::MatrixWorkspace_sptr cropWorkspace(const API::MatrixWorkspace_sptr &workspace, const double xMax)
API::MatrixWorkspace_sptr extractFFTSpectrum(const API::MatrixWorkspace_sptr &workspace)
API::MatrixWorkspace_sptr replaceSpecialValues(const API::MatrixWorkspace_sptr &workspace)
const std::vector< std::string > seeAlso() const override
Function to return all of the seeAlso algorithms related to this algorithm.
const std::string category() const override
function to return a category of the algorithm.
API::MatrixWorkspace_sptr monteCarloErrorCalculation(const API::MatrixWorkspace_sptr &sample, const API::MatrixWorkspace_sptr &resolution, const std::string &rebinParams, const int seed, const bool calculateErrors, const int nIterations)
API::MatrixWorkspace_sptr rebin(const API::MatrixWorkspace_sptr &workspace, const std::string ¶ms)
std::string rebinParamsAsString()
API::MatrixWorkspace_sptr convertToPointData(const API::MatrixWorkspace_sptr &workspace)
void init() override
Virtual method - must be overridden by concrete algorithm.
API::MatrixWorkspace_sptr divide(const API::MatrixWorkspace_sptr &lhsWorkspace, const API::MatrixWorkspace_sptr &rhsWorkspace)
API::MatrixWorkspace_sptr setErrorsToStandardDeviation(const std::vector< API::MatrixWorkspace_sptr > &simulatedWorkspaces)
API::MatrixWorkspace_sptr normalizedFourierTransform(API::MatrixWorkspace_sptr workspace, const std::string &rebinParams)
API::MatrixWorkspace_sptr integration(const API::MatrixWorkspace_sptr &workspace)
void exec() override
Virtual method - must be overridden by concrete algorithm.
API::MatrixWorkspace_sptr doSimulation(API::MatrixWorkspace_sptr sample, const API::MatrixWorkspace_sptr &resolution, const std::string &rebinParams, Kernel::MersenneTwister &mTwister)
API::MatrixWorkspace_sptr setErrorsToZero(const std::vector< API::MatrixWorkspace_sptr > &simulatedWorkspaces)
std::map< std::string, std::string > validateInputs() override
Method checking errors on ALL the inputs, before execution.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
This implements the the Mersenne Twister 19937 pseudo-random number generator algorithm as a specialz...
double nextValue() override
Generate the next random number in the sequence within the given range default range.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
@ Input
An input workspace.
@ Output
An output workspace.