Rebin.cpp 15.4 KB
Newer Older
1
2
3
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4
5
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6
// SPDX - License - Identifier: GPL - 3.0 +
7
#include "MantidAlgorithms/Rebin.h"
8
#include "MantidHistogramData/Exception.h"
9
#include "MantidHistogramData/Rebin.h"
10

11
#include "MantidAPI/Axis.h"
12
#include "MantidAPI/HistoWorkspace.h"
13
#include "MantidDataObjects/EventList.h"
14
#include "MantidDataObjects/EventWorkspace.h"
LamarMoore's avatar
LamarMoore committed
15
#include "MantidDataObjects/Workspace2D.h"
16
#include "MantidDataObjects/WorkspaceCreation.h"
17
#include "MantidKernel/ArrayProperty.h"
18
#include "MantidKernel/BoundedValidator.h"
19
#include "MantidKernel/RebinParamsValidator.h"
20
#include "MantidKernel/VectorHelper.h"
21

22
namespace Mantid::Algorithms {
23
24
25
26
27
28
29
30
31

// Register the class into the algorithm factory
DECLARE_ALGORITHM(Rebin)

using namespace Kernel;
using namespace API;
using DataObjects::EventList;
using DataObjects::EventWorkspace;
using DataObjects::EventWorkspace_const_sptr;
LamarMoore's avatar
LamarMoore committed
32
33
34
35
36
using DataObjects::EventWorkspace_sptr;
using HistogramData::BinEdges;
using HistogramData::Frequencies;
using HistogramData::FrequencyStandardDeviations;
using HistogramData::Histogram;
37
using HistogramData::Exception::InvalidBinEdgesError;
38
39
40
41
42
43

//---------------------------------------------------------------------------------------------
// Public static methods
//---------------------------------------------------------------------------------------------

/**
LamarMoore's avatar
LamarMoore committed
44
45
46
47
48
49
 * Return the rebin parameters from a user input
 * @param inParams Input vector from user
 * @param inputWS Input workspace from user
 * @param logger A reference to a logger
 * @returns A new vector containing the rebin parameters
 */
Samuel Jones's avatar
Samuel Jones committed
50
51
std::vector<double> Rebin::rebinParamsFromInput(const std::vector<double> &inParams,
                                                const API::MatrixWorkspace &inputWS, Kernel::Logger &logger) {
52
  std::vector<double> rbParams;
53
  // The validator only passes parameters with size 1, or 3xn. No need to check again here
54
55
56
57
58
59
60
61
  if (inParams.size() >= 3) {
    // Input are min, delta, max
    rbParams = inParams;
  } else if (inParams.size() == 1) {
    double xmin = 0.;
    double xmax = 0.;
    inputWS.getXMinMax(xmin, xmax);

Samuel Jones's avatar
Samuel Jones committed
62
    logger.information() << "Using the current min and max as default " << xmin << ", " << xmax << '\n';
63
64
65
66
    rbParams.resize(3);
    rbParams[0] = xmin;
    rbParams[1] = inParams[0];
    rbParams[2] = xmax;
67
68
    if ((rbParams[1] < 0.) && (xmin < 0.) && (xmax > 0.)) {
      std::stringstream msg;
69
      msg << "Cannot create logarithmic binning that changes sign (xmin=" << xmin << ", xmax=" << xmax << ")";
70
71
      throw std::runtime_error(msg.str());
    }
72
73
74
75
76
77
78
79
  }
  return rbParams;
}

//---------------------------------------------------------------------------------------------
// Public methods
//---------------------------------------------------------------------------------------------

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
/// Validate that the input properties are sane.
std::map<std::string, std::string> Rebin::validateInputs() {
  std::map<std::string, std::string> helpMessages;
  if (existsProperty("Power") && !isDefault("Power")) {
    const double power = getProperty("Power");

    // attempt to roughly guess how many bins these parameters imply
    double roughEstimate = 0;

    if (!isDefault("Params")) {
      const std::vector<double> params = getProperty("Params");

      // Five significant places of the Euler-Mascheroni constant is probably more than enough for our needs
      double eulerMascheroni = 0.57721;

      // Params is check by the validator first, so we can assume it is in a correct format
      for (size_t i = 0; i < params.size() - 2; i += 2) {
        double upperLimit = params[i + 2];
        double lowerLimit = params[i];
        double factor = params[i + 1];

        if (factor <= 0) {
          helpMessages["Params"] = "Provided width value cannot be negative for inverse power binning.";
          return helpMessages;
        }

        if (power == 1) {
          roughEstimate += std::exp((upperLimit - lowerLimit) / factor - eulerMascheroni);
        } else {
          roughEstimate += std::pow(((upperLimit - lowerLimit) / factor) * (1 - power) + 1, 1 / (1 - power));
        }
      }
    }

    // Prevent the user form creating too many bins
    if (roughEstimate > 10000) {
116
      helpMessages["Power"] = "This binning is expected to give more than 10000 bins.";
117
118
119
120
121
    }
  }
  return helpMessages;
}

122
/** Initialisation method. Declares properties to be used in algorithm.
LamarMoore's avatar
LamarMoore committed
123
124
 *
 */
125
void Rebin::init() {
Samuel Jones's avatar
Samuel Jones committed
126
  declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
Sam Jenkins's avatar
Sam Jenkins committed
127
                  "Workspace containing the input data");
Samuel Jones's avatar
Samuel Jones committed
128
  declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
129
                  "The name to give the output workspace");
130

Samuel Jones's avatar
Samuel Jones committed
131
132
  declareProperty(std::make_unique<ArrayProperty<double>>("Params", std::make_shared<RebinParamsValidator>()),
                  "A comma separated list of first bin boundary, width, last bin boundary. "
Mathieu Tillet's avatar
Mathieu Tillet committed
133
                  "Optionally this can be followed by a comma and more widths and last boundary pairs. "
134
135
                  "Optionally this can also be a single number, which is the bin width. In this case, the boundary of "
                  "binning will be determined by minimum and maximum TOF values among all events, or previous binning "
Mathieu Tillet's avatar
Mathieu Tillet committed
136
                  "boundary, in case of event Workspace, or non-event Workspace, respectively. "
137
                  "Negative width values indicate logarithmic binning.");
Samuel Jones's avatar
Samuel Jones committed
138
139

  declareProperty("PreserveEvents", true,
140
141
142
                  "Keep the output workspace as an EventWorkspace, if the input has events. If the input and output "
                  "EventWorkspace names are the same, only the X bins are set, which is very quick. If false, then the "
                  "workspace gets converted to a Workspace2D histogram.");
Samuel Jones's avatar
Samuel Jones committed
143

144
  declareProperty("FullBinsOnly", false, "Omit the final bin if its width is smaller than the step size");
145

146
  declareProperty("IgnoreBinErrors", false,
147
148
149
                  "Ignore errors related to zero/negative bin widths in input/output workspaces. When ignored, the "
                  "signal and errors are set to zero");

150
151
152
153
154
155
156
157
158
159
160
161
  declareProperty(
      "UseReverseLogarithmic", false,
      "For logarithmic intervals, the splitting starts from the end and goes back to the start, ie the bins are bigger "
      "at the start getting exponentially smaller until they reach the end. For these bins, the FullBinsOnly flag is "
      "ignored.");

  auto powerValidator = std::make_shared<Mantid::Kernel::BoundedValidator<double>>();
  powerValidator->setLower(0);
  powerValidator->setUpper(1);
  declareProperty("Power", 0., powerValidator,
                  "Splits the interval in bins which actual width is equal to requested width / (i ^ power); default "
                  "is linear. Power must be between 0 and 1.");
162
163
164
}

/** Executes the rebin algorithm
LamarMoore's avatar
LamarMoore committed
165
166
167
168
 *
 *  @throw runtime_error Thrown if the bin range does not intersect the range of
 *the input workspace
 */
169
170
171
172
173
174
175
176
177
178
179
void Rebin::exec() {
  // Get the input workspace
  MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
  MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");

  // Are we preserving event workspace-iness?
  bool PreserveEvents = getProperty("PreserveEvents");

  // Rebinning in-place
  bool inPlace = (inputWS == outputWS);

Samuel Jones's avatar
Samuel Jones committed
180
  std::vector<double> rbParams = rebinParamsFromInput(getProperty("Params"), *inputWS, g_log);
181
182
183
184
185

  const bool dist = inputWS->isDistribution();
  const bool isHist = inputWS->isHistogramData();

  // workspace independent determination of length
186
  const auto histnumber = static_cast<int>(inputWS->getNumberHistograms());
187
188

  bool fullBinsOnly = getProperty("FullBinsOnly");
189
  bool useReverseLog = getProperty("UseReverseLogarithmic");
190
  double power = getProperty("Power");
191
192
193
194

  double xmin = 0.;
  double xmax = 0.;
  inputWS->getXMinMax(xmin, xmax);
195

196
  HistogramData::BinEdges XValues_new(0);
197
  // create new output X axis
198
  static_cast<void>(VectorHelper::createAxisFromRebinParams(rbParams, XValues_new.mutableRawData(), true, fullBinsOnly,
199
                                                            xmin, xmax, useReverseLog, power));
200
201

  // Now, determine if the input workspace is actually an EventWorkspace
Samuel Jones's avatar
Samuel Jones committed
202
  EventWorkspace_const_sptr eventInputWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
203

204
  if (eventInputWS != nullptr) {
205
206
    //------- EventWorkspace as input -------------------------------------

207
208
    if (PreserveEvents) {
      if (!inPlace) {
209
        outputWS = inputWS->clone();
210
      }
211
      auto eventOutputWS = std::dynamic_pointer_cast<EventWorkspace>(outputWS);
212
213
214
215
216
      // This only sets the X axis. Actual rebinning will be done upon data
      // access.
      eventOutputWS->setAllX(XValues_new);
    } else {
      //--------- Different output, OR you're inplace but not preserving Events
Samuel Jones's avatar
Samuel Jones committed
217
218
      g_log.information() << "Creating a Workspace2D from the EventWorkspace " << eventInputWS->getName() << ".\n";
      outputWS = DataObjects::create<DataObjects::Workspace2D>(*inputWS, histnumber, XValues_new);
219
220
221
222
223

      // Initialize progress reporting.
      Progress prog(this, 0.0, 1.0, histnumber);

      // Go through all the histograms and set the data
224
      PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
225
226
227
      for (int i = 0; i < histnumber; ++i) {
        PARALLEL_START_INTERUPT_REGION
        // Get a const event list reference. eventInputWS->dataY() doesn't work.
228
        const EventList &el = eventInputWS->getSpectrum(i);
229
230
        MantidVec y_data, e_data;
        // The EventList takes care of histogramming.
231
        el.generateHistogram(XValues_new.rawData(), y_data, e_data);
232
233

        // Copy the data over.
234
235
        outputWS->mutableY(i) = y_data;
        outputWS->mutableE(i) = e_data;
236
237
238
239

        // Report progress
        prog.report(name());
        PARALLEL_END_INTERUPT_REGION
240
      }
241
      PARALLEL_CHECK_INTERUPT_REGION
242

243
244
      // Copy all the axes
      for (int i = 1; i < inputWS->axes(); i++) {
Samuel Jones's avatar
Samuel Jones committed
245
        outputWS->replaceAxis(i, std::unique_ptr<Axis>(inputWS->getAxis(i)->clone(outputWS.get())));
246
247
        outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit();
      }
248

249
250
251
252
253
254
      // Copy the units over too.
      for (int i = 0; i < outputWS->axes(); ++i)
        outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit();
      outputWS->setYUnit(eventInputWS->YUnit());
      outputWS->setYUnitLabel(eventInputWS->YUnitLabel());
    }
Nick Draper's avatar
re #100    
Nick Draper committed
255

256
257
258
    // Assign it to the output workspace property
    setProperty("OutputWorkspace", outputWS);

259
  } // END ---- EventWorkspace
260

261
  else
262

263
  { //------- Workspace2D or other MatrixWorkspace ---------------------------
264

265
    if (!isHist) {
Hahn, Steven's avatar
Hahn, Steven committed
266
      g_log.information() << "Rebin: Converting Data to Histogram.\n";
Samuel Jones's avatar
Samuel Jones committed
267
      Mantid::API::Algorithm_sptr ChildAlg = createChildAlgorithm("ConvertToHistogram");
268
269
270
271
272
      ChildAlg->initialize();
      ChildAlg->setProperty("InputWorkspace", inputWS);
      ChildAlg->execute();
      inputWS = ChildAlg->getProperty("OutputWorkspace");
    }
273

274
275
    // make output Workspace the same type is the input, but with new length of
    // signal array
Samuel Jones's avatar
Samuel Jones committed
276
    outputWS = DataObjects::create<API::HistoWorkspace>(*inputWS, histnumber, XValues_new);
277
278
279

    // Copy over the 'vertical' axis
    if (inputWS->axes() > 1)
Samuel Jones's avatar
Samuel Jones committed
280
      outputWS->replaceAxis(1, std::unique_ptr<Axis>(inputWS->getAxis(1)->clone(outputWS.get())));
281
    bool ignoreBinErrors = getProperty("IgnoreBinErrors");
282
283

    Progress prog(this, 0.0, 1.0, histnumber);
284
    PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
285
286
    for (int hist = 0; hist < histnumber; ++hist) {
      PARALLEL_START_INTERUPT_REGION
287
288

      try {
Samuel Jones's avatar
Samuel Jones committed
289
        outputWS->setHistogram(hist, HistogramData::rebin(inputWS->histogram(hist), XValues_new));
290
291
292
293
294
295
      } catch (InvalidBinEdgesError &) {
        if (ignoreBinErrors)
          outputWS->setBinEdges(hist, XValues_new);
        else
          throw;
      }
296
297
      prog.report(name());
      PARALLEL_END_INTERUPT_REGION
298
    }
299
    PARALLEL_CHECK_INTERUPT_REGION
300
    outputWS->setDistribution(dist);
301
302
303
304
305

    // Now propagate any masking correctly to the output workspace
    // More efficient to have this in a separate loop because
    // MatrixWorkspace::maskBins blocks multi-threading
    for (int i = 0; i < histnumber; ++i) {
Samuel Jones's avatar
Samuel Jones committed
306
      if (inputWS->hasMaskedBins(i)) // Does the current spectrum have any masked bins?
307
      {
308
        this->propagateMasks(inputWS, outputWS, i);
309
      }
310
311
312
313
314
    }
    // Copy the units over too.
    for (int i = 0; i < outputWS->axes(); ++i) {
      outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit();
    }
315

316
    if (!isHist) {
Hahn, Steven's avatar
Hahn, Steven committed
317
      g_log.information() << "Rebin: Converting Data back to Data Points.\n";
Samuel Jones's avatar
Samuel Jones committed
318
      Mantid::API::Algorithm_sptr ChildAlg = createChildAlgorithm("ConvertToPointData");
319
320
321
322
      ChildAlg->initialize();
      ChildAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", outputWS);
      ChildAlg->execute();
      outputWS = ChildAlg->getProperty("OutputWorkspace");
323
    }
324

325
326
327
328
329
330
331
    // Assign it to the output workspace property
    setProperty("OutputWorkspace", outputWS);

  } // END ---- Workspace2D
}

/** Takes the masks in the input workspace and apportions the weights into the
LamarMoore's avatar
LamarMoore committed
332
333
334
335
336
337
338
 *new bins that overlap
 *  with a masked bin. These bins are then masked with the calculated weight.
 *
 *  @param inputWS ::  The input workspace
 *  @param outputWS :: The output workspace
 *  @param hist ::    The index of the current histogram
 */
Samuel Jones's avatar
Samuel Jones committed
339
void Rebin::propagateMasks(const API::MatrixWorkspace_const_sptr &inputWS, const API::MatrixWorkspace_sptr &outputWS,
David Fairbrother's avatar
David Fairbrother committed
340
                           int hist) {
341
342
343
344
345
346
347
348
349
  // Not too happy with the efficiency of this way of doing it, but it's a lot
  // simpler to use the
  // existing rebin algorithm to distribute the weights than to re-implement it
  // for this

  MantidVec masked_bins, weights;
  // Get a reference to the list of masked bins for this spectrum
  const MatrixWorkspace::MaskList &mask = inputWS->maskedBins(hist);
  // Now iterate over the list, building up a vector of the masked bins
350
  auto it = mask.cbegin();
351
  auto &XValues = inputWS->x(hist);
352
353
354
  masked_bins.emplace_back(XValues[(*it).first]);
  weights.emplace_back((*it).second);
  masked_bins.emplace_back(XValues[(*it).first + 1]);
355
356
357
358
359
  for (++it; it != mask.end(); ++it) {
    const double currentX = XValues[(*it).first];
    // Add an intermediate bin with zero weight if masked bins aren't
    // consecutive
    if (masked_bins.back() != currentX) {
360
361
      weights.emplace_back(0.0);
      masked_bins.emplace_back(currentX);
362
    }
363
364
    weights.emplace_back((*it).second);
    masked_bins.emplace_back(XValues[(*it).first + 1]);
365
366
  }

Moore's avatar
Moore committed
367
368
  //// Create a zero vector for the errors because we don't care about them here
  auto errSize = weights.size();
Samuel Jones's avatar
Samuel Jones committed
369
  Histogram oldHist(BinEdges(std::move(masked_bins)), Frequencies(std::move(weights)),
Moore's avatar
Moore committed
370
                    FrequencyStandardDeviations(errSize, 0));
371
372
  // Use rebin function to redistribute the weights. Note that distribution flag
  // is set
373
  bool ignoreErrors = getProperty("IgnoreBinErrors");
374

375
376
  try {
    auto newHist = HistogramData::rebin(oldHist, outputWS->binEdges(hist));
Moore's avatar
Moore committed
377
    auto &newWeights = newHist.y();
378
379
380
381
382
383
384
385
386

    // Now process the output vector and fill the new masking list
    for (size_t index = 0; index < newWeights.size(); ++index) {
      if (newWeights[index] > 0.0)
        outputWS->flagMasked(hist, index, newWeights[index]);
    }
  } catch (InvalidBinEdgesError &) {
    if (!ignoreErrors)
      throw;
387
388
389
  }
}

390
} // namespace Mantid::Algorithms