From a66cbd201ce147ce30c7c946fe5e14211c6547b4 Mon Sep 17 00:00:00 2001 From: Thomas Hahn Date: Thu, 21 Nov 2024 14:32:15 -0500 Subject: [PATCH] Add tests for mc_generic and remove old ones --- c++/triqs/mc_tools/mc_generic.cpp | 12 +- test/c++/mc_tools/CMakeLists.txt | 8 +- test/c++/mc_tools/callback.cpp | 101 ------------ test/c++/mc_tools/different_moves_mc.cpp | 180 --------------------- test/c++/mc_tools/mctools_mc_generic.cpp | 191 +++++++++++++++++++++++ 5 files changed, 199 insertions(+), 293 deletions(-) delete mode 100644 test/c++/mc_tools/callback.cpp delete mode 100644 test/c++/mc_tools/different_moves_mc.cpp create mode 100644 test/c++/mc_tools/mctools_mc_generic.cpp diff --git a/c++/triqs/mc_tools/mc_generic.cpp b/c++/triqs/mc_tools/mc_generic.cpp index e26e5f246..542ab0bee 100644 --- a/c++/triqs/mc_tools/mc_generic.cpp +++ b/c++/triqs/mc_tools/mc_generic.cpp @@ -78,7 +78,7 @@ namespace triqs::mc_tools { for (; !stop_flag; ++cycle_counter) { try { // do cycle_length MC steps / cycle - for (std::int64_t i = 0; i < params.cycle_length; i++) { + for (std::int64_t i = 0; i < params.cycle_length; ++i) { if (triqs::signal_handler::received()) throw triqs::signal_handler::exception{}; // Metropolis step metropolis_step(); @@ -92,7 +92,7 @@ namespace triqs::mc_tools { } catch (std::exception const &err) { // log the exception and node number, either abort or report to other nodes std::cerr << fmt::format("[Rank {}] Error int mc_generic::run: Exception occured:\n{}\n", rank, err.what()); - if (params.propagate_exception) { + if (exception_monitor) { exception_monitor->report_local_event(); break; } else { @@ -101,7 +101,7 @@ namespace triqs::mc_tools { } // recompute fraction done and runtime so far - percentage_done_ = cycle_counter * 100.0 / params.ncycles; + percentage_done_ = static_cast(cycle_counter) * 100.0 / params.ncycles; double runtime = run_timer_; // is it time to print simulation info @@ -247,7 +247,7 @@ namespace triqs::mc_tools { info += fmt::format("[Rank {}] Measurement durations:\n{}", c.rank(), measures_.get_timings(fmt::format("[Rank {}] ", c.rank()))); info += fmt::format("[Rank {}] Move statistics:\n{}", c.rank(), moves_.get_statistics(fmt::format("[Rank {}] ", c.rank()))); - // gather all output string on rank 0 to print in order + // gather all output strings on rank 0 to print in order auto all_infos = mpi::gather(info, c); if (c.rank() == 0) { report_(3) << all_infos; @@ -262,7 +262,7 @@ namespace triqs::mc_tools { // current simulation parameters auto const rank = params.comm.rank(); double const runtime = run_timer_; - auto const cycles_per_sec = cycle_counter / runtime; + auto const cycles_per_sec = static_cast(cycle_counter) / runtime; // do the printing if (percentage_done_ < 0) { @@ -288,7 +288,7 @@ namespace triqs::mc_tools { params.after_cycle_duty(); if (params.enable_calibration) moves_.calibrate(params.comm); if (params.enable_measures) { - nmeasures_done_++; + ++nmeasures_done_; for (auto &m : measures_aux_) m(); measures_.accumulate(sign_); } diff --git a/test/c++/mc_tools/CMakeLists.txt b/test/c++/mc_tools/CMakeLists.txt index 633aac1cd..fff88234f 100644 --- a/test/c++/mc_tools/CMakeLists.txt +++ b/test/c++/mc_tools/CMakeLists.txt @@ -1,10 +1,6 @@ all_tests() -set(TEST_MPI_NUMPROC 2) -add_cpp_test(different_moves_mc) -add_cpp_test(callback) -set(TEST_MPI_NUMPROC 3) -add_cpp_test(different_moves_mc) set(TEST_MPI_NUMPROC 4) add_cpp_test(mctools_move_set) -add_cpp_test(different_moves_mc) +add_cpp_test(mctools_measure_set) +add_cpp_test(mctools_mc_generic) diff --git a/test/c++/mc_tools/callback.cpp b/test/c++/mc_tools/callback.cpp deleted file mode 100644 index 84d6ba951..000000000 --- a/test/c++/mc_tools/callback.cpp +++ /dev/null @@ -1,101 +0,0 @@ -// Copyright (c) 2023 Simons Foundation -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0.txt -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. -// -// Authors: Nils Wentzell - -#ifndef NDA_DEBUG -#define NDA_DEBUG -#endif - -#include - -#include -#include -#include - -// --------------- the configuration a spin, the inverse temperature, the external field --- -struct config_t { - double beta, h; - int spin = -1; -}; - -// --------------- a move: flip the spin --------------- -struct flip { - config_t &config; - double attempt() { return std::exp(-2 * config.spin * config.h * config.beta); } - double accept() { - config.spin *= -1; - return 1.0; - } - void reject() {} -}; - -// ----------------- a measurement: the magnetization ------------ -struct compute_m { - config_t &config; - double &Z; - double M = 0; - - void accumulate(double sign) { - Z += sign; - M += sign * config.spin; - } - - void collect_results(mpi::communicator c) { - auto Z_sum = mpi::reduce(Z, c); - auto M_sum = mpi::reduce(M, c); - if (c.rank() == 0) std::cout << "Magnetization: " << M_sum / Z_sum << std::endl; - } -}; - -TEST(mc_generic, callback) { - - // initialize mpi - mpi::communicator world; - - // prepare the MC parameters - int n_warmup_cycles = 0; - int n_cycles = 1000; - int length_cycle = 10; - std::string random_name = ""; - int random_seed = 374982 + world.rank() * 273894; - int verbosity = (world.rank() == 0 ? 2 : 0); - - // construct a Monte Carlo loop - triqs::mc_tools::mc_generic SpinMC(random_name, random_seed, verbosity); - - // parameters of the model - double beta = 0.3; - double field = 0.5; - - // construct config_t - config_t config{beta, field}; - - // add moves - SpinMC.add_move(flip{config}, "flip move"); - - // Add measures - double Z = 0.0; - SpinMC.add_measure(compute_m{config, Z}, "magnetization measure"); - - // Run and collect results - // We are using the callback to stop the simulation after n_cycles - auto call_back_cycles = [n_cycles, n = 0]() mutable { return ++n >= n_cycles; }; - SpinMC.warmup_and_accumulate(n_warmup_cycles, -1, length_cycle, call_back_cycles); - SpinMC.collect_results(world); - - EXPECT_EQ(Z, double(n_cycles)); -} - -MAKE_MAIN; diff --git a/test/c++/mc_tools/different_moves_mc.cpp b/test/c++/mc_tools/different_moves_mc.cpp deleted file mode 100644 index d804368c4..000000000 --- a/test/c++/mc_tools/different_moves_mc.cpp +++ /dev/null @@ -1,180 +0,0 @@ -// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) -// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS) -// Copyright (c) 2018-2021 Simons Foundation -// -// This program 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. -// -// This program 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 may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Laura Messio, Olivier Parcollet, Priyanka Seth, Nils Wentzell - -// The goal of this program is to give an example of the use of the probability for proposing a move (here, pl and pr) -// The walker try to go towards its left with probability 2*pl/(pl+pr) (pl is given to the move_left constructor), -// or towards the right with probability 2*pr/(pl+pr) (pr is given to the move_right constructor). -// The left move is accepted with probability proba = pr/pl and the right move with probability proba=pl/pr (in reality, each time with probability min(1,proba) ). -// Two h5 files are produced: -// - histo.h5 -// - params.h5 -// A python notebook is associated to this example, with the comparison to the theoretical gaussian histogram. - -#ifndef NDA_DEBUG -#define NDA_DEBUG -#endif - -#include -#include -#include -#include -nda::array, 1> make_array(std::complex c) { return {c}; }; - -struct configuration { - int x; //position of the walker - configuration() : x(0) {} //constructor -}; - -long node_that_will_raise_exception = 1; - -//------------------------------------------------------------ -// MOVES -//------------------------------------------------------------ -struct move_left { - configuration *config; - double proba = -9; - inline static long counter = 0; - mpi::communicator c; // for error test only - triqs::mc_tools::random_generator &RND; - move_left(configuration &config_, double pl, double pr, triqs::mc_tools::random_generator &RND_) - : config(&config_), proba(pr / pl), RND(RND_) {} //constructor - double attempt() { - if ((c.rank() == node_that_will_raise_exception) and counter++ == 2000) TRIQS_RUNTIME_ERROR << "A bad ERROR"; - return proba; - } - double accept() { - config->x += -1; - return 1; - } - void reject() {} - move_left(move_left const &x) : config(x.config), proba(x.proba), RND(x.RND) { std::cout << "copy move_left" << std::endl; } - move_left(move_left &&x) = default; //: config(x.config), proba(x.proba), RND(x.RND) { std::cout << "copy move_left"<x += 1; - return 1; - } - void reject() {} -}; - -//------------------------------------------------------------ -// MEASUREMENTS -//------------------------------------------------------------ -struct compute_histo { - configuration *config; - nda::array H; // an histogram of the positions of the walker - long xmax; //maximal position registered in the histo - mpi::communicator world; - long tot; //number of pointsregistered inthe histogram - compute_histo(configuration &config_, nda::array &H_, int xmax_) : config(&config_), H(H_), xmax(xmax_), tot(0) {} - void accumulate(double) { - auto x = config->x; - if (x + xmax >= 0 && x - xmax < 0) { - H(lround(floor((x + xmax)))) += 1; - tot += 1; - } - config->x = 0; // the walker returns to zero. - } - void collect_results(mpi::communicator) { - H /= tot; - if (world.rank() == 0) { - h5::file file("histo.h5", 'w'); - h5_write(file, "H", H); - } - } -}; - -mpi::communicator world; - -//------------------------------------------------------------ -// MAIN -//------------------------------------------------------------ - -bool test_mc(bool recover_from_exception) { - - // prepare the MC parameters - int N_Cycles = 100000; // number of different walks. - int Length_Cycle = 100; // the walker makes Length_Cycle steps during each walk - int N_Warmup_Cycles = 0; - std::string Random_Name = ""; - int Random_Seed = 374982 + world.rank() * 273894; - int Verbosity = (world.rank() == 0 ? 3 : 0); - int xmax = static_cast(floor(4 * sqrt(Length_Cycle))); // typical length of a walk - double pl = 2.5, pr = 1; //non normalized probabilities for proposing a left or right move - - if (world.rank() == 0) { - h5::file file("params.h5", 'w'); - h5_write(file, "pr", pr); - h5_write(file, "pl", pl); - h5_write(file, "xmax", xmax); - h5_write(file, "N_Cycles", N_Cycles); - h5_write(file, "Length_Cycle", Length_Cycle); - } - - // construct a Monte Carlo loop - triqs::mc_tools::mc_generic IntMC(Random_Name, Random_Seed, Verbosity); - - // construct configuration - configuration config; - nda::array H(2 * xmax); - H() = 0; - - // add moves and measures - IntMC.add_move(move_left(config, pl, pr, IntMC.get_rng()), "left move", pl); - IntMC.add_move(move_right(config, pl, pr, IntMC.get_rng()), "right move", pr); - IntMC.add_measure(compute_histo(config, H, xmax), "histogramn measure"); - - bool stopped_by_exception = false; - - try { - // Run and collect results - IntMC.warmup_and_accumulate(N_Warmup_Cycles, N_Cycles, Length_Cycle, triqs::utility::clock_callback(600)); - IntMC.collect_results(world); - } catch (std::exception const &e) { - std::cerr << "TEST : Exception occurred. Node " << world.rank() << " is stopping cleanly" << std::endl; - stopped_by_exception = true; - } - - return stopped_by_exception; -} - -TEST(mc_generic, WithException) { - node_that_will_raise_exception = 1; - bool stopped_by_exception = test_mc(true); - if (world.size() == 1) - EXPECT_FALSE(stopped_by_exception); // one node only. Should work - else - EXPECT_TRUE(stopped_by_exception); -} - -TEST(mc_generic, WithOutException) { - node_that_will_raise_exception = -1; - bool stopped_by_exception = test_mc(true); - EXPECT_FALSE(stopped_by_exception); // one node only. Should work -} - -MAKE_MAIN; diff --git a/test/c++/mc_tools/mctools_mc_generic.cpp b/test/c++/mc_tools/mctools_mc_generic.cpp new file mode 100644 index 000000000..727f4534c --- /dev/null +++ b/test/c++/mc_tools/mctools_mc_generic.cpp @@ -0,0 +1,191 @@ +// Copyright (c) 2015-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) +// Copyright (c) 2015-2018 Centre national de la recherche scientifique (CNRS) +// Copyright (c) 2018 Simons Foundation +// +// This program 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. +// +// This program 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 may obtain a copy of the License at +// https://www.gnu.org/licenses/gpl-3.0.txt +// +// Authors: Michel Ferrero, Olivier Parcollet, Nils Wentzell + +#include +#include +#include + +#include + +#include +#include +#include +#include + +// MC configuration for integrating a sin(x) using naive MC integration. +struct mc_config { + double a{0}; + double b{std::numbers::pi}; + double x{0}; + [[nodiscard]] double eval() const { return std::sin(x); } +}; + +// Draw a new x from a uniform distribution. +struct move_x { + using rng_type = triqs::mc_tools::random_generator; + std::reference_wrapper config; + std::reference_wrapper rng; + double new_x{0}; + move_x(mc_config &config, rng_type &rng) : config(config), rng(rng) {} + double attempt() { + new_x = rng.get()(config.get().a, config.get().b); + return 1.0; + } + double accept() { + config.get().x = new_x; + return 1.0; + } +}; + +// Measure the integral. +struct integral_1d { + std::reference_wrapper config; + double sum_fx{0.0}; + double sum_fx_sq{0.0}; + std::int64_t count{0}; + integral_1d(mc_config &config) : config(config) {} + void accumulate(double) { + auto fx = config.get().eval(); + sum_fx += fx; + sum_fx_sq += fx * fx; + ++count; + } + void collect_results(mpi::communicator c) { + auto sum_fx_red = mpi::reduce(sum_fx, c); + auto sum_fx_sq_red = mpi::reduce(sum_fx_sq, c); + auto count_red = static_cast(mpi::reduce(count, c)); + if (c.rank() == 0) { + auto ba = config.get().b - config.get().a; + auto mean = sum_fx_red / count_red * ba; + auto var = (sum_fx_sq_red - sum_fx_red / count_red) / (count_red - 1) * ba * ba; + fmt::print("Integral = {} +/- {}\n", mean, std::sqrt(var / count_red)); + } + } +}; + +// Fixture for mc_generic tests. +class TRIQSMCTools : public ::testing::Test { + protected: + void SetUp() override { + // add move + mc.add_move(move_x{config, mc.get_rng()}, "move_x"); + + // add measure + mc.add_measure(integral_1d{config}, "integral_1d"); + } + + mpi::communicator comm{}; + triqs::mc_tools::mc_generic mc{"mt19937", comm.rank() * 0x12af5988, 3}; + mc_config config{}; +}; + +// Test a basic MC simulation with specific parameters. +void test_mc_basic(triqs::mc_tools::mc_generic &mc, std::int64_t ncycles, std::int64_t cycle_length) { + auto params = mc.get_run_params(); + params.ncycles = ncycles; + params.cycle_length = cycle_length; + EXPECT_EQ(mc.accumulate(params), 0); + mc.collect_results(params.comm); + EXPECT_DOUBLE_EQ(mc.get_percent(), 100); + EXPECT_EQ(mc.get_current_cycle_number(), params.ncycles); + EXPECT_EQ(mc.get_nmeasures(), ncycles); + EXPECT_EQ(mc.get_config_id(), params.ncycles * params.cycle_length); + if (ncycles > 0) EXPECT_DOUBLE_EQ(mc.get_acceptance_rates().at("move_x"), 1); +} + +TEST_F(TRIQSMCTools, MCGenericBasicRunNCycles0) { test_mc_basic(mc, 0, 1); } + +TEST_F(TRIQSMCTools, MCGenericBasicRunCycleLength1) { test_mc_basic(mc, 1000000, 1); } + +TEST_F(TRIQSMCTools, MCGenericBasicRunCycleLength27) { test_mc_basic(mc, 10000, 27); } + +TEST_F(TRIQSMCTools, MCGenericStopCallback) { + auto params = mc.get_run_params(); + params.stop_callback = [this]() { + if (mc.get_config_id() == 100) return true; + return false; + }; + EXPECT_EQ(mc.accumulate(params), 1); + mc.collect_results(params.comm); + EXPECT_LT(mc.get_percent(), 0); + EXPECT_EQ(mc.get_current_cycle_number(), 100); + EXPECT_EQ(mc.get_nmeasures(), 100); + EXPECT_EQ(mc.get_config_id(), 100); + EXPECT_DOUBLE_EQ(mc.get_acceptance_rates().at("move_x"), 1); +} + +TEST_F(TRIQSMCTools, MCGenericAfterCycleDuty) { + int counter = 0; + auto params = mc.get_run_params(); + params.ncycles = 100; + params.after_cycle_duty = [&counter]() { ++counter; }; + EXPECT_EQ(mc.accumulate(params), 0); + mc.collect_results(params.comm); + EXPECT_DOUBLE_EQ(mc.get_percent(), 100); + EXPECT_EQ(mc.get_current_cycle_number(), params.ncycles); + EXPECT_EQ(mc.get_nmeasures(), params.ncycles); + EXPECT_EQ(mc.get_config_id(), params.ncycles * params.cycle_length); + EXPECT_DOUBLE_EQ(mc.get_acceptance_rates().at("move_x"), 1); +} + +TEST_F(TRIQSMCTools, MCGenericWarmuup) { + auto params = mc.get_run_params(); + params.ncycles = 100; + EXPECT_EQ(mc.warmup(params), 0); + EXPECT_DOUBLE_EQ(mc.get_percent(), 100); + EXPECT_EQ(mc.get_current_cycle_number(), params.ncycles); + EXPECT_EQ(mc.get_nmeasures(), 0); + EXPECT_EQ(mc.get_config_id(), params.ncycles * params.cycle_length); +} + +TEST_F(TRIQSMCTools, MCGenericContinueAfterNCyclesDone) { + auto params = mc.get_run_params(); + params.ncycles = static_cast(1000000) * (comm.rank() + 1); + params.continue_after_ncycles_done = true; + params.check_cycles_interval = 0.01; + EXPECT_EQ(mc.accumulate(params), 0); + mc.collect_results(params.comm); + EXPECT_GE(mc.get_percent(), 100); + EXPECT_GE(mc.get_current_cycle_number(), params.ncycles); + EXPECT_GE(mc.get_nmeasures(), params.ncycles); + EXPECT_GE(mc.get_config_id(), params.ncycles * params.cycle_length); + EXPECT_DOUBLE_EQ(mc.get_acceptance_rates().at("move_x"), 1); +} + +TEST_F(TRIQSMCTools, MCGenericPropagateException) { + if (comm.size() != 1) { + auto params = mc.get_run_params(); + params.ncycles = -1; + params.propagate_exception = true; + params.check_exception_interval = 0.1; + params.after_cycle_duty = [this]() { + if (comm.rank() == 0 && mc.get_config_id() == 10) throw std::runtime_error("Some exception on rank 0!"); + }; + bool exc_caught = false; + try { + mc.accumulate(params); + } catch (std::exception const &e) { + fmt::print("[Rank {}] Exception caught: {}\n", comm.rank(), e.what()); + exc_caught = true; + } + EXPECT_TRUE(exc_caught); + } +} + +MAKE_MAIN;