BinnedAxis.h

Go to the documentation of this file.

// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2023 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_BINNEDAXIS_H
#define YODA_BINNEDAXIS_H
 
#include "YODA/Utils/BinEstimators.h"
#include "YODA/Utils/BinningUtils.h"
#include "YODA/Utils/MathUtils.h"
#include "YODA/Utils/MetaUtils.h"
#include <cmath>
#include <cstdlib>
#include <limits>
#include <memory>
#include <algorithm>
#include <set>
#include <string>
#include <stdexcept>
#include <type_traits>
#include <vector>
#include <iostream>
#include <iomanip>
 
namespace YODA {
 
  const size_t SEARCH_SIZElc = 16;
  const size_t BISECT_LINEAR_THRESHOLDlc = 32;
 
  namespace YODAConcepts {
 
      using MetaUtils::conjunction;
 
      template <typename T>
      struct AxisImpl {
 
 
        using index_sig        = size_t (T::*)(const typename T::EdgeT&) const;
        using edge_sig         = typename T::EdgeT (T::*)(const size_t) const;
        using edges_sig        = const std::vector<std::reference_wrapper<const typename T::EdgeT>> (T::*)() const noexcept;
        using size_sig         = size_t (T::*)() const noexcept;
        using same_edges_sig   = bool (T::*)(const T&) const noexcept;
        using shared_edges_sig = std::vector<typename T::EdgeT> (T::*)(const T&) const noexcept;
 
        using checkResult = conjunction<
            std::is_same<index_sig,         decltype(&T::index)>,
            std::is_same<edge_sig,          decltype(&T::edge)>,
            //std::is_same<edges_sig,         decltype(&T::edges)>,
            std::is_same<size_sig,          decltype(&T::size)>,
            std::is_same<same_edges_sig,    decltype(&T::hasSameEdges)>,
            std::is_same<shared_edges_sig,  decltype(&T::sharedEdges)>
            >;
      };
  }
 
 
  namespace {
    template <typename EdgeT>
    using isCAxis = std::enable_if_t<std::is_floating_point<EdgeT>::value>;
 
    template <typename T>
    struct hasWidth : std::false_type {};
 
    template <>
    struct hasWidth<std::string> : std::true_type {};
  }
 
 
  template <typename T, typename>
  class Axis;
 
 
  template <typename T, typename = void>
  class Axis {
    public:
    using EdgeT = T;
    using ContainerT = std::vector<T>;
    using const_iterator = typename ContainerT::const_iterator;
 
 
 
    Axis() { }
 
    Axis(const std::vector<T>& edges);
 
    Axis(std::vector<T>&& edges);
 
    Axis(std::initializer_list<T>&& edges);
 
    Axis(Axis<T>&& other) : _edges(std::move(other._edges)) {}
 
    Axis(const Axis<T>& other) : _edges(other._edges) {}
 
    Axis& operator=(const Axis& other) {
      if (this != &other)  _edges = other._edges;
      return *this;
    }
 
    Axis& operator=(Axis&& other) {
      if (this != &other)  _edges = std::move(other._edges);
      return *this;
    }
 
 
    // @{
 
    void _renderYODA(std::ostream& os) const noexcept {
      os << "[";
      for (size_t i = 0; i < _edges.size(); ++i) {
        if (i)  os << ", ";
        if constexpr(std::is_same<T, std::string>::value) {
          os << std::quoted(_edges[i]);
        }
        else {
          os << _edges[i];
        }
      }
      os << "]";
    }
 
    std::string type() const noexcept { return TypeID<EdgeT>::name(); }
 
    int maxEdgeWidth() const noexcept {
      int maxwidth = 0;
      if constexpr (hasWidth<EdgeT>::value) {
        auto it = std::max_element(_edges.begin(), _edges.end(),
                                   [](const auto& a, const auto& b) {
                                      return a.size() < b.size();
        });
        maxwidth = (*it).size();
      }
      return maxwidth;
    }
 
 
    size_t index(const T& x) const;
 
    EdgeT edge(const size_t i) const;
 
    std::vector<EdgeT> edges() const noexcept;
 
    const_iterator begin() const { return _edges.cbegin(); }
 
    const_iterator end() const { return _edges.cend(); }
 
    size_t size() const noexcept;
 
    size_t numBins(const bool includeOverflows = false) const noexcept;
 
    bool hasSameEdges(const Axis<T>& other) const noexcept;
 
    std::vector<T> sharedEdges(const Axis<T>& other) const noexcept;
 
    bool isSubsetEdges(const Axis<T>& other) const noexcept;
 
    protected:
 
    // @{
 
    void fillEdges(std::vector<EdgeT>&& edges) noexcept;
 
    // @}
 
    std::vector<T> _edges;
 
  };
 
 
  // @todo Document!
  template <typename T, typename U>
  Axis<T, U>::Axis(const std::vector<T>& edges) : Axis(std::vector<T>(edges)) {
      static_assert(MetaUtils::checkConcept<Axis<EdgeT>, YODAConcepts::AxisImpl>(),
        "Axis<T> should implement Axis concept.");
  }
 
  template <typename T, typename U>
  Axis<T, U>::Axis(std::vector<T>&& edges) {
    fillEdges(std::move(edges));
  }
 
  template <typename T, typename U>
  Axis<T, U>::Axis(std::initializer_list<T>&& edges) {
    fillEdges(std::vector<T>{edges});
  }
 
  template <typename T, typename U>
  size_t Axis<T, U>::index(const T& x) const {
    auto it = std::find(this->_edges.begin(), this->_edges.end(), x);
    if (it == this->_edges.end())  return 0; // otherflow
    return it - this->_edges.begin() + 1;
  }
 
  template <typename T, typename U>
  typename Axis<T, U>::EdgeT Axis<T, U>::edge(const size_t i) const {
    if (this->_edges.empty()) {
      throw RangeError("Axis has no edges!");
     }
    if (!i || i > this->_edges.size()) {
      throw RangeError("Invalid index, must be in range 1.." + std::to_string(this->_edges.size()));
    }
    return this->_edges.at(i-1);
  }
 
  template <typename T, typename U>
  std::vector<typename Axis<T, U>::EdgeT> Axis<T, U>::edges() const noexcept {
    return this->_edges;
  }
 
  template <typename T, typename U>
  size_t Axis<T, U>::size() const noexcept {
    return numBins(true);
  }
 
  template <typename T, typename U>
  size_t Axis<T, U>::numBins(const bool includeOverflows) const noexcept {
    return _edges.size() + (includeOverflows? 1 : 0);
  }
 
  template <typename T, typename U>
  bool Axis<T, U>::hasSameEdges(const Axis<T>& other) const noexcept {
    return _edges.size() == other._edges.size() &&
           std::equal(_edges.begin(), _edges.end(), other._edges.begin());
  }
 
  template <typename T, typename U>
  std::vector<T> Axis<T, U>::sharedEdges(const Axis<T>& other) const noexcept {
    std::vector<EdgeT> res;
    const auto& otherBegin = other._edges.begin();
    const auto& otherEnd = other._edges.end();
    for (const T& edge : this->_edges) {
      if (std::find(otherBegin, otherEnd, edge) != otherEnd)
        res.emplace_back(std::move(edge));
    }
    return res;
  }
 
  template <typename T, typename U>
  bool Axis<T, U>::isSubsetEdges(const Axis<T>& other) const noexcept {
    if (_edges.size() == other._edges.size()) return hasSameEdges(other);
 
    std::vector<T> edgesLhs(edges());
    std::vector<T> edgesRhs(other.edges());
 
    std::sort(edgesLhs.begin(), edgesLhs.end());
    std::sort(edgesRhs.begin(), edgesRhs.end());
 
    return std::includes(edgesLhs.begin(), edgesLhs.end(),
                         edgesRhs.begin(), edgesRhs.end());
  }
 
  template <typename T, typename U>
  void Axis<T, U>::fillEdges(std::vector<EdgeT>&& edges) noexcept {
    for (auto& edge : edges) {
      if (std::find(this->_edges.begin(),
                    this->_edges.end(), edge) == this->_edges.end()) // no duplicate edges allowed
        _edges.emplace_back(std::move(edge));
    }
  }
 
 
 
 
 
 
 
 
  template <typename T>
  class Axis<T, isCAxis<T>> {
    public:
    using EdgeT = T;
    using ContainerT = std::vector<T>;
    using const_iterator = typename ContainerT::const_iterator;
    using CAxisT = isCAxis<T>;
 
    Axis() {
      _updateEdges({});
      _setEstimator();
    }
 
    Axis(const Axis<T, CAxisT>& other);
 
    Axis(Axis<T, CAxisT>&& other)
      : _est(other._est),
        _maskedBins(std::move(other._maskedBins)),
        _edges(std::move(other._edges)) {}
 
    Axis(std::vector<EdgeT> edges);
 
    Axis(std::initializer_list<T>&& edges);
 
    Axis(std::vector<std::pair<EdgeT, EdgeT>> binsEdges);
 
    Axis(size_t nBins, EdgeT lower, EdgeT upper);
 
    Axis& operator=(const Axis& other) {
      if (this != &other) {
        _est = other._est;
        _maskedBins = other._maskedBins;
        _edges = other._edges;
      }
      return *this;
    }
 
    Axis& operator=(Axis&& other) {
      if (this != &other) {
        _est = std::move(other._est);
        _maskedBins = std::move(other._maskedBins);
        _edges = std::move(other._edges);
      }
      return *this;
    }
 
    // /// Explicit constructor, specifying the edges and estimation strategy
    // Axis(const std::vector<double>& edges, bool log) {
    //   _updateEdges(edges);
    //   // Internally use a log or linear estimator as requested
    //   if (log) {
    //     _est.reset(new LogBinEstimator(edges.size()-1, edges.front(), edges.back()));
    //   } else {
    //     _est.reset(new LinBinEstimator(edges.size()-1, edges.front(), edges.back()));
    //   }
    // }
 
 
    size_t index(const EdgeT& x) const;
 
    size_t size() const noexcept;
 
    size_t numBins(const bool includeOverflows = false) const noexcept;
 
    EdgeT edge(const size_t i) const;
 
    std::vector<EdgeT> edges() const noexcept;
 
    const_iterator begin() const { return _edges.cbegin(); }
 
    const_iterator end() const { return _edges.cend(); }
 
    bool hasSameEdges(const Axis<EdgeT, CAxisT>& other) const noexcept;
 
    std::vector<T> sharedEdges(const Axis<EdgeT, CAxisT>& other) const noexcept;
 
    bool isSubsetEdges(const Axis<EdgeT, CAxisT>& other) const noexcept;
 
    std::vector<size_t> maskedBins() const noexcept {  return _maskedBins; }
 
    // ssize_t index_inrange(double x) const {
    //   const size_t i = index(x);
    //   /// Change to <= and >=
    //   if (i == 0 || i == _edges.size()-1) return -1;
    //   return i;
    // }
 
    // @{
 
    void _renderYODA(std::ostream& os) const noexcept {
      os << "[";
      size_t nEdges = _edges.size() - 2; // exclude under-/overflow
      for (size_t i = 0; i < nEdges; ++i) {
        if (i) {
          os << ", ";
        }
        os << _edges[i+1];
      }
      os << "]";
    }
 
    std::string type() const noexcept { return TypeID<EdgeT>::name(); }
 
    // @}
 
    // @{
 
    void mergeBins(std::pair<size_t, size_t>);
 
    // @}
 
    // @{
 
    std::vector<T> widths(const bool includeOverflows = false) const noexcept {
      // number of widths = number of edges - 1
      // unless you exclude under-/overflow
      size_t offset = includeOverflows? 1 : 3;
      std::vector<T> ret(_edges.size()-offset);
      size_t start = 1 + !includeOverflows;
      size_t end = _edges.size() - !includeOverflows;
      for (size_t i = start; i < end; ++i) {
        ret[i-start] = _edges[i] - _edges[i-1];
      }
      return ret;
    }
 
    EdgeT width(size_t binNum) const noexcept {
      return _edges[binNum+1] - _edges[binNum];
    }
 
    EdgeT max(size_t binNum) const noexcept {
      return _edges[binNum+1];
    }
 
 
    EdgeT min(size_t binNum) const noexcept {
      return _edges[binNum];
    }
 
    EdgeT mid(size_t binNum) const noexcept {
      if(binNum == 0)
        return std::numeric_limits<EdgeT>::min();
      if (binNum == (numBins(true) - 1))
        return std::numeric_limits<EdgeT>::max();
 
      EdgeT minVal = min(binNum);
      return (max(binNum) - minVal)/2 + minVal;
    }
    // @}
 
    protected:
 
    void _updateEdges(std::vector<EdgeT>&& edges) noexcept;
 
    void _setEstimator() noexcept;
 
    ssize_t _linsearch_forward(size_t istart, double x, size_t nmax) const noexcept;
 
    ssize_t _linsearch_backward(size_t istart, double x, size_t nmax) const noexcept;
 
    size_t _bisect(double x, size_t imin, size_t imax) const noexcept;
 
    std::shared_ptr<BinEstimator> _est;
 
    std::vector<size_t> _maskedBins;
 
    std::vector<T> _edges;
 
  };
 
  template <typename T>
  Axis<T, isCAxis<T>>::Axis(const Axis<T, CAxisT>& other) {
    static_assert(MetaUtils::checkConcept<Axis<EdgeT>, YODAConcepts::AxisImpl>(),
      "Axis<T> should implement Axis concept.");
 
    _est = other._est;
    _edges = other._edges;
    _maskedBins = other._maskedBins;
  }
 
  template <typename T>
  Axis<T, isCAxis<T>>::Axis(const size_t nBins, const EdgeT lower, const EdgeT upper) {
    if(upper <= lower)
      throw(std::logic_error("Upper bound should be larger than lower."));
    _edges.resize(nBins+1+2);
    double step = (upper - lower) / nBins;
 
    _edges[0] = -std::numeric_limits<double>::infinity();
 
    _edges[1] = lower;
 
    for(size_t i = 2; i < _edges.size()-1; i++) {
      _edges[i] = _edges[i-1] + step;
    }
 
    _edges[_edges.size()-1] = std::numeric_limits<double>::infinity();
 
    _setEstimator();
  }
 
  template <typename T>
  Axis<T, isCAxis<T>>::Axis(std::vector<std::pair<EdgeT, EdgeT>> binsEdges) {
    if (binsEdges.empty()) throw BinningError("Expected at least one discrete edge");
 
    std::sort(binsEdges.begin(), binsEdges.end(), [](auto &left, auto &right) {
      return left.first < right.first;
    });
 
    _edges.resize(binsEdges.size()*2+2);
 
    _edges[0] = -std::numeric_limits<double>::infinity();
 
 
    /*             Edges pairs from binsEdges vector
    */
 
    size_t i = 1;
    for (const auto& pair : binsEdges) {
      if (!fuzzyGtrEquals(pair.first, _edges[i-1])) throw BinningError("Bin edges shouldn't intersect");
      if (i == 1 && std::isinf(pair.first) && pair.first < 0) {
        _edges[i++] = pair.second;
        continue;
      }
      if (fuzzyEquals(pair.first, _edges[i-1])) {
        _edges[i++] = pair.second; 
        continue;
      }
      if (i != 1 && pair.first > _edges[i-1]) {
        _maskedBins.emplace_back(i-1);
      }
      _edges[i++] = pair.first;
      _edges[i++] = pair.second;
    }
 
    _edges[i] = std::numeric_limits<double>::infinity();
 
    _edges.resize(i+1); 
 
    _setEstimator();
  }
 
  template <typename T>
  Axis<T, isCAxis<T>>::Axis(std::vector<EdgeT> edges) {
    std::sort(edges.begin(), edges.end());
    edges.erase( std::unique(edges.begin(), edges.end()), edges.end() );
 
    _updateEdges(std::move(edges));
 
    _setEstimator();
  }
 
  template <typename T>
  Axis<T, isCAxis<T>>::Axis(std::initializer_list<T>&& edgeList) {
    std::vector<T> edges{edgeList};
    std::sort(edges.begin(), edges.end());
    edges.erase( std::unique(edges.begin(), edges.end()), edges.end() );
 
    _updateEdges(std::move(edges));
 
    _setEstimator();
  }
 
  template <typename T>
  size_t Axis<T, isCAxis<T>>::index(const EdgeT& x) const {
      if (_edges.size() <= 2) throw BinningError("Axis initialised without specifying edges");
      // Check overflows
      if (std::isinf(x)) { return x < 0? 0 : _edges.size() - 2; }
      // Get initial estimate
      size_t index = std::min(_est->estindex(x), _edges.size()-2);
      // Return now if this is the correct bin
      if (x >= this->_edges[index] && x < this->_edges[index+1]) return index;
 
      // Otherwise refine the estimate, if x is not exactly on a bin edge
      if (x > this->_edges[index]) {
        const ssize_t newindex = _linsearch_forward(index, x, SEARCH_SIZElc);
        index = (newindex > 0) ? newindex : _bisect(x, index, this->_edges.size()-1);
      } else if (x < this->_edges[index]) {
        const ssize_t newindex = _linsearch_backward(index, x, SEARCH_SIZElc);
        index = (newindex > 0) ? newindex : _bisect(x, 0, index+1);
      }
 
      assert(x >= this->_edges[index] && (x < this->_edges[index+1] || std::isinf(x)));
      return index;
  }
 
  template <typename T>
  size_t Axis<T, isCAxis<T>>::size() const noexcept {
    return numBins(true); // number of visible + 2 for +-inf
  }
 
  template <typename T>
  size_t Axis<T, isCAxis<T>>::numBins(const bool includeOverflows) const noexcept {
    if (_edges.size() < 3)  return 0; // no visible bin edge
    return this->_edges.size() - (includeOverflows? 1 : 3); // has 2 extra edges for +-inf
  }
 
  template <typename T>
  typename Axis<T, isCAxis<T>>::EdgeT Axis<T, isCAxis<T>>::edge(const size_t i) const {
    return this->_edges.at(i);
  }
 
  template <typename T>
  std::vector<typename Axis<T, isCAxis<T>>::EdgeT> Axis<T, isCAxis<T>>::edges() const noexcept {
    return this->_edges;
  }
 
  template <typename T>
  bool Axis<T, isCAxis<T>>::hasSameEdges(const Axis<T, CAxisT>& other) const noexcept{
    if (this->numBins(true) != other.numBins(true)) return false;
    for (size_t i = 1; i < this->numBins(true) - 1; i++) {
      if (!fuzzyEquals(edge(i), other.edge(i))) return false;
    }
    return true;
  }
 
  template <typename T>
  std::vector<T> Axis<T, isCAxis<T>>::sharedEdges(const Axis<T, CAxisT>& other) const noexcept {
    std::vector<T> intersection;
 
    if (_edges.size() > 2 && other._edges.size() > 2) {
      std::set_intersection(std::next(_edges.begin()), std::prev(_edges.end()),
                            std::next(other._edges.begin()), std::prev(other._edges.end()),
                            std::back_inserter(intersection));
    }
 
    std::vector<T> rtn;
    rtn.reserve(intersection.size()+2);
 
    rtn.emplace_back(-std::numeric_limits<Axis<T, isCAxis<T>>::EdgeT>::infinity());
    rtn.insert(std::next(rtn.begin()),
                  std::make_move_iterator(intersection.begin()),
                  std::make_move_iterator(intersection.end()));
    rtn.emplace_back(std::numeric_limits<Axis<T, isCAxis<T>>::EdgeT>::infinity());
 
    return rtn;
  }
 
  template <typename T>
  bool Axis<T, isCAxis<T>>::isSubsetEdges(const Axis<T, CAxisT>& other) const noexcept {
    if (_edges.size() == other._edges.size()) return hasSameEdges(other);
 
    if (_edges.size() > 2 && other._edges.size() > 2) {
      return std::includes(std::next(_edges.begin()), std::prev(_edges.end()),
                            std::next(other._edges.begin()), std::prev(other._edges.end()));
    }
 
    return true;
  }
 
  template <typename T>
  void Axis<T, isCAxis<T>>::mergeBins(std::pair<size_t, size_t> mergeRange) {
    if (_edges.size() <= 2) throw BinningError("Axis initialised without specifying edges");
    if (mergeRange.second < mergeRange.first) throw RangeError("Upper index comes before lower index");
    if (mergeRange.second >= numBins(true)) throw RangeError("Upper index exceeds last visible bin");
    _edges.erase(_edges.begin()+mergeRange.first+1, _edges.begin() + mergeRange.second + 1);
    // masked bins above the merge range need to be re-indexed
    // masked bins within the merge range need to be unmasked
    std::vector<size_t> toRemove;
    size_t mrange = mergeRange.second - mergeRange.first;
    for (size_t i = 0; i < _maskedBins.size(); ++i) {
      if (_maskedBins[i] > mergeRange.second)  _maskedBins[i] -= mrange;
      else if (_maskedBins[i] >= mergeRange.first) toRemove.push_back(_maskedBins[i]);
    }
    if (toRemove.size()) {
      _maskedBins.erase(
        std::remove_if(_maskedBins.begin(), _maskedBins.end(), [&](const auto& idx) {
          return std::find(toRemove.begin(), toRemove.end(), idx) != toRemove.end();
      }), _maskedBins.end());
    }
  }
 
 
  template <typename T>
  void Axis<T, isCAxis<T>>::_updateEdges(std::vector<EdgeT>&& edges) noexcept {
    // Array of in-range edges, plus underflow and overflow sentinels
    _edges.clear();
 
    // Move vector and insert -+inf at ends
    _edges.emplace_back(-std::numeric_limits<Axis<T, isCAxis<T>>::EdgeT>::infinity());
    _edges.insert(std::next(_edges.begin()),
                  std::make_move_iterator(edges.begin()),
                  std::make_move_iterator(edges.end()));
    _edges.emplace_back(std::numeric_limits<Axis<T, isCAxis<T>>::EdgeT>::infinity());
  }
 
  template <typename T>
  void Axis<T, isCAxis<T>>::_setEstimator() noexcept {
    if (_edges.empty()) {
        _est = std::make_shared<LinBinEstimator>(0, 0, 1);
      } else if (_edges.front() <= 0.0) {
        _est = std::make_shared<LinBinEstimator>(_edges.size()-1, _edges.front(), _edges.back());
      } else {
        LinBinEstimator linEst(_edges.size()-1, _edges.front(), _edges.back());
        LogBinEstimator logEst(_edges.size()-1, _edges.front(), _edges.back());
 
        // Calculate mean index estimate deviations from the correct answers (for bin edges)
        double logsum = 0, linsum = 0;
        for (size_t i = 0; i < _edges.size(); i++) {
          logsum += logEst(_edges[i]) - i;
          linsum += linEst(_edges[i]) - i;
        }
        const double log_avg = logsum / _edges.size();
        const double lin_avg = linsum / _edges.size();
 
        // This also implicitly works for NaN returned from the log There is a
        // subtle bug here if the if statement is the other way around, as
        // (nan < linsum) -> false always.  But (nan > linsum) -> false also.
        if (log_avg < lin_avg) { //< Use log estimator if its avg performance is better than lin
          _est = std::make_shared<LogBinEstimator>(logEst);
        } else { // Else use linear estimation
          _est = std::make_shared<LinBinEstimator>(linEst);
        }
    }
  }
 
 
  template <typename T>
  ssize_t Axis<T, isCAxis<T>>::_linsearch_forward(size_t istart, double x, size_t nmax) const noexcept {
    assert(x >= this->_edges[istart]); // assumption that x >= start is wrong
    for (size_t i = 0; i < nmax; i++) {
      const size_t j = istart + i + 1; // index of _next_ edge
      if (j > this->_edges.size()-1) return -1;
      if (x < this->_edges[j]) {
        assert(x >= this->_edges[j-1] && (x < this->_edges[j] || std::isinf(x)));
        return j-1; // note one more iteration needed if x is on an edge
      }
    }
    return -1;
  }
 
  template <typename T>
  ssize_t Axis<T, isCAxis<T>>::_linsearch_backward(size_t istart, double x, size_t nmax) const noexcept {
    assert(x < this->_edges[istart]); // assumption that x < start is wrong
    for (size_t i = 0; i < nmax; i++) {
      const int j = istart - i - 1; // index of _next_ edge (working backwards)
      if (j < 0) return -1;
      if (x >= this->_edges[j]) {
        assert(x >= this->_edges[j] && (x < this->_edges[j+1] || std::isinf(x)));
        return (ssize_t) j; // note one more iteration needed if x is on an edge
      }
    }
    return -1;
  }
 
  template <typename T>
  size_t Axis<T, isCAxis<T>>::_bisect(double x, size_t imin, size_t imax) const noexcept {
    size_t len = imax - imin;
    while (len >= BISECT_LINEAR_THRESHOLDlc) {
      const size_t half = len >> 1;
      const size_t imid = imin + half;
      if (x >= this->_edges[imid]) {
        if (x < this->_edges[imid+1]) return imid; // Might as well return directly if we get lucky!
        imin = imid;
      } else {
        imax = imid;
      }
      len = imax - imin;
    }
    assert(x >= this->_edges[imin] && (x < this->_edges[imax] || std::isinf(x)));
    return _linsearch_forward(imin, x, BISECT_LINEAR_THRESHOLDlc);
  }
 
}
 
#endif