even better flow bin handling

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;

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");

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']