BinFinder.cpp

Go to the documentation of this file.

// 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 "MantidKernel/BinFinder.h"
#include <cmath>
#include <cstddef>
#include <stdexcept>
 
using std::size_t;
 
namespace Mantid::Kernel {
 
BinFinder::BinFinder(const std::vector<double> &binParams) {
  boundaries.clear();
  stepSizes.clear();
 
  size_t n = binParams.size();
  if (n < 3)
    throw std::invalid_argument("BinFinder: not enough bin parameters.");
  if (n % 2 == 0)
    throw std::invalid_argument("BinFinder: the number of bin parameters should be odd.");
 
  for (size_t i = 0; i < n / 2; i++) {
    // The boundaries
    double min = binParams[i * 2];
    double max = binParams[i * 2 + 2];
    // Only the first bin needs the min boundary
    if (i == 0)
      boundaries.emplace_back(min);
    boundaries.emplace_back(max);
    // The step
    double step = binParams[i * 2 + 1];
    stepSizes.emplace_back(step);
    if (step == 0)
      throw std::invalid_argument("BinFinder: step size of 0.");
    if ((step < 0) && (min <= 0))
      throw std::invalid_argument("BinFinder: logarithmic binning with 0.0 starting bin.");
    if (max <= min)
      throw std::invalid_argument("BinFinder: final bin must be > starting bin boundary.");
 
    int numBins = 0;
 
    // Pre-do some calculations for log binning.
    if (step < 0) {
      double log_step = log1p(fabs(step));
      logSteps.emplace_back(log_step);
      if (i == 0)
        logBoundaries.emplace_back(log(min));
      logBoundaries.emplace_back(log(max));
      // How many bins is that?
      numBins = static_cast<int>(ceil((log(max) - log(min)) / log_step));
 
      // Check that the last bin is at least .25 x the previous step
      // This is because the VectorHelper removes that final bin. Annoying!
      double nextToLastValue = min * pow(1.0 + fabs(step), numBins - 1);
      double nextToNextToLastValue = min * pow(1.0 + fabs(step), numBins - 2);
      double lastBinSize = max - nextToLastValue;
      double nextToLastBinSize = nextToLastValue - nextToNextToLastValue;
      if (lastBinSize < nextToLastBinSize * 0.25)
        numBins--;
      if (numBins < 1)
        numBins = 1;
 
    } else {
      // Empty log values; these won't be used
      logSteps.emplace_back(0);
      if (i == 0)
        logBoundaries.emplace_back(0);
      logBoundaries.emplace_back(0);
      //# of linear bins
      numBins = static_cast<int>(ceil((max - min) / step));
 
      // Check that the last bin is at least .25 x the previous step
      // This is because the VectorHelper removes that final bin. Annoying!
      double lastBinSize = max - ((numBins - 1) * step + min);
      if (lastBinSize < step * 0.25)
        numBins--;
      if (numBins < 1)
        numBins = 1;
    }
 
    // Find the end bin index
    int startBinIndex = 0;
    if (i > 0)
      startBinIndex = this->endBinIndex[i - 1];
    endBinIndex.emplace_back(numBins + startBinIndex);
  }
  // How many binning regions?
  numRegions = static_cast<int>(stepSizes.size());
}
 
int BinFinder::lastBinIndex() {
  if (!endBinIndex.empty())
    return endBinIndex.back();
  else
    return -1;
}
 
int BinFinder::bin(double x) {
  int index;
  double min(boundaries[0]);
 
  // Too small?
  if (x < boundaries[0])
    return -1;
 
  // Find which binning region to use
  int i = -1;
  for (i = 0; i < numRegions; i++) {
    min = boundaries[i];
    double max = boundaries[i + 1];
    if ((x >= min) && (x < max))
      break;
  }
  // Didn't find it?
  if (i >= numRegions)
    return -1;
 
  // Step size in this region
  double step = stepSizes[i];
  if (step > 0) {
    // Linear binning. Truncate when you divide by the step size
    index = static_cast<int>((x - min) / step);
    // Add bin index offset if not in the first region
    if (i > 0)
      index += endBinIndex[i - 1];
    // In the event that a final bin was skipped, cap to the max
    if (index >= endBinIndex[i])
      index = endBinIndex[i] - 1;
    return index;
  } else {
    double log_x = log(x); // Just one log to call per event!
    double log_step = logSteps[i];
    double log_min = logBoundaries[i];
    index = static_cast<int>((log_x - log_min) / log_step);
    // Add bin index offset if not in the first region
    if (i > 0)
      index += endBinIndex[i - 1];
    // In the event that a final bin was skipped, cap to the max
    if (index >= endBinIndex[i])
      index = endBinIndex[i] - 1;
    return index;
  }
}
} // namespace Mantid::Kernel