Mantid
Loading...
Searching...
No Matches
BinMD.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 +
22#include "MantidKernel/System.h"
23#include "MantidKernel/Utils.h"
24#include <boost/algorithm/string.hpp>
25
26namespace Mantid::MDAlgorithms {
27
28// Register the algorithm into the AlgorithmFactory
30
31using namespace Mantid::Kernel;
32using namespace Mantid::API;
33using namespace Mantid::Geometry;
34using namespace Mantid::DataObjects;
35
36//----------------------------------------------------------------------------------------------
40 : SlicingAlgorithm(), outWS(), implicitFunction(), indexMultiplier(), signals(nullptr), errors(nullptr),
41 numEvents(nullptr) {}
42
43//----------------------------------------------------------------------------------------------
47 declareProperty(std::make_unique<WorkspaceProperty<IMDWorkspace>>("InputWorkspace", "", Direction::Input),
48 "An input MDWorkspace.");
49
50 // Properties for specifying the slice to perform.
51 this->initSlicingProps();
52
53 // --------------- Processing methods and options
54 // ---------------------------------------
55 std::string grp = "Methods";
56 declareProperty(std::make_unique<PropertyWithValue<std::string>>("ImplicitFunctionXML", "", Direction::Input),
57 "XML string describing the implicit function determining "
58 "which bins to use.");
59 setPropertyGroup("ImplicitFunctionXML", grp);
60
61 declareProperty(std::make_unique<PropertyWithValue<bool>>("IterateEvents", true, Direction::Input),
62 "Alternative binning method where you iterate through every event, "
63 "placing them in the proper bin.\n"
64 "This may be faster for workspaces with few events and lots of output "
65 "bins.");
66 setPropertyGroup("IterateEvents", grp);
67
68 declareProperty(std::make_unique<PropertyWithValue<bool>>("Parallel", false, Direction::Input),
69 "Temporary parameter: true to run in parallel. This is ignored for "
70 "file-backed workspaces, where running in parallel makes things slower "
71 "due to disk thrashing.");
72 setPropertyGroup("Parallel", grp);
73
74 declareProperty(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("TemporaryDataWorkspace", "", Direction::Input,
76 "An input MDHistoWorkspace used to accumulate results from "
77 "multiple MDEventWorkspaces. If unspecified a blank "
78 "MDHistoWorkspace will be created.");
79
80 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Direction::Output),
81 "A name for the output MDHistoWorkspace.");
82}
83
84//----------------------------------------------------------------------------------------------
93template <typename MDE, size_t nd>
94inline void BinMD::binMDBox(MDBox<MDE, nd> *box, const size_t *const chunkMin, const size_t *const chunkMax) {
95 // An array to hold the rotated/transformed coordinates
96 auto outCenter = std::vector<coord_t>(m_outD);
97
98 // Evaluate whether the entire box is in the same bin
99 if (box->getNPoints() > (1 << nd) * 2) {
100 // There is a check that the number of events is enough for it to make sense
101 // to do all this processing.
102 size_t numVertexes = 0;
103 auto vertexes = box->getVertexesArray(numVertexes);
104
105 // All vertexes have to be within THE SAME BIN = have the same linear index.
106 size_t lastLinearIndex = 0;
107 bool badOne = false;
108
109 for (size_t i = 0; i < numVertexes; i++) {
110 // Cache the center of the event (again for speed)
111 const coord_t *inCenter = vertexes.get() + i * nd;
112
113 // Now transform to the output dimensions
114 m_transform->apply(inCenter, outCenter.data());
115
116 // To build up the linear index
117 size_t linearIndex = 0;
118 // To mark VERTEXES outside range
119 badOne = false;
120
122 for (size_t bd = 0; bd < m_outD; bd++) {
123 // What is the bin index in that dimension
124 coord_t x = outCenter[bd];
125 auto ix = size_t(x);
126 // Within range (for this chunk)?
127 if ((x >= 0) && (ix >= chunkMin[bd]) && (ix < chunkMax[bd])) {
128 // Build up the linear index
129 linearIndex += indexMultiplier[bd] * ix;
130 } else {
131 // Outside the range
132 badOne = true;
133 break;
134 }
135 } // (for each dim in MDHisto)
136
137 // Is the vertex at the same place as the last one?
138 if (!badOne) {
139 if ((i > 0) && (linearIndex != lastLinearIndex)) {
140 // Change of index
141 badOne = true;
142 break;
143 }
144 lastLinearIndex = linearIndex;
145 }
146
147 // Was the vertex completely outside the range?
148 if (badOne)
149 break;
150 } // (for each vertex)
151
152 if (!badOne) {
153 // Yes, the entire box is within a single bin
154 // std::cout << "Box at " << box->getExtentsStr() << " is within a
155 // single bin.\n";
156 // Add the CACHED signal from the entire box
157 signals[lastLinearIndex] += box->getSignal();
158 errors[lastLinearIndex] += box->getErrorSquared();
159 // TODO: If DataObjects get a weight, this would need to get the summed
160 // weight.
161 numEvents[lastLinearIndex] += static_cast<signal_t>(box->getNPoints());
162
163 // And don't bother looking at each event. This may save lots of time
164 // loading from disk.
165 return;
166 }
167 }
168
169 // If you get here, you could not determine that the entire box was in the
170 // same bin.
171 // So you need to iterate through events.
172 const std::vector<MDE> &events = box->getConstEvents();
173 for (auto it = events.begin(); it != events.end(); ++it) {
174 // Cache the center of the event (again for speed)
175 const coord_t *inCenter = it->getCenter();
176
177 // Now transform to the output dimensions
178 m_transform->apply(inCenter, outCenter.data());
179
180 // To build up the linear index
181 size_t linearIndex = 0;
182 // To mark events outside range
183 bool badOne = false;
184
186 for (size_t bd = 0; bd < m_outD; bd++) {
187 // What is the bin index in that dimension
188 coord_t x = outCenter[bd];
189 auto ix = size_t(x);
190 // Within range (for this chunk)?
191 if ((x >= 0) && (ix >= chunkMin[bd]) && (ix < chunkMax[bd])) {
192 // Build up the linear index
193 linearIndex += indexMultiplier[bd] * ix;
194 } else {
195 // Outside the range
196 badOne = true;
197 break;
198 }
199 } // (for each dim in MDHisto)
200
201 if (!badOne) {
202 // Sum the signals as doubles to preserve precision
203 signals[linearIndex] += static_cast<signal_t>(it->getSignal());
204 errors[linearIndex] += static_cast<signal_t>(it->getErrorSquared());
205 // TODO: If DataObjects get a weight, this would need to get the summed
206 // weight.
207 numEvents[linearIndex] += 1.0;
208 }
209 }
210 // Done with the events list
211 box->releaseEvents();
212}
213
214//----------------------------------------------------------------------------------------------
220template <typename MDE, size_t nd> void BinMD::binByIterating(typename MDEventWorkspace<MDE, nd>::sptr ws) {
222 // store exisiting write buffer size for the future
223 // uint64_t writeBufSize =bc->getDiskBuffer().getWriteBufferSize();
224 // and disable write buffer (if any) for input MD Events for this algorithm
225 // purposes;
226 // bc->setCacheParameters(1,0);
227
228 // Cache some data to speed up accessing them a bit
229 indexMultiplier.resize(m_outD);
230 for (size_t d = 0; d < m_outD; d++) {
231 if (d > 0)
232 indexMultiplier[d] = outWS->getIndexMultiplier()[d - 1];
233 else
234 indexMultiplier[d] = 1;
235 }
236 signals = outWS->mutableSignalArray();
237 errors = outWS->mutableErrorSquaredArray();
238 numEvents = outWS->mutableNumEventsArray();
239
240 if (!m_accumulate) {
241 // Start with signal/error/numEvents at 0.0
242 outWS->setTo(0.0, 0.0, 0.0);
243 }
244
245 // The dimension (in the output workspace) along which we chunk for parallel
246 // processing
247 // TODO: Find the smartest dimension to chunk against
248 size_t chunkDimension = 0;
249
250 // How many bins (in that dimension) per chunk.
251 // Try to split it so each core will get 2 tasks:
252 auto chunkNumBins = int(m_binDimensions[chunkDimension]->getNBins() / (PARALLEL_GET_MAX_THREADS * 2));
253 if (chunkNumBins < 1)
254 chunkNumBins = 1;
255
256 // Do we actually do it in parallel?
257 bool doParallel = getProperty("Parallel");
258 // Not if file-backed!
259 if (bc->isFileBacked())
260 doParallel = false;
261 if (!doParallel)
262 chunkNumBins = int(m_binDimensions[chunkDimension]->getNBins());
263
264 // Total number of steps
265 size_t progNumSteps = 0;
266 if (prog) {
267 prog->setNotifyStep(0.1);
268 prog->resetNumSteps(100, 0.00, 1.0);
269 }
270
271 // Run the chunks in parallel. There is no overlap in the output workspace so
272 // it is thread safe to write to it..
273
274 PRAGMA_OMP( parallel for schedule(dynamic,1) if (doParallel) )
275 for (int chunk = 0; chunk < int(m_binDimensions[chunkDimension]->getNBins()); chunk += chunkNumBins) {
277 // Region of interest for this chunk.
278 std::vector<size_t> chunkMin(m_outD);
279 std::vector<size_t> chunkMax(m_outD);
280 for (size_t bd = 0; bd < m_outD; bd++) {
281 // Same limits in the other dimensions
282 chunkMin[bd] = 0;
283 chunkMax[bd] = m_binDimensions[bd]->getNBins();
284 }
285 // Parcel out a chunk in that single dimension dimension
286 chunkMin[chunkDimension] = size_t(chunk);
287 if (size_t(chunk + chunkNumBins) > m_binDimensions[chunkDimension]->getNBins())
288 chunkMax[chunkDimension] = m_binDimensions[chunkDimension]->getNBins();
289 else
290 chunkMax[chunkDimension] = size_t(chunk + chunkNumBins);
291
292 // Build an implicit function (it needs to be in the space of the
293 // MDEventWorkspace)
294 auto function = this->getImplicitFunctionForChunk(chunkMin.data(), chunkMax.data());
295
296 // Use getBoxes() to get an array with a pointer to each box
297 std::vector<API::IMDNode *> boxes;
298 // Leaf-only; no depth limit; with the implicit function passed to it.
299 ws->getBox()->getBoxes(boxes, 1000, true, function.get());
300
301 // Sort boxes by file position IF file backed. This reduces seeking time,
302 // hopefully.
303 if (bc->isFileBacked())
305
306 // For progress reporting, the # of boxes
307 if (prog) {
308 PARALLEL_CRITICAL(BinMD_progress) {
309 g_log.debug() << "Chunk " << chunk << ": found " << boxes.size() << " boxes within the implicit function.\n";
310 progNumSteps += boxes.size();
311 prog->setNumSteps(progNumSteps);
312 }
313 }
314
315 // Go through every box for this chunk.
316 for (auto &boxe : boxes) {
317 auto *box = dynamic_cast<MDBox<MDE, nd> *>(boxe);
318 // Perform the binning in this separate method.
319 if (box && !box->getIsMasked())
320 this->binMDBox(box, chunkMin.data(), chunkMax.data());
321
322 // Progress reporting
323 if (prog)
324 prog->report();
325 // For early cancelling of the loop
326 if (this->m_cancel)
327 break;
328 } // for each box in the vector
330 } // for each chunk in parallel
332
333 // Now the implicit function
334 if (implicitFunction) {
335 if (prog)
336 prog->report("Applying implicit function.");
337 signal_t nan = std::numeric_limits<signal_t>::quiet_NaN();
338 outWS->applyImplicitFunction(implicitFunction.get(), nan, nan);
339 }
340}
341
342//----------------------------------------------------------------------------------------------
346 // Input MDEventWorkspace/MDHistoWorkspace
347 m_inWS = getProperty("InputWorkspace");
348 // Look at properties, create either axis-aligned or general transform.
349 // This (can) change m_inWS
350 this->createTransform();
351
352 // De serialize the implicit function
353 std::string ImplicitFunctionXML = getPropertyValue("ImplicitFunctionXML");
354 implicitFunction = nullptr;
355 if (!ImplicitFunctionXML.empty())
356 implicitFunction = std::unique_ptr<MDImplicitFunction>(
357 Mantid::API::ImplicitFunctionFactory::Instance().createUnwrapped(ImplicitFunctionXML));
358
359 // This gets deleted by the thread pool; don't delete it in here.
360 prog = std::make_unique<Progress>(this, 0.0, 1.0, 1);
361
362 // Create the dense histogram. This allocates the memory
363 std::shared_ptr<IMDHistoWorkspace> tmp = this->getProperty("TemporaryDataWorkspace");
364 outWS = std::dynamic_pointer_cast<MDHistoWorkspace>(tmp);
365 if (!outWS) {
366 outWS = std::make_shared<MDHistoWorkspace>(m_binDimensions);
367 } else {
368 m_accumulate = true;
369 }
370
371 // Saves the geometry transformation from original to binned in the workspace
372 outWS->setTransformFromOriginal(this->m_transformFromOriginal.release(), 0);
373 outWS->setTransformToOriginal(this->m_transformToOriginal.release(), 0);
374 for (size_t i = 0; i < m_bases.size(); i++)
375 outWS->setBasisVector(i, m_bases[i]);
376 outWS->setOrigin(this->m_translation);
377 outWS->setOriginalWorkspace(m_inWS, 0);
378
379 // And the intermediate WS one too, if any
380 if (m_intermediateWS) {
381 outWS->setOriginalWorkspace(m_intermediateWS, 1);
382 outWS->setTransformFromOriginal(m_transformFromIntermediate.release(), 1);
383 outWS->setTransformToOriginal(m_transformToIntermediate.release(), 1);
384 }
385
386 // Wrapper to cast to MDEventWorkspace then call the function
387 bool IterateEvents = getProperty("IterateEvents");
388 if (!IterateEvents) {
389 g_log.warning() << "IterateEvents=False is no longer supported. Setting "
390 "IterateEvents=True.\n";
391 IterateEvents = true;
392 }
393
394 /*
395 We should fail noisily here. CALL_MDEVENT_FUNCTION will silently allow
396 IMDHistoWorkspaces to cascade through to the end
397 and result in an empty output. The only way we allow InputWorkspaces to be
398 IMDHistoWorkspaces is if they also happen to contain original workspaces
399 that are MDEventWorkspaces.
400 */
401 if (m_inWS->isMDHistoWorkspace()) {
402 throw std::runtime_error("Cannot rebin a workspace that is histogrammed and has no original "
403 "workspace that is an MDEventWorkspace. "
404 "Reprocess the input so that it contains full MDEvents.");
405 }
406
408
409 // Copy the coordinate system & experiment infos to the output
410 IMDEventWorkspace_sptr inEWS = std::dynamic_pointer_cast<IMDEventWorkspace>(m_inWS);
411 if (inEWS) {
412 outWS->setCoordinateSystem(inEWS->getSpecialCoordinateSystem());
413 try {
414 outWS->copyExperimentInfos(*inEWS);
415 } catch (std::runtime_error &) {
416 g_log.warning() << this->name() << " was not able to copy experiment info to output workspace "
417 << outWS->getName() << '\n';
418 }
419 }
420
421 // Pass on the display normalization from the input workspace
422 outWS->setDisplayNormalization(m_inWS->displayNormalizationHisto());
423
424 outWS->updateSum();
425 // Save the output
426 setProperty("OutputWorkspace", std::dynamic_pointer_cast<Workspace>(outWS));
427}
428
429} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
gsl_vector * tmp
#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_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 PRAGMA_OMP(expression)
#define PARALLEL_GET_MAX_THREADS
#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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
std::atomic< bool > m_cancel
Set to true to stop execution.
Definition: Algorithm.h:433
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
static void sortObjByID(std::vector< IMDNode * > &boxes)
Static method for sorting a list of MDBoxBase pointers by their file position, ascending.
Definition: IMDNode.h:292
A property class for workspaces.
std::unique_ptr< coord_t[]> getVertexesArray(size_t &numVertices) const override
signal_t getSignal() const override
Returns the integrated signal from all points within.
Definition: MDBoxBase.h:232
signal_t getErrorSquared() const override
Returns the integrated error squared from all points within.
Definition: MDBoxBase.h:242
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.
uint64_t getNPoints() const override
std::shared_ptr< MDEventWorkspace< MDE, nd > > sptr
Typedef for a shared pointer of this kind of event workspace.
Mantid::API::BoxController_sptr getBoxController() override
Returns the BoxController used in this workspace.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
void binMDBox(DataObjects::MDBox< MDE, nd > *box, const size_t *const chunkMin, const size_t *const chunkMax)
Method to bin a single MDBox.
Definition: BinMD.cpp:94
std::unique_ptr< Mantid::Geometry::MDImplicitFunction > implicitFunction
ImplicitFunction used.
Definition: BinMD.h:73
signal_t * numEvents
Definition: BinMD.h:79
std::vector< size_t > indexMultiplier
Cached values for speed up.
Definition: BinMD.h:76
void exec() override
Run the algorithm.
Definition: BinMD.cpp:345
const std::string name() const override
Algorithm's name for identification.
Definition: BinMD.h:42
std::unique_ptr< Mantid::API::Progress > prog
Progress reporting.
Definition: BinMD.h:71
void binByIterating(typename DataObjects::MDEventWorkspace< MDE, nd >::sptr ws)
Helper method.
Definition: BinMD.cpp:220
BinMD()
Constructor.
Definition: BinMD.cpp:39
void init() override
Initialise the properties.
Definition: BinMD.cpp:46
Mantid::DataObjects::MDHistoWorkspace_sptr outWS
The output MDHistoWorkspace.
Definition: BinMD.h:69
Abstract Algorithm class that will be used by: BinMD and SliceMD and shares code for getting a slice ...
std::unique_ptr< API::CoordTransform > m_transformFromOriginal
Coordinate transformation to save in the output workspace (original->binned)
Mantid::API::IMDWorkspace_sptr m_intermediateWS
Intermediate original workspace.
std::unique_ptr< Mantid::Geometry::MDImplicitFunction > getImplicitFunctionForChunk(const size_t *const chunkMin, const size_t *const chunkMax)
Create an implicit function for picking boxes, based on the indexes in the output MDHistoWorkspace.
std::unique_ptr< API::CoordTransform > m_transform
Coordinate transformation to apply.
void createTransform()
Read the algorithm properties and creates the appropriate transforms for slicing the MDEventWorkspace...
Mantid::API::IMDWorkspace_sptr m_inWS
Input workspace.
std::unique_ptr< API::CoordTransform > m_transformToOriginal
Coordinate transformation to save in the output workspace (binned->original)
std::unique_ptr< DataObjects::CoordTransformAffine > m_transformFromIntermediate
Coordinate transformation to save in the output WS, from the intermediate WS.
size_t m_outD
Number of dimensions in the output (binned) workspace.
std::vector< Mantid::Kernel::VMD > m_bases
Basis vectors of the output dimensions, normalized to unity length.
void initSlicingProps()
Initialise the properties.
Mantid::Kernel::VMD m_translation
Translation from the OUTPUT to the INPUT workspace i.e.
std::unique_ptr< DataObjects::CoordTransformAffine > m_transformToIntermediate
Coordinate transformation to save in the intermediate WS.
std::vector< Mantid::Geometry::MDHistoDimension_sptr > m_binDimensions
Bin dimensions to actually use.
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< BoxController > BoxController_sptr
Shared ptr to BoxController.
std::size_t numEvents(::NeXus::File &file, bool &hasTotalCounts, bool &oldNeXusFileNames, const std::string &prefix, const NexusHDF5Descriptor &descriptor)
Get the number of events in the currently opened group.
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
Definition: MDTypes.h:27
double signal_t
Typedef for the signal recorded in a MDBox, etc.
Definition: MDTypes.h:36
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54