EIC code displayed by LXR master/​include/​YODA/​BinnedAxis.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-09-11 09:31:34

 // -*- C++ -*-
 //
 // This file is part of YODA -- Yet more Objects for Data Analysis
 // Copyright (C) 2008-2025 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;
 
       /// @brief Axis concept
       template <typename T>
       struct AxisImpl {
 
         /// @note const/volatile function parameter qualifiers are not enforced:
         ///  C++ Standard (C++ 17, 16.1 Overloadable declarations):
         ///   (3.4) — Parameter declarations that differ only in the presence or
         ///   absence of const and/or volatile are equivalent. That is, the const and
         ///   volatile type-specifiers for each parameter type are ignored when
         ///   determining which function is being declared, defined, or called.
         ///
         /// noexcept qualifier is not enforced too since it's not part of type.
 
         /// @brief Function signatures
         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)>
             >;
       };
   }
 
 
   /// @note Anonymous namespace to limit visibility to this file
   namespace {
     /// Checks if edge types is continuous
     template <typename EdgeT>
     using isCAxis = std::enable_if_t<std::is_floating_point<EdgeT>::value>;
 
     /// Checks if edge type has width measure
     template <typename T>
     struct hasWidth : std::false_type {};
 
     template <>
     struct hasWidth<std::string> : std::true_type {};
   }
 
 
   template <typename T, typename>
   class Axis;
 
 
   /// @brief Discrete axis with edges of non-floating-point-type
   ///
   /// @note Based on *unsorted* std::vector<T>
   template <typename T, typename = void>
   class Axis {
     public:
     using EdgeT = T;
     using ContainerT = std::vector<T>;
     using const_iterator = typename ContainerT::const_iterator;
 
     /// @name Constructors
     /// @{
 
     /// @brief Nullary constructor for unique pointers etc.
     Axis() { }
 
     /// @brief Constructs discrete Axis from edges vector.
     ///
     /// @note Vector is not sorted by constructor.
     Axis(const std::vector<T>& edges);
 
     /// @brief Constructs discrete Axis from edges vector (rvalue).
     ///
     /// @note Vector is not sorted by constructor.
     Axis(std::vector<T>&& edges);
 
     /// @brief Constructs discrete Axis from an initializer list.
     ///
     /// @note Vector is not sorted by constructor.
     Axis(std::initializer_list<T>&& edges);
 
     /// @brief Move constructs Axis
     Axis(Axis<T>&& other) : _edges(std::move(other._edges)) {}
 
     /// @brief Copy constructs Axis
     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;
     }
 
     /// @}
 
     /// @name I/O
     /// @{
 
     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 << "]";
     }
 
     /// @brief Returns a string representation of EdgeT
     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;
     }
 
     /// @}
 
     /// @brief Returns index of edge x
     ///
     /// @note Returns 0 if there is no @a x on this axis.
     size_t index(const T& x) const;
 
     /// @brief Returns edge corresponding to index @a i
     EdgeT edge(const size_t i) const;
 
     /// @brief Returns a copy of the container of edges.
     std::vector<EdgeT> edges() const noexcept;
 
     /// @brief Returns the const begin iterator for the edges container
     const_iterator begin() const { return _edges.cbegin(); }
 
     /// @brief Returns the const end iterator for the edges container
     const_iterator end() const { return _edges.cend(); }
 
     /// @brief Returns number of edges + 1 for this Axis
     size_t size() const noexcept;
 
     /// @brief Returns number of bins of this axis
     size_t numBins(const bool includeOverflows=false) const noexcept;
 
     /// @brief Checks if two axes have exactly the same edges
     bool hasSameEdges(const Axis<T>& other) const noexcept;
 
     /// @brief Finds shared between Axes edges.
     ///
     /// @note Not sorted, order similar to initial (in constructor) not guaranteed.
     std::vector<T> sharedEdges(const Axis<T>& other) const noexcept;
 
     /// @brief Check if other axis edges are a subset of edges of this one
     /// @remark Should it be symmetrical? E.g. ax1.subsetEdges(ax2) == ax2.subsetEdges(ax1)?
     bool isSubsetEdges(const Axis<T>& other) const noexcept;
 
     protected:
 
     /// @name Utility
     /// @{
 
     /// @brief Fills edge storage. Used in constructors.
     void fillEdges(std::vector<EdgeT>&& edges) noexcept;
 
     /// @}
 
     /// @brief Axis edges
     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)) {
       /// @brief Concept check shall appear inside body of type's member function
       /// or outside of type, since otherwise type is considered incomplete.
       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;
   }
 
   /// Includes +1 for the otherflow bin
   template <typename T, typename U>
   size_t Axis<T, U>::size() const noexcept {
     return numBins(true);
   }
 
   /// Includes +1 for the otherflow bin
   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());
 
     /// @note std::includes demands sorted ranges
     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));
     }
   }
 
 
 
 
 
 
 
 
   /// @brief Continuous axis with floating-point-type edges
   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>;
 
     /// @brief Nullary constructor for unique pointers etc.
     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)) {}
 
     /// @brief Constructs continuous Axis from a vector.
     ///
     /// @note Edges are sorted on construction stage.
     Axis(std::vector<EdgeT> edges);
 
     /// @brief Constructs continuous Axis from an initializer list.
     ///
     /// @note Edges are sorted on construction stage
     Axis(std::initializer_list<T>&& edges);
 
     /// @note Vector shouldn't contain any intersecting pairs,
     /// e.g. pair1.second > pair2.first. Order of pairs does not
     /// matter. BinsEdges is sorted on construction stage.
     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()));
     //   }
     // }
 
 
     /// @brief Bin searcher
     ///
     /// @author David Mallows
     /// @author Andy Buckley
     ///
     /// Handles low-level bin lookups using a hybrid algorithm that is
     /// considerably faster for regular (logarithmic or linear) and near-regular
     /// binnings. Comparable performance for irregular binnings.
     ///
     /// The reason this works is that linear search is faster than bisection
     /// search up to about 32-64 elements. So we make a guess, and we then do a
     /// linear search. If that fails, then we bisect on the remainder,
     /// terminating once bisection search has got the range down to about 32. So
     /// we actually pay for the fanciness of predicting the bin out of speeding
     /// up the bisection search by finishing it with a linear search. So in most
     /// cases, we get constant-time lookups regardless of the space.
     ///
     size_t index(const EdgeT& x) const;
 
     /// @brief Returns number of edges + 2 for this Axis
     size_t size() const noexcept;
 
     /// @brief Returns number of bins of this axis
     size_t numBins(const bool includeOverflows=false) const noexcept;
 
     /// @brief Returns edge corresponding to index @a i
     EdgeT edge(const size_t i) const;
 
     /// @brief Returns a copy of all axis edges. Includes -inf and +inf edges.
     std::vector<EdgeT> edges() const noexcept;
 
     /// @brief Returns the const begin iterator for the edges container
     const_iterator begin() const { return _edges.cbegin(); }
 
     /// @brief Returns the const end iterator for the edges container
     const_iterator end() const { return _edges.cend(); }
 
     /// @brief Check if two BinnedAxis objects have the same edges
     bool hasSameEdges(const Axis<EdgeT, CAxisT>& other) const noexcept;
 
     /// @brief Find edges which are shared between BinnedAxis objects, within numeric tolerance
     /// @note The return vector is sorted and includes -inf and inf
     std::vector<T> sharedEdges(const Axis<EdgeT, CAxisT>& other) const noexcept;
 
     /// @brief Check if other axis edges are a subset of edges of this one
     /// @remark Should it be symmetrical? E.g. ax1.subsetEdges(ax2) == ax2.subsetEdges(ax1)?
     bool isSubsetEdges(const Axis<EdgeT, CAxisT>& other) const noexcept;
 
     /// @brief Returns the masked indices
     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;
     // }
 
     /// @name I/O
     /// @{
 
     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 << "]";
     }
 
     /// @brief Returns a string representation of EdgeT
     std::string type() const noexcept { return TypeID<EdgeT>::name(); }
 
     /// @}
 
     /// @name Transformations
     /// @{
 
     /// @brief Merges bins, e.g. removes edges in between.
     ///
     /// @note At least 1 non-overflow bin must exist after merging.
     void mergeBins(std::pair<size_t, size_t>);
 
     /// @}
 
     /// @name Space characteristics
     /// @{
 
     /// @brief List of bin widths on this axis.
     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;
     }
 
     /// @brief List of max values for each bin on this axis.
     std::vector<T> maxs(const bool includeOverflows=false) const noexcept {
       // number of maxs = 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];
       }
       return ret;
     }
 
     /// @brief List of min values for each bin on this axis.
     std::vector<T> mins(const bool includeOverflows=false) const noexcept {
       // number of mins = 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-1];
       }
       return ret;
     }
 
     /// @brief List of mid values for each bin on this axis.
     std::vector<T> mids(const bool includeOverflows=false) const noexcept {
       // number of mins = 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] = 0.5*(_edges[i-1] + _edges[i]);
       }
       return ret;
     }
 
     /// @brief Width of bin side on this axis.
     EdgeT width(size_t binNum) const noexcept {
       return _edges[binNum+1] - _edges[binNum];
     }
 
     /// @brief Max of bin side on this axis.
     EdgeT max(size_t binNum) const noexcept {
       return _edges[binNum+1];
     }
 
 
     /// @brief Min of bin side on this axis.
     EdgeT min(size_t binNum) const noexcept {
       return _edges[binNum];
     }
 
     /// @brief Mid of bin side on this axis.
     EdgeT mid(size_t binNum) const noexcept {
       /// @note Corner bins are overflow bins, thus infinite limit values.
       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:
 
     /// @brief Set the edges array and related member variables
     void _updateEdges(std::vector<EdgeT>&& edges) noexcept;
 
     /// @brief Set the estimator.
     /// @note Used in constructors.
     void _setEstimator() noexcept;
 
     /// @brief Linear search in the forward direction
     ///
     /// Return bin index or -1 if not found within linear search range. Assumes that edges[istart] <= x
     ssize_t _linsearch_forward(size_t istart, double x, size_t nmax) const noexcept;
 
     /// @brief Linear search in the backward direction
     ///
     /// Return bin index or -1 if not found within linear search range. Assumes that edges[istart] > x
     ssize_t _linsearch_backward(size_t istart, double x, size_t nmax) const noexcept;
 
     /// Truncated bisection search, adapted from C++ std lib implementation
     size_t _bisect(double x, size_t imin, size_t imax) const noexcept;
 
     /// BinEstimator object to be used for making fast bin index guesses
     std::shared_ptr<BinEstimator> _est;
 
     /// @brief masked bins indices.
     std::vector<size_t> _maskedBins;
 
     /// @brief Axis edges
     std::vector<T> _edges;
 
   };
 
   template <typename T>
   Axis<T, isCAxis<T>>::Axis(const Axis<T, CAxisT>& other) {
     /// @brief Concept check shall appear inside body of type's member function
     /// or outside of type, since otherwise type is considered incomplete.
     /// @note Concept check appears once to check whether the type Axis satisfies
     /// the concept.
     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 pair of edges");
 
     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
     ///                      ____|____
     ///            __{1,3}__/         \_{5,6}__
     ///           /        |    GAP    |       \
     ///   -inf    1        3   MASKED  5       6   +inf
     ///     | BIN |  BIN   |    BIN    |  BIN  | BIN |
     ///    i=0   i=1      i=2         i=3     i=4   i=5
     */
 
     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; /// Merge if equal
         continue;
       }
       if (i != 1 && pair.first > _edges[i-1]) {
         /// @note If there is a gap, mark bin as masked.
         _maskedBins.emplace_back(i-1);
       }
       _edges[i++] = pair.first;
       _edges[i++] = pair.second;
     }
 
     _edges[i] = std::numeric_limits<double>::infinity();
 
     _edges.resize(i+1); /// In case there have been merges. +1 to account for +inf.
 
     _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");
       // Only one edge (i.e. axis has under- and overflow bin, but no visible bin)
       if (_edges.size() == 3) return x >= _edges[1];
       // 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 includeOverflows? 1 : 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++) {
       /// @todo Be careful about using fuzzyEquals... should be an exact comparison?
       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;
 
     /// Skip if any of axes only has two limit edges
     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);
 
     /// Skip if any of axes only has two limit edges
     if (_edges.size() > 2 && other._edges.size() > 2) {
       /// @note std::includes demands sorted ranges
       return std::includes(std::next(_edges.begin()), std::prev(_edges.end()),
                             std::next(other._edges.begin()), std::prev(other._edges.end()));
     }
 
     /// Since one of the axes consists only of limits (-inf, +inf), it's a
     /// subset of the other one
     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 {
 
     // Empty set for nullary constructor
     if (_edges.size() <= 2) {
       _est = std::make_shared<LinBinEstimator>(0, 0, 1);
       return;
     }
 
     // There is at least one visible edge
     int front = 1, back = (int)_edges.size()-2;
     if (_edges[front] <= 0.0) {
       _est = std::make_shared<LinBinEstimator>(back - front, _edges[front], _edges[back]);
     }
     else {
       LinBinEstimator linEst(back - front, _edges[front], _edges[back]);
       LogBinEstimator logEst(back - front, _edges[front], _edges[back]);
 
       // Calculate mean index estimate deviations from the correct answers (for bin edges)
       double logsum = 0, linsum = 0;
       for (int i = front; i <= back; ++i) {
         logsum += std::abs(logEst(_edges[i]) - double(i-1));
         linsum += std::abs(linEst(_edges[i]) - double(i-1));
       }
       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

This page was automatically generated by the 2.3.7 LXR engine. The LXR team