diff --git a/Framework/Utilities/SaveBoostHistogram.hpp b/Framework/Utilities/SaveBoostHistogram.hpp
index ea02dcaba8ef02530554de90c041b6898da98479..8418457f0ab8c9fe86412a77352cca278e919f35 100644
--- a/Framework/Utilities/SaveBoostHistogram.hpp
+++ b/Framework/Utilities/SaveBoostHistogram.hpp
@@ -23,30 +23,34 @@ namespace corsika {
void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
std::string const& filename) {
int const rank = h.rank();
- int const size = h.size();
- using value_type = typename boost::histogram::histogram<Axes, Storage>::value_type;
std::vector<size_t> axes_dims;
axes_dims.reserve(rank);
+ auto overflow = std::make_unique<bool[]>(rank);
+ auto underflow = std::make_unique<bool[]>(rank);
+
std::vector<char> axis_types;
axis_types.reserve(rank);
for (int i = 0; i < rank; ++i) {
auto const& ax = h.axis(i);
- int const has_underflow = (ax.options() & 0x01) ? 1 : 0;
- int const has_overflow = (ax.options() & 0x02) ? 1 : 0;
+ unsigned const has_underflow = (ax.options() & 0x01) ? 1 : 0;
+ unsigned const has_overflow = (ax.options() & 0x02) ? 1 : 0;
+
+ underflow[i] = has_underflow;
+ overflow[i] = has_overflow;
axes_dims.emplace_back(ax.size() + has_underflow + has_overflow);
if (ax.continuous()) {
axis_types.push_back('c');
std::vector<double> ax_edges;
- ax_edges.reserve(ax.size());
+ ax_edges.reserve(ax.size() + 1);
+
+ for (int j = 0; j < ax.size(); ++j) { ax_edges.push_back(ax.bin(j).lower()); }
- for (int j = 0; j < ax.size() + has_overflow; ++j) {
- ax_edges.push_back(ax.bin(j).lower());
- }
+ ax_edges.push_back(ax.bin(ax.size() - 1).upper());
cnpy::npz_save(filename, std::string{"binedges_"} + std::to_string(i),
ax_edges.data(), {ax_edges.size()}, "a");
@@ -63,29 +67,30 @@ namespace corsika {
cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(),
{axis_types.size()}, "a");
+
+ cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(),
+ {axis_types.size()}, "a");
+ cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(),
+ {axis_types.size()}, "a");
}
- auto const prod_axis_size =
- std::accumulate(axes_dims.cbegin(), axes_dims.cend(), 1, std::multiplies<>());
- auto temp = std::make_unique<value_type[]>(prod_axis_size);
+ auto const prod_axis_size = std::accumulate(axes_dims.cbegin(), axes_dims.cend(),
+ unsigned{1}, std::multiplies<>());
+ auto temp = std::make_unique<float[]>(prod_axis_size);
- std::cout << "rank = " << rank << std::endl;
- std::cout << "size = " << size << std::endl;
- std::cout << "prod_axis_size = " << prod_axis_size << std::endl;
+ assert(prod_axis_size == h.size());
// reduce multi-dim. to 1-dim, row-major (i.e., last axis index is contiguous in
- // memory) take special care of underflow bins, which have -1 as index and thus need
+ // memory) take special care of underflow bins which have -1 as index and thus need
// to be shifted by +1
for (auto&& x : indexed(h, boost::histogram::coverage::all)) {
- int const offset_underflow = (h.axis(0).options() & 0x01) ? 1 : 0;
- auto p = x.index(0) + offset_underflow; // 1-d-index
+ int p = 0;
- for (int axis_index = 1; axis_index < rank; ++axis_index) {
+ for (int axis_index = 0; axis_index < rank; ++axis_index) {
int const offset_underflow = (h.axis(axis_index).options() & 0x01) ? 1 : 0;
auto k = x.index(axis_index) + offset_underflow;
p = k + p * axes_dims.at(axis_index);
- std::cout << " " << x.index(axis_index);
}
temp[p] = *x;
diff --git a/Framework/Utilities/testSaveBoostHistogram.cc b/Framework/Utilities/testSaveBoostHistogram.cc
index b00777480d769ec330e6cc528ae7fe211ff17d68..aa05a1947a9c989724dee5244f41859fa634f59f 100644
--- a/Framework/Utilities/testSaveBoostHistogram.cc
+++ b/Framework/Utilities/testSaveBoostHistogram.cc
@@ -11,23 +11,28 @@
#include <random>
+namespace bh = boost::histogram;
+
TEST_CASE("SaveHistogram") {
std::mt19937 rng;
std::normal_distribution n1{2., 6.};
std::exponential_distribution e{1.};
std::uniform_int_distribution u{1, 10};
+ std::uniform_real_distribution<double> r{-3, 3};
- auto h =
- boost::histogram::make_histogram(boost::histogram::axis::regular{5, 0, 10, "x"},
- boost::histogram::axis::regular{3, 0, 4, "y"},
- boost::histogram::axis::category{1, 4, 8});
+ auto h = bh::make_histogram(
+ bh::axis::regular{5, 0, 10, "normal"}, bh::axis::regular{3, 0, 4, "exponential"},
+ bh::axis::category<int>{{2, 3, 5, 7}, "integer category"},
+ bh::axis::regular<double, bh::use_default, bh::use_default,
+ bh::axis::option::growth_t>{10, -1, 1, "integer category"});
for (int i{0}; i < 100'000; ++i) {
auto const a = n1(rng);
auto const b = e(rng);
auto const c = u(rng);
+ auto const d = r(rng);
- h(a, b, c);
+ h(a, b, c, d);
}
corsika::utl::save_hist(h, "hist.npz");
diff --git a/Tools/read_hist.py b/Tools/read_hist.py
index 1d3a69490ba6ccea7dde99fc5d3075a35582ed60..4348e118bdcb4fe4bc3e11fcd50b626425238788 100755
--- a/Tools/read_hist.py
+++ b/Tools/read_hist.py
@@ -11,14 +11,16 @@ def read_hist(filename):
"""
d = np.load(filename)
- axistypes = d['axistypes'].view('c')
+ axistypes = d['axistypes'].view('c')
+ overflow = d['overflow']
+ underflow = d['underflow']
axes = []
- for i, at in enumerate(axistypes):
+ for i, (at, has_overflow, has_underflow) in enumerate(zip(axistypes, overflow, underflow)):
if at == b'c':
- axes.append(bh.axis.Variable(d[f'binedges_{i}'], overflow=True, underflow=True))
+ axes.append(bh.axis.Variable(d[f'binedges_{i}'], overflow=has_overflow, underflow=has_underflow))
elif at == b'd':
- axes.append(bh.axis.IntCategory(d[f'binedges_{i}']))
+ axes.append(bh.axis.IntCategory(d[f'bins_{i}'], growth=(not has_overflow)))
h = bh.Histogram(*axes)
h.view(flow=True)[:] = d['data']