better tests

Framework/Utilities/SaveBoostHistogram.hpp

@@ -15,10 +15,14 @@
 #include <numeric>
 #include <utility>
 #include <vector>
+#include <string>
 
 #pragma once
 
 namespace corsika::utl {
+
+  enum class SaveMode { overwrite, append };
+
   /**
    * This functions saves a boost::histogram into a numpy file. Only rather basic axis
    * types are supported: regular, variable, integer, category<int>. Only "ordinary" bin
@@ -29,9 +33,12 @@ namespace corsika::utl {
    */
   template <class Axes, class Storage>
   inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
-                        std::string const& filename) {
+                        std::string const& filename, SaveMode mode = SaveMode::append) {
     unsigned const rank = h.rank();
 
+    // append vs. overwrite
+    const std::string mode_str = (mode == SaveMode::append ? "a" : "w");
+
     std::vector<size_t> axes_dims;
     axes_dims.reserve(rank);
 
@@ -61,7 +68,7 @@ namespace corsika::utl {
         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");
+                       ax_edges.data(), {ax_edges.size()}, mode_str);
       } else {
         axis_types.push_back('d');
         std::vector<int64_t> bins; // we assume that discrete axes have integer bins
@@ -70,13 +77,14 @@ namespace corsika::utl {
         for (int j = 0; j < ax.size(); ++j) { bins.push_back(ax.bin(j).lower()); }
 
         cnpy::npz_save(filename, std::string{"bins_"} + std::to_string(i), bins.data(),
-                       {bins.size()}, "a");
+                       {bins.size()}, mode_str);
       }
     }
 
-    cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(), {rank}, "a");
-    cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(), {rank}, "a");
-    cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(), {rank}, "a");
+    cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(), {rank},
+                   mode_str);
+    cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(), {rank}, mode_str);
+    cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(), {rank}, mode_str);
 
     auto const prod_axis_size = std::accumulate(axes_dims.cbegin(), axes_dims.cend(),
                                                 unsigned{1}, std::multiplies<>());
@@ -100,7 +108,7 @@ namespace corsika::utl {
       temp[p] = *x;
     }
 
-    cnpy::npz_save(filename, "data", temp.get(), axes_dims, "a");
+    cnpy::npz_save(filename, "data", temp.get(), axes_dims, mode_str);
     // In Python this array can directly be assigned to a histogram view if that
     // histogram has its axes correspondingly: hist.view(flow=True)[:] = file['data']
   }

Processes/InteractionCounter/CMakeLists.txt

@@ -50,6 +50,11 @@ CORSIKA_ADD_TEST(testInteractionCounter SOURCES
   ${MODEL_HEADERS}
 )
 
+target_compile_definitions (
+  testInteractionCounter
+  PRIVATE
+  REFDATADIR="${CMAKE_CURRENT_SOURCE_DIR}"
+  )
 target_link_libraries (
   testInteractionCounter
   ProcessInteractionCounter

Processes/InteractionCounter/InteractionHistogram.cc

@@ -41,12 +41,14 @@ void InteractionHistogram::fill(particles::Code projectile_id,
   }
 }
 
-void InteractionHistogram::saveLab(std::string const& filename) const {
-  corsika::utl::save_hist(inthist_lab_, filename);
+void InteractionHistogram::saveLab(std::string const& filename,
+                                   utl::SaveMode mode) const {
+  corsika::utl::save_hist(inthist_lab_, filename, mode);
 }
 
-void InteractionHistogram::saveCMS(std::string const& filename) const {
-  corsika::utl::save_hist(inthist_cms_, filename);
+void InteractionHistogram::saveCMS(std::string const& filename,
+                                   utl::SaveMode mode) const {
+  corsika::utl::save_hist(inthist_cms_, filename, mode);
 }
 
 InteractionHistogram& InteractionHistogram::operator+=(

Processes/InteractionCounter/InteractionHistogram.hpp

@@ -54,9 +54,11 @@ namespace corsika::process::interaction_counter {
 
     hist_type const& labHist() const { return inthist_lab_; }
 
-    void saveLab(std::string const& filename) const;
+    void saveLab(std::string const& filename,
+                 utl::SaveMode mode = utl::SaveMode::append) const;
 
-    void saveCMS(std::string const& filename) const;
+    void saveCMS(std::string const& filename,
+                 utl::SaveMode mode = utl::SaveMode::append) const;
 
     InteractionHistogram& operator+=(InteractionHistogram const& other);
 

Processes/InteractionCounter/testInteractionCounter.cc

@@ -27,6 +27,8 @@ using namespace corsika::process::interaction_counter;
 using namespace corsika::units;
 using namespace corsika::units::si;
 
+const std::string refDataDir = std::string(REFDATADIR);
+
 struct DummyProcess {
   template <typename TParticle>
   GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) {
@@ -39,7 +41,7 @@ struct DummyProcess {
   }
 };
 
-TEST_CASE("InteractionCounter") {
+TEST_CASE("InteractionCounter", "[process]") {
 
   logging::SetLevel(logging::level::debug);
 
@@ -69,15 +71,14 @@ TEST_CASE("InteractionCounter") {
 
     auto const& h2 = countedProcess.GetHistogram().CMSHist();
     REQUIRE(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
-    // REQUIRE(h2.at(1'000'070'140, 92) == 1); // bin 1.584 .. 1.995 TeV √s
     REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
 
-    countedProcess.GetHistogram().saveLab("testInteractionCounter_file1.npz");
-    countedProcess.GetHistogram().saveCMS("testInteractionCounter_file2.npz");
+    countedProcess.GetHistogram().saveLab("testInteractionCounter_file1.npz",
+                                          utl::SaveMode::overwrite);
+    countedProcess.GetHistogram().saveCMS("testInteractionCounter_file2.npz",
+                                          utl::SaveMode::overwrite);
   }
 
-  SECTION("check saving") {}
-
   SECTION("DoInteraction Lambda") {
     auto constexpr code = particles::Code::Lambda0;
     auto [stackPtr, secViewPtr] =
@@ -97,3 +98,38 @@ TEST_CASE("InteractionCounter") {
     REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
   }
 }
+
+#include <algorithm>
+#include <iterator>
+#include <string>
+#include <fstream>
+
+TEST_CASE("InteractionCounterOutput", "[output validation]") {
+
+  auto file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
+                       "testInteractionCounter_file2");
+
+  SECTION(std::string("check saved data, ") + file + ".npz") {
+
+    std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;
+
+    // compare to binary reference data
+    std::ifstream file1(file + ".npz");
+    std::ifstream file1ref(refDataDir + "/" + file + "_REF.npz");
+
+    std::istreambuf_iterator<char> begin1(file1);
+    std::istreambuf_iterator<char> begin1ref(file1ref);
+
+    std::istreambuf_iterator<char> end;
+
+    while (begin1 != end && begin1ref != end) {
+      CHECK(*begin1 == *begin1ref);
+      ++begin1;
+      ++begin1ref;
+    }
+    CHECK(begin1 == end);
+    CHECK(begin1ref == end);
+    file1.close();
+    file1ref.close();
+  }
+}