Mantid
Loading...
Searching...
No Matches
CreateWorkspace.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 +
8
13#include "MantidAPI/TextAxis.h"
16#include "MantidHistogramData/HistogramBuilder.h"
17#include "MantidIndexing/GlobalSpectrumIndex.h"
18#include "MantidIndexing/IndexInfo.h"
19#include "MantidIndexing/SpectrumIndexSet.h"
26
27namespace Mantid::Algorithms {
28
29using namespace Kernel;
30using namespace API;
31using namespace DataObjects;
32using namespace HistogramData;
33using namespace Indexing;
34
35DECLARE_ALGORITHM(CreateWorkspace)
36
37
38void CreateWorkspace::init() {
39
40 std::vector<std::string> unitOptions = UnitFactory::Instance().getKeys();
41 unitOptions.emplace_back("SpectraNumber");
42 unitOptions.emplace_back("Text");
43
44 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
45 "Name to be given to the created workspace.");
46
47 auto required = std::make_shared<MandatoryValidator<std::vector<double>>>();
48 declareProperty(std::make_unique<ArrayProperty<double>>("DataX", required), "X-axis data values for workspace.");
49 declareProperty(std::make_unique<ArrayProperty<double>>("DataY", required),
50 "Y-axis data values for workspace (measures).");
51 declareProperty(std::make_unique<ArrayProperty<double>>("DataE"), "Error values for workspace.");
52 declareProperty(std::make_unique<PropertyWithValue<int>>("NSpec", 1), "Number of spectra to divide data into.");
53 declareProperty("UnitX", "", "The unit to assign to the XAxis");
54
55 declareProperty("VerticalAxisUnit", "SpectraNumber", std::make_shared<StringListValidator>(unitOptions),
56 "The unit to assign to the second Axis (leave blank for "
57 "default Spectra number)");
58 declareProperty(std::make_unique<ArrayProperty<std::string>>("VerticalAxisValues"), "Values for the VerticalAxis.");
59
60 declareProperty(std::make_unique<PropertyWithValue<bool>>("Distribution", false),
61 "Whether OutputWorkspace should be marked as a distribution.");
62 declareProperty("YUnitLabel", "", "Label for Y Axis");
63
64 declareProperty("WorkspaceTitle", "", "Title for Workspace");
65
66 declareProperty(
67 std::make_unique<WorkspaceProperty<>>("ParentWorkspace", "", Direction::Input, PropertyMode::Optional),
68 "Name of a parent workspace.");
69 declareProperty(std::make_unique<ArrayProperty<double>>("Dx"), "X error values for workspace (optional).");
70 std::vector<std::string> propOptions{Parallel::toString(Parallel::StorageMode::Cloned),
71 Parallel::toString(Parallel::StorageMode::Distributed),
72 Parallel::toString(Parallel::StorageMode::MasterOnly)};
73 declareProperty("ParallelStorageMode", Parallel::toString(Parallel::StorageMode::Cloned),
74 std::make_shared<StringListValidator>(propOptions),
75 "The parallel storage mode of the output workspace for MPI builds");
76 setPropertySettings("ParallelStorageMode", std::make_unique<InvisibleProperty>());
77}
78
80std::map<std::string, std::string> CreateWorkspace::validateInputs() {
81 std::map<std::string, std::string> issues;
82
83 const std::string vUnit = getProperty("VerticalAxisUnit");
84 const std::vector<std::string> vAxis = getProperty("VerticalAxisValues");
85
86 if (vUnit == "SpectraNumber" && !vAxis.empty())
87 issues["VerticalAxisValues"] = "Axis values cannot be provided when using a spectra axis";
88
89 return issues;
90}
91
94 // Contortions to get at the vector in the property without copying it
95 const Property *const dataXprop = getProperty("DataX");
96 const Property *const dataYprop = getProperty("DataY");
97 const Property *const dataEprop = getProperty("DataE");
98 const Property *const errorDxprop = getProperty("Dx");
99
100 const auto xCheck = dynamic_cast<const ArrayProperty<double> *>(dataXprop);
101 if (!xCheck)
102 throw std::invalid_argument("DataX cannot be casted to a double vector");
103 const std::vector<double> &dataX = *xCheck;
104
105 const auto yCheck = dynamic_cast<const ArrayProperty<double> *>(dataYprop);
106 if (!yCheck)
107 throw std::invalid_argument("DataY cannot be casted to a double vector");
108 const std::vector<double> &dataY = *yCheck;
109
110 const auto eCheck = dynamic_cast<const ArrayProperty<double> *>(dataEprop);
111 if (!eCheck)
112 throw std::invalid_argument("DataE cannot be casted to a double vector");
113 const std::vector<double> &dataE = *eCheck;
114
115 const auto dxCheck = dynamic_cast<const ArrayProperty<double> *>(errorDxprop);
116 if (!dxCheck)
117 throw std::invalid_argument("Dx cannot be casted to a double vector");
118 const std::vector<double> &dX = *dxCheck;
119
120 const int nSpec = getProperty("NSpec");
121 const std::string xUnit = getProperty("UnitX");
122 const std::string vUnit = getProperty("VerticalAxisUnit");
123 const std::vector<std::string> vAxis = getProperty("VerticalAxisValues");
124
125 // Verify the size of the vertical axis.
126 const auto vAxisSize = static_cast<int>(vAxis.size());
127 if (vUnit != "SpectraNumber") {
128 // In the case of numerical axis, the vertical axis can represent either
129 // point data or bin edges.
130 if ((vUnit == "Text" && vAxisSize != nSpec) || (vAxisSize != nSpec && vAxisSize != nSpec + 1)) {
131 throw std::invalid_argument("The number of vertical axis values doesn't "
132 "match the number of histograms.");
133 }
134 }
135
136 // Verify length of vectors makes sense with NSpec
137 if ((dataY.size() % nSpec) != 0) {
138 throw std::invalid_argument("Length of DataY must be divisible by NSpec");
139 }
140 const std::size_t ySize = dataY.size() / nSpec;
141
142 // Check whether the X values provided are to be re-used for (are common to)
143 // every spectrum
144 const bool commonX(dataX.size() == ySize || dataX.size() == ySize + 1);
145
146 std::size_t xSize{dataX.size()};
147 HistogramBuilder histogramBuilder;
148 if (commonX) {
149 histogramBuilder.setX(dataX);
150 } else {
151 if (xSize % nSpec != 0) {
152 throw std::invalid_argument("Length of DataX must be divisible by NSpec");
153 }
154 xSize /= nSpec;
155 histogramBuilder.setX(xSize);
156 }
157 histogramBuilder.setY(ySize);
158
159 if (!dX.empty()) {
160 if (dX.size() != dataY.size())
161 throw std::runtime_error("Dx must have the same size as DataY");
162 histogramBuilder.setDx(ySize);
163 }
164
165 histogramBuilder.setDistribution(getProperty("Distribution"));
166 auto histogram = histogramBuilder.build();
167
168 const bool dataE_provided = !dataE.empty();
169 if (dataE_provided && dataY.size() != dataE.size()) {
170 throw std::runtime_error("DataE (if provided) must be the same size as DataY");
171 }
172
173 // Create the OutputWorkspace
174 MatrixWorkspace_const_sptr parentWS = getProperty("ParentWorkspace");
175 MatrixWorkspace_sptr outputWS;
176 if (parentWS) {
177 outputWS = create<HistoWorkspace>(*parentWS, nSpec, histogram);
178 } else {
179 auto storageMode = Parallel::fromString(getProperty("ParallelStorageMode"));
180 IndexInfo indexInfo(nSpec, storageMode, communicator());
181 outputWS = create<Workspace2D>(indexInfo, histogram);
182 }
183
184 Progress progress(this, 0.0, 1.0, nSpec);
185 const auto &indexInfo = outputWS->indexInfo();
186
188 for (int i = 0; i < nSpec; i++) {
190
191 // In an MPI run the global index i is not necessarily on this rank, i.e.,
192 // there might not be a corrsponding workspace index.
193 const auto localIndices = indexInfo.makeIndexSet({static_cast<GlobalSpectrumIndex>(i)});
194 if (localIndices.empty())
195 continue;
196
197 const std::vector<double>::difference_type xStart = i * xSize;
198 const std::vector<double>::difference_type xEnd = xStart + xSize;
199 const std::vector<double>::difference_type yStart = i * ySize;
200 const std::vector<double>::difference_type yEnd = yStart + ySize;
201 auto local_i = localIndices[0];
202
203 // Just set the pointer if common X bins. Otherwise, copy in the right chunk
204 // (as we do for Y).
205 if (!commonX)
206 outputWS->mutableX(local_i).assign(dataX.begin() + xStart, dataX.begin() + xEnd);
207
208 outputWS->mutableY(local_i).assign(dataY.begin() + yStart, dataY.begin() + yEnd);
209
210 if (dataE_provided)
211 outputWS->mutableE(local_i).assign(dataE.begin() + yStart, dataE.begin() + yEnd);
212
213 if (!dX.empty())
214 outputWS->mutableDx(local_i).assign(dX.begin() + yStart, dX.begin() + yEnd);
215
216 progress.report();
218 }
220
221 // Set the Unit of the X Axis
222 try {
223 outputWS->getAxis(0)->unit() = UnitFactory::Instance().create(xUnit);
224 } catch (Exception::NotFoundError &) {
225 outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("Label");
226 Unit_sptr unit = outputWS->getAxis(0)->unit();
227 std::shared_ptr<Units::Label> label = std::dynamic_pointer_cast<Units::Label>(unit);
228 label->setLabel(xUnit, xUnit);
229 }
230
231 // Populate the VerticalAxis. A spectra one is there by default with a 1->N
232 // mapping
233 if (vUnit != "SpectraNumber") {
234 if (vUnit == "Text") {
235 auto newAxis = std::make_unique<TextAxis>(vAxis.size());
236 auto newAxisRaw = newAxis.get();
237 outputWS->replaceAxis(1, std::move(newAxis));
238 for (size_t i = 0; i < vAxis.size(); i++) {
239 newAxisRaw->setLabel(i, vAxis[i]);
240 }
241 } else {
242 std::unique_ptr<NumericAxis> newAxis(nullptr);
243 if (vAxisSize == nSpec)
244 newAxis = std::make_unique<NumericAxis>(vAxisSize); // treat as points
245 else if (vAxisSize == nSpec + 1)
246 newAxis = std::make_unique<BinEdgeAxis>(vAxisSize); // treat as bin edges
247 else
248 throw std::range_error("Invalid vertical axis length. It must be the "
249 "same length as NSpec or 1 longer.");
250 auto newAxisRaw = newAxis.get();
251 newAxis->unit() = UnitFactory::Instance().create(vUnit);
252 outputWS->replaceAxis(1, std::move(newAxis));
253 for (size_t i = 0; i < vAxis.size(); i++) {
254 try {
255 newAxisRaw->setValue(i, boost::lexical_cast<double, std::string>(vAxis[i]));
256 } catch (boost::bad_lexical_cast &) {
257 throw std::invalid_argument("CreateWorkspace - YAxisValues property "
258 "could not be converted to a double.");
259 }
260 }
261 }
262 }
263
264 // Set Y Unit label
265 if (!parentWS || !getPropertyValue("YUnitLabel").empty()) {
266 outputWS->setYUnitLabel(getProperty("YUnitLabel"));
267 }
268
269 // Set Workspace Title
270 if (!parentWS || !getPropertyValue("WorkspaceTitle").empty()) {
271 outputWS->setTitle(getProperty("WorkspaceTitle"));
272 }
273
274 // Set OutputWorkspace property
275 setProperty("OutputWorkspace", outputWS);
276}
277
278Parallel::ExecutionMode
279CreateWorkspace::getParallelExecutionMode(const std::map<std::string, Parallel::StorageMode> &storageModes) const {
280 const auto storageMode = Parallel::fromString(getProperty("ParallelStorageMode"));
281 if (!storageModes.empty())
282 if (storageModes.begin()->second != storageMode)
283 throw std::invalid_argument("Input workspace storage mode differs from "
284 "requested output workspace storage mode.");
285 return Parallel::getCorrespondingExecutionMode(storageMode);
286}
287
288} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
const Parallel::Communicator & communicator() const
Returns a const reference to the (MPI) communicator of the algorithm.
Definition: Algorithm.cpp:1870
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
CreateWorkspace Algorithm.
std::map< std::string, std::string > validateInputs() override
Input validation.
void exec() override
Execute the Algorithm.
Parallel::ExecutionMode getParallelExecutionMode(const std::map< std::string, Parallel::StorageMode > &storageModes) const override
Get correct execution mode based on input storage modes for an MPI run.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Exception for when an item is not found in a collection.
Definition: Exception.h:145
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
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::shared_ptr< Unit > Unit_sptr
Shared pointer to the Unit base class.
Definition: Unit.h:229
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