Mantid
Loading...
Searching...
No Matches
SetUncertainties.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 +
17
18#include <algorithm>
19#include <vector>
20
21namespace Mantid::Algorithms {
22
23// Register the algorithm into the AlgorithmFactory
24DECLARE_ALGORITHM(SetUncertainties)
25
26using namespace Kernel;
27using namespace API;
28
29namespace {
31const double TOLERANCE = 1.e-10;
32
33const std::string ZERO("zero");
34const std::string SQRT("sqrt");
35const std::string ONE_IF_ZERO("oneIfZero");
36const std::string SQRT_OR_ONE("sqrtOrOne");
37const std::string CUSTOM("custom");
38
39struct SetError {
40 explicit SetError(const double setTo, const double ifEqualTo, const double tolerance)
41 : valueToSet(setTo), valueToCompare(ifEqualTo), tolerance(tolerance) {}
42
43 double operator()(const double error) {
44 double deviation = error - valueToCompare;
45 if (deviation < tolerance && deviation >= 0) {
46 return valueToSet;
47 } else {
48 return error;
49 }
50 }
51
52 double valueToSet;
54 double tolerance;
55};
56
57struct SqrtError {
58 explicit SqrtError(const double constant) : zeroSqrtValue(constant) {}
59
60 double operator()(const double intensity) {
61 const double localIntensity = fabs(intensity);
62 if (localIntensity > TOLERANCE) {
63 return sqrt(localIntensity);
64 } else {
65 return zeroSqrtValue;
66 }
67 }
68
70};
71} // namespace
72
74const std::string SetUncertainties::name() const { return "SetUncertainties"; }
75
77int SetUncertainties::version() const { return (1); }
78
80 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
81 auto mustBePositiveInt = std::make_shared<BoundedValidator<int>>();
82 mustBePositive->setLower(0);
83 mustBePositiveInt->setLower(0);
84 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Direction::Input));
85 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Direction::Output));
86 std::vector<std::string> errorTypes = {ZERO, SQRT, SQRT_OR_ONE, ONE_IF_ZERO, CUSTOM};
87 declareProperty("SetError", ZERO, std::make_shared<StringListValidator>(errorTypes),
88 "How to reset the uncertainties");
89 declareProperty("SetErrorTo", 1.000, mustBePositive, "The error value to set when using custom mode");
90 setPropertySettings("SetErrorTo", std::make_unique<VisibleWhenProperty>("SetError", IS_EQUAL_TO, "custom"));
91
92 declareProperty("IfEqualTo", 0.000, mustBePositive,
93 "Which error values in the input workspace should be "
94 "replaced when using custom mode");
95 setPropertySettings("IfEqualTo", std::make_unique<VisibleWhenProperty>("SetError", IS_EQUAL_TO, "custom"));
96
97 declareProperty("Precision", 3, mustBePositiveInt,
98 "How many decimal places of ``IfEqualTo`` are taken into "
99 "account for matching when using custom mode");
100 setPropertySettings("Precision", std::make_unique<VisibleWhenProperty>("SetError", IS_EQUAL_TO, "custom"));
101}
102
104 MatrixWorkspace_const_sptr inputWorkspace = getProperty("InputWorkspace");
105 auto inputEventWorkspace = std::dynamic_pointer_cast<const DataObjects::EventWorkspace>(inputWorkspace);
106 MatrixWorkspace_sptr outputWorkspace = getProperty("OutputWorkspace");
107 std::string errorType = getProperty("SetError");
108 bool zeroError = (errorType == ZERO);
109 bool takeSqrt = ((errorType == SQRT) || (errorType == SQRT_OR_ONE));
110 bool resetOne = ((errorType == ONE_IF_ZERO) || (errorType == SQRT_OR_ONE));
111 bool customError = (errorType == CUSTOM);
112
113 double valueToSet = resetOne ? 1.0 : getProperty("SetErrorTo");
114 double valueToCompare = resetOne ? 0.0 : getProperty("IfEqualTo");
115 int precision = getProperty("Precision");
116 double tolerance = resetOne ? 1E-10 : std::pow(10.0, -1. * precision);
117
118 // Create the output workspace. This will copy many aspects from the input
119 // one.
120
121 const int64_t numHists = static_cast<int64_t>(inputWorkspace->getNumberHistograms());
122 if (inputEventWorkspace) {
123 if (inputWorkspace == outputWorkspace && errorType == SQRT &&
124 inputEventWorkspace->getEventType() == Mantid::API::EventType::TOF) {
125 // Uncertainty is already square root of the y value
126 return;
127 }
128 // Copy the histogram representation over to a Workspace2D
129 outputWorkspace = WorkspaceFactory::Instance().create(inputWorkspace);
130 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace))
131 for (int64_t i = 0; i < numHists; ++i) {
133 // copy X/Y/E
134 outputWorkspace->setSharedX(i, inputWorkspace->sharedX(i));
135 outputWorkspace->setSharedY(i, inputWorkspace->sharedY(i));
136 outputWorkspace->setSharedE(i, inputWorkspace->sharedE(i));
138 }
140 } else if (inputWorkspace != outputWorkspace) {
141 outputWorkspace = inputWorkspace->clone();
142 }
143
144 const auto &spectrumInfo = inputWorkspace->spectrumInfo();
145 Progress prog(this, 0.0, 1.0, numHists);
146 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace))
147 for (int64_t i = 0; i < numHists; ++i) {
149 // copy the E or set to zero depending on the mode
150 if (errorType == ONE_IF_ZERO || customError) {
151 } else {
152 outputWorkspace->mutableE(i) = 0.0;
153 }
154 // ZERO mode doesn't calculate anything further
155 if ((!zeroError) && (!(spectrumInfo.hasDetectors(i) && spectrumInfo.isMasked(i)))) {
156 auto &E = outputWorkspace->mutableE(i);
157 if (takeSqrt) {
158 const auto &Y = outputWorkspace->y(i);
159 std::transform(Y.begin(), Y.end(), E.begin(), SqrtError(resetOne ? 1. : 0.));
160 } else {
161 std::transform(E.begin(), E.end(), E.begin(), SetError(valueToSet, valueToCompare, tolerance));
162 }
163 }
164 prog.report();
166 }
168 setProperty("OutputWorkspace", outputWorkspace);
169}
170
171} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const int64_t TOLERANCE(1000000)
double error
Definition: IndexPeaks.cpp:133
#define fabs(x)
Definition: Matrix.cpp:22
#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_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.
double zeroSqrtValue
double valueToSet
double tolerance
double valueToCompare
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void exec() override
Execution code.
int version() const override
Algorithm's version.
const std::string name() const override
Algorithm's name.
void init() override
Initialisation code.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
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
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54