Mantid
Loading...
Searching...
No Matches
CalculateIqt.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
11#include "MantidAPI/Progress.h"
12#include "MantidHistogramData/HistogramY.h"
14
15#include <boost/numeric/conversion/cast.hpp>
16#include <cmath>
17#include <functional>
18#include <utility>
19
20using namespace Mantid::Algorithms;
21using namespace Mantid::API;
22using namespace Mantid::Kernel;
23using namespace Mantid::HistogramData;
24
25namespace {
26constexpr int DEFAULT_ITERATIONS = 100;
27constexpr int DEFAULT_SEED = 89631139;
28
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();
34}
35
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]);
40}
41
42MatrixWorkspace_sptr randomizeWorkspaceWithinError(MatrixWorkspace_sptr workspace, MersenneTwister &mTwister) {
43 auto randomNumberGenerator = [&mTwister](const double error) { return mTwister.nextValue(-error, error); };
44 for (auto i = 0u; i < workspace->getNumberHistograms(); ++i)
45 randomizeHistogramWithinError(workspace->mutableY(i), workspace->e(i), randomNumberGenerator);
46 return workspace;
47}
48
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));
56}
57
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;
63}
64
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());
73 for (auto &&workspace : workspaces) {
74 const auto values = workspace->y(index).rawData();
75 for (auto j = 0u; j < values.size(); ++j)
76 yValues[j].emplace_back(values[j]);
77 }
78 return yValues;
79}
80
81int getWorkspaceNumberOfHistograms(const MatrixWorkspace_sptr &workspace) {
82 return boost::numeric_cast<int>(workspace->getNumberHistograms());
83}
84
85} // namespace
86
87namespace Mantid::Algorithms {
88
90
91const std::string CalculateIqt::name() const { return "CalculateIqt"; }
92
93int CalculateIqt::version() const { return 1; }
94
95const std::vector<std::string> CalculateIqt::seeAlso() const { return {"TransformToIqt"}; }
96
97const std::string CalculateIqt::category() const { return "Inelastic\\Indirect"; }
98
99const std::string CalculateIqt::summary() const {
100 return "Calculates I(Q,t) from S(Q,w) and computes the errors using a "
101 "monte-carlo routine.";
102}
103
105
106 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
107 "The name of the sample workspace.");
108 declareProperty(std::make_unique<WorkspaceProperty<>>("ResolutionWorkspace", "", Direction::Input),
109 "The name of the resolution workspace.");
110
111 declareProperty("EnergyMin", -0.5, "Minimum energy for fit. Default = -0.5.");
112 declareProperty("EnergyMax", 0.5, "Maximum energy for fit. Default = 0.5.");
113 declareProperty("EnergyWidth", 0.1, "Width of energy bins for fit.");
114
115 auto positiveInt = std::make_shared<Kernel::BoundedValidator<int>>();
116 positiveInt->setLower(1);
117
118 declareProperty("NumberOfIterations", DEFAULT_ITERATIONS, positiveInt,
119 "Number of randomised simulations within error to run.");
120 declareProperty("SeedValue", DEFAULT_SEED, positiveInt,
121 "Seed the random number generator for monte-carlo error calculation.");
122
123 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
124 "The name to use for the output workspace.");
125
126 declareProperty("CalculateErrors", true, "Calculate monte-carlo errors.");
127}
128
130 const auto rebinParams = rebinParamsAsString();
131 const MatrixWorkspace_sptr sampleWorkspace = getProperty("InputWorkspace");
132 MatrixWorkspace_sptr resolution = getProperty("ResolutionWorkspace");
133 bool calculateErrors = getProperty("CalculateErrors");
134 const int nIterations = getProperty("NumberOfIterations");
135 const int seed = getProperty("SeedValue");
136 resolution = normalizedFourierTransform(resolution, rebinParams);
137
138 auto outputWorkspace =
139 monteCarloErrorCalculation(sampleWorkspace, resolution, rebinParams, seed, calculateErrors, nIterations);
140
141 outputWorkspace = replaceSpecialValues(outputWorkspace);
142 setProperty("OutputWorkspace", outputWorkspace);
143}
144
146 const double e_min = getProperty("EnergyMin");
147 const double e_max = getProperty("EnergyMax");
148 const double e_width = getProperty("EnergyWidth");
149 return createRebinString(e_min, e_max, e_width);
150}
151
153 const MatrixWorkspace_sptr &resolution,
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);
160
161 MersenneTwister mTwister(seed);
162 if (calculateErrors) {
163 Progress errorCalculationProg(this, 0.0, 1.0, nIterations);
164 PARALLEL_FOR_IF(Kernel::threadSafe(*sample, *resolution))
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);
169 PARALLEL_CRITICAL(emplace_back)
170 simulatedWorkspaces.emplace_back(simulated);
172 }
174 return setErrorsToStandardDeviation(simulatedWorkspaces);
175 }
176 return setErrorsToZero(simulatedWorkspaces);
177}
178
179std::map<std::string, std::string> CalculateIqt::validateInputs() {
180 std::map<std::string, std::string> issues;
181 const double eMin = getProperty("EnergyMin");
182 const double eMax = getProperty("EnergyMax");
183 if (eMin > eMax) {
184 auto energy_swapped = "EnergyMin is greater than EnergyMax";
185 issues["EnergyMin"] = energy_swapped;
186 issues["EnergyMax"] = energy_swapped;
187 }
188 return issues;
189}
190
192 auto rebinAlgorithm = createChildAlgorithm("Rebin");
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");
199}
200
202 auto integrationAlgorithm = createChildAlgorithm("Integration");
203 integrationAlgorithm->initialize();
204 integrationAlgorithm->setProperty("InputWorkspace", workspace);
205 integrationAlgorithm->setProperty("OutputWorkspace", "_");
206 integrationAlgorithm->execute();
207 return integrationAlgorithm->getProperty("OutputWorkspace");
208}
209
211 auto pointDataAlgorithm = createChildAlgorithm("ConvertToPointData");
212 pointDataAlgorithm->initialize();
213 pointDataAlgorithm->setProperty("InputWorkspace", workspace);
214 pointDataAlgorithm->setProperty("OutputWorkspace", "_");
215 pointDataAlgorithm->execute();
216 return pointDataAlgorithm->getProperty("OutputWorkspace");
217}
218
220 auto FFTAlgorithm = createChildAlgorithm("ExtractFFTSpectrum");
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");
227}
228
230 const MatrixWorkspace_sptr &rhsWorkspace) {
231 auto divideAlgorithm = createChildAlgorithm("Divide");
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");
238}
239
241 auto cropAlgorithm = createChildAlgorithm("CropWorkspace");
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");
248}
249
251 auto specialValuesAlgorithm = createChildAlgorithm("ReplaceSpecialValues");
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");
260}
261
263 const std::string &rebinParams) {
264 workspace = rebin(workspace, rebinParams);
265 auto workspace_int = integration(workspace);
268 return divide(workspace, workspace_int);
269}
270
272 const MatrixWorkspace_sptr &resolutionWorkspace,
273 const std::string &rebinParams) {
275 return divide(workspace, resolutionWorkspace);
276}
277
279 const std::string &rebinParams, MersenneTwister &mTwister) {
280 auto simulatedWorkspace = randomizeWorkspaceWithinError(std::move(sample), mTwister);
281 return calculateIqt(simulatedWorkspace, resolution, rebinParams);
282}
283
285CalculateIqt::setErrorsToStandardDeviation(const std::vector<MatrixWorkspace_sptr> &simulatedWorkspaces) {
286 auto outputWorkspace = simulatedWorkspaces.front();
288 for (auto i = 0; i < getWorkspaceNumberOfHistograms(outputWorkspace); ++i) {
290 outputWorkspace->mutableE(i) = standardDeviationArray(allYValuesAtIndex(simulatedWorkspaces, i));
292 }
294 return outputWorkspace;
295}
296
297MatrixWorkspace_sptr CalculateIqt::setErrorsToZero(const std::vector<MatrixWorkspace_sptr> &simulatedWorkspaces) {
298 auto outputWorkspace = simulatedWorkspaces.front();
300 for (auto i = 0; i < getWorkspaceNumberOfHistograms(outputWorkspace); ++i) {
302 outputWorkspace->mutableE(i) = 0;
304 }
306 return outputWorkspace;
307}
308
309} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double error
Definition: IndexPeaks.cpp:133
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#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...
Definition: MultiThreaded.h:94
#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.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
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.
Definition: Algorithm.cpp:842
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
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 &params)
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.
Definition: ProgressBase.h:51
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.
Definition: MultiThreaded.h:22
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54