Mantid
Loading...
Searching...
No Matches
MergeMD.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 +
15
16using namespace Mantid::Kernel;
17using namespace Mantid::API;
18using namespace Mantid::Geometry;
19using namespace Mantid::DataObjects;
20
21namespace Mantid::MDAlgorithms {
22
23// Register the algorithm into the AlgorithmFactory
24DECLARE_ALGORITHM(MergeMD)
25
26
27const std::string MergeMD::name() const { return "MergeMD"; }
28
30int MergeMD::version() const { return 1; }
31
33const std::string MergeMD::category() const { return "MDAlgorithms\\Creation"; }
34
38 // declare arbitrary number of input m_workspaces as a list of strings at the
39 // moment
41 "InputWorkspaces", std::make_shared<MandatoryValidator<std::vector<std::string>>>()),
42 "The names of the input MDWorkspaces as a comma-separated list");
43
44 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("OutputWorkspace", "", Direction::Output),
45 "Name of the output MDWorkspace.");
46
47 // Set the box controller properties
48 this->initBoxControllerProps("2", 500, 16);
49}
50
55void MergeMD::createOutputWorkspace(std::vector<std::string> &inputs) {
56 auto it = inputs.begin();
57 for (; it != inputs.end(); it++) {
59 std::dynamic_pointer_cast<IMDEventWorkspace>(AnalysisDataService::Instance().retrieve(*it));
60 if (!ws)
61 throw std::invalid_argument("Workspace " + *it + " is not a MDEventWorkspace. Cannot merge this workspace.");
62 else
63 m_workspaces.emplace_back(ws);
64 }
65 if (m_workspaces.empty())
66 throw std::invalid_argument("No valid m_workspaces specified.");
67
68 // Number of dimensions, start with the 0th one. They all have to match # of
69 // dims anyway
71 size_t numDims = ws0->getNumDims();
72 auto displNorm = ws0->displayNormalization();
73 auto displNormH = ws0->displayNormalizationHisto();
74
75 // Extents to create.
76 std::vector<coord_t> dimMin(numDims, +1e30f);
77 std::vector<coord_t> dimMax(numDims, -1e30f);
78
79 // Validate each workspace
80 for (auto &ws : m_workspaces) {
81 if (ws->getNumDims() != numDims)
82 throw std::invalid_argument(
83 "Workspace " + ws->getName() + " does not match the number of dimensions of the others (" +
84 Strings::toString(ws->getNumDims()) + ", expected " + Strings::toString(numDims) + ")");
85
86 if (ws->getEventTypeName() != ws0->getEventTypeName())
87 throw std::invalid_argument("Workspace " + ws->getName() + " does not match the MDEvent type of the others (" +
88 ws->getEventTypeName() + ", expected " + ws0->getEventTypeName() + ")");
89
90 for (size_t d = 0; d < numDims; d++) {
91 IMDDimension_const_sptr dim = ws->getDimension(d);
92 IMDDimension_const_sptr dim0 = ws0->getDimension(d);
93 if (dim->getName() != dim0->getName())
94 throw std::invalid_argument("Workspace " + ws->getName() + " does not have the same dimension " +
95 Strings::toString(d) + " as the others (" + dim->getName() + ", expected " +
96 dim0->getName() + ")");
97
98 // Find the extents
99 if (dim->getMaximum() > dimMax[d])
100 dimMax[d] = dim->getMaximum();
101 if (dim->getMinimum() < dimMin[d])
102 dimMin[d] = dim->getMinimum();
103 }
104 }
105
106 // OK, now create the blank MDWorkspace
107
108 // Have the factory create it
109 out = MDEventFactory::CreateMDWorkspace(numDims, ws0->getEventTypeName());
110 out->setDisplayNormalization(displNorm);
111 out->setDisplayNormalizationHisto(displNormH);
112
113 // Give all the dimensions
114 for (size_t d = 0; d < numDims; d++) {
115 IMDDimension_const_sptr dim0 = ws0->getDimension(d);
116 MDHistoDimension *dim = new MDHistoDimension(dim0->getName(), dim0->getDimensionId(), dim0->getMDFrame(), dimMin[d],
117 dimMax[d], dim0->getNBins());
118 out->addDimension(MDHistoDimension_sptr(dim));
119 }
120
121 // Initialize it using the dimension
122 out->initialize();
123
124 // Set the box controller settings from the properties
125 this->setBoxController(out->getBoxController());
126
127 // Perform the initial box splitting
128 out->splitBox();
129
130 // copy experiment infos
131 uint16_t nExperiments(0);
132 if (m_workspaces.size() > std::numeric_limits<uint16_t>::max())
133 throw std::invalid_argument("currently we can not combine more then 65535 experiments");
134 else
135 nExperiments = static_cast<uint16_t>(m_workspaces.size());
136
137 for (uint16_t i = 0; i < nExperiments; i++) {
138 uint16_t nWSexperiments = m_workspaces[i]->getNumExperimentInfo();
139 experimentInfoNo.emplace_back(nWSexperiments);
140 for (uint16_t j = 0; j < nWSexperiments; j++) {
142 API::ExperimentInfo_sptr(m_workspaces[i]->getExperimentInfo(j)->cloneExperimentInfo());
143 out->addExperimentInfo(ei);
144 }
145 }
146
147 // Cumulative sum of number of experimentInfo and reverse order
148 std::partial_sum(experimentInfoNo.begin(), experimentInfoNo.end(), experimentInfoNo.begin());
149 std::reverse(std::begin(experimentInfoNo), std::end(experimentInfoNo));
150}
151
152//----------------------------------------------------------------------------------------------
160template <size_t nd, size_t ond>
161inline void copyEvent(const MDLeanEvent<nd> &srcEvent, MDLeanEvent<ond> &newEvent, const uint16_t expInfoIndexOffset) {
162 // Nothing extra copy - this is no-op
163 UNUSED_ARG(srcEvent);
164 UNUSED_ARG(newEvent);
165 UNUSED_ARG(expInfoIndexOffset);
166}
167
168//----------------------------------------------------------------------------------------------
176template <size_t nd, size_t ond>
177inline void copyEvent(const MDEvent<nd> &srcEvent, MDEvent<ond> &newEvent, const uint16_t expInfoIndexOffset) {
178 newEvent.setDetectorId(srcEvent.getDetectorID());
179 newEvent.setExpInfoIndex(static_cast<uint16_t>(srcEvent.getExpInfoIndex() + expInfoIndexOffset));
180}
181
182//----------------------------------------------------------------------------------------------
188template <typename MDE, size_t nd> void MergeMD::doPlus(typename MDEventWorkspace<MDE, nd>::sptr ws2) {
189 // CPUTimer tim;
190 typename MDEventWorkspace<MDE, nd>::sptr ws1 = std::dynamic_pointer_cast<MDEventWorkspace<MDE, nd>>(out);
191 if (!ws1 || !ws2)
192 throw std::runtime_error("Incompatible workspace types passed to MergeMD.");
193
194 MDBoxBase<MDE, nd> *box1 = ws1->getBox();
195 MDBoxBase<MDE, nd> *box2 = ws2->getBox();
196
197 uint16_t expInfoIndexOffset = experimentInfoNo.back();
198 experimentInfoNo.pop_back();
199
200 // How many events you started with
201 size_t initial_numEvents = ws1->getNPoints();
202
203 // Make a leaf-only iterator through all boxes with events in the RHS
204 // workspace
205 std::vector<API::IMDNode *> boxes;
206 box2->getBoxes(boxes, 1000, true);
207 auto numBoxes = int(boxes.size());
208
209 bool fileBasedSource(false);
210 if (ws2->isFileBacked())
211 fileBasedSource = true;
212
213 // Add the boxes in parallel. They should be spread out enough on each
214 // core to avoid stepping on each other.
215
216 PRAGMA_OMP( parallel for if (!ws2->isFileBacked()) )
217 for (int i = 0; i < numBoxes; i++) {
219 auto *box = dynamic_cast<MDBox<MDE, nd> *>(boxes[i]);
220 if (box && !box->getIsMasked()) {
221 // Copy the events from WS2 and add them into WS1
222 const std::vector<MDE> &events = box->getConstEvents();
223 // Add events, with bounds checking
224
225 for (auto it = events.cbegin(); it != events.cend(); ++it) {
226 // Create the event
227 MDE newEvent(it->getSignal(), it->getErrorSquared(), it->getCenter());
228 // Copy extra data, if any
229 copyEvent(*it, newEvent, expInfoIndexOffset);
230 // Add it to the workspace
231 box1->addEvent(newEvent);
232 }
233 if (fileBasedSource)
234 box->clear();
235 else
236 box->releaseEvents();
237 }
239 }
241
242 // Progress * prog2 = new Progress(this, 0.4, 0.9, 100);
243 Progress *prog2 = nullptr;
245 ThreadPool tp(ts, 0, prog2);
246 ws1->splitAllIfNeeded(ts);
247 // prog2->resetNumSteps( ts->size(), 0.4, 0.6);
248 tp.joinAll();
249
250 // Set a marker that the file-back-end needs updating if the # of events
251 // changed.
252 if (ws1->getNPoints() != initial_numEvents)
253 ws1->setFileNeedsUpdating(true);
254 //
255 // std::cout << tim << " to add workspace " << ws2->name() << '\n';
256}
257
258//----------------------------------------------------------------------------------------------
262 CPUTimer tim;
263 // Check that all input workspaces exist and match in certain important ways
264 const std::vector<std::string> inputs_orig = getProperty("InputWorkspaces");
265
266 // This will hold the inputs, with the groups separated off
267 std::vector<std::string> inputs;
268 for (const auto &input : inputs_orig) {
270 if (wsgroup) { // Workspace group
271 std::vector<std::string> group = wsgroup->getNames();
272 inputs.insert(inputs.end(), group.begin(), group.end());
273 } else {
274 // Single workspace
275 inputs.emplace_back(input);
276 }
277 }
278
279 if (inputs.size() == 1) {
280 g_log.error("Only one input workspace specified");
281 throw std::invalid_argument("Only one input workspace specified");
282 }
283
284 // Create a blank output workspace
285 this->createOutputWorkspace(inputs);
286
287 // Run PlusMD on each of the input workspaces, in order.
288 double progStep = 1.0 / double(m_workspaces.size());
289 for (size_t i = 0; i < m_workspaces.size(); i++) {
290 g_log.information() << "Adding workspace " << m_workspaces[i]->getName() << '\n';
291 progress(double(i) * progStep, m_workspaces[i]->getName());
293 }
294
295 this->progress(0.95, "Refreshing cache");
296 out->refreshCache();
297
298 this->setProperty("OutputWorkspace", out);
299
300 g_log.debug() << tim << " to merge all workspaces.\n";
301}
302
303} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
static std::unique_ptr< QThreadPool > tp
#define CALL_MDEVENT_FUNCTION(funcname, workspace)
Macro that makes it possible to call a templated method for a MDEventWorkspace using a IMDEventWorksp...
#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 PRAGMA_OMP(expression)
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
std::string getName(const IMDDimension &self)
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
void setBoxController(const Mantid::API::BoxController_sptr &bc, const Mantid::Geometry::Instrument_const_sptr &instrument)
Set the settings in the given box controller.
void initBoxControllerProps(const std::string &SplitInto="5", int SplitThreshold=1000, int MaxRecursionDepth=5)
Initialise the properties.
void setFileNeedsUpdating(bool value)
Sets the marker set to true when a file-backed workspace needs its back-end file updated (by calling ...
virtual void getBoxes(std::vector< IMDNode * > &boxes, size_t maxDepth, bool leafOnly)=0
Fill a vector with all the boxes who are the childred of this one up to a certain depth.
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
Class to hold a set of workspaces.
std::vector< std::string > getNames() const
Returns the names of workspaces that make up this group.
A property class for workspaces.
Templated super-class of a multi-dimensional event "box".
Definition: MDBoxBase.h:50
virtual size_t addEvent(const MDE &point)=0
Add a single event.
Templated class for a multi-dimensional event "box".
Definition: MDBox.h:45
const std::vector< MDE > & getConstEvents() const
Get vector of constant events to use.
static API::IMDEventWorkspace_sptr CreateMDWorkspace(size_t nd, const std::string &eventType="MDLeanEvent", const Mantid::API::MDNormalization &preferredNormalization=Mantid::API::MDNormalization::VolumeNormalization, const Mantid::API::MDNormalization &preferredNormalizationHisto=Mantid::API::MDNormalization::VolumeNormalization)
Create a MDEventWorkspace of the given type.
uint64_t getNPoints() const override
std::shared_ptr< MDEventWorkspace< MDE, nd > > sptr
Typedef for a shared pointer of this kind of event workspace.
void splitAllIfNeeded(Kernel::ThreadScheduler *ts) override
Split all boxes that exceed the split threshold.
Templated class holding data about a neutron detection event in N-dimensions (for example,...
Definition: MDEvent.h:36
uint16_t getExpInfoIndex() const
Definition: MDEvent.h:175
void setDetectorId(int32_t id)
Sets the detectorId of this event.
Definition: MDEvent.h:195
void setExpInfoIndex(uint16_t index)
Sets the expInfoIndex of this event.
Definition: MDEvent.h:179
int32_t getDetectorID() const
Definition: MDEvent.h:191
Templated class holding data about a neutron detection event in N-dimensions (for example,...
Definition: MDLeanEvent.h:60
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
CPUTimer : Timer that uses the CPU time, rather than wall-clock time to measure execution time.
Definition: CPUTimer.h:24
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Validator to check that a property is not left empty.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A Thread Pool implementation that keeps a certain number of threads running (normally,...
Definition: ThreadPool.h:36
A First-In-First-Out Thread Scheduler.
The ThreadScheduler object defines how tasks are allocated to threads and in what order.
Merge several MDWorkspaces into one.
Definition: MergeMD.h:22
void doPlus(typename Mantid::DataObjects::MDEventWorkspace< MDE, nd >::sptr ws)
Perform the adding.
Definition: MergeMD.cpp:188
void init() override
Initialize the algorithm's properties.
Definition: MergeMD.cpp:37
int version() const override
Algorithm's version for identification.
Definition: MergeMD.cpp:30
void createOutputWorkspace(std::vector< std::string > &inputs)
Create the output MDWorkspace from a list of input.
Definition: MergeMD.cpp:55
std::vector< Mantid::API::IMDEventWorkspace_sptr > m_workspaces
Vector of input MDWorkspaces.
Definition: MergeMD.h:40
void exec() override
Execute the algorithm.
Definition: MergeMD.cpp:261
Mantid::API::IMDEventWorkspace_sptr out
Output MDEventWorkspace.
Definition: MergeMD.h:46
std::vector< uint16_t > experimentInfoNo
Vector of number of experimentalInfos for each input workspace.
Definition: MergeMD.h:43
const std::string category() const override
Algorithm's category for identification.
Definition: MergeMD.cpp:33
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< const IMDEventWorkspace > IMDEventWorkspace_const_sptr
Shared pointer to Mantid::API::IMDEventWorkspace (const version)
std::shared_ptr< WorkspaceGroup > WorkspaceGroup_sptr
shared pointer to Mantid::API::WorkspaceGroup
std::shared_ptr< ExperimentInfo > ExperimentInfo_sptr
Shared pointer to ExperimentInfo.
std::shared_ptr< MDHistoDimension > MDHistoDimension_sptr
Shared pointer to a MDHistoDimension.
std::shared_ptr< const IMDDimension > IMDDimension_const_sptr
Shared Pointer to const IMDDimension.
Definition: IMDDimension.h:101
std::string toString(const T &value)
Convert a number to a string.
Definition: Strings.cpp:703
void copyEvent(const MDLeanEvent< nd > &srcEvent, MDLeanEvent< ond > &newEvent, const uint16_t expInfoIndexOffset)
Copy the extra data (not signal, error or coordinates) from one event to another with different numbe...
Definition: MergeMD.cpp:161
STL namespace.
@ Output
An output workspace.
Definition: Property.h:54