histogram.icc

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
 
#include "H5Cpp.h"
 
#include <vector>
#include <cassert>
#include <array>
#include <cmath>
#include <cstdint>
#include <numeric>
#include <optional>
#include <set>
#include <string>
#include <utility>
 
#include <boost/histogram.hpp>
 
 
namespace H5Utils::hist::detail {
 
  // various stuff to inspect the boost histograms
 
  // Concepts to inspect which members the accumulator exposes.
  template <typename T>
  concept HasValue = requires(T t) { t.begin()->value(); };
  template <typename T>
  concept HasVariance = requires(T t) { t.begin()->variance(); };
  template <typename T>
  concept HasCount = requires(T t) { t.begin()->count(); };
  template <typename T>
  concept HasSumOfWeights = requires(T t) { t.begin()->sum_of_weights(); };
  template <typename T>
  concept HasSumOfWeightsSquared = requires(T t) {
    t.begin()->sum_of_weights_squared();
  };
 
  // figure out the output array type
  template <typename T>
  struct array_type
  {
    using type = typename T::value_type;
  };
  // this needs further specialization in some cases
  template <typename T>
  using iter_value_t = decltype(
    std::declval<T>().begin()->value());
  template <HasValue T>
  struct array_type<T>
  {
    using type = std::remove_cv_t<std::remove_reference_t<iter_value_t<T>>>;
  };
  template <typename T>
  using array_t = typename array_type<T>::type;
 
  // H5 typing magic
  template <typename T>
  const H5::DataType hdf5_t;
  template<>
  const H5::DataType hdf5_t<double> = H5::PredType::NATIVE_DOUBLE;
  template<>
  const H5::DataType hdf5_t<float> = H5::PredType::NATIVE_FLOAT;
  template<>
  const H5::DataType hdf5_t<int> = H5::PredType::NATIVE_INT;
  template<>
  const H5::DataType hdf5_t<long> = H5::PredType::NATIVE_LONG;
  template<>
  const H5::DataType hdf5_t<long long> = H5::PredType::NATIVE_LLONG;
  template<>
  const H5::DataType hdf5_t<unsigned long> = H5::PredType::NATIVE_ULONG;
  template<>
  const H5::DataType hdf5_t<unsigned long long> = H5::PredType::NATIVE_ULLONG;
 
 
  // function to flatten the histogram array
  template <typename T, typename F>
  std::vector<array_t<T>> get_flat_hist_array(const T& hist,
                                              const F& value_getter)
  {
    using out_t = std::vector<array_t<T>>;
    namespace bh = boost::histogram;
 
    // calculate a global index for the flat output array
    auto global_index = [&hist](const auto& hist_idx)
    {
      // This is unrolling the index to a c-style array. We calculate
      // a global index, starting with the last dimension. The last
      // dimension is just added to the global index. Each previous
      // dimension needs to multiply the index by the cumulative
      // product of the maximum index of subsequent dimensions.
      size_t global = 0;
      size_t stride = 1;
      for (size_t dim_p1 = hist.rank(); dim_p1 != 0; dim_p1--) {
        size_t dim = dim_p1 - 1;
        auto ops = hist.axis(dim).options();
        size_t underflow = (ops & bh::axis::option::underflow) ? 1 : 0;
        size_t overflow = (ops & bh::axis::option::overflow) ? 1 : 0;
        // ugly check to decide if we have to add underflow offset
        int h5idx_signed = hist_idx.index(dim) + underflow;
        if (h5idx_signed < 0) throw std::logic_error("neg signed index");
        size_t h5idx = static_cast<size_t>(h5idx_signed);
        global += h5idx * stride;
        stride *= hist.axis(dim).size() + underflow + overflow;
      }
      return global;
    };
 
    out_t out(hist.size());
    std::set<size_t> used_indices;
    for (const auto& x : bh::indexed(hist, bh::coverage::all)) {
      size_t global = global_index(x);
      if (!used_indices.insert(global).second) {
        throw std::logic_error(
          "repeat global index: " + std::to_string(global));
      }
      out.at(global) = value_getter(x);
    }
    return out;
  }
 
  // Returns true when all bin widths are equal and positive (regular axis).
  inline bool is_regular_axis(const std::vector<float>& edges) {
    if (edges.size() < 2) return false;
    const double width0 = edges.at(1) - edges.at(0);
    if (width0 <= 0.0) return false;
    for (size_t i = 2; i < edges.size(); i++) {
      const double width = edges.at(i) - edges.at(i-1);
      if (std::abs(width - width0) > 1e-5 * std::abs(width0)) return false;
    }
    return true;
  }
 
  struct Axis
  {
    std::vector<float> edges{};
    std::string name{};
    bool underflow{false};
    bool overflow{false};
    bool is_int_or_category{false};
    // set iff the axis has equal-width bins (regular axis)
    std::optional<int> n_regular_bins{};
  };
 
  template <typename T>
  std::vector<Axis> get_axes(const T& hist) {
    using out_t = std::vector<Axis>;
    namespace opts = boost::histogram::axis::option;
    out_t axes;
    for (size_t ax_n = 0; ax_n < hist.rank(); ax_n++) {
      const auto& axis = hist.axis(ax_n);
      out_t::value_type ax;
      bool underflow = (axis.options() & opts::underflow);
      bool overflow = (axis.options() & opts::overflow);
      ax.underflow = underflow;
      ax.overflow = overflow;
      for (const auto& bin: axis) {
        ax.edges.push_back(bin.lower());
      }
      // add the upper edge too, unless this is an integer/category axis
      // where every bin has the same lower and upper value — in that case
      // we have exactly N edges for N bins rather than N+1
      if (axis.size() > 0) {
        const auto& last = *axis.rbegin();
        if (last.upper() == last.lower()) {
          ax.is_int_or_category = true;
        } else {
          ax.edges.push_back(last.upper());
        }
      }
      ax.name = axis.metadata();
 
      if (is_regular_axis(ax.edges)) {
        ax.n_regular_bins = static_cast<int>(ax.edges.size()) - 1;
      }
 
      axes.push_back(std::move(ax));
    }
    return axes;
  }
 
  inline void chkerr(herr_t code, const std::string& error) {
    if (code < 0) throw std::runtime_error("error setting " + error);
  }
 
  template <typename T>
  auto product(const std::vector<T>& v) {
    return std::accumulate(
      v.begin(), v.end(), T(1), std::multiplies<hsize_t>{});
  }
 
  // Attribute helpers
 
  inline void write_str_attr(H5::H5Object& obj,
                             const std::string& key,
                             const std::string& val)
  {
    H5::StrType strtype(H5::PredType::C_S1, H5T_VARIABLE);
    H5::DataSpace scalar(H5S_SCALAR);
    H5::Attribute attr = obj.createAttribute(key, strtype, scalar);
    const char* cstr = val.c_str();
    attr.write(strtype, &cstr);
  }
 
  inline void write_bool_attr(H5::H5Object& obj,
                              const std::string& key,
                              bool val)
  {
    H5::DataSpace scalar(H5S_SCALAR);
    uint8_t v = val ? 1 : 0;
    H5::Attribute attr = obj.createAttribute(
      key, H5::PredType::NATIVE_UINT8, scalar);
    attr.write(H5::PredType::NATIVE_UINT8, &v);
  }
 
  inline void write_int_attr(H5::H5Object& obj,
                             const std::string& key,
                             int64_t val)
  {
    H5::DataSpace scalar(H5S_SCALAR);
    H5::Attribute attr = obj.createAttribute(
      key, H5::PredType::NATIVE_INT64, scalar);
    attr.write(H5::PredType::NATIVE_INT64, &val);
  }
 
  inline void write_double_attr(H5::H5Object& obj,
                                const std::string& key,
                                double val)
  {
    H5::DataSpace scalar(H5S_SCALAR);
    H5::Attribute attr = obj.createAttribute(
      key, H5::PredType::NATIVE_DOUBLE, scalar);
    attr.write(H5::PredType::NATIVE_DOUBLE, &val);
  }
 
  // Return the UHI storage type string for histogram type T.
  template <typename T>
  std::string storage_type_name() {
    if constexpr (HasSumOfWeightsSquared<T>) {
      return "weighted_mean";
    } else if constexpr (HasCount<T>) {
      return "mean";
    } else if constexpr (HasVariance<T>) {
      return "weighted";
    } else if constexpr (std::is_integral_v<array_t<T>>) {
      return "int";
    } else {
      return "double";
    }
  }
}
 
namespace H5Utils::hist {
 
  // main hist writing function, this should be all anyone needs to
  // call
  template <typename T>
  void write_hist_to_group(
    H5::Group& parent_group,
    const T& hist,
    const std::string& name)
  {
    using namespace detail;
 
    // build the data space
    auto axes = get_axes(hist);
    std::vector<hsize_t> ds_dims;
    for (const auto& axis: axes) {
      int total_size = static_cast<int>(axis.edges.size())
                     - (axis.is_int_or_category ? 0 : 1)
                     + (axis.underflow ? 1 : 0)
                     + (axis.overflow  ? 1 : 0);
      if (total_size < 0) throw std::logic_error("negative bin count");
      ds_dims.push_back(static_cast<hsize_t>(total_size));
    }
    H5::DataSpace space(ds_dims.size(), ds_dims.data());
 
    // Create histogram group and write UHI schema version
    H5::Group group = parent_group.createGroup(name);
    write_int_attr(group, "uhi_schema", 1);
    group.createGroup("metadata");
    group.createGroup("writer_info");
 
    // Dataset creation properties for storage datasets
    H5::DSetCreatPropList props;
    props.setChunk(ds_dims.size(), ds_dims.data());
    props.setDeflate(7);
    chkerr(
      H5Pset_dset_no_attrs_hint(props.getId(), true),
      "no attribute hint");
    H5::DataType type = hdf5_t<array_t<T>>;
 
    // Write ref_axes
    H5::Group ref_axes_group = group.createGroup("ref_axes");
    for (size_t ax_n = 0; ax_n < axes.size(); ax_n++) {
      const auto& ax = axes.at(ax_n);
      H5::Group ax_group = ref_axes_group.createGroup(
        "axis_" + std::to_string(ax_n));
 
      if (ax.n_regular_bins) {
        write_str_attr(ax_group, "type", "regular");
        write_double_attr(ax_group, "lower",
                          static_cast<double>(ax.edges.front()));
        write_double_attr(ax_group, "upper",
                          static_cast<double>(ax.edges.back()));
        write_int_attr(ax_group, "bins", *ax.n_regular_bins);
      } else {
        write_str_attr(ax_group, "type", "variable");
        std::array<hsize_t,1> edgedims{ax.edges.size()};
        H5::DataSpace edge_space(1, edgedims.data());
        H5::DSetCreatPropList edge_props;
        edge_props.setChunk(1, edgedims.data());
        edge_props.setDeflate(7);
        H5::DataSet edge_ds = ax_group.createDataSet(
          "edges", hdf5_t<float>, edge_space, edge_props);
        edge_ds.write(ax.edges.data(), hdf5_t<float>);
      }
 
      write_bool_attr(ax_group, "underflow", ax.underflow);
      write_bool_attr(ax_group, "overflow", ax.overflow);
      write_bool_attr(ax_group, "circular", false);
 
      // Write axis name into a metadata subgroup (key matches Python output)
      H5::Group ax_meta = ax_group.createGroup("metadata");
      if (!ax.name.empty()) {
        write_str_attr(ax_meta, "metadata", ax.name);
      }
      ax_group.createGroup("writer_info");
    }
 
    // Write "axes" object-reference dataset (required by the UHI reader)
    {
      std::vector<hobj_ref_t> axis_refs(axes.size());
      for (size_t ax_n = 0; ax_n < axes.size(); ax_n++) {
        std::string ax_path = "ref_axes/axis_" + std::to_string(ax_n);
        chkerr(H5Rcreate(&axis_refs.at(ax_n), group.getId(),
                         ax_path.c_str(), H5R_OBJECT, -1),
               "axis object reference");
      }
      hsize_t ref_dims[1] = {axes.size()};
      H5::DataSpace ref_space(1, ref_dims);
      group.createDataSet("axes", H5::PredType::STD_REF_OBJ, ref_space)
           .write(axis_refs.data(), H5::PredType::STD_REF_OBJ);
    }
 
    // Write storage group
    H5::Group storage = group.createGroup("storage");
    write_str_attr(storage, "type", storage_type_name<T>());
 
    // Helper to write a named dataset into the storage group
    auto add_dataset = [&hist, &type, &space, &props, &storage](
      const std::string& ds_name,
      const auto& func) {
      auto flat = get_flat_hist_array<T>(hist, func);
      assert(std::cmp_equal(space.getSelectNpoints(),flat.size()));
      storage.createDataSet(ds_name, type, space, props)
             .write(flat.data(), type);
    };
 
    // values (the primary bin contents)
    if constexpr (HasValue<T>) {
      add_dataset("values",
        [](const auto& x) { return x->value(); });
    } else {
      add_dataset("values",
        [](const auto& x) { return *x; });
    }
    // variances
    if constexpr (HasVariance<T>) {
      add_dataset("variances",
        [](const auto& x) { return x->variance(); });
    }
    // count (mean / weighted_mean accumulators)
    if constexpr (HasCount<T>) {
      add_dataset("count",
        [](const auto& x) { return x->count(); });
    }
    // sum_of_weights (weighted_mean accumulator)
    if constexpr (HasSumOfWeights<T>) {
      add_dataset("sum_of_weights",
        [](const auto& x) { return x->sum_of_weights(); });
    }
    // sum_of_weights_squared (weighted_mean accumulator)
    if constexpr (HasSumOfWeightsSquared<T>) {
      add_dataset("sum_of_weights_squared",
        [](const auto& x) { return x->sum_of_weights_squared(); });
    }
 
  }
}