LookupTable.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
#pragma once
 
#include <functional>
#include <algorithm>
#include <vector>
#include <cmath>
#include <cassert>
 
namespace Belle2 {
  namespace TrackFindingCDC {
    template <class AResultType = double>
    std::vector<AResultType> linspace(double start,
                                      double final,
                                      std::size_t n,
                                      const std::function<AResultType(double)>& map)
    {
      std::vector<AResultType> result;
      result.reserve(n);
      result.push_back(map(start));
      for (std::size_t i = 1; i < n - 1; ++i) {
        result.push_back(map((start * (n - 1 - i) + final * i) / (n - 1)));
      }
      result.push_back(map(final));
      assert(result.size() == n);
      return result;
    }
 
    template <class AResultType = double>
    std::vector<AResultType> linspace(double start, double final, std::size_t n)
    {
      auto map = [](const double in) -> AResultType {return AResultType(in);};
      return linspace<AResultType>(start, final, n, map);
    }
 
    template<class T = double>
    class LookupTable {
 
    public:
      LookupTable(const std::function<T(double)>& map,
                  std::size_t nBins,
                  double lowerBound,
                  double upperBound)
        : m_lowerBound(lowerBound)
        , m_upperBound(upperBound)
        , m_binWidth((m_upperBound - m_lowerBound) / nBins)
        , m_values(linspace(lowerBound, upperBound, nBins + 1, map))
      {
        // Add a sentinel at the back.
        m_values.push_back(map(NAN));
      }
 
      T operator()(double x) const
      {
        if (not std::isfinite(x)) return m_values.back(); // Return sentinel
        const int iMax =  m_values.size() - 2; // Subtracting sentinel index
        double iBin = 0;
        double delta = std::modf((x - m_lowerBound) / m_binWidth, &iBin);
        int i = std::min(std::max(0, static_cast<int>(iBin)), iMax);
        return m_values[i] * (1 - delta) + m_values[i + 1] * delta;
      }
 
      const T& nearest(double x) const
      {
        if (not std::isfinite(x)) return m_values.back(); // Return sentinel value
        const int iMax =  m_values.size() - 2; // Subtracting sentinel index
        int iBin = std::round((x - m_lowerBound) / m_binWidth);
        int i = std::min(std::max(0, iBin), iMax);
        return m_values[i];
      }
 
      const T& at(int i) const
      {
        const int iMax =  m_values.size() - 2; // Subtracting sentinel index
        i = std::min(std::max(0, i), iMax);
        return m_values[i];
      }
 
      int getNBins() const
      {
        return m_values.size() - 2;
      }
 
      int getNPoints() const
      {
        return m_values.size() - 1;
      }
 
    private:
      double m_lowerBound;
 
      double m_upperBound;
 
      double m_binWidth;
 
      std::vector<T> m_values;
    };
  }
}