From c80eb49e033a9965f754fabfe73d5754f07c014e Mon Sep 17 00:00:00 2001 From: ec147 <134504305+ec147@users.noreply.github.com> Date: Thu, 25 Jul 2024 14:07:59 +0200 Subject: [PATCH] Measurement of the average update time MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Update impurity_trace.cpp Update density_matrix.cpp Update density_matrix.hpp Delete c++/triqs_cthyb/measures/update_time.hppZone.Identifier Update update_time.hpp Update double_remove.cpp Update insert.cpp Update remove.cpp Update shift.cpp Update qmc_data.hpp Update solver_core.cpp Update solver_core.cpp Update solver_core.hpp test test2 --- c++/triqs_cthyb/impurity_trace.cpp | 24 ++++--- c++/triqs_cthyb/impurity_trace.hpp | 2 +- c++/triqs_cthyb/measures/density_matrix.cpp | 34 ++++++---- c++/triqs_cthyb/measures/density_matrix.hpp | 3 + c++/triqs_cthyb/measures/update_time.hpp | 71 +++++++++++++++++++++ c++/triqs_cthyb/moves/double_insert.cpp | 2 + c++/triqs_cthyb/moves/double_remove.cpp | 2 + c++/triqs_cthyb/moves/global.cpp | 2 + c++/triqs_cthyb/moves/insert.cpp | 2 + c++/triqs_cthyb/moves/remove.cpp | 2 + c++/triqs_cthyb/moves/shift.cpp | 2 + c++/triqs_cthyb/qmc_data.hpp | 4 ++ c++/triqs_cthyb/solver_core.cpp | 3 + c++/triqs_cthyb/solver_core.hpp | 6 ++ python/triqs_cthyb/solver_core_desc.py | 4 ++ 15 files changed, 140 insertions(+), 23 deletions(-) create mode 100644 c++/triqs_cthyb/measures/update_time.hpp diff --git a/c++/triqs_cthyb/impurity_trace.cpp b/c++/triqs_cthyb/impurity_trace.cpp index ec84b14c..32b408d4 100644 --- a/c++/triqs_cthyb/impurity_trace.cpp +++ b/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 impurity_trace::compute(double p_yee, double u_yee) { + std::pair 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>{}; //FIXME complex -- can histos handle this? @@ -382,13 +382,21 @@ 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; } } @@ -396,7 +404,7 @@ namespace triqs_cthyb { // 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; } diff --git a/c++/triqs_cthyb/impurity_trace.hpp b/c++/triqs_cthyb/impurity_trace.hpp index 9463f77e..772881da 100644 --- a/c++/triqs_cthyb/impurity_trace.hpp +++ b/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 compute(double p_yee = -1, double u_yee = 0); + std::pair compute(double p_yee = -1, double u_yee = 0, bool meas_den = false); // ------- Configuration and h_loc data ---------------- diff --git a/c++/triqs_cthyb/measures/density_matrix.cpp b/c++/triqs_cthyb/measures/density_matrix.cpp index 5b9a1289..766102b7 100644 --- a/c++/triqs_cthyb/measures/density_matrix.cpp +++ b/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){ diff --git a/c++/triqs_cthyb/measures/density_matrix.hpp b/c++/triqs_cthyb/measures/density_matrix.hpp index 95616c64..9f244dba 100644 --- a/c++/triqs_cthyb/measures/density_matrix.hpp +++ b/c++/triqs_cthyb/measures/density_matrix.hpp @@ -27,6 +27,9 @@ namespace triqs_cthyb { qmc_data const &data; std::vector &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 &density_matrix); void accumulate(mc_weight_t s); diff --git a/c++/triqs_cthyb/measures/update_time.hpp b/c++/triqs_cthyb/measures/update_time.hpp new file mode 100644 index 00000000..cc0e5406 --- /dev/null +++ b/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 . + * + ******************************************************************************/ +#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 diff --git a/c++/triqs_cthyb/moves/double_insert.cpp b/c++/triqs_cthyb/moves/double_insert.cpp index a78ab91d..b4f8c9ae 100644 --- a/c++/triqs_cthyb/moves/double_insert.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/double_remove.cpp b/c++/triqs_cthyb/moves/double_remove.cpp index 4988f5ca..b6e3a4c1 100644 --- a/c++/triqs_cthyb/moves/double_remove.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/global.cpp b/c++/triqs_cthyb/moves/global.cpp index c0de2cfe..cf022875 100644 --- a/c++/triqs_cthyb/moves/global.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/insert.cpp b/c++/triqs_cthyb/moves/insert.cpp index 9e2acf6e..627961f0 100644 --- a/c++/triqs_cthyb/moves/insert.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/remove.cpp b/c++/triqs_cthyb/moves/remove.cpp index 3c262d3d..06d3f03d 100644 --- a/c++/triqs_cthyb/moves/remove.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/shift.cpp b/c++/triqs_cthyb/moves/shift.cpp index 81604e9d..ed62881e 100644 --- a/c++/triqs_cthyb/moves/shift.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/qmc_data.hpp b/c++/triqs_cthyb/qmc_data.hpp index 12223822..826def97 100644 --- a/c++/triqs_cthyb/qmc_data.hpp +++ b/c++/triqs_cthyb/qmc_data.hpp @@ -41,6 +41,8 @@ namespace triqs_cthyb { mutable impurity_trace imp_trace; // Calculator of the trace std::vector n_inner; block_gf 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 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(); diff --git a/c++/triqs_cthyb/solver_core.cpp b/c++/triqs_cthyb/solver_core.cpp index e4d951f1..ca2a3d28 100644 --- a/c++/triqs_cthyb/solver_core.cpp +++ b/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 diff --git a/c++/triqs_cthyb/solver_core.hpp b/c++/triqs_cthyb/solver_core.hpp index 12995864..309b07af 100644 --- a/c++/triqs_cthyb/solver_core.hpp +++ b/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); diff --git a/python/triqs_cthyb/solver_core_desc.py b/python/triqs_cthyb/solver_core_desc.py index 1f199ce6..ea132ca4 100644 --- a/python/triqs_cthyb/solver_core_desc.py +++ b/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.""")