Store and expose the DLR node fine grid indices and the symmetrize flag

CMakeLists.txt

@@ -24,7 +24,7 @@ cmake_policy(VERSION ${CMAKE_VERSION})
 
 # ############
 # Define Project
-project(cppdlr VERSION 1.3.0 LANGUAGES CXX)
+project(cppdlr VERSION 1.4.0 LANGUAGES CXX)
 get_directory_property(IS_SUBPROJECT PARENT_DIRECTORY)
 
 message(STATUS "-------- cppdlr version and git hash detection -------------")

c++/cppdlr/dlr_build.cpp

@@ -313,4 +313,30 @@ namespace cppdlr {
 
   nda::vector<double> build_dlr_rf(double lambda, double eps) { return build_dlr_rf(lambda, eps, NONSYM); }
 
+  nda::vector<int> recover_itnode_idx(double lambda, nda::vector_const_view<double> dlr_it) {
+
+    int r      = dlr_it.size();
+    auto t     = std::get<0>(build_it_fine(fineparams(lambda), NONSYM));
+    auto igrid = nda::range(t.size());
+    auto idx   = nda::vector<int>(r);
+
+    // Nearest-neighbor match: stored nodes may differ from the grid in the last few digits
+    for (int l = 0; l < r; ++l) {
+      auto dist = [&](long i) { return std::abs(t(i) - dlr_it(l)); };
+      idx(l)    = *std::ranges::min_element(igrid, {}, dist);
+      if (dist(idx(l)) > 1e-12) throw std::runtime_error("Could not match the DLR imaginary time nodes to the fine grid for the given lambda.");
+    }
+
+    return idx;
+  }
+
+  void check_unsymmetrized(nda::vector_const_view<double> dlr_rf) {
+
+    int r            = dlr_rf.size();
+    auto is_mirrored = [&](long i) { return dlr_rf(i) == -dlr_rf(r - 1 - i); };
+
+    if (std::ranges::all_of(range(r), is_mirrored))
+      throw std::runtime_error("Mirror-symmetric DLR frequency grid: symmetrized grids from cppdlr <= 1.3.0 are no longer supported.");
+  }
+
 } // namespace cppdlr

c++/cppdlr/dlr_build.hpp

@@ -173,4 +173,28 @@ namespace cppdlr {
   */
   nda::vector<double> build_dlr_rf(double lambda, double eps);
 
+  /**
+  * @brief Recover the indices of the DLR imaginary time nodes in the unsymmetrized
+  * fine imaginary time grid
+  *
+  * @param[in] lambda DLR cutoff parameter
+  * @param[in] dlr_it DLR imaginary time nodes
+  *
+  * @return Indices of the DLR imaginary time nodes in the fine grid
+  *
+  * @note Nodes are matched to their nearest fine grid point; throws if a node is
+  * further than 1e-12 from the grid.
+  */
+  nda::vector<int> recover_itnode_idx(double lambda, nda::vector_const_view<double> dlr_it);
+
+  /**
+  * @brief Throw if the given DLR frequency grid is mirror-symmetric about omega=0
+  *
+  * Guards the code paths predating the stored symmetrize flag against symmetrized
+  * grids, including the pair-only grids of cppdlr <= 1.3.0.
+  *
+  * @param[in] dlr_rf DLR frequencies
+  */
+  void check_unsymmetrized(nda::vector_const_view<double> dlr_rf);
+
 } // namespace cppdlr

c++/cppdlr/dlr_imfreq.cpp

@@ -23,7 +23,7 @@ using namespace nda;
 namespace cppdlr {
 
   imfreq_ops::imfreq_ops(double lambda, nda::vector_const_view<double> dlr_rf, statistic_t statistic, bool symmetrize)
-     : lambda_(lambda), statistic(statistic), r(dlr_rf.size()), dlr_rf(dlr_rf) {
+     : lambda_(lambda), symmetrize_(symmetrize), statistic(statistic), r(dlr_rf.size()), dlr_rf(dlr_rf) {
 
     // # DLR imaginary frequency nodes. The symmetrized grid is mirror-symmetric
     // about i*nu=0. The bosonic grid contains the self-paired node n=0, so niom = r,

c++/cppdlr/dlr_imfreq.hpp

@@ -57,8 +57,10 @@ namespace cppdlr {
                nda::vector_const_view<int> dlr_if,                                          //
                nda::matrix_const_view<nda::dcomplex> cf2if,                                 //
                nda::matrix_const_view<nda::dcomplex> if2cf_lu,                              //
-               nda::vector_const_view<int> if2cf_piv)
+               nda::vector_const_view<int> if2cf_piv,                                       //
+               bool symmetrize)
        : lambda_(lambda),
+         symmetrize_(symmetrize),
          statistic(statistic),
          r(cf2if.extent(1)),
          niom(dlr_if.size()),
@@ -67,6 +69,15 @@ namespace cppdlr {
          cf2if(cf2if),
          if2cf{if2cf_lu, if2cf_piv} {};
 
+    imfreq_ops(double lambda, nda::vector_const_view<double> dlr_rf, statistic_t statistic, //
+               nda::vector_const_view<int> dlr_if,                                          //
+               nda::matrix_const_view<nda::dcomplex> cf2if,                                 //
+               nda::matrix_const_view<nda::dcomplex> if2cf_lu,                              //
+               nda::vector_const_view<int> if2cf_piv)
+       : imfreq_ops(lambda, dlr_rf, statistic, dlr_if, cf2if, if2cf_lu, if2cf_piv, NONSYM) {
+      check_unsymmetrized(dlr_rf); // this signature predates the symmetrize flag
+    };
+
     imfreq_ops() = default;
 
     /** 
@@ -242,8 +253,12 @@ namespace cppdlr {
     double lambda() const { return lambda_; }
     statistic_t get_statistic() const { return statistic; }
 
+    /** Whether this object was built with symmetrized DLR grids, SYM/true or NONSYM/false */
+    bool is_symmetrized() const { return symmetrize_; }
+
     private:
     double lambda_;                   ///< Energy cutoff divided by temperature
+    bool symmetrize_ = false;         ///< Whether the DLR grids are symmetrized
     statistic_t statistic;            ///< Particle statistic: Fermion or Boson
     int r;                            ///< DLR rank
     int niom;                         ///< # DLR imaginary freq nodes (= r + 1 in the symmetrized fermionic case, else = r)
@@ -268,7 +283,7 @@ namespace cppdlr {
      *
      * @param[in] ar Archive to serialize into
      */
-    void serialize(auto &ar) const { ar & lambda_ & statistic & r & niom & dlr_rf & dlr_if & cf2if & if2cf.lu & if2cf.piv; }
+    void serialize(auto &ar) const { ar & lambda_ & symmetrize_ & statistic & r & niom & dlr_rf & dlr_if & cf2if & if2cf.lu & if2cf.piv; }
 
     /**
      * Deserialize an object from the archive. This will initialize all members.
@@ -277,7 +292,7 @@ namespace cppdlr {
      *
      * @param[in] ar Archive to deserialize from
      */
-    void deserialize(auto &ar) { ar & lambda_ & statistic & r & niom & dlr_rf & dlr_if & cf2if & if2cf.lu & if2cf.piv; }
+    void deserialize(auto &ar) { ar & lambda_ & symmetrize_ & statistic & r & niom & dlr_rf & dlr_if & cf2if & if2cf.lu & if2cf.piv; }
 
     // -------------------- hdf5 -------------------
 
@@ -295,6 +310,7 @@ namespace cppdlr {
       h5::write(gr, "cf2if", m.get_cf2if());
       h5::write(gr, "if2cf_lu", m.get_if2cf_lu());
       h5::write(gr, "if2cf_piv", m.get_if2cf_piv());
+      h5::write(gr, "symmetrize", m.is_symmetrized());
     }
 
     friend void h5_read(h5::group fg, std::string const &subgroup_name, imfreq_ops &m) {
@@ -318,7 +334,11 @@ namespace cppdlr {
       h5::read(gr, "if2cf_lu", if2cf_lu);
       h5::read(gr, "if2cf_piv", if2cf_piv);
 
-      m = imfreq_ops(lambda, rf, statistic_, if_, cf2if_, if2cf_lu, if2cf_piv);
+      // Archives written before the symmetrization option was stored hold an unsymmetrized grid
+      bool symmetrize = NONSYM;
+      if (not h5::try_read(gr, "symmetrize", symmetrize)) check_unsymmetrized(rf);
+
+      m = imfreq_ops(lambda, rf, statistic_, if_, cf2if_, if2cf_lu, if2cf_piv, symmetrize);
     }
   };
 

c++/cppdlr/dlr_imtime.cpp

@@ -25,7 +25,8 @@ using namespace nda;
 
 namespace cppdlr {
 
-  imtime_ops::imtime_ops(double lambda, nda::vector_const_view<double> dlr_rf, bool symmetrize) : lambda_(lambda), r(dlr_rf.size()), dlr_rf(dlr_rf) {
+  imtime_ops::imtime_ops(double lambda, nda::vector_const_view<double> dlr_rf, bool symmetrize)
+     : lambda_(lambda), symmetrize_(symmetrize), r(dlr_rf.size()), dlr_rf(dlr_rf) {
 
     // The symmetrized selection always takes the self-paired tau=beta/2 node plus
     // mirror pairs, so it can only produce an odd number of nodes.
@@ -45,6 +46,7 @@ namespace cppdlr {
     // Pivoted Gram-Schmidt to obtain DLR imaginary time nodes
     auto [q, norms, piv] = (symmetrize ? pivrgs_sym(kmat, 1e-100) : pivrgs(kmat, 1e-100));
     std::sort(piv.begin(), piv.end()); // Sort pivots in ascending order
+    dlr_it_idx = piv;                  // Integer identity of the node set, see get_itnodes_idx
     for (int i = 0; i < r; ++i) { dlr_it(i) = t(piv(i)); }
 
     // Obtain coefficients to imaginary time values transformation matrix

c++/cppdlr/dlr_imtime.hpp

@@ -66,9 +66,23 @@ namespace cppdlr {
     */
     imtime_ops(double lambda, nda::vector_const_view<double> dlr_rf);
 
+    imtime_ops(double lambda, nda::vector_const_view<double> dlr_rf, nda::vector_const_view<double> dlr_it, nda::matrix_const_view<double> cf2it,
+               nda::matrix_const_view<double> it2cf_lu, nda::vector_const_view<int> it2cf_piv, nda::vector_const_view<int> dlr_it_idx,
+               bool symmetrize)
+       : lambda_(lambda),
+         symmetrize_(symmetrize),
+         r(dlr_rf.size()),
+         dlr_rf(dlr_rf),
+         dlr_it(dlr_it),
+         dlr_it_idx(dlr_it_idx),
+         cf2it(cf2it),
+         it2cf{it2cf_lu, it2cf_lu, it2cf_piv} {};
+
     imtime_ops(double lambda, nda::vector_const_view<double> dlr_rf, nda::vector_const_view<double> dlr_it, nda::matrix_const_view<double> cf2it,
                nda::matrix_const_view<double> it2cf_lu, nda::vector_const_view<int> it2cf_piv)
-       : lambda_(lambda), r(dlr_rf.size()), dlr_rf(dlr_rf), dlr_it(dlr_it), cf2it(cf2it), it2cf{it2cf_lu, it2cf_lu, it2cf_piv} {};
+       : imtime_ops(lambda, dlr_rf, dlr_it, cf2it, it2cf_lu, it2cf_piv, recover_itnode_idx(lambda, dlr_it), NONSYM) {
+      check_unsymmetrized(dlr_rf); // this signature predates the symmetrize flag
+    };
 
     imtime_ops() = default;
 
@@ -721,6 +735,19 @@ namespace cppdlr {
     nda::vector_const_view<double> get_itnodes() const { return dlr_it; };
     double get_itnodes(int i) const { return dlr_it(i); };
 
+    /**
+    * @brief Get indices of the DLR imaginary time nodes in the fine imaginary time
+    * grid built with this object's symmetrize option
+    *
+    * @return Indices of the DLR imaginary time nodes in the fine grid
+    *
+    * @note The nodes are exact copies of the fine grid points they index, so the
+    * indices, together with is_symmetrized, identify the node set exactly, e.g. for
+    * hashing or comparing grids across toolchains.
+    */
+    nda::vector_const_view<int> get_itnodes_idx() const { return dlr_it_idx; };
+    int get_itnodes_idx(int i) const { return dlr_it_idx(i); };
+
     /** Access DLR imaginary real frequency nodes*/
     /**
     * @brief Get DLR real frequency nodes
@@ -767,6 +794,9 @@ namespace cppdlr {
     int rank() const { return r; }
     double lambda() const { return lambda_; }
 
+    /** Whether this object was built with symmetrized DLR grids, SYM/true or NONSYM/false */
+    bool is_symmetrized() const { return symmetrize_; }
+
     /**
     * @brief Get inner product matrix
     *
@@ -914,10 +944,12 @@ namespace cppdlr {
 
     private:
     double lambda_;
-    int r;                      ///< DLR rank
-    nda::vector<double> dlr_rf; ///< DLR frequencies
-    nda::vector<double> dlr_it; ///< DLR imaginary time nodes
-    nda::matrix<double> cf2it;  ///< Transformation matrix from DLR coefficients to values at DLR imaginary time nodes
+    bool symmetrize_ = false;    ///< Whether the DLR grids are symmetrized
+    int r;                       ///< DLR rank
+    nda::vector<double> dlr_rf;  ///< DLR frequencies
+    nda::vector<double> dlr_it;  ///< DLR imaginary time nodes
+    nda::vector<int> dlr_it_idx; ///< Indices of the DLR imaginary time nodes in the fine time grid
+    nda::matrix<double> cf2it;   ///< Transformation matrix from DLR coefficients to values at DLR imaginary time nodes
 
     /**
     * @brief Struct for transformation from DLR imaginary time values to coefficients
@@ -952,7 +984,8 @@ namespace cppdlr {
      * @param[in] ar Archive to serialize into
      */
     void serialize(auto &ar) const {
-      ar & lambda_ & r & dlr_rf & dlr_it & cf2it & it2cf.lu & it2cf.zlu & it2cf.piv & hilb & tcf2it & thilb & ttcf2it & ipmat & refl;
+      ar & lambda_ & symmetrize_ & r & dlr_rf & dlr_it & dlr_it_idx & cf2it & it2cf.lu & it2cf.zlu & it2cf.piv & hilb & tcf2it & thilb & ttcf2it
+         & ipmat & refl;
     }
 
     /**
@@ -963,7 +996,8 @@ namespace cppdlr {
      * @param[in] ar Archive to deserialize from
      */
     void deserialize(auto &ar) {
-      ar & lambda_ & r & dlr_rf & dlr_it & cf2it & it2cf.lu & it2cf.zlu & it2cf.piv & hilb & tcf2it & thilb & ttcf2it & ipmat & refl;
+      ar & lambda_ & symmetrize_ & r & dlr_rf & dlr_it & dlr_it_idx & cf2it & it2cf.lu & it2cf.zlu & it2cf.piv & hilb & tcf2it & thilb & ttcf2it
+         & ipmat & refl;
     }
 
     // -------------------- hdf5 -------------------
@@ -981,6 +1015,8 @@ namespace cppdlr {
       h5::write(gr, "cf2it", m.get_cf2it());
       h5::write(gr, "it2cf_lu", m.get_it2cf_lu());
       h5::write(gr, "it2cf_piv", m.get_it2cf_piv());
+      h5::write(gr, "it_idx", m.get_itnodes_idx());
+      h5::write(gr, "symmetrize", m.is_symmetrized());
     }
 
     friend void h5_read(h5::group fg, std::string const &subgroup_name, imtime_ops &m) {
@@ -995,7 +1031,17 @@ namespace cppdlr {
       auto it2cf_lu  = h5::read<nda::matrix<double>>(gr, "it2cf_lu");
       auto it2cf_piv = h5::read<nda::vector<int>>(gr, "it2cf_piv");
 
-      m = imtime_ops(lambda, rf, it, cf2it_, it2cf_lu, it2cf_piv);
+      // Both datasets were introduced in cppdlr 1.4.0; an archive without them holds an unsymmetrized grid
+      bool symmetrize = NONSYM;
+      auto it_idx     = nda::vector<int>{};
+      if (h5::try_read(gr, "symmetrize", symmetrize)) {
+        h5::read(gr, "it_idx", it_idx);
+      } else {
+        check_unsymmetrized(rf);
+        it_idx = recover_itnode_idx(lambda, it);
+      }
+
+      m = imtime_ops(lambda, rf, it, cf2it_, it2cf_lu, it2cf_piv, it_idx, symmetrize);
     }
   };
 

doc/ChangeLog.md

@@ -2,6 +2,20 @@
 
 # Changelog
 
+## Version 1.4.0
+
+### New features
+* Add `imtime_ops::get_itnodes_idx` accessor exposing the indices of the DLR imaginary time nodes in the fine time grid. Being integers, they identify the node set exactly and are therefore suited for fingerprinting the grid across toolchains, e.g. for hashing
+* Add free function `recover_itnode_idx` reconstructing those indices from the nodes alone
+* `imtime_ops` and `imfreq_ops` retain the `symmetrize` option they were constructed with, exposed via the new `is_symmetrized` accessor
+* Add free function `check_unsymmetrized`, rejecting a mirror-symmetric DLR frequency grid
+
+### Other changes
+* `imtime_ops` h5 archives carry the new `it_idx` dataset. Archives written without it are still read, with the indices recovered via `recover_itnode_idx`
+* `imtime_ops` and `imfreq_ops` h5 archives carry the new `symmetrize` dataset. Archives written without it are read as unsymmetrized
+* Symmetrized DLR grids from cppdlr <= 1.3.0 are rejected. Their pair-only selection is mirror-symmetric about `omega=0` without containing it, so its imaginary time nodes lie on today's unsymmetrized fine grid, and reporting their indices would fingerprint the grid against one the `symmetrize` flag does not name. Reading such an archive, or passing such a grid to the constructors predating the flag, now throws via `check_unsymmetrized`
+
+
 ## Version 1.3.0
 
 This update contains additional functionality, bug fixes, and build system improvements.

test/c++/dlr_build.cpp

@@ -101,3 +101,26 @@ TEST(dlr_build, symmetric_fixed_points) {
   EXPECT_EQ(t_non.size(), fine.nt);
   EXPECT_EQ(w_non.size(), fine.nt);
 }
+
+/**
+* @brief Test the unsymmetrized check guarding the paths that predate the stored
+* symmetrize flag
+*/
+TEST(dlr_build, check_unsymmetrized) {
+
+  double eps = 1e-10;
+
+  for (double lambda : {10.0, 100.0, 1000.0, 10000.0}) {
+    EXPECT_NO_THROW(check_unsymmetrized(build_dlr_rf(lambda, eps)));
+    EXPECT_THROW(check_unsymmetrized(build_dlr_rf(lambda, eps, SYM)), std::runtime_error);
+  }
+
+  // Symmetrized grids from cppdlr <= 1.3.0 contain mirror pairs only, hence have even
+  // rank. Mimic one by dropping omega=0 from a symmetrized grid.
+  auto rf = build_dlr_rf(1000.0, eps, SYM);
+  int r   = rf.size();
+  EXPECT_EQ(rf(r / 2), 0.0);
+
+  auto rf_legacy = nda::concatenate(rf(nda::range(0, r / 2)), rf(nda::range(r / 2 + 1, r)));
+  EXPECT_THROW(check_unsymmetrized(rf_legacy), std::runtime_error);
+}

test/c++/imfreq_ops.cpp

@@ -430,6 +430,7 @@ static void check_h5_roundtrip(imfreq_ops const &ifops, std::string const &filen
 
   EXPECT_EQ(ifops.lambda(), ifops_ref.lambda());
   EXPECT_EQ(ifops.rank(), ifops_ref.rank());
+  EXPECT_EQ(ifops.is_symmetrized(), ifops_ref.is_symmetrized());
   EXPECT_EQ_ARRAY(ifops.get_rfnodes(), ifops_ref.get_rfnodes());
   EXPECT_EQ_ARRAY(ifops.get_ifnodes(), ifops_ref.get_ifnodes());
   EXPECT_EQ_ARRAY(ifops.get_cf2if(), ifops_ref.get_cf2if());