[mesh] Make the dlr_imtime mesh hash blind to the particle statistic

c++/triqs/mesh/dlr_imtime.hpp

@@ -71,6 +71,11 @@ namespace triqs::mesh {
    * time points \f$ \tau_l \f$, i.e. \f$ f_l = f(\tau_l) \f$. In contrast to DLR and imaginary-time meshes, the GF
    * container cannot evaluate the function at an arbitrary imaginary time \f$ \tau \in [0, \beta] \f$
    * (evaluation at arbitrary points is intentionally unsupported).
+   *
+   * The imaginary-time node layout does not depend on the particle statistic, so the mesh hash ignores it:
+   * mesh points of a fermionic and a bosonic mesh with otherwise identical parameters are interchangeable,
+   * e.g. to accumulate a bosonic susceptibility while iterating a fermionic mesh. Whole-mesh equality still
+   * compares the statistic.
    */
   class C2PY_RENAME(MeshDLRImTime) dlr_imtime {
     public:
@@ -145,7 +150,7 @@ namespace triqs::mesh {
          w_max_(w_max),
          eps_(eps),
          symmetrize_(symmetrize),
-         mesh_hash_(hash(beta, statistic, w_max, eps, symmetrize, hash_bytes(ops.imf.get_ifnodes()), std::string_view{"dlr_imtime"})),
+         mesh_hash_(hash(beta, w_max, eps, symmetrize, hash_bytes(ops.imt.get_itnodes_idx()), std::string_view{"dlr_imtime"})),
          dlr_{std::make_shared<detail::dlr_ops>(std::move(ops))} {}
 
     public:
@@ -182,10 +187,10 @@ namespace triqs::mesh {
          w_max_(m.w_max_),
          eps_(m.eps_),
          symmetrize_(m.symmetrize_),
-         mesh_hash_(hash(beta_, stat_, w_max_, eps_, symmetrize_, hash_bytes(m.dlr_->imf.get_ifnodes()), std::string_view{"dlr_imtime"})),
+         mesh_hash_(hash(beta_, w_max_, eps_, symmetrize_, hash_bytes(m.dlr_->imt.get_itnodes_idx()), std::string_view{"dlr_imtime"})),
          dlr_(m.dlr_) {}
 
-    /// Equal-to comparison operator compares the hash values
+    /// Equal-to comparison operator compares the hash value and the particle statistic
     bool operator==(dlr_imtime const &m) const { return mesh_hash_ == m.mesh_hash_ and stat_ == m.stat_; }
 
     /**
@@ -276,7 +281,7 @@ namespace triqs::mesh {
     /// Get the Matsubara frequency DLR operations object (see also `cppdlr::imfreq_ops`).
     [[nodiscard]] C2PY_IGNORE auto const &dlr_if() const { return dlr_->imf; }
 
-    /// Get the hash value of the mesh.
+    /// Get the hash value of the mesh. It identifies the node layout only and ignores the particle statistic.
     [[nodiscard]] C2PY_PROPERTY_GET(mesh_hash) uint64_t mesh_hash() const noexcept { return mesh_hash_; }
 
     /// Get the size \f$ N \f$ of the mesh, i.e. the DLR rank \f$ r \f$.

c++/triqs/mesh/imtime.hpp

@@ -64,6 +64,11 @@ namespace triqs::mesh {
    * Green's function containers that are based on an imaginary time mesh store the function values at the discrete time
    * points \f$ \tau(n) \f$, i.e. \f$ f_n = f(\tau(n)) \f$, and use linear interpolation to evaluate the function at an
    * arbitrary imaginary time \f$ \tau \in [0, \beta] \f$.
+   *
+   * The point layout depends only on \f$ N \f$ and \f$ \beta \f$, not on the particle statistic, so the mesh
+   * hash ignores it: mesh points of a fermionic and a bosonic mesh with matching \f$ N \f$ and \f$ \beta \f$
+   * are interchangeable, e.g. to accumulate a bosonic susceptibility while iterating a fermionic mesh.
+   * Whole-mesh equality still compares the statistic. triqs::mesh::dlr_imtime follows the same contract.
    */
   class C2PY_RENAME(MeshImTime) imtime : public detail::linear<imtime, double> {
     public:

deps/CMakeLists.txt

@@ -130,6 +130,6 @@ external_dependency(nda
 # -- cppdlr --
 external_dependency(cppdlr
   GIT_REPO https://github.com/flatironinstitute/cppdlr
-  VERSION 1.1
+  VERSION 1.4
   GIT_TAG main
 )

test/c++/gfs/gf_bubble_time.cpp

@@ -51,6 +51,16 @@ TEST(Gf, LocalChi2) {
   EXPECT_NEAR(real(chi(beta)), std::exp(-2 * beta), precision);
 }
 
+// A particle-particle bubble: the imtime mesh hash is statistic-blind, so the bosonic mesh points may index the fermionic g.
+TEST(Gf, LocalChiPP) {
+
+  for (auto tau : g.mesh()) g[tau] = std::exp(-2 * tau);
+
+  for (auto tau : chi.mesh()) chi[tau] = g[tau] * g[tau];
+
+  for (auto tau : chi.mesh()) EXPECT_NEAR(real(chi[tau]), std::exp(-4 * tau), precision);
+}
+
 TEST(Gf, GtauEvalBoundary) {
 
   for (auto tau : g.mesh()) g[tau] = std::exp(-2 * tau);

test/c++/gfs/gf_dlr.cpp

@@ -339,6 +339,35 @@ TEST(Gf, DLR_imtime_interpolation) {
 
 // ----------------------------------------------------------------
 
+// A particle-particle bubble: the dlr_imtime mesh hash is statistic-blind, so the bosonic mesh points may index the fermionic g.
+TEST(Gf, dlr_pp_bubble) {
+
+  double beta     = 2.0;
+  double w_max    = 5.0;
+  double eps      = 1e-12;
+  bool symmetrize = true;
+
+  double e = 1.5; // 2 * e < w_max, so the squared propagator is representable in the basis
+
+  auto g   = gf<dlr_imtime, scalar_valued>{{beta, Fermion, w_max, eps, symmetrize}};
+  auto chi = gf<dlr_imtime, scalar_valued>{{beta, Boson, w_max, eps, symmetrize}};
+
+  EXPECT_EQ(g.mesh().mesh_hash(), chi.mesh().mesh_hash());
+
+  for (auto tau : g.mesh()) g[tau] = onefermion(tau, e, beta);
+  for (auto tau : chi.mesh()) chi[tau] = g[tau] * g[tau];
+
+  // evaluate(dlr_imtime, ...) is deleted, so compare through the coefficient mesh
+  auto chi_c = make_gf_dlr(chi);
+  for (auto x : {0.1, 0.3, 0.5, 0.7, 0.9}) {
+    double tau = x * beta;
+    auto g_tau = onefermion(tau, e, beta);
+    EXPECT_COMPLEX_NEAR(chi_c(tau), g_tau * g_tau, 1e-9);
+  }
+}
+
+// ----------------------------------------------------------------
+
 TEST(Gf, DLR_imfreq_interpolation) {
 
   double beta  = 2.0;

test/c++/meshes/mesh_dlr.cpp

@@ -128,4 +128,23 @@ TEST(TRIQSMesh, DLRConversions) {
   EXPECT_EQ(m_dlr_imtime1, m_dlr_imtime4);
 }
 
+// The tau nodes are statistic-independent: fermionic and bosonic dlr_imtime meshes share the hash
+// and the mesh points, but compare unequal.
+TEST(TRIQSMesh, DLRImtimeStatisticBlindHash) {
+  auto m_f = triqs::mesh::dlr_imtime{10, triqs::mesh::Fermion, 1, 1e-6};
+  auto m_b = triqs::mesh::dlr_imtime{10, triqs::mesh::Boson, 1, 1e-6};
+
+  EXPECT_EQ(m_f.mesh_hash(), m_b.mesh_hash());
+  EXPECT_NE(m_f, m_b);
+
+  EXPECT_EQ(m_f.size(), m_b.size());
+  for (auto [mp_f, mp_b] : itertools::zip(m_f, m_b)) {
+    EXPECT_EQ(mp_f.index(), mp_b.index());
+    EXPECT_EQ(mp_f.value(), mp_b.value());
+  }
+
+  // The imaginary frequency grids genuinely differ between statistics
+  EXPECT_NE(triqs::mesh::dlr_imfreq(10, triqs::mesh::Fermion, 1, 1e-6), triqs::mesh::dlr_imfreq(10, triqs::mesh::Boson, 1, 1e-6));
+}
+
 MAKE_MAIN;

test/c++/meshes/mesh_dlr_stability.cpp

@@ -23,7 +23,8 @@
  *
  * For each parameter set we rebuild the meshes and compare their characterizing quantities: the DLR
  * frequencies and imaginary-time nodes (floating point, compared with a tight relative tolerance, as
- * compilers/BLAS backends differ in the last ulps) and the integer Matsubara indices (compared exactly).
+ * compilers/BLAS backends differ in the last ulps) and the integer Matsubara and imaginary-time
+ * fine-grid indices (compared exactly).
  * The mesh hashes are intentionally not pinned: they depend on the standard-library hashing
  * implementation and so are not portable across toolchains. Failures report the parameter set via
  * SCOPED_TRACE.
@@ -100,6 +101,7 @@ namespace {
     h5::write(g, "dlr_freq", m_imfreq.dlr_freq());
     h5::write(g, "ifnodes", m_imfreq.dlr_if().get_ifnodes());
     h5::write(g, "itnodes", m_imtime.dlr_it().get_itnodes());
+    h5::write(g, "itnodes_idx", m_imtime.dlr_it().get_itnodes_idx());
   }
 
 } // namespace
@@ -136,6 +138,8 @@ TEST(MeshDLRRef, CompareAgainstReference) {
     // Matsubara index nodes (dlr_imfreq, exact) and imaginary-time nodes (dlr_imtime, up to rounding).
     EXPECT_EQ_ARRAY(m_imfreq.dlr_if().get_ifnodes(), h5::read<nda::vector<int>>(g, "ifnodes"));
     EXPECT_ARRAY_REL_NEAR(m_imtime.dlr_it().get_itnodes(), h5::read<nda::vector<double>>(g, "itnodes"), node_rtol);
+    // Fine-grid indices of the imaginary-time nodes (exact).
+    EXPECT_EQ_ARRAY(m_imtime.dlr_it().get_itnodes_idx(), h5::read<nda::vector<int>>(g, "itnodes_idx"));
     ++i;
   }
 }

test/python/base/gf_dlr.py

@@ -230,6 +230,26 @@ def test_dlr_symmetrized(self):
 
         assert_gfs_are_close(make_gf_dlr_imtime(make_gf_dlr(giw)), gtau)
 
+    def test_dlr_pp_bubble(self):
+        # The dlr_imtime mesh hash is statistic-blind, so the bosonic mesh points may index the fermionic g.
+        beta, eps, w_max = 2.0, 1e-12, 5.
+        e = 1.5 # 2 * e < w_max, so the squared propagator is representable in the basis
+
+        g   = Gf(mesh = MeshDLRImTime(beta, 'Fermion', w_max, eps), target_shape = [])
+        chi = Gf(mesh = MeshDLRImTime(beta, 'Boson', w_max, eps), target_shape = [])
+
+        assert g.mesh.mesh_hash == chi.mesh.mesh_hash
+
+        for tau in g.mesh:
+            g[tau] = onefermion(tau, e, beta)
+        for tau in chi.mesh:
+            chi[tau] = g[tau] * g[tau]
+
+        chi_c = make_gf_dlr(chi)
+        taus = np.array([0.1, 0.3, 0.5, 0.7, 0.9]) * beta
+        ref = [onefermion(tau, e, beta)**2 for tau in taus]
+        assert np.allclose([chi_c(tau) for tau in taus], ref, atol = 1e-9)
+
     def test_dlr_l2_norm(self):
 
         beta, eps, w_max = 40.1, 1e-10, 5.