Print of the average update time

c++/triqs_cthyb/impurity_trace.cpp

@@ -254,7 +254,7 @@ namespace triqs_cthyb {
 
   //-------- Compute the full trace ------------------------------------------
   // Returns MC atomic weight and reweighting = trace/(atomic weight)
-  std::pair<h_scalar_t, h_scalar_t> impurity_trace::compute(double p_yee, double u_yee) {
+  std::pair<h_scalar_t, h_scalar_t> impurity_trace::compute(double p_yee, double u_yee, bool meas_den) {
 
     double epsilon         = 1.e-15; // Machine precision
     auto log_epsilon0      = -std::log(1.e-15);
@@ -264,7 +264,7 @@ namespace triqs_cthyb {
     // simplifies later code
     if (tree_size == 0) {
       if (use_norm_as_weight) {
-        density_matrix = atomic_rho;
+        if (meas_den) density_matrix = atomic_rho;
         return {atomic_norm, atomic_z / atomic_norm};
       } else
         return {atomic_z, 1};
@@ -329,7 +329,7 @@ namespace triqs_cthyb {
     double norm_trace_sq = 0, trace_abs = 0;
 
     // Put density_matrix to "not recomputed"
-    for (int bl = 0; bl < n_blocks; ++bl) density_matrix[bl].is_valid = false;
+    if (meas_den) for (int bl = 0; bl < n_blocks; ++bl) density_matrix[bl].is_valid = false;
 
     auto trace_contrib_block = std::vector<std::pair<double, int>>{}; //FIXME complex -- can histos handle this?
 
@@ -382,21 +382,29 @@ namespace triqs_cthyb {
 
       if (use_norm_as_weight) { // else we are not allowed to compute this matrix, may make no sense
         // recompute the density matrix
-        density_matrix[block_index].is_valid = true;
         double norm_trace_sq_partial         = 0;
-        auto &mat                            = density_matrix[block_index].mat;
+        matrix_t M = {};
+        matrix_t *mat;
+        if (meas_den) {
+          mat = &density_matrix[block_index].mat;
+          density_matrix[block_index].is_valid = true;
+        }
+        else {
+          M = matrix_t(dim,dim);
+          mat = &M;
+        }
         for (int u = 0; u < dim; ++u) {
           for (int v = 0; v < dim; ++v) {
-            mat(u, v) = b_mat.second(u, v) * std::exp(-dtau_beta * get_block_eigenval(block_index, u) - dtau_0 * get_block_eigenval(block_index, v));
-            double xx = std::abs(mat(u, v));
+            (*mat)(u, v) = b_mat.second(u, v) * std::exp(-dtau_beta * get_block_eigenval(block_index, u) - dtau_0 * get_block_eigenval(block_index, v));
+            double xx = std::abs((*mat)(u, v));
             norm_trace_sq_partial += xx * xx;
           }
         }
         norm_trace_sq += norm_trace_sq_partial;
         // internal check
         if (std::abs(trace_partial) - 1.0000001 * std::sqrt(norm_trace_sq_partial) * get_block_dim(block_index) > 1.e-15)
           TRIQS_RUNTIME_ERROR << "|trace| > dim * norm" << trace_partial << " " << std::sqrt(norm_trace_sq_partial) << "  " << trace_abs;
-        auto dev = std::abs(trace_partial - trace(mat));
+        auto dev = std::abs(trace_partial - trace(*mat));
         if (dev > 1.e-14) TRIQS_RUNTIME_ERROR << "Internal error : trace and density mismatch. Deviation: " << dev;
       }
 

c++/triqs_cthyb/impurity_trace.hpp

@@ -53,7 +53,7 @@ namespace triqs_cthyb {
       cancel_insert_impl(); // in case of an exception, we need to remove any trial nodes before cleaning the tree!
     }
 
-    std::pair<h_scalar_t, h_scalar_t> compute(double p_yee = -1, double u_yee = 0);
+    std::pair<h_scalar_t, h_scalar_t> compute(double p_yee = -1, double u_yee = 0, bool meas_den = false);
 
     // ------- Configuration and h_loc data ----------------
 

c++/triqs_cthyb/measures/density_matrix.cpp

@@ -35,25 +35,31 @@ namespace triqs_cthyb {
   // --------------------
 
   void measure_density_matrix::accumulate(mc_weight_t s) {
-    // we assume here that we are in "Norm" mode, i.e. qmc weight is norm, not trace
-
-    // We need to recompute since the density_matrix in the trace is changed at each computatation,
-    // in particular at the last failed attempt.
-    // So we need to compute it, without any Yee threshold.
-    data.imp_trace.compute();
-    z += s * data.atomic_reweighting;
-    s /= data.atomic_weight; // accumulate matrix / norm since weight is norm * det
-
-    // Careful: there is no reweighting factor here!
-    int size = block_dm.size();
-    for (int i = 0; i < size; ++i)
-      if (data.imp_trace.get_density_matrix()[i].is_valid) { block_dm[i] += s * data.imp_trace.get_density_matrix()[i].mat; }
+    if (data.updated) {
+      if (!flag) {
+        z += data.n_acc * old_z;
+        mc_weight_t s_temp = data.n_acc * old_s;
+        int size = block_dm.size();
+        for (int i = 0; i < size; ++i) 
+          if (data.imp_trace.get_density_matrix()[i].is_valid) { block_dm[i] += s_temp * data.imp_trace.get_density_matrix()[i].mat; }
+      }	
+      data.imp_trace.compute(-1,0,true);
+      old_z = s * data.atomic_reweighting;
+      old_s = s / data.atomic_weight;
+    }
+    flag = false;
   }
 
   // ---------------------------------------------
 
   void measure_density_matrix::collect_results(mpi::communicator const &c) {
-
+
+    z += data.n_acc * old_z;
+    mc_weight_t s_temp = data.n_acc * old_s;
+    int size = block_dm.size();
+    for (int i = 0; i < size; ++i)
+      if (data.imp_trace.get_density_matrix()[i].is_valid) { block_dm[i] += s_temp * data.imp_trace.get_density_matrix()[i].mat; }
+
     z                          = mpi::all_reduce(z, c);
     block_dm                   = mpi::all_reduce(block_dm, c);
     for (auto &b : block_dm){

c++/triqs_cthyb/measures/density_matrix.hpp

@@ -27,6 +27,9 @@ namespace triqs_cthyb {
     qmc_data const &data;
     std::vector<matrix_t> &block_dm; // density matrix of each block
     mc_weight_t z = 0;
+    bool flag = true;  
+    mc_weight_t old_z = 0; 
+    mc_weight_t old_s = 0; 
 
     measure_density_matrix(qmc_data const &data, std::vector<matrix_t> &density_matrix);
     void accumulate(mc_weight_t s);

c++/triqs_cthyb/measures/update_time.hpp

@@ -0,0 +1,71 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2024, Simons Foundation
+ *
+ * TRIQS is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include "../qmc_data.hpp"
+
+namespace triqs_cthyb {
+
+  /// Measure of the average update time
+  struct measure_update_time {
+
+    measure_update_time(qmc_data &_data, double &_update_time) : data(_data), update_time(_update_time) { update_time = 0.0; }
+
+    void accumulate(mc_weight_t) {
+      if (data.updated && !flag) {
+        update_time += data.n_acc;
+        data.n_acc = 1.;
+        ++N;
+        data.updated = false;
+      }
+      else
+        data.n_acc += 1.;
+      if (flag) data.updated = false; // need to set it back to false at the first measurement since it has been set to true during warmup
+      flag = false;
+    }
+
+    void collect_results(mpi::communicator const &comm) {
+      if (N == 0) {
+        N = 1;
+        update_time = data.n_acc;
+      }
+
+      N = mpi::all_reduce(N, comm);
+
+      // Reduce and normalize
+      update_time = mpi::all_reduce(update_time, comm);
+      update_time = update_time / N;
+    }
+
+    private:
+    // The Monte-Carlo configuration
+    qmc_data &data;
+
+    // Flag to indicate whether or not this is the first measure
+    bool flag = true;
+
+    // Reference to double for accumulation
+    double &update_time;
+
+    // Accumulation counter
+    long long N = 0;
+  };
+
+} // namespace triqs_cthyb

c++/triqs_cthyb/moves/double_insert.cpp

@@ -189,6 +189,8 @@ namespace triqs_cthyb {
 
   mc_weight_t move_insert_c_c_cdag_cdag::accept() {
 
+    data.updated = true;
+
     // insert in the tree
     data.imp_trace.confirm_insert();
 

c++/triqs_cthyb/moves/double_remove.cpp

@@ -144,6 +144,8 @@ namespace triqs_cthyb {
   }
 
   mc_weight_t move_remove_c_c_cdag_cdag::accept() {
+
+    data.updated = true;
 
     // remove from the tree
     data.imp_trace.confirm_delete();

c++/triqs_cthyb/moves/global.cpp

@@ -167,6 +167,8 @@ namespace triqs_cthyb {
 
   mc_weight_t move_global::accept() {
 
+    data.updated = true;
+
     for (auto const &o : updated_ops) data.config.replace(o.first, o.second);
     config.finalize();
 

c++/triqs_cthyb/moves/insert.cpp

@@ -141,6 +141,8 @@ namespace triqs_cthyb {
   }
 
   mc_weight_t move_insert_c_cdag::accept() {
+
+    data.updated = true;
 
     // insert in the tree
     data.imp_trace.confirm_insert();

c++/triqs_cthyb/moves/remove.cpp

@@ -117,6 +117,8 @@ namespace triqs_cthyb {
   }
 
   mc_weight_t move_remove_c_cdag::accept() {
+
+    data.updated = true;
 
     // remove from the tree
     data.imp_trace.confirm_delete();

c++/triqs_cthyb/moves/shift.cpp

@@ -205,6 +205,8 @@ namespace triqs_cthyb {
   }
 
   mc_weight_t move_shift_operator::accept() {
+
+    data.updated = true;
 
     // Update the tree
     data.imp_trace.confirm_shift();

c++/triqs_cthyb/qmc_data.hpp

@@ -41,6 +41,8 @@ namespace triqs_cthyb {
     mutable impurity_trace imp_trace;            // Calculator of the trace
     std::vector<int> n_inner;
     block_gf<imtime, delta_target_t> delta; // Hybridization function
+    bool updated;
+    double n_acc;
 
     /// This callable object adapts the Delta function for the call of the det.
     struct delta_block_adaptor {
@@ -76,6 +78,8 @@ namespace triqs_cthyb {
          imp_trace(beta, h_diag, histo_map, p.use_norm_as_weight, p.measure_density_matrix, p.performance_analysis),
          n_inner(n_inner),
          delta(map([](gf_const_view<imtime> d) { return real(d); }, delta)),
+         updated(false),
+         n_acc(0.),
          current_sign(1),
          old_sign(1) {
       std::tie(atomic_weight, atomic_reweighting) = imp_trace.compute();

c++/triqs_cthyb/solver_core.cpp

@@ -44,6 +44,7 @@
 #include "./measures/average_sign.hpp"
 #include "./measures/average_order.hpp"
 #include "./measures/auto_corr_time.hpp"
+#include "./measures/update_time.hpp" 
 #ifdef CTHYB_G2_NFFT
 #include "./measures/G2_tau.hpp"
 #include "./measures/G2_iw.hpp"
@@ -428,6 +429,7 @@ namespace triqs_cthyb {
     qmc.add_measure(measure_average_sign{data, _average_sign}, "Average sign");
     qmc.add_measure(measure_average_order{data, _average_order}, "Average order");
     qmc.add_measure(measure_auto_corr_time{data, _auto_corr_time}, "Auto-correlation time");
+    qmc.add_measure(measure_update_time{data, _update_time}, "Update time");  // Careful, this needs to be added last
 
     // --------------------------------------------------------------------------
 
@@ -441,6 +443,7 @@ namespace triqs_cthyb {
       std::cout << "Average sign: " << _average_sign << std::endl;
       std::cout << "Average order: " << _average_order << std::endl;
       std::cout << "Auto-correlation time: " << _auto_corr_time << std::endl;
+      std::cout << "Average update time: " << _update_time << std::endl; 
     }
 
     // Copy local (real or complex) G_tau back to complex G_tau

c++/triqs_cthyb/solver_core.hpp

@@ -52,6 +52,7 @@ namespace triqs_cthyb {
     mc_weight_t _average_sign;             // average sign of the QMC
     double _average_order;                 // average perturbation order
     double _auto_corr_time;                // Auto-correlation time
+    double _update_time;                   // average update time
     int _solve_status;                     // Status of the solve upon exit: 0 for clean termination, > 0 otherwise.
 
     // Single-particle Green's function containers
@@ -149,6 +150,9 @@ namespace triqs_cthyb {
 
     /// Auto-correlation time
     double auto_corr_time() const { return _auto_corr_time; }
+
+	/// Average update time
+    double update_time() const { return _update_time; } 
 
     /// Status of the ``solve()`` on exit.
     int solve_status() const { return _solve_status; }
@@ -192,6 +196,7 @@ namespace triqs_cthyb {
       h5_write(grp, "average_sign", s._average_sign);
       h5_write(grp, "average_order", s._average_order);
       h5_write(grp, "auto_corr_time", s._auto_corr_time);
+      h5_write(grp, "update_time", s._update_time);
       h5_write(grp, "solve_status", s._solve_status);
       h5_write(grp, "Delta_infty_vec", s.Delta_infty_vec);
     }
@@ -213,6 +218,7 @@ namespace triqs_cthyb {
       h5::try_read(grp, "average_sign", s._average_sign);
       h5::try_read(grp, "average_order", s._average_order);
       h5::try_read(grp, "auto_corr_time", s._auto_corr_time);
+      h5::try_read(grp, "update_time", s._update_time);
       h5::try_read(grp, "solve_status", s._solve_status);
       h5::try_read(grp, "Delta_infty_vec", s.Delta_infty_vec);
 

python/triqs_cthyb/solver_core_desc.py

@@ -309,6 +309,10 @@
                getter = cfunction("double auto_corr_time ()"),
                doc = r"""Auto-correlation time""")
 
+c.add_property(name = "update_time",
+               getter = cfunction("double update_time ()"),
+               doc = r"""Average update time""")			   
+
 c.add_property(name = "solve_status",
                getter = cfunction("int solve_status ()"),
                doc = r"""status of the ``solve()`` on exit.""")