Mantid
Loading...
Searching...
No Matches
CropWorkspaceRagged.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 +
13
14#include <algorithm>
15
16namespace {
17std::vector<double> getSubVector(Mantid::MantidVec &data, const int64_t &lowerIndex, const int64_t &upperIndex) {
18 auto low = std::next(data.begin(), lowerIndex);
19 auto up = std::next(data.begin(), upperIndex);
20 // get new vectors
21 std::vector<double> newData(low, up);
22 return newData;
23}
24
25} // namespace
26
27namespace Mantid::Algorithms {
28
29using namespace Kernel;
30using namespace API;
31
32DECLARE_ALGORITHM(CropWorkspaceRagged)
33
34
35void CropWorkspaceRagged::init() {
36 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
37 "The input workspace");
38 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
39 "Name to be given to the cropped workspace.");
40
41 auto required = std::make_shared<MandatoryValidator<std::vector<double>>>();
42 declareProperty(std::make_unique<ArrayProperty<double>>("XMin", required),
43 "The value(s) to start the cropping from. Should be either a "
44 "single value or a list.");
45 declareProperty(std::make_unique<ArrayProperty<double>>("XMax", required),
46 "The value(s) to end the cropping at. Should be either a "
47 "single value or a list.");
48}
49
51std::map<std::string, std::string> CropWorkspaceRagged::validateInputs() {
52 std::map<std::string, std::string> issues;
53 MatrixWorkspace_sptr ws = getProperty("InputWorkspace");
54 auto numSpectra = ws->getNumberHistograms();
55 std::vector<double> xMin = getProperty("XMin");
56 std::vector<double> xMax = getProperty("XMax");
57 if (xMin.size() == 0 || (xMin.size() != numSpectra && xMin.size() > 1)) {
58 issues["XMin"] = "XMin must be a single value or one value per sepctrum.";
59 }
60 if (xMax.size() == 0 || (xMax.size() > 1 && xMax.size() != numSpectra)) {
61 issues["XMax"] = "XMax must be a single value or one value per sepctrum.";
62 }
63 if (xMin.size() == 1 && xMax.size() == 1 && xMin[0] > xMax[0]) {
64 issues["XMax"] = "XMax must be greater than XMin.";
65 } else if (xMin.size() == 1 && xMax.size() > 1) {
66 auto it = std::find_if(xMax.cbegin(), xMax.cend(), [&xMin](auto max) { return max < xMin[0]; });
67 if (it != xMax.cend()) {
68 issues["XMax"] = "XMax must be greater than XMin.";
69 return issues;
70 }
71 } else if (xMin.size() > 1 && xMax.size() == 1) {
72 auto it = std::find_if(xMin.cbegin(), xMin.cend(), [&xMax](auto min) { return min > xMax[0]; });
73 if (it != xMin.cend()) {
74 issues["XMin"] = "XMin must be less than XMax.";
75 return issues;
76 }
77 } else if (xMin.size() > 1 && xMax.size() > 1) {
78 for (size_t k = 0; k < xMin.size(); k++) {
79 if (xMin[k] > xMax[k]) {
80 issues["XMin"] = "XMin must be less than XMax.";
81 return issues;
82 }
83 }
84 }
85 return issues;
86} // namespace Algorithms
87
90 MatrixWorkspace_sptr ws = getProperty("InputWorkspace");
91 auto numSpectra = ws->getNumberHistograms();
92 // clone ws to copy logs etc.
93 MatrixWorkspace_sptr outputWS = ws->clone();
94
95 std::vector<double> xMin = getProperty("XMin");
96 std::vector<double> xMax = getProperty("XMax");
97 if (xMin.size() == 1) {
98 auto value = xMin[0];
99 xMin.assign(numSpectra, value);
100 }
101 if (xMax.size() == 1) {
102 auto value = xMax[0];
103 xMax.assign(numSpectra, value);
104 }
105
106 // Its easier to work with point data -> index is same for x, y, E
107 MatrixWorkspace_sptr tmp = outputWS;
108 bool histogram = false;
109 if (outputWS->isHistogramData()) {
110 auto alg = createChildAlgorithm("ConvertToPointData");
111 alg->initialize();
112 alg->setRethrows(true);
113 alg->setProperty("InputWorkspace", outputWS);
114 alg->setProperty("OutputWorkspace", outputWS);
115 alg->execute();
116 tmp = alg->getProperty("OutputWorkspace");
117 histogram = true;
118 }
120 for (int64_t i = 0; i < int64_t(numSpectra); ++i) {
122 auto points = tmp->points(i);
123 auto &dataX = outputWS->dataX(i);
124 auto &dataY = outputWS->dataY(i);
125 auto &dataE = outputWS->dataE(i);
126
127 // get iterators for cropped region using points
128 auto low = std::lower_bound(points.begin(), points.end(), xMin[i]);
129 auto up = std::upper_bound(points.begin(), points.end(), xMax[i]);
130 // convert to index
131 int64_t lowerIndex = std::distance(points.begin(), low);
132 int64_t upperIndex = std::distance(points.begin(), up);
133
134 // get new vectors
135 std::vector<double> newY = getSubVector(dataY, lowerIndex, upperIndex);
136 std::vector<double> newE = getSubVector(dataE, lowerIndex, upperIndex);
137 if (histogram && upperIndex + (size_t)1 <= dataX.size()) {
138 // the offset adds one to the upper index for histograms
139 // only use the offset if the end is cropped
140 upperIndex += 1;
141 }
142 std::vector<double> newX = getSubVector(dataX, lowerIndex, upperIndex);
143
144 // resize the data
145 dataX.resize(newX.size());
146 dataY.resize(newY.size());
147 dataE.resize(newE.size());
148
149 // update the data
150 outputWS->mutableX(i) = newX;
151 outputWS->mutableY(i) = newY;
152 outputWS->mutableE(i) = newE;
154 }
156
157 setProperty("OutputWorkspace", outputWS);
158}
159
160} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
gsl_vector * tmp
double value
The value of the point.
Definition: FitMW.cpp:51
#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.
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
A property class for workspaces.
Extracts a 'block' from a workspace and places it in a new workspace.
std::map< std::string, std::string > validateInputs() override
Input validation.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54