A simple binned storage implementation with bisection lookup

// This is a pretty simple implementation of a "binning", which allows you to store objects of any
// type TY in an axis denominated in any variable TX, with lookup of the bin at position x using the
// bisection std::upper_bound function. It allows irregular binnings, but not *gaps* in the binning.
// This can be made faster, e.g. by guessing at uniform binnings first, using linear search below a // certain problem size, etc., but it gets quite a bit more fiddly than I want to write on a Saturday // afternoon! See https://yoda.hepforge.org/trac/browser/include/YODA/Utils/BinSearcher.h for that.

// Definition... normally in Binning.h header file
// #pragma once

#include <cstddef>
#include <vector>
#include <algorithm>
#include <stdexcept>
#include <cassert>

template <typename TY, typename TX=double>
class Binning {
public:

  /// Constructor taking a list of bin edges
  Binning(const std::vector<TX>& binedges)
    : edges(binedges)
  {
    reset();
  }

  /// Constructor taking lists of bin edges and values
  Binning(const std::vector<TX>& binedges, const std::vector<TY>& binvalues)
    : edges(binedges), values(binvalues)
  {
    check();
  }

  /// Get the bin index enclosing position @a x
  size_t get_index(const TX& x) const {
    check();
    if (x < edges.front()) throw std::out_of_range("x value underflow");
    if (x > edges.back()) throw std::out_of_range("x value overflow");
    // Find the closest knot below the requested value
    size_t i = upper_bound(edges.begin(), edges.end(), x) - edges.begin();
    if (i == edges.size()) i -= 1; // can't return the last knot index
    i -= 1; // have to step back to get the knot <= x behaviour
    return i;
  }

  /// Get the value in bin number @a ix
  const TY& get_at_index(size_t ix) const {
    check();
    return values[ix];
  }

  /// Get the value in the bin at position @a x
  const TY& get_at(const TX& x) const {
    return get_at_index(get_index(x));
  }

  /// Get the value in bin number @a ix
  void set_at_index(size_t ix, const TY& val) {
    check();
    values[ix] = val;
  }

  /// Get the value in the bin at position @a x
  void set_at(const TX& x, const TY& val) {
    set_at_index(get_index(x), val);
  }

  /// Get the number of bins
  size_t num_bins() {
    check();
    return values.size();
  }

  /// Clear the bin contents (but leave the binning intact)
  void reset() {
    assert(edges.size() > 1);
    values.clear();
    values.resize(edges.size()-1);
    check();
  }

  /// Check consistency of the edges and values vectors
  void check() const {
    if (edges.size() <= 1) throw std::length_error("There must be 2 or more bin edges");
    if (values.size() < 1) throw std::length_error("There must be 1 or more bin values");
    if (edges.size()-1 != values.size()) throw std::length_error("There must be one more bin edge than there are bin values");
  }

  /// The list of bin edges
  std::vector<TX> edges;

  /// The list of values
  std::vector<TY> values;

};

//////////////////////////////////////
// And usage... normally from a separate .cc file

// #include "Binning.h"
#include <iostream>
using namespace std;

struct F {
  F(int a) : i(a) {}
  int i;
};
ostream& operator << (ostream& os, const F& f) { os << "F" << f.i; return os; }

int main() {
  Binning<float> b1({{0,1,2,3,4,5}}, {{10,20,30,40,50}});
  cout << b1.get_index(3.3) << endl;
  cout << b1.get_at(3.3) << endl;
  cout << b1.get_at_index(3) << endl;

  Binning<F> b2({{0,1,2,3,4,5}}, {{10,20,30,40,50}});
  cout << b2.get_index(2.3) << endl;
  cout << b2.get_at(2.3) << endl;
  cout << b2.get_at_index(2) << endl;

  return 0;
}