/
SliceMD.cpp
324 lines (277 loc) · 12.5 KB
/
SliceMD.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidMDAlgorithms/SliceMD.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidGeometry/MDGeometry/MDHistoDimension.h"
#include "MantidGeometry/MDGeometry/MDImplicitFunction.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidKernel/ThreadPool.h"
#include "MantidKernel/ThreadScheduler.h"
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace Mantid::DataObjects;
namespace Mantid::MDAlgorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SliceMD)
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void SliceMD::init() {
declareProperty(std::make_unique<WorkspaceProperty<IMDWorkspace>>("InputWorkspace", "", Direction::Input),
"An input MDWorkspace.");
// Properties for specifying the slice to perform.
this->initSlicingProps();
declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Direction::Output),
"Name of the output MDEventWorkspace.");
declareProperty(std::make_unique<FileProperty>("OutputFilename", "", FileProperty::OptionalSave, ".nxs"),
"Optional: Specify a NeXus file to write if you want the output "
"workspace to be file-backed.");
declareProperty(std::make_unique<PropertyWithValue<int>>("Memory", -1),
"If OutputFilename is specified to use a file back end:\n"
" The amount of memory (in MB) to allocate to the in-memory cache.\n"
" If not specified, a default of 40% of free physical memory is used.");
declareProperty("TakeMaxRecursionDepthFromInput", true, "Copy the maximum recursion depth from the input workspace.");
auto mustBePositiveInteger = std::make_shared<BoundedValidator<int>>();
mustBePositiveInteger->setLower(0);
declareProperty("MaxRecursionDepth", 1000, mustBePositiveInteger,
"Sets the maximum recursion depth to use. Can be used to "
"constrain the workspaces internal structure");
setPropertySettings("MaxRecursionDepth",
std::make_unique<EnabledWhenProperty>("TakeMaxRecursionDepthFromInput", IS_EQUAL_TO, "0"));
setPropertyGroup("OutputFilename", "File Back-End");
setPropertyGroup("Memory", "File Back-End");
}
//----------------------------------------------------------------------------------------------
/** Copy the extra data (not signal, error or coordinates) from one event to
*another
* with different numbers of dimensions
*
* @param srcEvent :: the source event, being copied
* @param newEvent :: the destination event
*/
template <size_t nd, size_t ond> inline void copyEvent(const MDLeanEvent<nd> &srcEvent, MDLeanEvent<ond> &newEvent) {
// Nothing extra copy - this is no-op
UNUSED_ARG(srcEvent);
UNUSED_ARG(newEvent);
}
//----------------------------------------------------------------------------------------------
/** Copy the extra data (not signal, error or coordinates) from one event to
*another
* with different numbers of dimensions
*
* @param srcEvent :: the source event, being copied
* @param newEvent :: the destination event
*/
template <size_t nd, size_t ond> inline void copyEvent(const MDEvent<nd> &srcEvent, MDEvent<ond> &newEvent) {
newEvent.setDetectorId(srcEvent.getDetectorID());
newEvent.setExpInfoIndex(srcEvent.getExpInfoIndex());
}
//----------------------------------------------------------------------------------------------
/** Perform the slice from nd input dimensions to ond output dimensions
*
* @param ws :: input workspace with nd dimensions
* @tparam OMDE :: MDEvent type for the OUTPUT workspace
* @tparam ond :: number of dimensions in the OUTPUT workspace
*/
template <typename MDE, size_t nd, typename OMDE, size_t ond>
void SliceMD::slice(typename MDEventWorkspace<MDE, nd>::sptr ws) {
// Create the ouput workspace
typename MDEventWorkspace<OMDE, ond>::sptr outWS(new MDEventWorkspace<OMDE, ond>());
for (auto &binDimension : m_binDimensions) {
outWS->addDimension(binDimension);
}
outWS->setCoordinateSystem(ws->getSpecialCoordinateSystem());
outWS->initialize();
// Copy settings from the original box controller
BoxController_sptr bc = ws->getBoxController();
// store wrute buffer size for the future
// uint64_t writeBufSize =
// bc->getFileIO()getDiskBuffer().getWriteBufferSize();
// and disable write buffer (if any) for input MD Events for this algorithm
// purposes;
// bc->setCacheParameters(1,0);
BoxController_sptr obc = outWS->getBoxController();
// Use the "number of bins" as the "split into" parameter for the top level
for (size_t od = 0; od < m_binDimensions.size(); od++) {
obc->setSplitTopInto(od, m_binDimensions[od]->getNBins());
obc->setSplitInto(od, bc->getSplitInto(od));
}
obc->setSplitThreshold(bc->getSplitThreshold());
bool bTakeDepthFromInputWorkspace = getProperty("TakeMaxRecursionDepthFromInput");
int tempDepth = getProperty("MaxRecursionDepth");
size_t maxDepth = bTakeDepthFromInputWorkspace ? bc->getMaxDepth() : size_t(tempDepth);
obc->setMaxDepth(maxDepth);
// size_t outputSize = writeBufSize;
// obc->setCacheParameters(sizeof(OMDE),outputSize);
obc->resetNumBoxes();
// Perform the first box splitting
outWS->splitBox();
size_t lastNumBoxes = obc->getTotalNumMDBoxes();
// --- File back end ? ----------------
std::string filename = getProperty("OutputFilename");
if (!filename.empty()) {
// First save to the NXS file
g_log.notice() << "Running SaveMD to create file back-end\n";
auto alg = createChildAlgorithm("SaveMD");
alg->setPropertyValue("Filename", filename);
alg->setProperty("InputWorkspace", outWS);
alg->setProperty("MakeFileBacked", true);
alg->executeAsChildAlg();
if (!obc->isFileBacked())
throw std::runtime_error("SliceMD with file-backed output: Can not set "
"up file-backed output workspace ");
auto IOptr = obc->getFileIO();
size_t outBufSize = IOptr->getWriteBufferSize();
// the buffer size for resulting workspace; reasonable size is at least 10
// data chunk sizes (nice to verify)
if (outBufSize < 10 * IOptr->getDataChunk()) {
outBufSize = 10 * IOptr->getDataChunk();
IOptr->setWriteBufferSize(outBufSize);
}
}
// Function defining which events (in the input dimensions) to place in the
// output
auto function = this->getImplicitFunctionForChunk(nullptr, nullptr);
std::vector<API::IMDNode *> boxes;
// Leaf-only; no depth limit; with the implicit function passed to it.
ws->getBox()->getBoxes(boxes, 1000, true, function.get());
// Sort boxes by file position IF file backed. This reduces seeking time,
// hopefully.
bool fileBackedWS = bc->isFileBacked();
if (fileBackedWS)
API::IMDNode::sortObjByID(boxes);
auto prog = std::make_unique<Progress>(this, 0.0, 1.0, boxes.size());
// The root of the output workspace
MDBoxBase<OMDE, ond> *outRootBox = outWS->getBox();
// if target workspace has events, we should count them as added
uint64_t totalAdded = outWS->getNEvents();
uint64_t numSinceSplit = 0;
// Go through every box for this chunk.
// PARALLEL_FOR_IF( !bc->isFileBacked() )
for (int i = 0; i < int(boxes.size()); i++) {
auto *box = dynamic_cast<MDBox<MDE, nd> *>(boxes[i]);
// Perform the binning in this separate method.
if (box && !box->getIsMasked()) {
// An array to hold the rotated/transformed coordinates
coord_t outCenter[ond];
const std::vector<MDE> &events = box->getConstEvents();
for (auto it = events.cbegin(); it != events.cend(); ++it) {
// Cache the center of the event (again for speed)
const coord_t *inCenter = it->getCenter();
if (function->isPointContained(inCenter)) {
// Now transform to the output dimensions
m_transformFromOriginal->apply(inCenter, outCenter);
// Create the event
OMDE newEvent(it->getSignal(), it->getErrorSquared(), outCenter);
// Copy extra data, if any
copyEvent(*it, newEvent);
// Add it to the workspace
if (outRootBox->addEvent(newEvent))
numSinceSplit++;
}
}
box->releaseEvents();
// Ask BC if one needs to split boxes
if (obc->shouldSplitBoxes(totalAdded, numSinceSplit, lastNumBoxes))
// if (numSinceSplit > 20000000 || (i == int(boxes.size()-1)))
{
// This splits up all the boxes according to split thresholds and sizes.
Kernel::ThreadScheduler *ts = new ThreadSchedulerFIFO();
ThreadPool tp(ts);
outWS->splitAllIfNeeded(ts);
tp.joinAll();
// Accumulate stats
totalAdded += numSinceSplit;
numSinceSplit = 0;
lastNumBoxes = obc->getTotalNumMDBoxes();
// Progress reporting
if (!fileBackedWS)
prog->report(i);
}
if (fileBackedWS) {
if (!(i % 10))
prog->report(i);
}
} // is box
} // for each box in the vector
prog->report();
outWS->splitAllIfNeeded(nullptr);
// Refresh all cache.
outWS->refreshCache();
// Account for events that were added after the last split
totalAdded += numSinceSplit;
g_log.notice() << totalAdded << " " << OMDE::getTypeName() << "s added to the output workspace.\n";
if (outWS->isFileBacked()) {
// Update the file-back-end
g_log.notice() << "Running SaveMD\n";
auto alg = createChildAlgorithm("SaveMD");
alg->setProperty("UpdateFileBackEnd", true);
alg->setProperty("InputWorkspace", outWS);
alg->executeAsChildAlg();
}
try {
outWS->copyExperimentInfos(*ws);
} catch (std::runtime_error &) {
g_log.warning() << this->name() << " was not able to copy experiment info to output workspace " << outWS->getName()
<< '\n';
}
// Pass on the display normalization from the input event workspace to the
// output event workspace
IMDEventWorkspace_sptr outEvent = std::dynamic_pointer_cast<IMDEventWorkspace>(outWS);
outEvent->setDisplayNormalization(ws->displayNormalization());
outEvent->setDisplayNormalizationHisto(ws->displayNormalizationHisto());
// return the size of the input workspace write buffer to its initial value
// bc->setCacheParameters(sizeof(MDE),writeBufSize);
this->setProperty("OutputWorkspace", outEvent);
}
//----------------------------------------------------------------------------------------------
/// Helper method
template <typename MDE, size_t nd> void SliceMD::doExec(typename MDEventWorkspace<MDE, nd>::sptr ws) {
if (m_outD == 0)
throw std::runtime_error("No output dimensions specified!");
// Templated method needs to call another templated method depending on the #
// of output dimensions.
if (MDE::getTypeName() == "MDLeanEvent") {
if (m_outD == 1)
this->slice<MDE, nd, MDLeanEvent<1>, 1>(ws);
else if (m_outD == 2)
this->slice<MDE, nd, MDLeanEvent<2>, 2>(ws);
else if (m_outD == 3)
this->slice<MDE, nd, MDLeanEvent<3>, 3>(ws);
else if (m_outD == 4)
this->slice<MDE, nd, MDLeanEvent<4>, 4>(ws);
else
throw std::runtime_error("Number of output dimensions > 4. This is not currently handled.");
} else if (MDE::getTypeName() == "MDEvent") {
if (m_outD == 1)
this->slice<MDE, nd, MDEvent<1>, 1>(ws);
else if (m_outD == 2)
this->slice<MDE, nd, MDEvent<2>, 2>(ws);
else if (m_outD == 3)
this->slice<MDE, nd, MDEvent<3>, 3>(ws);
else if (m_outD == 4)
this->slice<MDE, nd, MDEvent<4>, 4>(ws);
else
throw std::runtime_error("Number of output dimensions > 4. This is not currently handled.");
} else
throw std::runtime_error("Unexpected MDEvent type '" + MDE::getTypeName() + "'. This is not currently handled.");
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void SliceMD::exec() {
// Input MDEventWorkspace
m_inWS = getProperty("InputWorkspace");
// Run through the properties to create the transform you need
createTransform();
CALL_MDEVENT_FUNCTION(this->doExec, m_inWS);
}
} // namespace Mantid::MDAlgorithms