Merge branch '341-remove-save_hist-append-write-argument' into 'master'

Resolve "Remove save_hist append/write argument"

Closes #341

See merge request AirShowerPhysics/corsika!307

corsika/detail/framework/process/InteractionHistogram.inl

@@ -42,14 +42,6 @@ namespace corsika {
     }
   }
 
-  void InteractionHistogram::saveLab(std::string const& filename, SaveMode mode) const {
-    corsika::save_hist(inthist_lab_, filename, mode);
-  }
-
-  void InteractionHistogram::saveCMS(std::string const& filename, SaveMode mode) const {
-    corsika::save_hist(inthist_cms_, filename, mode);
-  }
-
   InteractionHistogram& InteractionHistogram::operator+=(
       InteractionHistogram const& other) {
     inthist_lab_ += other.inthist_lab_;
@@ -64,4 +56,4 @@ namespace corsika {
 
     return other;
   }
-} // namespace corsika
+} // namespace corsika
\ No newline at end of file

corsika/detail/framework/utility/SaveBoostHistogram.inl

@@ -11,6 +11,7 @@
 #include <cnpy.hpp>
 
 #include <boost/histogram.hpp>
+#include <boost/filesystem.hpp> // can be changed to std::filesystem if compiler supports it
 
 #include <functional>
 #include <memory>
@@ -23,11 +24,19 @@ namespace corsika {
 
   template <class Axes, class Storage>
   inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
-                        std::string const& filename, SaveMode mode) {
-    unsigned const rank = h.rank();
+                        std::string const& filename, bool overwrite) {
+    if (boost::filesystem::status(filename).type() !=
+        boost::filesystem::file_type::file_not_found) {
+      if (overwrite) {
+        boost::filesystem::remove(filename);
+      } else {
+        using namespace std::literals;
+        throw std::runtime_error(
+            ("save_hist(): "s + filename + " already exists"s).c_str());
+      }
+    }
 
-    // append vs. overwrite
-    const std::string mode_str = (mode == SaveMode::append ? "a" : "w");
+    unsigned const rank = h.rank();
 
     std::vector<size_t> axes_dims;
     axes_dims.reserve(rank);
@@ -58,7 +67,7 @@ namespace corsika {
         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()}, mode_str);
+                       ax_edges.data(), {ax_edges.size()}, "a");
       } else {
         axis_types.push_back('d');
         std::vector<int64_t> bins; // we assume that discrete axes have integer bins
@@ -67,14 +76,13 @@ namespace corsika {
         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()}, mode_str);
+                       {bins.size()}, "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);
+    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");
 
     auto const prod_axis_size = std::accumulate(axes_dims.cbegin(), axes_dims.cend(),
                                                 unsigned{1}, std::multiplies<>());
@@ -98,9 +106,9 @@ namespace corsika {
       temp[p] = *x;
     }
 
-    cnpy::npz_save(filename, "data", temp.get(), axes_dims, mode_str);
+    cnpy::npz_save(filename, "data", temp.get(), axes_dims, "a");
     // 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']
   } // namespace corsika
 
-} // namespace corsika
+} // namespace corsika
\ No newline at end of file

corsika/framework/process/InteractionHistogram.hpp

@@ -43,9 +43,6 @@ namespace corsika {
     hist_type const& CMSHist() const { return inthist_cms_; }
     hist_type const& labHist() const { return inthist_lab_; }
 
-    void saveLab(std::string const& filename, SaveMode mode = SaveMode::append) const;
-    void saveCMS(std::string const& filename, SaveMode mode = SaveMode::append) const;
-
     InteractionHistogram& operator+=(InteractionHistogram const& other);
     InteractionHistogram operator+(InteractionHistogram other) const;
   };

corsika/framework/utility/SaveBoostHistogram.hpp

@@ -12,8 +12,6 @@
 
 namespace corsika {
 
-  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
@@ -21,10 +19,13 @@ namespace corsika {
    *
    * Note that this function makes a temporary, dense copy of the histogram, which could
    * be an issue for huge sizes (e.g. for high dimensions)
+   *
+   * @param overwrite silently overwrite existing files if true, otherwise throw
+   * runtime_error
    */
   template <class Axes, class Storage>
   inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
-                        std::string const& filename, SaveMode mode = SaveMode::append);
+                        std::string const& filename, bool overwrite = true);
 } // namespace corsika
 
 #include <corsika/detail/framework/utility/SaveBoostHistogram.inl>

examples/em_shower.cpp

@@ -7,6 +7,7 @@
  */
 
 #include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
 #include <corsika/framework/geometry/Plane.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
 #include <corsika/framework/geometry/PhysicalGeometry.hpp>
@@ -180,7 +181,7 @@ int main(int argc, char** argv) {
   em_continuous.reset();
 
   auto const hists = proposalCounted.getHistogram();
-  hists.saveLab("inthist_lab_emShower.npz");
-  hists.saveCMS("inthist_cms_emShower.npz");
+  save_hist(hists.labHist(), "inthist_lab_emShower.npz", true);
+  save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true);
   longprof.save("longprof_emShower.txt");
 }

examples/hybrid_MC.cpp

@@ -13,6 +13,7 @@
 #include <corsika/framework/process/InteractionCounter.hpp>
 /* clang-format on */
 #include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
 #include <corsika/framework/core/Logging.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
@@ -267,11 +268,8 @@ int main(int argc, char** argv) {
   auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
                      urqmdCounted.getHistogram();
 
-  hists.saveLab("inthist_lab.txt");
-  hists.saveCMS("inthist_cms.txt");
-
-  hists.saveLab("inthist_lab.txt");
-  hists.saveCMS("inthist_cms.txt");
+  save_hist(hists.labHist(), "inthist_lab_hybrid.npz", true);
+  save_hist(hists.CMSHist(), "inthist_cms_hybrid.npz", true);
 
   longprof.save("longprof.txt");
 

examples/vertical_EAS.cpp

@@ -15,6 +15,7 @@
 #include <corsika/framework/geometry/Plane.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
 #include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
 #include <corsika/framework/process/SwitchProcessSequence.hpp>
 #include <corsika/framework/process/InteractionCounter.hpp>
@@ -281,7 +282,7 @@ int main(int argc, char** argv) {
   auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
                      urqmdCounted.getHistogram() + proposalCounted.getHistogram();
 
-  hists.saveLab("inthist_lab_verticalEAS.npz");
-  hists.saveCMS("inthist_cms_verticalEAS.npz");
+  save_hist(hists.labHist(), "inthist_lab_verticalEAS.npz", true);
+  save_hist(hists.CMSHist(), "inthist_cms_verticalEAS.npz", true);
   longprof.save("longprof_verticalEAS.txt");
 }

tests/framework/testInteractionCounter.cpp

@@ -21,6 +21,11 @@
 #include <catch2/catch.hpp>
 
 #include <numeric>
+#include <algorithm>
+#include <iterator>
+#include <string>
+#include <fstream>
+#include <cstdio>
 
 using namespace corsika;
 
@@ -69,10 +74,31 @@ TEST_CASE("InteractionCounter", "[process]") {
     CHECK(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
     CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
 
-    countedProcess.getHistogram().saveLab("testInteractionCounter_file1.npz",
-                                          SaveMode::overwrite);
-    countedProcess.getHistogram().saveCMS("testInteractionCounter_file2.npz",
-                                          SaveMode::overwrite);
+    save_hist(countedProcess.getHistogram().labHist(), "testInteractionCounter_file1.npz",
+              true);
+    save_hist(countedProcess.getHistogram().CMSHist(), "testInteractionCounter_file2.npz",
+              true);
+
+    SECTION("output validation") {
+      auto const file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
+                                 "testInteractionCounter_file2");
+
+      std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;
+
+      // compare to binary reference data
+      // note that this currenly compares the whole files byte by byte. If the new
+      // or the reference file are compressed this would be a false-negative outcome
+      // of this test
+      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;
+
+      CHECK(std::equal(begin1, end, begin1ref));
+    }
   }
 
   SECTION("DoInteraction Lambda") {
@@ -93,41 +119,3 @@ TEST_CASE("InteractionCounter", "[process]") {
     CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
   }
 }
-
-#include <algorithm>
-#include <iterator>
-#include <string>
-#include <fstream>
-
-TEST_CASE("InteractionCounterOutput", "[output validation]") {
-
-  logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
-  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();
-  }
-}

tests/framework/testSaveBoostHistogram.cpp

@@ -40,5 +40,6 @@ TEST_CASE("SaveHistogram") {
     h(a, b, c, d);
   }
 
-  corsika::save_hist(h, "hist.npz");
+  REQUIRE_NOTHROW(corsika::save_hist(h, "hist.npz", true));
+  REQUIRE_THROWS(corsika::save_hist(h, "hist.npz", false));
 }

tests/modules/testInteractionCounter.cpp

-/*
- * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/modules/InteractionCounter.hpp>
-
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/HomogeneousMedium.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-#include <corsika/setup/SetupStack.hpp>
-
-#include <catch2/catch.hpp>
-
-#include <numeric>
-
-using namespace corsika;
-
-const std::string refDataDir = std::string(REFDATADIR); // from cmake
-
-struct DummyProcess {
-  template <typename TParticle>
-  GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) {
-    return 100_g / 1_cm / 1_cm;
-  }
-
-  template <typename TParticle>
-  auto DoInteraction([[maybe_unused]] TParticle& projectile) {
-    return nullptr;
-  }
-};
-
-TEST_CASE("InteractionCounter", "[process]") {
-
-  logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
-  DummyProcess d;
-  InteractionCounter countedProcess(d);
-
-  SECTION("GetInteractionLength") {
-    CHECK(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
-  }
-
-  auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
-  [[maybe_unused]] auto& env_dummy = env;
-
-  SECTION("DoInteraction nucleus") {
-    unsigned short constexpr A = 14, Z = 7;
-    auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
-                                                             Z, 105_TeV, nodePtr, *csPtr);
-    CHECK(stackPtr->getEntries() == 1);
-    CHECK(secViewPtr->getEntries() == 0);
-
-    auto const ret = countedProcess.DoInteraction(*secViewPtr);
-    CHECK(ret == nullptr);
-
-    auto const& h = countedProcess.GetHistogram().labHist();
-    CHECK(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
-    CHECK(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
-
-    auto const& h2 = countedProcess.GetHistogram().CMSHist();
-    CHECK(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
-    CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
-
-    countedProcess.GetHistogram().saveLab("testInteractionCounter_file1.npz",
-                                          utl::SaveMode::overwrite);
-    countedProcess.GetHistogram().saveCMS("testInteractionCounter_file2.npz",
-                                          utl::SaveMode::overwrite);
-  }
-
-  SECTION("DoInteraction Lambda") {
-    auto constexpr code = particles::Code::Lambda0;
-    auto [stackPtr, secViewPtr] =
-        setup::testing::setupStack(code, 0, 0, 105_TeV, nodePtr, *csPtr);
-    CHECK(stackPtr->getEntries() == 1);
-    CHECK(secViewPtr->getEntries() == 0);
-
-    auto const ret = countedProcess.DoInteraction(*secViewPtr);
-    CHECK(ret == nullptr);
-
-    auto const& h = countedProcess.GetHistogram().labHist();
-    CHECK(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);
-    CHECK(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
-
-    auto const& h2 = countedProcess.GetHistogram().CMSHist();
-    CHECK(h2.at(h2.axis(0).index(3122), h2.axis(1).index(1.6e12)) == 1);
-    CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
-  }
-}
-
-#include <algorithm>
-#include <iterator>
-#include <string>
-#include <fstream>
-
-TEST_CASE("InteractionCounterOutput", "[output validation]") {
-
-  logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
-  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();
-  }
-}

tools/read_hist.py

+import numpy as np
+import matplotlib.pyplot as plt
+import boost_histogram as bh
+
+def read_hist(filename):
+    """
+    read numpy file produced with CORSIKA 8's save_hist() function into
+    boost-histogram object.
+    """
+
+    d = np.load(filename)
+    axistypes = d['axistypes'].view('c')
+    overflow = d['overflow']
+    underflow = d['underflow']
+
+    axes = []
+    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=has_overflow, underflow=has_underflow))
+        elif at == b'd':
+            axes.append(bh.axis.IntCategory(d[f'bins_{i}'], growth=(not has_overflow)))
+ 
+    h = bh.Histogram(*axes)
+    h.view(flow=True)[:] = d['data']
+    
+    return h
+