BinMD.cpp 16 KB
Newer Older
1
#include "MantidAPI/ImplicitFunctionFactory.h"
2
#include "MantidGeometry/MDGeometry/MDBoxImplicitFunction.h"
3
#include "MantidGeometry/MDGeometry/MDHistoDimension.h"
4
#include "MantidKernel/CPUTimer.h"
5
6
#include "MantidKernel/Strings.h"
#include "MantidKernel/System.h"
7
#include "MantidKernel/Utils.h"
8
9
10
11
12
13
14
#include "MantidDataObjects/CoordTransformAffineParser.h"
#include "MantidDataObjects/CoordTransformAligned.h"
#include "MantidDataObjects/MDBoxBase.h"
#include "MantidDataObjects/MDBox.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidDataObjects/MDEventWorkspace.h"
#include "MantidDataObjects/MDHistoWorkspace.h"
15
#include "MantidMDAlgorithms/BinMD.h"
16
#include <boost/algorithm/string.hpp>
17
#include "MantidKernel/EnabledWhenProperty.h"
18
#include "MantidDataObjects/CoordTransformAffine.h"
19
20

using Mantid::Kernel::CPUTimer;
21
using Mantid::Kernel::EnabledWhenProperty;
22

23
24
25
26
27
28
29
30
31
namespace Mantid {
namespace MDAlgorithms {

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(BinMD)

using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
32
using namespace Mantid::DataObjects;
33
34
35
36

//----------------------------------------------------------------------------------------------
/** Constructor
 */
37
BinMD::BinMD()
Hahn, Steven's avatar
Hahn, Steven committed
38
39
40
    : outWS(), prog(nullptr), implicitFunction(nullptr),
      indexMultiplier(nullptr), signals(nullptr), errors(nullptr),
      numEvents(nullptr) {}
41
42
43
44
45

//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
 */
void BinMD::init() {
46
47
  declareProperty(make_unique<WorkspaceProperty<IMDWorkspace>>(
                      "InputWorkspace", "", Direction::Input),
48
49
50
51
52
53
54
55
                  "An input MDWorkspace.");

  // Properties for specifying the slice to perform.
  this->initSlicingProps();

  // --------------- Processing methods and options
  // ---------------------------------------
  std::string grp = "Methods";
56
57
  declareProperty(make_unique<PropertyWithValue<std::string>>(
                      "ImplicitFunctionXML", "", Direction::Input),
58
59
60
61
62
                  "XML string describing the implicit function determining "
                  "which bins to use.");
  setPropertyGroup("ImplicitFunctionXML", grp);

  declareProperty(
63
64
      make_unique<PropertyWithValue<bool>>("IterateEvents", true,
                                           Direction::Input),
65
66
67
68
69
70
71
      "Alternative binning method where you iterate through every event, "
      "placing them in the proper bin.\n"
      "This may be faster for workspaces with few events and lots of output "
      "bins.");
  setPropertyGroup("IterateEvents", grp);

  declareProperty(
72
      make_unique<PropertyWithValue<bool>>("Parallel", false, Direction::Input),
73
74
75
76
77
      "Temporary parameter: true to run in parallel. This is ignored for "
      "file-backed workspaces, where running in parallel makes things slower "
      "due to disk thrashing.");
  setPropertyGroup("Parallel", grp);

78
79
  declareProperty(make_unique<WorkspaceProperty<Workspace>>(
                      "OutputWorkspace", "", Direction::Output),
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
                  "A name for the output MDHistoWorkspace.");
}

//----------------------------------------------------------------------------------------------
/** Bin the contents of a MDBox
 *
 * @param box :: pointer to the MDBox to bin
 * @param chunkMin :: the minimum index in each dimension to consider "valid"
 *(inclusive)
 * @param chunkMax :: the maximum index in each dimension to consider "valid"
 *(exclusive)
 */
template <typename MDE, size_t nd>
inline void BinMD::binMDBox(MDBox<MDE, nd> *box, const size_t *const chunkMin,
                            const size_t *const chunkMax) {
  // An array to hold the rotated/transformed coordinates
96
  auto outCenter = new coord_t[m_outD];
97
98
99
100
101
102
103
104
105
106
107
108
109

  // Evaluate whether the entire box is in the same bin
  if (box->getNPoints() > (1 << nd) * 2) {
    // There is a check that the number of events is enough for it to make sense
    // to do all this processing.
    size_t numVertexes = 0;
    coord_t *vertexes = box->getVertexesArray(numVertexes);

    // All vertexes have to be within THE SAME BIN = have the same linear index.
    size_t lastLinearIndex = 0;
    bool badOne = false;

    for (size_t i = 0; i < numVertexes; i++) {
110
      // Cache the center of the event (again for speed)
111
      const coord_t *inCenter = vertexes + i * nd;
112
113
114

      // Now transform to the output dimensions
      m_transform->apply(inCenter, outCenter);
115
116
      // std::cout << "Input coord " << VMD(nd,inCenter) << " transformed to "
      // <<  VMD(nd,outCenter) << std::endl;
117
118
119

      // To build up the linear index
      size_t linearIndex = 0;
120
121
      // To mark VERTEXES outside range
      badOne = false;
122
123

      /// Loop through the dimensions on which we bin
124
      for (size_t bd = 0; bd < m_outD; bd++) {
125
126
127
        // What is the bin index in that dimension
        coord_t x = outCenter[bd];
        size_t ix = size_t(x);
128
        // Within range (for this chunk)?
129
        if ((x >= 0) && (ix >= chunkMin[bd]) && (ix < chunkMax[bd])) {
130
          // Build up the linear index
131
          linearIndex += indexMultiplier[bd] * ix;
132
        } else {
133
          // Outside the range
134
          badOne = true;
135
136
137
138
          break;
        }
      } // (for each dim in MDHisto)

139
140
141
142
143
144
145
146
      // Is the vertex at the same place as the last one?
      if (!badOne) {
        if ((i > 0) && (linearIndex != lastLinearIndex)) {
          // Change of index
          badOne = true;
          break;
        }
        lastLinearIndex = linearIndex;
147
      }
148

149
150
151
152
153
154
155
156
157
158
159
160
161
162
      // Was the vertex completely outside the range?
      if (badOne)
        break;
    } // (for each vertex)

    delete[] vertexes;

    if (!badOne) {
      // Yes, the entire box is within a single bin
      //        std::cout << "Box at " << box->getExtentsStr() << " is within a
      //        single bin.\n";
      // Add the CACHED signal from the entire box
      signals[lastLinearIndex] += box->getSignal();
      errors[lastLinearIndex] += box->getErrorSquared();
163
      // TODO: If DataObjects get a weight, this would need to get the summed
164
165
166
167
168
169
170
171
      // weight.
      numEvents[lastLinearIndex] += static_cast<signal_t>(box->getNPoints());

      // And don't bother looking at each event. This may save lots of time
      // loading from disk.
      delete[] outCenter;
      return;
    }
172
173
  }

174
175
176
177
  // If you get here, you could not determine that the entire box was in the
  // same bin.
  // So you need to iterate through events.
  const std::vector<MDE> &events = box->getConstEvents();
Hahn, Steven's avatar
Hahn, Steven committed
178
  for (auto it = events.begin(); it != events.end(); ++it) {
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
    // Cache the center of the event (again for speed)
    const coord_t *inCenter = it->getCenter();

    // Now transform to the output dimensions
    m_transform->apply(inCenter, outCenter);

    // To build up the linear index
    size_t linearIndex = 0;
    // To mark events outside range
    bool badOne = false;

    /// Loop through the dimensions on which we bin
    for (size_t bd = 0; bd < m_outD; bd++) {
      // What is the bin index in that dimension
      coord_t x = outCenter[bd];
      size_t ix = size_t(x);
      // Within range (for this chunk)?
      if ((x >= 0) && (ix >= chunkMin[bd]) && (ix < chunkMax[bd])) {
        // Build up the linear index
        linearIndex += indexMultiplier[bd] * ix;
      } else {
        // Outside the range
        badOne = true;
        break;
      }
    } // (for each dim in MDHisto)

    if (!badOne) {
      // Sum the signals as doubles to preserve precision
      signals[linearIndex] += static_cast<signal_t>(it->getSignal());
      errors[linearIndex] += static_cast<signal_t>(it->getErrorSquared());
210
      // TODO: If DataObjects get a weight, this would need to get the summed
211
212
      // weight.
      numEvents[linearIndex] += 1.0;
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
  }
  // Done with the events list
  box->releaseEvents();

  delete[] outCenter;
}

//----------------------------------------------------------------------------------------------
/** Perform binning by iterating through every event and placing them in the
 *output workspace
 *
 * @param ws :: MDEventWorkspace of the given type.
 */
template <typename MDE, size_t nd>
void BinMD::binByIterating(typename MDEventWorkspace<MDE, nd>::sptr ws) {
  BoxController_sptr bc = ws->getBoxController();
  // store exisiting write buffer size for the future
  // uint64_t writeBufSize =bc->getDiskBuffer().getWriteBufferSize();
  // and disable write buffer (if any) for input MD Events for this algorithm
  // purposes;
  // bc->setCacheParameters(1,0);

  // Cache some data to speed up accessing them a bit
  indexMultiplier = new size_t[m_outD];
  for (size_t d = 0; d < m_outD; d++) {
    if (d > 0)
      indexMultiplier[d] = outWS->getIndexMultiplier()[d - 1];
    else
      indexMultiplier[d] = 1;
  }
  signals = outWS->getSignalArray();
  errors = outWS->getErrorSquaredArray();
  numEvents = outWS->getNumEventsArray();

  // Start with signal/error/numEvents at 0.0
  outWS->setTo(0.0, 0.0, 0.0);

  // The dimension (in the output workspace) along which we chunk for parallel
  // processing
  // TODO: Find the smartest dimension to chunk against
  size_t chunkDimension = 0;

  // How many bins (in that dimension) per chunk.
  // Try to split it so each core will get 2 tasks:
  int chunkNumBins = int(m_binDimensions[chunkDimension]->getNBins() /
                         (PARALLEL_GET_MAX_THREADS * 2));
  if (chunkNumBins < 1)
    chunkNumBins = 1;

  // Do we actually do it in parallel?
  bool doParallel = getProperty("Parallel");
  // Not if file-backed!
  if (bc->isFileBacked())
    doParallel = false;
  if (!doParallel)
    chunkNumBins = int(m_binDimensions[chunkDimension]->getNBins());

  // Total number of steps
  size_t progNumSteps = 0;
  if (prog)
    prog->setNotifyStep(0.1);
  if (prog)
    prog->resetNumSteps(100, 0.00, 1.0);

  // Run the chunks in parallel. There is no overlap in the output workspace so
279
  // it is thread safe to write to it..
280
  // cppcheck-suppress syntaxError
Janik Zikovsky's avatar
Janik Zikovsky committed
281
    PRAGMA_OMP( parallel for schedule(dynamic,1) if (doParallel) )
282
283
284
    for (int chunk = 0;
         chunk < int(m_binDimensions[chunkDimension]->getNBins());
         chunk += chunkNumBins) {
285
      PARALLEL_START_INTERUPT_REGION
286
      // Region of interest for this chunk.
287
288
      std::vector<size_t> chunkMin(m_outD);
      std::vector<size_t> chunkMax(m_outD);
289
      for (size_t bd = 0; bd < m_outD; bd++) {
290
        // Same limits in the other dimensions
291
        chunkMin[bd] = 0;
292
        chunkMax[bd] = m_binDimensions[bd]->getNBins();
293
294
      }
      // Parcel out a chunk in that single dimension dimension
295
      chunkMin[chunkDimension] = size_t(chunk);
296
297
      if (size_t(chunk + chunkNumBins) >
          m_binDimensions[chunkDimension]->getNBins())
298
        chunkMax[chunkDimension] = m_binDimensions[chunkDimension]->getNBins();
299
      else
300
        chunkMax[chunkDimension] = size_t(chunk + chunkNumBins);
301

302
303
304
305
      // Build an implicit function (it needs to be in the space of the
      // MDEventWorkspace)
      MDImplicitFunction *function =
          this->getImplicitFunctionForChunk(chunkMin.data(), chunkMax.data());
306

307
      // Use getBoxes() to get an array with a pointer to each box
308
      std::vector<API::IMDNode *> boxes;
309
310
      // Leaf-only; no depth limit; with the implicit function passed to it.
      ws->getBox()->getBoxes(boxes, 1000, true, function);
311

312
313
      // Sort boxes by file position IF file backed. This reduces seeking time,
      // hopefully.
314
      if (bc->isFileBacked())
315
        API::IMDNode::sortObjByID(boxes);
316

317
      // For progress reporting, the # of boxes
318
319
320
321
      if (prog) {
        PARALLEL_CRITICAL(BinMD_progress) {
          g_log.debug() << "Chunk " << chunk << ": found " << boxes.size()
                        << " boxes within the implicit function." << std::endl;
322
          progNumSteps += boxes.size();
323
          prog->setNumSteps(progNumSteps);
324
325
        }
      }
326

327
      // Go through every box for this chunk.
328
329
      for (auto &boxe : boxes) {
        MDBox<MDE, nd> *box = dynamic_cast<MDBox<MDE, nd> *>(boxe);
330
        // Perform the binning in this separate method.
331
        if (box && !box->getIsMasked())
332
          this->binMDBox(box, chunkMin.data(), chunkMax.data());
333
334

        // Progress reporting
335
336
        if (prog)
          prog->report();
337
338
339
        // For early cancelling of the loop
        if (this->m_cancel)
          break;
340
      } // for each box in the vector
341
      PARALLEL_END_INTERUPT_REGION
342
    } // for each chunk in parallel
343
    PARALLEL_CHECK_INTERUPT_REGION
344

345
    // Now the implicit function
346
    if (implicitFunction) {
347
348
      if (prog)
        prog->report("Applying implicit function.");
349
350
351
      signal_t nan = std::numeric_limits<signal_t>::quiet_NaN();
      outWS->applyImplicitFunction(implicitFunction, nan, nan);
    }
Alex Buts's avatar
Alex Buts committed
352
353

    // return the size of the input workspace write buffer to its initial value
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
    // bc->setCacheParameters(sizeof(MDE),writeBufSize);
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
 */
void BinMD::exec() {
  // Input MDEventWorkspace/MDHistoWorkspace
  m_inWS = getProperty("InputWorkspace");
  // Look at properties, create either axis-aligned or general transform.
  // This (can) change m_inWS
  this->createTransform();

  // De serialize the implicit function
  std::string ImplicitFunctionXML = getPropertyValue("ImplicitFunctionXML");
369
  implicitFunction = nullptr;
370
371
372
373
374
  if (!ImplicitFunctionXML.empty())
    implicitFunction =
        Mantid::API::ImplicitFunctionFactory::Instance().createUnwrapped(
            ImplicitFunctionXML);

375
376
  // This gets deleted by the thread pool; don't delete it in here.
  prog = new Progress(this, 0, 1.0, 1);
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393

  // Create the dense histogram. This allocates the memory
  outWS = MDHistoWorkspace_sptr(new MDHistoWorkspace(m_binDimensions));

  // Saves the geometry transformation from original to binned in the workspace
  outWS->setTransformFromOriginal(this->m_transformFromOriginal, 0);
  outWS->setTransformToOriginal(this->m_transformToOriginal, 0);
  for (size_t i = 0; i < m_bases.size(); i++)
    outWS->setBasisVector(i, m_bases[i]);
  outWS->setOrigin(this->m_translation);
  outWS->setOriginalWorkspace(m_inWS, 0);

  // And the intermediate WS one too, if any
  if (m_intermediateWS) {
    outWS->setOriginalWorkspace(m_intermediateWS, 1);
    outWS->setTransformFromOriginal(m_transformFromIntermediate, 1);
    outWS->setTransformToOriginal(m_transformToIntermediate, 1);
394
  }
395

396
397
398
399
400
401
  // Wrapper to cast to MDEventWorkspace then call the function
  bool IterateEvents = getProperty("IterateEvents");
  if (!IterateEvents) {
    g_log.warning() << "IterateEvents=False is no longer supported. Setting "
                       "IterateEvents=True." << std::endl;
    IterateEvents = true;
402
403
  }

404
  /*
405
406
407
408
  We should fail noisily here. CALL_MDEVENT_FUNCTION will silently allow
  IMDHistoWorkspaces to cascade through to the end
  and result in an empty output. The only way we allow InputWorkspaces to be
  IMDHistoWorkspaces is if they also happen to contain original workspaces
409
410
  that are MDEventWorkspaces.
  */
411
412
413
414
415
  if (boost::dynamic_pointer_cast<IMDHistoWorkspace>(m_inWS)) {
    throw std::runtime_error(
        "Cannot rebin a workspace that is histogrammed and has no original "
        "workspace that is an MDEventWorkspace. "
        "Reprocess the input so that it contains full MDEvents.");
416
417
  }

418
419
  CALL_MDEVENT_FUNCTION(this->binByIterating, m_inWS);

420
  // Copy the coordinate system & experiment infos to the output
421
422
  IMDEventWorkspace_sptr inEWS =
      boost::dynamic_pointer_cast<IMDEventWorkspace>(m_inWS);
423
424
  if (inEWS) {
    outWS->setCoordinateSystem(inEWS->getSpecialCoordinateSystem());
425
426
427
428
429
430
431
432
    try {
      outWS->copyExperimentInfos(*inEWS);
    } catch (std::runtime_error &) {
      g_log.warning()
          << this->name()
          << " was not able to copy experiment info to output workspace "
          << outWS->getName() << std::endl;
    }
433
  }
434

435
  // Pass on the display normalization from the input workspace
436
437
  outWS->setDisplayNormalization(m_inWS->displayNormalizationHisto());

438
439
440
441
  outWS->updateSum();
  // Save the output
  setProperty("OutputWorkspace", boost::dynamic_pointer_cast<Workspace>(outWS));
}
442
443

} // namespace Mantid
444
} // namespace DataObjects