From 3b673bbef831dab4186b5989bc00afd0a185e69f Mon Sep 17 00:00:00 2001 From: Igor Krivenko Date: Tue, 22 Feb 2022 12:07:17 +0100 Subject: [PATCH] Draft implementation of `som_core::compute_final_solution_cc()` This method computes the final solution using the consistent constraints approach as described in Section II.B of O. Goulko et al. Phys. Rev. B 95, 014102 (2017). The new method is also exposed in Python, although I could not mark it as `release_GIL_and_enable_signal = True` because of cpp2py issue TRIQS/cpp2py#44. --- c++/som/kernels/mesh_traits.hpp | 5 + c++/som/som_core/accumulate.cpp | 19 +- c++/som/som_core/adjust_f.cpp | 16 +- c++/som/som_core/common.hxx | 22 +- c++/som/som_core/final_solution.cpp | 468 +++++++++++++++++++- c++/som/som_core/parameters.hpp | 70 +++ c++/som/som_core/som_core.hpp | 17 +- doc/notes/final_solution_cc.nb | 638 ++++++++++++++++++++++++++++ python/som/som_core_desc.py | 128 +++++- test/c++/CMakeLists.txt | 2 +- 10 files changed, 1337 insertions(+), 48 deletions(-) create mode 100644 doc/notes/final_solution_cc.nb diff --git a/c++/som/kernels/mesh_traits.hpp b/c++/som/kernels/mesh_traits.hpp index b7521dea4..5fba30ef1 100644 --- a/c++/som/kernels/mesh_traits.hpp +++ b/c++/som/kernels/mesh_traits.hpp @@ -21,6 +21,7 @@ #pragma once #include +#include #include @@ -49,4 +50,8 @@ template <> struct mesh_traits { inline static std::string name() { return "real frequency"; } }; +using mesh_variant_t = std::variant, + triqs::gfs::gf_mesh, + triqs::gfs::gf_mesh>; + } // namespace som diff --git a/c++/som/som_core/accumulate.cpp b/c++/som/som_core/accumulate.cpp index 8ee541d17..d7d6d0893 100644 --- a/c++/som/som_core/accumulate.cpp +++ b/c++/som/som_core/accumulate.cpp @@ -26,7 +26,8 @@ #include #include -#include +#include +#include #include @@ -226,16 +227,12 @@ void som_core::accumulate(accumulate_parameters_t const& p) { triqs::signal_handler::start(); accumulate_status = 0; try { -#define RUN_IMPL_CASE(r, okmk) \ - case(int(BOOST_PP_SEQ_ELEM(0, okmk)) + \ - n_observable_kinds * mesh_traits::index): \ - accumulate_impl>(); \ - break; - switch(int(kind) + n_observable_kinds * mesh.index()) { - BOOST_PP_SEQ_FOR_EACH_PRODUCT(RUN_IMPL_CASE, - (ALL_OBSERVABLES)(ALL_INPUT_MESHES)) - } -#undef RUN_IMPL_CASE +#define IMPL_CASE(r, okmk) \ + case(kernel_id(BOOST_PP_SEQ_ELEM(0, okmk), BOOST_PP_SEQ_ELEM(1, okmk){})): \ + accumulate_impl>(); break; + + SELECT_KERNEL(IMPL_CASE, som_core::accumulate()) +#undef IMPL_CASE } catch(stopped& e) { accumulate_status = e.code; triqs::signal_handler::received(true); diff --git a/c++/som/som_core/adjust_f.cpp b/c++/som/som_core/adjust_f.cpp index 5a746d5a4..a5648c9db 100644 --- a/c++/som/som_core/adjust_f.cpp +++ b/c++/som/som_core/adjust_f.cpp @@ -168,18 +168,12 @@ int som_core::adjust_f(adjust_f_parameters_t const& p) { << " particular solutions ..." << std::endl; } -#define RUN_IMPL_CASE(r, okmk) \ - case(int(BOOST_PP_SEQ_ELEM(0, okmk)) + \ - n_observable_kinds * mesh_traits::index): \ +#define IMPL_CASE(r, okmk) \ + case(kernel_id(BOOST_PP_SEQ_ELEM(0, okmk), BOOST_PP_SEQ_ELEM(1, okmk){})): \ return adjust_f_impl>(p); - switch(int(kind) + n_observable_kinds * mesh.index()) { - BOOST_PP_SEQ_FOR_EACH_PRODUCT(RUN_IMPL_CASE, - (ALL_OBSERVABLES)(ALL_INPUT_MESHES)) - default: - fatal_error("Unknown observable kind " + to_string(kind) + - " in adjust_f()"); - } -#undef RUN_IMPL_CASE + + SELECT_KERNEL(IMPL_CASE, som_core::adjust_f()) +#undef IMPL_CASE } } // namespace som diff --git a/c++/som/som_core/common.hxx b/c++/som/som_core/common.hxx index 5f7df1312..4bb2a9ace 100644 --- a/c++/som/som_core/common.hxx +++ b/c++/som/som_core/common.hxx @@ -44,7 +44,7 @@ inline void warning(std::string const& message) { std::cout << "WARNING: " << message << std::endl; } -inline std::ostream & mpi_cout(mpi::communicator const& comm) { +inline std::ostream& mpi_cout(mpi::communicator const& comm) { return (std::cout << "[Rank " << comm.rank() << "] "); } @@ -74,4 +74,24 @@ void check_gf_stat(gf_view g, statistic_enum expected_stat) { return check_gf_stat(make_const_view(g), expected_stat); } +template +inline constexpr int kernel_id(observable_kind kind, MeshType) { + return int(kind) + n_observable_kinds * mesh_traits::index; +} + +inline int kernel_id(observable_kind kind, mesh_variant_t const& mesh) { + return int(kind) + n_observable_kinds * mesh.index(); +} + +#define FOR_EACH_KERNEL(F) \ + BOOST_PP_SEQ_FOR_EACH_PRODUCT(F, (ALL_OBSERVABLES)(ALL_INPUT_MESHES)) + +#define SELECT_KERNEL(F, F_NAME) \ + switch(kernel_id(kind, mesh)) { \ + FOR_EACH_KERNEL(F) \ + default: \ + fatal_error("Unknown observable kind " + std::to_string(kind) + \ + " in " #F_NAME); \ + } + } // namespace som diff --git a/c++/som/som_core/final_solution.cpp b/c++/som/som_core/final_solution.cpp index 549511777..3e753a69a 100644 --- a/c++/som/som_core/final_solution.cpp +++ b/c++/som/som_core/final_solution.cpp @@ -19,11 +19,20 @@ * ******************************************************************************/ +#include +#include #include +#include +#include -#include +#include +#include +#include + +#include #include +#include #include #include "common.hxx" @@ -31,6 +40,9 @@ namespace som { +using namespace triqs::arrays; +using namespace itertools; + //////////////////////////////////////////// // som_core::data_t::compute_objf_final() // //////////////////////////////////////////// @@ -38,8 +50,8 @@ namespace som { template void som_core::data_t::compute_objf_final(KernelType const& kernel) { using mesh_t = typename KernelType::mesh_type; - objective_function of(kernel, get_rhs(), - get_error_bars()); + objective_function of( + kernel, get_rhs(), get_error_bars()); objf_final = of(final_solution); } @@ -55,7 +67,7 @@ std::vector som_core::compute_objf_final_impl() { KernelType kernel(m); std::vector objf_final(data.size()); - for(int n = 0; n < data.size(); ++n) { + for(int n : range(data.size())) { auto& d = data[n]; d.compute_objf_final(kernel); objf_final[n] = d.objf_final; @@ -65,19 +77,12 @@ std::vector som_core::compute_objf_final_impl() { } std::vector som_core::compute_objf_final() { -#define RUN_IMPL_CASE(r, okmk) \ - case(int(BOOST_PP_SEQ_ELEM(0, okmk)) + \ - n_observable_kinds * mesh_traits::index): \ +#define IMPL_CASE(r, okmk) \ + case(kernel_id(BOOST_PP_SEQ_ELEM(0, okmk), BOOST_PP_SEQ_ELEM(1, okmk){})): \ return compute_objf_final_impl>(); - switch(int(kind) + n_observable_kinds * mesh.index()) { - BOOST_PP_SEQ_FOR_EACH_PRODUCT(RUN_IMPL_CASE, - (ALL_OBSERVABLES)(ALL_INPUT_MESHES)) - default: - fatal_error("Unknown observable kind " + to_string(kind) + - " in compute_objf_final()"); - } -#undef RUN_IMPL_CASE + SELECT_KERNEL(IMPL_CASE, som_core::compute_objf_final()) +#undef IMPL_CASE } //////////////////////////////////////// @@ -86,7 +91,7 @@ std::vector som_core::compute_objf_final() { std::vector som_core::compute_final_solution(double good_chi_rel, double good_chi_abs) { - for(int n = 0; n < data.size(); ++n) { + for(int n : range(data.size())) { auto& d = data[n]; configuration sol_sum(ci); @@ -104,14 +109,441 @@ std::vector som_core::compute_final_solution(double good_chi_rel, } // Sum over all processes - n_good_solutions = mpi::all_reduce(n_good_solutions); - d.final_solution = mpi::all_reduce(sol_sum); + n_good_solutions = mpi::all_reduce(n_good_solutions, comm); + d.final_solution = mpi::all_reduce(sol_sum, comm); d.final_solution *= 1.0 / double(n_good_solutions); } return compute_objf_final(); } +/////////////////////////////////////////// +// som_core::compute_final_solution_cc() // +/////////////////////////////////////////// + +void final_solution_cc_parameters_t::validate() const { + using std::to_string; + + struct { + int mesh_size; + void operator()(triqs::gfs::gf_mesh const& m) { + mesh_size = m.size(); + if(mesh_size < 3) + fatal_error("There must be at least 3 energy points in refreq_mesh"); + } + void operator()(array const& m) { + mesh_size = m.size(); + if(mesh_size < 3) + fatal_error("There must be at least 3 energy points in refreq_mesh"); + if(!std::is_sorted(m.begin(), m.end(), std::less_equal())) + fatal_error( + "Energy points in refreq_mesh must be listed in strictly ascending " + "order"); + } + } validate_mesh; + std::visit(validate_mesh, refreq_mesh); + + if(good_chi_rel < 1.0) + fatal_error("Parameter good_chi_rel must be >= 1.0 (got " + + to_string(good_chi_rel) + ")"); + + if(good_chi_abs < 0) + fatal_error("Parameter good_chi_abs must be non-negative (got " + + to_string(good_chi_abs) + ")"); + + if((default_model.size() != 0) && // TODO: Use nda::array::empty() + (default_model.size() != validate_mesh.mesh_size)) + fatal_error("default_model must have the same size as refreq_mesh (got " + + to_string(default_model.size()) + " vs " + + to_string(validate_mesh.mesh_size) + ")"); + + if((default_model_weights.size() != 0) && // TODO: Use nda::array::empty() + (default_model_weights.size() != validate_mesh.mesh_size)) + fatal_error( + "default_model_weights must have the same size as refreq_mesh (got " + + to_string(default_model_weights.size()) + " vs " + + to_string(validate_mesh.mesh_size) + ")"); + + if(max_iter <= 0) + fatal_error("Number of iterations max_iter must be positive (got " + + to_string(max_iter) + ")"); + + if(unity_sum_coeff < 0) + fatal_error("Parameter unity_sum_coeff must positive (got " + + to_string(unity_sum_coeff) + ")"); + + if(amp_penalty_max <= 0) + fatal_error("Parameter amp_penalty_max must be positive (got " + + to_string(amp_penalty_max) + ")"); + + if(amp_penalty_divisor <= 1.0) + fatal_error("Parameter amp_penalty_divisor must exceed 1.0 (got " + + to_string(amp_penalty_divisor) + ")"); + + if(der_penalty_init <= 0) + fatal_error("Parameter der_penalty_init must be positive (got " + + to_string(der_penalty_init) + ")"); + + if(der_penalty_coeff <= 1.0) + fatal_error("Parameter der_penalty_coeff must exceed 1.0 (got " + + to_string(der_penalty_coeff) + ")"); +} + +vector make_rhs_vector(array const& A_j_block, + double U, + array const& T_k, + array const& A_T_k, + global_index_map const& index_map, + mpi::communicator const& comm) { + int J = index_map.size(); + vector b(J); + b() = 0; + + // RHS generated by O_5 + if(T_k.size() != 0) { // TODO: Use nda::array::empty() + for(int j : range(first_dim(A_j_block))) + b(j) += sum(T_k * A_T_k * A_j_block(j, range())); + b = mpi::all_reduce(b, comm); + } + + // RHS generated by O_2 and O_3 + b() += U + 1.0 / double(J); + + return b; +} + +void update_O_mat(array& A_j_local_block, + array& Ap_j_local_block, + array& App_j_local_block, + double U, + array const& Q_k, + array const& D_k, + array const& B_k, + array const& T_k, + global_index_map const& index_map, + mpi::communicator const& comm, + bool verbose, + matrix& O_mat) { + int const j_block1 = comm.rank(); + int const J1 = first_dim(A_j_local_block); + assert(first_dim(Ap_j_local_block) == J1); + assert(first_dim(App_j_local_block) == J1); + + if(verbose) + mpi_cout(comm) << " Updating matrix of quadratic form O" << std::endl; + + // Terms O_1, O_4 and O_5 + + // A_j(e_k) block received from other ranks + array A_j_remote_block; + // A'_j(e_k) block received from other ranks + array Ap_j_remote_block; + // A''_j(e_k) block received from other ranks + array App_j_remote_block; + + auto const& A_j_block1 = A_j_local_block; + auto const& Ap_j_block1 = Ap_j_local_block; + auto const& App_j_block1 = App_j_local_block; + + O_mat() = 0; + + for(int j_block2 : range(comm.size())) { + bool diag_block = j_block2 == j_block1; + + auto& A_j_block2 = diag_block ? A_j_local_block : A_j_remote_block; + auto& Ap_j_block2 = diag_block ? Ap_j_local_block : Ap_j_remote_block; + auto& App_j_block2 = diag_block ? App_j_local_block : App_j_remote_block; + + mpi::broadcast(A_j_block2, comm, j_block2); + mpi::broadcast(Ap_j_block2, comm, j_block2); + mpi::broadcast(App_j_block2, comm, j_block2); + + int const J2 = first_dim(A_j_block2); + assert(first_dim(Ap_j_block2) == J2); + assert(first_dim(App_j_block2) == J2); + + if(verbose) { + mpi_cout(comm) << " Combining " << J1 << " solutions from rank " + << j_block1 << " with " << J2 << " solutions from rank " + << j_block2 << std::endl; + } + + for(int j1_local : range(J1)) { + int j1 = index_map(j_block1, j1_local); + for(int j2_local : range(J2)) { + int j2 = index_map(j_block2, j2_local); + + O_mat(j1, j2) += sum(Q_k * A_j_block1(j1_local, range()) * + A_j_block2(j2_local, range())); + O_mat(j1, j2) += sum(D_k * D_k * Ap_j_block1(j1_local, range()) * + Ap_j_block2(j2_local, range())); + O_mat(j1, j2) += sum(B_k * B_k * App_j_block1(j1_local, range()) * + App_j_block2(j2_local, range())); + + if(T_k.size() != 0) // TODO: Use nda::array::empty() + O_mat(j1, j2) += sum(T_k * A_j_block1(j1_local, range()) * + A_j_block2(j2_local, range())); + } + } + } + + O_mat = mpi::all_reduce(O_mat, comm); + + // Unity sum constraint O_2 + O_mat.as_array_view() += U; + + // Penalty for large deviations from the equal-weight superposition, O_3 + O_mat() += 1.0; +} + +vector cc_protocol( + std::vector> const& particular_solutions, + std::vector const& good_solution_indices, + array const e_k_list, + double convergence_tol, + final_solution_cc_parameters_t const& p, + global_index_map const& index_map, + mpi::communicator const& comm) { + int n_good_solutions = good_solution_indices.size(); + + int const my_rank = comm.rank(); + auto my_j_range = + range(index_map.range_start(my_rank), + index_map.range_start(my_rank) + index_map.range_size(my_rank)); + assert(my_j_range.size() == n_good_solutions); + + // Rank-local portion of A_j(e_k) array (partition over the first index) + array A_block(n_good_solutions, e_k_list.size()); + // Rank-local portion of A'_j(e_k) array (partition over the first index) + array Ap_block(n_good_solutions, e_k_list.size() - 1); + // Rank-local portion of A''_j(e_k) array (partition over the first index) + array App_block(n_good_solutions, e_k_list.size() - 2); + + for(auto const& [j, s_index] : enumerate(good_solution_indices)) { + // Fill A + for(auto [k, e_k] : enumerate(e_k_list)) + A_block(j, k) = particular_solutions[s_index].first(e_k); + auto A_view = make_const_view(A_block(j, range())); + + // Fill Ap + finite_diff_forward(A_view, e_k_list, Ap_block(j, range())); + // Fill App + finite_diff_2_symm(A_view, e_k_list, App_block(j, range())); + } + + // Sum of A_j(e_k) over j + array A_k(second_dim(A_block)); + // Sum of A'_j(e_k) over j + array Ap_k(second_dim(Ap_block)); + // Sum of A''_j(e_k) over j + array App_k(second_dim(App_block)); + + double DD = p.der_penalty_init; + + // Amplitude regularization parameters Q_k + array Q_k(A_k.size()); + Q_k() = 0; // TODO: Use nda::zeros() + // First derivative regularization parameters D_k + array D_k(Ap_k.size()); + D_k() = DD; // TODO: Use nda::ones() + // Second derivative regularization parameters B_k + array B_k(App_k.size()); + B_k() = DD; // TODO: Use nda::ones() + + // Total number of selected particular solutions (summed over all ranks) + int J = index_map.size(); + + // Coefficients c_j + vector c(J); + c() = 1.0 / double(J); // TODO: Use nda::ones() + + // Prepare RHS vector + auto b = make_rhs_vector(A_block, + p.unity_sum_coeff, + p.default_model_weights, + p.default_model, + index_map, + comm); + + // Matrix of the quadratic form + matrix O_mat(J, J); + + int iter; + for(iter = 0; iter < p.max_iter; ++iter) { + if(p.verbosity >= 1) + mpi_cout(comm) << " Iteration " << (iter + 1) << " out of " << p.max_iter + << std::endl; + + // Update matrix of the quadratic form + update_O_mat(A_block, + Ap_block, + App_block, + p.unity_sum_coeff, + Q_k, + D_k, + B_k, + p.default_model_weights, + index_map, + comm, + p.verbosity >= 2, + O_mat); + + // Minimize the functional + c = inverse(O_mat) * b; + + // + // Update regularization parameters + // + + DD *= p.der_penalty_coeff; + D_k *= p.der_penalty_coeff; + B_k *= p.der_penalty_coeff; + + triqs::clef::placeholder<0> j_; + triqs::clef::placeholder<1> k_; + + // Linear combination of particular solutions: Rank-local stage + if(n_good_solutions != 0) { + A_k(k_) << triqs::clef::sum(c(j_) * A_block(j_, k_), j_ = my_j_range); + Ap_k(k_) << triqs::clef::sum(c(j_) * Ap_block(j_, k_), j_ = my_j_range); + App_k(k_) << triqs::clef::sum(c(j_) * App_block(j_, k_), j_ = my_j_range); + } else { + A_k() = 0; + Ap_k() = 0; + App_k() = 0; + } + + // Linear combination of particular solutions: Inter-rank stage + A_k = mpi::all_reduce(A_k, comm); + Ap_k = mpi::all_reduce(Ap_k, comm); + App_k = mpi::all_reduce(App_k, comm); + + // Update Q + for(int k : range(A_k.size())) { + if(A_k(k) < 0) + Q_k(k) = p.amp_penalty_max; + else + Q_k(k) /= p.amp_penalty_divisor; + } + + // Update D + for(int k : range(Ap_k.size())) { + double val = std::abs(Ap_k(k)); + if(D_k(k) * val > DD) D_k(k) = DD / val; + } + + // Update B + for(int k : range(App_k.size())) { + double val = std::abs(App_k(k)); + if(B_k(k) * val > DD) B_k(k) = DD / val; + } + + double diff = std::abs(sum(c) - 1); + if(p.verbosity >= 1) { + double sum_abs = sum(abs(c)); + mpi_cout(comm) << " End of iteration " << (iter + 1) + << ": |sum(c) - 1| = " << diff + << ", sum(abs(c)) = " << sum_abs << std::endl; + } + + if(diff < convergence_tol) { + if(p.verbosity >= 1) mpi_cout(comm) << "Convergence reached" << std::endl; + break; + } + } + + if(p.verbosity >= 1 && iter == p.max_iter) + mpi_cout(comm) << "Maximum number of iterations reached" << std::endl; + + return c; +} + +template +std::vector som_core::compute_final_solution_cc_impl( + final_solution_cc_parameters_t const& p) { + + // List of real frequency points + auto e_k_list = std::visit( + [](auto const& m) { + array res(m.size()); + for(auto const& [k, e] : enumerate(m)) res(k) = double(e); + return res; + }, + p.refreq_mesh); + + for(int n : range(data.size())) { + auto& d = data[n]; + + if(p.verbosity > 0) + mpi_cout(comm) << "Using CC protocol to construct the final solution " + "for observable component [" + << n << "," << n << "]" << std::endl; + + // Select good solutions + std::vector good_solution_indices; + good_solution_indices.reserve(d.particular_solutions.size()); + double chi_min = std::sqrt(d.objf_min); + for(auto const& [s_index, s] : enumerate(d.particular_solutions)) { + double chi = std::sqrt(s.second); + if(chi / chi_min <= p.good_chi_rel && chi <= p.good_chi_abs) + good_solution_indices.emplace_back(s_index); + } + + if(p.verbosity > 1) + mpi_cout(comm) << "Selected " << good_solution_indices.size() + << " good solutions" << std::endl; + + global_index_map index_map(comm, good_solution_indices.size()); + + // Convergence tolerance for CC iterations + using mesh_t = typename KernelType::mesh_type; + double convergence_tol = + min_element(abs(d.get_error_bars() / d.get_rhs())); + if(p.verbosity > 1) + mpi_cout(comm) << "Convergence tolerance is " << convergence_tol + << std::endl; + + // Compute coefficients c_j using CC optimization protocol + auto c = cc_protocol(d.particular_solutions, + good_solution_indices, + e_k_list, + convergence_tol, + p, + index_map, + comm); + + if(p.verbosity >= 1) { + mpi_cout(comm) + << "Forming the resulting linear combination of particular solutions" + << std::endl; + if(p.verbosity >= 2) + mpi_cout(comm) << "Coefficients of the linear combination: " << c + << std::endl; + } + + configuration sol_sum(ci); + int const j_range_start = index_map.range_start(comm.rank()); + for(auto const& [j, s_index] : enumerate(good_solution_indices)) { + sol_sum += c(j_range_start + j) * d.particular_solutions[s_index].first; + } + d.final_solution = mpi::all_reduce(sol_sum, comm); + } + + return compute_objf_final(); +} + +std::vector +som_core::compute_final_solution_cc(final_solution_cc_parameters_t const& p) { + p.validate(); + +#define IMPL_CASE(r, okmk) \ + case(kernel_id(BOOST_PP_SEQ_ELEM(0, okmk), BOOST_PP_SEQ_ELEM(1, okmk){})): \ + return compute_final_solution_cc_impl>(p); + + SELECT_KERNEL(IMPL_CASE, som_core::compute_final_solution_cc()) +#undef IMPL_CASE +} + ///////////////////// // som_core::run() // ///////////////////// diff --git a/c++/som/som_core/parameters.hpp b/c++/som/som_core/parameters.hpp index 491723ae4..a92c1ed3c 100644 --- a/c++/som/som_core/parameters.hpp +++ b/c++/som/som_core/parameters.hpp @@ -20,7 +20,13 @@ ******************************************************************************/ #pragma once +#include #include +#include + +#include + +#include #include @@ -98,4 +104,68 @@ struct accumulate_parameters_t : public worker_parameters_t { void validate(observable_kind kind) const; }; +struct final_solution_cc_parameters_t { + + /// Grid of energy points used in derivative regularization procedure. + std::variant, + triqs::arrays::array> + refreq_mesh; + + /// Verbosity level (max level - 3). + /// default: 2 on MPI rank 0, 0 otherwise. + int verbosity = + ((mpi::communicator().rank() == 0) ? 2 : 0); // silence the slave nodes + + /// Maximal ratio :math:`\chi/\chi_\mathrm{min}` for a particular solution to + /// be selected. This criterion must be fulfilled together with the one set + /// by `good_chi_abs`. + double good_chi_rel = 2.0; + + /// Maximal value of :math:`\chi` for a particular solution to be selected. + /// This criterion must be fulfilled together with the one set + /// by `good_chi_rel`. + double good_chi_abs = HUGE_VAL; + + /// Default model of the spectral function evaluated at + /// energy points of `refreq_mesh`. + triqs::arrays::array default_model; + + /// Weights determining how much deviations from `default_model` are penalized + /// at each energy point of `refreq_mesh`. + triqs::arrays::array default_model_weights; + + /// Maximum allowed number of parameter adjustment iterations. + int max_iter = 20; + + /// Coefficient of the term that enforces the unity sum constraint. + double unity_sum_coeff = 1e6; + + /// Maximum value of the regularization parameter that penalizes + /// negative values of the spectral function. + double amp_penalty_max = 1e6; + + /// Divisor used to reduce the regularization parameter that penalizes + /// negative values of the spectral function. + double amp_penalty_divisor = 100; + + /// Initial value of the regularization parameters that penalize large + /// derivatives of the solution. + double der_penalty_init = 1.0; + + /// Coefficient used to increase the regularization parameters that penalize + /// large derivatives of the solution. + double der_penalty_coeff = 2.0; + + final_solution_cc_parameters_t() = default; + explicit final_solution_cc_parameters_t( + triqs::gfs::gf_mesh refreq_mesh) + : refreq_mesh(std::move(refreq_mesh)) {} + explicit final_solution_cc_parameters_t( + triqs::arrays::array refreq_mesh) + : refreq_mesh(std::move(refreq_mesh)) {} + + // Validate values of parameters + void validate() const; +}; + } // namespace som diff --git a/c++/som/som_core/som_core.hpp b/c++/som/som_core/som_core.hpp index 8db082cc3..53f6fb4c0 100644 --- a/c++/som/som_core/som_core.hpp +++ b/c++/som/som_core/som_core.hpp @@ -36,6 +36,7 @@ #include "parameters.hpp" #include +#include #include namespace som { @@ -56,10 +57,7 @@ class som_core { observable_kind kind; // Mesh of the input functions - std::variant, - triqs::gfs::gf_mesh, - triqs::gfs::gf_mesh> - mesh; + mesh_variant_t mesh; using histogram_t = triqs::statistics::histogram; @@ -140,6 +138,11 @@ class som_core { std::vector compute_objf_final(); template std::vector compute_objf_final_impl(); + // Implementation details of compute_final_solution_cc_impl() + template + std::vector compute_final_solution_cc_impl( + final_solution_cc_parameters_t const& params); + public: /// Construct on imaginary-time quantities som_core(triqs::gfs::gf_const_view g_tau, @@ -170,6 +173,12 @@ class som_core { std::vector compute_final_solution(double good_chi_rel = 2.0, double good_chi_abs = HUGE_VAL); + /// Compute the final solution using the SOCC protocol + TRIQS_WRAP_ARG_AS_DICT + std::vector compute_final_solution_cc( + final_solution_cc_parameters_t const& p + ); + /// Set of parameters used in the last call to accumulate() accumulate_parameters_t get_last_accumulate_parameters() const { return params; diff --git a/doc/notes/final_solution_cc.nb b/doc/notes/final_solution_cc.nb new file mode 100644 index 000000000..ed26f7c06 --- /dev/null +++ b/doc/notes/final_solution_cc.nb @@ -0,0 +1,638 @@ +(* Content-type: application/vnd.wolfram.mathematica *) + +(*** Wolfram Notebook File ***) +(* http://www.wolfram.com/nb *) + +(* CreatedBy='Mathematica 12.1' *) + +(*CacheID: 234*) +(* Internal cache information: +NotebookFileLineBreakTest +NotebookFileLineBreakTest +NotebookDataPosition[ 158, 7] +NotebookDataLength[ 21363, 630] +NotebookOptionsPosition[ 18775, 571] +NotebookOutlinePosition[ 19170, 587] +CellTagsIndexPosition[ 19127, 584] +WindowFrame->Normal*) + +(* Beginning of Notebook Content *) +Notebook[{ + +Cell[CellGroupData[{ +Cell["Construction of the final solution in SOCC", "Title", + CellChangeTimes->{{3.854528504995222*^9, + 3.8545285219388123`*^9}},ExpressionUUID->"ef5e19c9-dbcf-417f-960f-\ +3b951d31dd8a"], + +Cell[BoxData[{ + RowBox[{ + RowBox[{"J", "=", "5"}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{"Afin", "[", "z_", "]"}], ":=", + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + RowBox[{"c", "[", "j", "]"}], + RowBox[{"A", "[", + RowBox[{"j", ",", "z"}], "]"}]}], ",", + RowBox[{"{", + RowBox[{"j", ",", "1", ",", "J"}], "}"}]}], "]"}]}], ";"}], + "\[IndentingNewLine]"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{"Contrib", "[", + RowBox[{"O_", ",", "b_", ",", "const_"}], "]"}], ":=", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + RowBox[{"c", "[", "j1", "]"}], + RowBox[{"O", "[", + RowBox[{"[", + RowBox[{"j1", ",", "j2"}], "]"}], "]"}], + RowBox[{"c", "[", "j2", "]"}]}], ",", + RowBox[{"{", + RowBox[{"j1", ",", "1", ",", "J"}], "}"}], ",", + RowBox[{"{", + RowBox[{"j2", ",", "1", ",", "J"}], "}"}]}], "]"}], "-", + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + RowBox[{"b", "[", + RowBox[{"[", "j", "]"}], "]"}], + RowBox[{"c", "[", "j", "]"}]}], ",", + RowBox[{"{", + RowBox[{"j", ",", "1", ",", "J"}], "}"}]}], "]"}], "+", "const"}]}], + ";"}]}], "Input", + CellChangeTimes->{{3.8545282403819*^9, 3.8545283739082623`*^9}, + 3.854528409288911*^9, 3.8545289454930983`*^9, {3.854529528268736*^9, + 3.854529567782641*^9}, {3.854529751002017*^9, 3.854529751381201*^9}}, + CellLabel->"In[1]:=",ExpressionUUID->"915fc75b-64ed-43a5-bc09-92c53b9cbae5"], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"Unity", "-", + RowBox[{"sum", " ", "constraint", " ", + SubscriptBox["O", "2"]}]}]], "Section", + CellChangeTimes->{{3.8545285450250473`*^9, 3.854528556762639*^9}, + 3.854528754445278*^9},ExpressionUUID->"e0b1030e-1f40-436e-a2a7-\ +4ea282a5b8db"], + +Cell[CellGroupData[{ + +Cell[BoxData[{ + RowBox[{ + RowBox[{"O2", "=", + RowBox[{"ConstantArray", "[", + RowBox[{"U", ",", + RowBox[{"{", + RowBox[{"J", ",", "J"}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"b2", "=", + RowBox[{"2", + RowBox[{"ConstantArray", "[", + RowBox[{"U", ",", + RowBox[{"{", "J", "}"}]}], "]"}]}]}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{"const2", "=", "U"}], ";"}], + "\[IndentingNewLine]"}], "\[IndentingNewLine]", + RowBox[{"FullSimplify", "[", + RowBox[{ + RowBox[{"U", + SuperscriptBox[ + RowBox[{"(", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{"c", "[", "j", "]"}], ",", + RowBox[{"{", + RowBox[{"j", ",", "1", ",", "J"}], "}"}]}], "]"}], "-", "1"}], + ")"}], "2"]}], "\[Equal]", + RowBox[{"Contrib", "[", + RowBox[{"O2", ",", "b2", ",", "const2"}], "]"}]}], "]"}]}], "Input", + CellChangeTimes->{{3.8545280948075533`*^9, 3.854528223008564*^9}, { + 3.854528258861807*^9, 3.85452828348925*^9}, {3.854528329507614*^9, + 3.854528361580317*^9}, {3.854528417693653*^9, 3.854528446459695*^9}, { + 3.854528559595565*^9, 3.854528560071701*^9}, 3.854528801390521*^9, { + 3.854530704100696*^9, 3.854530705850648*^9}}, + CellLabel->"In[4]:=",ExpressionUUID->"cce08e03-98a1-4309-8648-d5e9e91e9743"], + +Cell[BoxData["True"], "Output", + CellChangeTimes->{ + 3.854528115843223*^9, 3.854528155722949*^9, 3.85452836302817*^9, { + 3.854528414615176*^9, 3.854528436347947*^9}, 3.854528565906988*^9, + 3.85452888707019*^9, 3.854528957517851*^9, 3.854529843685512*^9, + 3.854530057879571*^9, 3.854530890464913*^9}, + CellLabel->"Out[7]=",ExpressionUUID->"dd62c5ee-5f27-4442-a231-c1e183977e22"] +}, Open ]] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[TextData[{ + "Penalty for large deviations from the equal weight superposition ", + Cell[BoxData[ + SubscriptBox["O", "3"]], "Section", + CellChangeTimes->{{3.8545285450250473`*^9, 3.854528556762639*^9}, + 3.854528754445278*^9},ExpressionUUID-> + "45df1657-1a08-41a1-9bba-ddb6ec4c84b4"] +}], "Section", + CellChangeTimes->{{3.854528730812148*^9, 3.8545287464124937`*^9}, { + 3.855401384301261*^9, 3.85540138467835*^9}, {3.85540150073554*^9, + 3.8554015035687943`*^9}},ExpressionUUID->"f9ae31b2-466b-474b-a277-\ +71a7bf590b84"], + +Cell[CellGroupData[{ + +Cell[BoxData[{ + RowBox[{ + RowBox[{"O3", "=", + RowBox[{"IdentityMatrix", "[", "J", "]"}]}], + ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"b3", "=", + RowBox[{"2", + RowBox[{"ConstantArray", "[", + RowBox[{ + RowBox[{"1", "/", "J"}], ",", + RowBox[{"{", "J", "}"}]}], "]"}]}]}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{"const3", "=", + RowBox[{"1", "/", "J"}]}], ";"}], + "\[IndentingNewLine]"}], "\[IndentingNewLine]", + RowBox[{"FullSimplify", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + SuperscriptBox[ + RowBox[{"(", + RowBox[{ + RowBox[{"c", "[", "j", "]"}], "-", + RowBox[{"1", "/", "J"}]}], ")"}], "2"], ",", + RowBox[{"{", + RowBox[{"j", ",", "1", ",", "J"}], "}"}]}], "]"}], "\[Equal]", + RowBox[{"Contrib", "[", + RowBox[{"O3", ",", "b3", ",", "const3"}], "]"}]}], "]"}]}], "Input", + CellChangeTimes->{{3.854528767387372*^9, 3.854528841361432*^9}, { + 3.8545306993258266`*^9, 3.854530701361587*^9}}, + CellLabel->"In[8]:=",ExpressionUUID->"6dfc5ea2-74ea-43a9-8754-b3fb6ada0e31"], + +Cell[BoxData["True"], "Output", + CellChangeTimes->{{3.85452884235756*^9, 3.854528887245104*^9}, + 3.854528957539507*^9, 3.8545298436978407`*^9, 3.854530057895117*^9, + 3.8545308904791117`*^9}, + CellLabel->"Out[11]=",ExpressionUUID->"a221b550-1076-41c2-b490-42af5d291fff"] +}, Open ]] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[TextData[{ + "Amplitude penalty ", + Cell[BoxData[ + SubscriptBox["O", "4"]], "Section", + CellChangeTimes->{{3.8545285450250473`*^9, 3.854528556762639*^9}, + 3.854528754445278*^9},ExpressionUUID-> + "b0fa5e1b-6661-4f30-a5d7-f48b49bf8678"] +}], "Section", + CellChangeTimes->{{3.8545294526826143`*^9, 3.85452946239362*^9}, { + 3.855401513270331*^9, + 3.85540151595057*^9}},ExpressionUUID->"a609d4dc-6646-4e0f-8b76-\ +164552d62ee1"], + +Cell[CellGroupData[{ + +Cell[BoxData[{ + RowBox[{ + RowBox[{"K", "=", "10"}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"O4", "=", + RowBox[{"Table", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + RowBox[{"Q", "[", "k", "]"}], + RowBox[{"A", "[", + RowBox[{"j1", ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}], + RowBox[{"A", "[", + RowBox[{"j2", ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}], ",", + RowBox[{"{", + RowBox[{"j1", ",", "1", ",", "J"}], "}"}], ",", + RowBox[{"{", + RowBox[{"j2", ",", "1", ",", "J"}], "}"}]}], "]"}]}], + ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"b4", "=", + RowBox[{"2", + RowBox[{"ConstantArray", "[", + RowBox[{"0", ",", + RowBox[{"{", "J", "}"}]}], "]"}]}]}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{"const4", "=", "0"}], ";"}], + "\[IndentingNewLine]"}], "\[IndentingNewLine]", + RowBox[{"FullSimplify", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + RowBox[{"Q", "[", "k", "]"}], + SuperscriptBox[ + RowBox[{"Afin", "[", + RowBox[{"z", "[", "k", "]"}], "]"}], "2"]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}], "\[Equal]", + RowBox[{"Contrib", "[", + RowBox[{"O4", ",", "b4", ",", "const4"}], "]"}]}], + "]"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"Unset", "[", "K", "]"}], ";"}]}], "Input", + CellChangeTimes->{{3.8545294657323008`*^9, 3.854529523940433*^9}, { + 3.854529576105775*^9, 3.854529714619088*^9}, {3.85452976439645*^9, + 3.8545298364566*^9}, {3.854529871932619*^9, 3.8545298765514383`*^9}, { + 3.854530029393412*^9, 3.8545300641326313`*^9}, {3.8545306299090967`*^9, + 3.8545306320916348`*^9}, 3.8545307165053043`*^9, {3.854530809123551*^9, + 3.854530812099329*^9}}, + CellLabel->"In[12]:=",ExpressionUUID->"268d40dc-ed6f-4397-9764-3b26d6d7bd38"], + +Cell[BoxData["True"], "Output", + CellChangeTimes->{ + 3.854529580881155*^9, {3.854529699095257*^9, 3.854529722260611*^9}, { + 3.8545297612857227`*^9, 3.8545297651077127`*^9}, {3.8545297979115553`*^9, + 3.854529832153233*^9}, {3.854529862601556*^9, 3.854529876806319*^9}, { + 3.854530036078031*^9, 3.854530058334656*^9}, 3.854530890817486*^9}, + CellLabel->"Out[16]=",ExpressionUUID->"86af3eb0-9459-437d-8118-54cc26c18484"] +}, Open ]] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[TextData[{ + "Deviation from a default model ", + Cell[BoxData[ + SubscriptBox["O", "5"]], "Section", + CellChangeTimes->{{3.8545285450250473`*^9, 3.854528556762639*^9}, + 3.854528754445278*^9},ExpressionUUID-> + "660835d1-f96d-47ef-833c-72f3042f960b"] +}], "Section", + CellChangeTimes->{{3.8545302902819653`*^9, 3.854530308610468*^9}, { + 3.855401520855886*^9, + 3.855401526412497*^9}},ExpressionUUID->"9dc645a4-9545-48e3-b4e4-\ +81c0f81c4231"], + +Cell[CellGroupData[{ + +Cell[BoxData[{ + RowBox[{ + RowBox[{"K", "=", "10"}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"O5", "=", + RowBox[{"Table", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + RowBox[{"T", "[", "k", "]"}], + RowBox[{"A", "[", + RowBox[{"j1", ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}], + RowBox[{"A", "[", + RowBox[{"j2", ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}], ",", + RowBox[{"{", + RowBox[{"j1", ",", "1", ",", "J"}], "}"}], ",", + RowBox[{"{", + RowBox[{"j2", ",", "1", ",", "J"}], "}"}]}], "]"}]}], + ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"b5", "=", + RowBox[{"2", + RowBox[{"Table", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + RowBox[{"T", "[", "k", "]"}], + RowBox[{"AT", "[", + RowBox[{"z", "[", "k", "]"}], "]"}], + RowBox[{"A", "[", + RowBox[{"j", ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}], ",", + RowBox[{"{", + RowBox[{"j", ",", "1", ",", "J"}], "}"}]}], "]"}]}]}], + ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{"const5", "=", + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + RowBox[{"T", "[", "k", "]"}], + SuperscriptBox[ + RowBox[{"AT", "[", + RowBox[{"z", "[", "k", "]"}], "]"}], "2"]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}]}], ";"}], + "\[IndentingNewLine]"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"FullSimplify", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + RowBox[{"T", "[", "k", "]"}], + SuperscriptBox[ + RowBox[{"(", + RowBox[{ + RowBox[{"Afin", "[", + RowBox[{"z", "[", "k", "]"}], "]"}], "-", + RowBox[{"AT", "[", + RowBox[{"z", "[", "k", "]"}], "]"}]}], ")"}], "2"]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}], "\[Equal]", + RowBox[{"Contrib", "[", + RowBox[{"O5", ",", "b5", ",", "const5"}], "]"}]}], "]"}], + "\[IndentingNewLine]"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"Unset", "[", "K", "]"}], ";"}]}], "Input", + CellChangeTimes->{{3.85453059595428*^9, 3.85453083806176*^9}, { + 3.854530871017926*^9, 3.854530881495779*^9}}, + CellLabel->"In[18]:=",ExpressionUUID->"0da8566a-98bf-4acd-9b03-702e21f2ecd5"], + +Cell[BoxData["True"], "Output", + CellChangeTimes->{{3.854530835486794*^9, 3.854530891526266*^9}}, + CellLabel->"Out[22]=",ExpressionUUID->"100c9a79-4a6b-4f67-a933-10c595ddaf45"] +}, Open ]] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[TextData[{ + "First derivative penalty ", + Cell[BoxData[ + SubscriptBox["O", "1"]], "Section", + CellChangeTimes->{{3.8545285450250473`*^9, 3.854528556762639*^9}, + 3.854528754445278*^9},ExpressionUUID-> + "3985c5fc-9588-4e18-888f-1a8c55454f6a"] +}], "Section", + CellChangeTimes->{{3.854531307927421*^9, 3.854531334498663*^9}, { + 3.85540153499741*^9, + 3.855401538836618*^9}},ExpressionUUID->"2fe3ee8d-0645-4ccd-b9c9-\ +aa22085adbd9"], + +Cell[CellGroupData[{ + +Cell[BoxData[{ + RowBox[{ + RowBox[{"K", "=", "10"}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"O11", "=", + RowBox[{"Table", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + SuperscriptBox[ + RowBox[{"DD", "[", "k", "]"}], "2"], + RowBox[{"D", "[", + RowBox[{ + RowBox[{"A", "[", + RowBox[{"j1", ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}], ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}], + RowBox[{"D", "[", + RowBox[{ + RowBox[{"A", "[", + RowBox[{"j2", ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}], ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}], ",", + RowBox[{"{", + RowBox[{"j1", ",", "1", ",", "J"}], "}"}], ",", + RowBox[{"{", + RowBox[{"j2", ",", "1", ",", "J"}], "}"}]}], "]"}]}], + ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"b11", "=", + RowBox[{"2", + RowBox[{"ConstantArray", "[", + RowBox[{"0", ",", + RowBox[{"{", "J", "}"}]}], "]"}]}]}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{"const11", "=", "0"}], ";"}], + "\[IndentingNewLine]"}], "\[IndentingNewLine]", + RowBox[{"FullSimplify", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + SuperscriptBox[ + RowBox[{"DD", "[", "k", "]"}], "2"], + SuperscriptBox[ + RowBox[{"(", + RowBox[{"D", "[", + RowBox[{ + RowBox[{"Afin", "[", + RowBox[{"z", "[", "k", "]"}], "]"}], ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}], ")"}], "2"]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}], "\[Equal]", + RowBox[{"Contrib", "[", + RowBox[{"O11", ",", "b11", ",", "const11"}], "]"}]}], + "]"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"Unset", "[", "K", "]"}], ";"}]}], "Input", + CellChangeTimes->{{3.854531353047722*^9, 3.8545313858509827`*^9}, { + 3.854531442554275*^9, 3.8545315180976954`*^9}}, + CellLabel->"In[40]:=",ExpressionUUID->"1fad85e4-a4a1-4024-9813-da8b4a079727"], + +Cell[BoxData["True"], "Output", + CellChangeTimes->{3.85453147611681*^9, 3.8545315196841908`*^9}, + CellLabel->"Out[44]=",ExpressionUUID->"efa97a37-8a21-4e29-bd7f-127a87bb370e"] +}, Open ]] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[TextData[{ + "Second derivative penalty ", + Cell[BoxData[ + SubscriptBox["O", "1"]], "Section", + CellChangeTimes->{{3.8545285450250473`*^9, 3.854528556762639*^9}, + 3.854528754445278*^9},ExpressionUUID-> + "3dad49c2-3022-4776-8d66-8edbde5a5183"] +}], "Section", + CellChangeTimes->{{3.854531307927421*^9, 3.854531334498663*^9}, { + 3.8545315258834457`*^9, 3.854531527073267*^9}, {3.855401541433428*^9, + 3.855401544882609*^9}},ExpressionUUID->"fe5580d4-3c81-436a-aa57-\ +44a779ccb341"], + +Cell[CellGroupData[{ + +Cell[BoxData[{ + RowBox[{ + RowBox[{"K", "=", "10"}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"O12", "=", + RowBox[{"Table", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + SuperscriptBox[ + RowBox[{"B", "[", "k", "]"}], "2"], + RowBox[{"D", "[", + RowBox[{ + RowBox[{"A", "[", + RowBox[{"j1", ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}], ",", + RowBox[{"{", + RowBox[{ + RowBox[{"z", "[", "k", "]"}], ",", "2"}], "}"}]}], "]"}], + RowBox[{"D", "[", + RowBox[{ + RowBox[{"A", "[", + RowBox[{"j2", ",", + RowBox[{"z", "[", "k", "]"}]}], "]"}], ",", + RowBox[{"{", + RowBox[{ + RowBox[{"z", "[", "k", "]"}], ",", "2"}], "}"}]}], "]"}]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}], ",", + RowBox[{"{", + RowBox[{"j1", ",", "1", ",", "J"}], "}"}], ",", + RowBox[{"{", + RowBox[{"j2", ",", "1", ",", "J"}], "}"}]}], "]"}]}], + ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"b12", "=", + RowBox[{"2", + RowBox[{"ConstantArray", "[", + RowBox[{"0", ",", + RowBox[{"{", "J", "}"}]}], "]"}]}]}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{"const12", "=", "0"}], ";"}], + "\[IndentingNewLine]"}], "\[IndentingNewLine]", + RowBox[{"FullSimplify", "[", + RowBox[{ + RowBox[{"Sum", "[", + RowBox[{ + RowBox[{ + SuperscriptBox[ + RowBox[{"B", "[", "k", "]"}], "2"], + SuperscriptBox[ + RowBox[{"(", + RowBox[{"D", "[", + RowBox[{ + RowBox[{"Afin", "[", + RowBox[{"z", "[", "k", "]"}], "]"}], ",", + RowBox[{"{", + RowBox[{ + RowBox[{"z", "[", "k", "]"}], ",", "2"}], "}"}]}], "]"}], ")"}], + "2"]}], ",", + RowBox[{"{", + RowBox[{"k", ",", "1", ",", "K"}], "}"}]}], "]"}], "\[Equal]", + RowBox[{"Contrib", "[", + RowBox[{"O12", ",", "b12", ",", "const12"}], "]"}]}], + "]"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"Unset", "[", "K", "]"}], ";"}]}], "Input", + CellChangeTimes->{{3.854531353047722*^9, 3.8545313858509827`*^9}, { + 3.854531442554275*^9, 3.854531622853261*^9}, 3.854531723563972*^9, { + 3.854531753816245*^9, 3.8545317697801027`*^9}}, + CellLabel->"In[71]:=",ExpressionUUID->"297ce869-3ff4-401f-87ba-f185ce565bef"], + +Cell[BoxData["True"], "Output", + CellChangeTimes->{ + 3.85453162392544*^9, {3.854531720938303*^9, 3.854531724268659*^9}, { + 3.854531766408155*^9, 3.854531770883284*^9}}, + CellLabel->"Out[75]=",ExpressionUUID->"dff29bc5-4545-4ecb-81b6-9238f0f8c2fe"] +}, Open ]] +}, Open ]] +}, Open ]] +}, +WindowSize->{904.5, 658.5}, +WindowMargins->{{195, Automatic}, {30, Automatic}}, +FrontEndVersion->"12.1 for Linux x86 (64-bit) (June 19, 2020)", +StyleDefinitions->"Default.nb", +ExpressionUUID->"6b574963-b77f-42f2-84cd-8f8ddc710d81" +] +(* End of Notebook Content *) + +(* Internal cache information *) +(*CellTagsOutline +CellTagsIndex->{} +*) +(*CellTagsIndex +CellTagsIndex->{} +*) +(*NotebookFileOutline +Notebook[{ +Cell[CellGroupData[{ +Cell[580, 22, 187, 3, 98, "Title",ExpressionUUID->"ef5e19c9-dbcf-417f-960f-3b951d31dd8a"], +Cell[770, 27, 1558, 45, 113, "Input",ExpressionUUID->"915fc75b-64ed-43a5-bc09-92c53b9cbae5"], +Cell[CellGroupData[{ +Cell[2353, 76, 273, 6, 67, "Section",ExpressionUUID->"e0b1030e-1f40-436e-a2a7-4ea282a5b8db"], +Cell[CellGroupData[{ +Cell[2651, 86, 1346, 36, 114, "Input",ExpressionUUID->"cce08e03-98a1-4309-8648-d5e9e91e9743"], +Cell[4000, 124, 386, 6, 33, "Output",ExpressionUUID->"dd62c5ee-5f27-4442-a231-c1e183977e22"] +}, Open ]] +}, Open ]], +Cell[CellGroupData[{ +Cell[4435, 136, 530, 11, 72, "Section",ExpressionUUID->"f9ae31b2-466b-474b-a277-71a7bf590b84"], +Cell[CellGroupData[{ +Cell[4990, 151, 1095, 32, 114, "Input",ExpressionUUID->"6dfc5ea2-74ea-43a9-8754-b3fb6ada0e31"], +Cell[6088, 185, 275, 4, 33, "Output",ExpressionUUID->"a221b550-1076-41c2-b490-42af5d291fff"] +}, Open ]] +}, Open ]], +Cell[CellGroupData[{ +Cell[6412, 195, 435, 11, 72, "Section",ExpressionUUID->"a609d4dc-6646-4e0f-8b76-164552d62ee1"], +Cell[CellGroupData[{ +Cell[6872, 210, 2020, 56, 156, "Input",ExpressionUUID->"268d40dc-ed6f-4397-9764-3b26d6d7bd38"], +Cell[8895, 268, 425, 6, 33, "Output",ExpressionUUID->"86af3eb0-9459-437d-8118-54cc26c18484"] +}, Open ]] +}, Open ]], +Cell[CellGroupData[{ +Cell[9369, 280, 450, 11, 72, "Section",ExpressionUUID->"9dc645a4-9545-48e3-b4e4-81c0f81c4231"], +Cell[CellGroupData[{ +Cell[9844, 295, 2605, 79, 178, "Input",ExpressionUUID->"0da8566a-98bf-4acd-9b03-702e21f2ecd5"], +Cell[12452, 376, 176, 2, 33, "Output",ExpressionUUID->"100c9a79-4a6b-4f67-a933-10c595ddaf45"] +}, Open ]] +}, Open ]], +Cell[CellGroupData[{ +Cell[12677, 384, 442, 11, 72, "Section",ExpressionUUID->"2fe3ee8d-0645-4ccd-b9c9-aa22085adbd9"], +Cell[CellGroupData[{ +Cell[13144, 399, 2160, 64, 158, "Input",ExpressionUUID->"1fad85e4-a4a1-4024-9813-da8b4a079727"], +Cell[15307, 465, 175, 2, 33, "Output",ExpressionUUID->"efa97a37-8a21-4e29-bd7f-127a87bb370e"] +}, Open ]] +}, Open ]], +Cell[CellGroupData[{ +Cell[15531, 473, 492, 11, 72, "Section",ExpressionUUID->"fe5580d4-3c81-436a-aa57-44a779ccb341"], +Cell[CellGroupData[{ +Cell[16048, 488, 2434, 72, 180, "Input",ExpressionUUID->"297ce869-3ff4-401f-87ba-f185ce565bef"], +Cell[18485, 562, 250, 4, 33, "Output",ExpressionUUID->"dff29bc5-4545-4ecb-81b6-9238f0f8c2fe"] +}, Open ]] +}, Open ]] +}, Open ]] +} +] +*) + diff --git a/python/som/som_core_desc.py b/python/som/som_core_desc.py index ecb09ea55..d6b95d728 100644 --- a/python/som/som_core_desc.py +++ b/python/som/som_core_desc.py @@ -36,6 +36,7 @@ #include #include #include +#include #include #include #include @@ -355,8 +356,6 @@ def add_worker_parameters(conv): module.add_converter(accumulate_params_conv) - - # # SomCore.accumulate() # @@ -410,6 +409,131 @@ def add_worker_parameters(conv): c.add_method("std::vector compute_final_solution(double good_chi_rel = 2.0, double good_chi_abs = HUGE_VAL)", doc = """Select particular solutions according to the standard SOM criterion and compute the final solution""") +# +# Converter for final_solution_cc_parameters_t +# + +compute_final_solution_cc_params_conv = converter_( + c_type = "som::final_solution_cc_parameters_t", + doc = """Arguments of SomCore.compute_final_solution_cc()""", +) + +compute_final_solution_cc_params_conv.add_member(c_name = "refreq_mesh", + c_type = "std::variant, triqs::arrays::array>", + initializer = """ """, + doc = """Grid of energy points used in derivative regularization procedure.""") + +compute_final_solution_cc_params_conv.add_member(c_name = "verbosity", + c_type = "int", + initializer = """ ((mpi::communicator().rank() == 0) ? 1 : 0) """, + doc = """Verbosity level (max level - 2).""") + +compute_final_solution_cc_params_conv.add_member(c_name = "good_chi_rel", + c_type = "double", + initializer = "2.0", + doc = """Maximal ratio :math:`\chi/\chi_\mathrm{min}` for a particular solution to be selected. +This criterion must be fulfilled together with the one set by `good_chi_abs`.""") + +compute_final_solution_cc_params_conv.add_member(c_name = "good_chi_abs", + c_type = "double", + initializer = "HUGE_VAL", + doc = """Maximal value of :math:`\chi` for a particular solution to be selected. +This criterion must be fulfilled together with the one set by `good_chi_rel`.""") + +compute_final_solution_cc_params_conv.add_member(c_name = "default_model", + c_type = "triqs::arrays::array", + initializer = """{}""", + doc = """Default model of the spectral function evaluated at energy points of `refreq_mesh`.""") + +compute_final_solution_cc_params_conv.add_member(c_name = "default_model_weights", + c_type = "triqs::arrays::array", + initializer = """{}""", + doc = """Weights determining how much deviations from `default_model` are penalized at each energy point of `refreq_mesh`.""") + +compute_final_solution_cc_params_conv.add_member(c_name = "max_iter", + c_type = "int", + initializer = """20""", + doc = """Maximum allowed number of parameter adjustment iterations.""") +compute_final_solution_cc_params_conv.add_member(c_name = "unity_sum_coeff", + c_type = "double", + initializer = """1e6""", + doc = """Coefficient of the term that enforces the unity sum constraint.""") + +compute_final_solution_cc_params_conv.add_member(c_name = "amp_penalty_max", + c_type = "double", + initializer = """1e6""", + doc = """Maximum value of the regularization parameter that penalizes negative values of the spectral function.""") + +compute_final_solution_cc_params_conv.add_member(c_name = "amp_penalty_divisor", + c_type = "double", + initializer = """100""", + doc = """Divisor used to reduce the regularization parameter that penalizes negative values of the spectral function.""") + +compute_final_solution_cc_params_conv.add_member(c_name = "der_penalty_init", + c_type = "double", + initializer = """1.0""", + doc = """Initial value of the regularization parameters that penalize large derivatives of the solution.""") + +compute_final_solution_cc_params_conv.add_member(c_name = "der_penalty_coeff", + c_type = "double", + initializer = """2.0""", + doc = """Coefficient used to increase the regularization parameters that penalize large derivatives of the solution.""") + +module.add_converter(compute_final_solution_cc_params_conv) + +# +# SomCore.compute_final_solution_cc() +# + +c.add_method("std::vector compute_final_solution_cc(**som::final_solution_cc_parameters_t)", + #release_GIL_and_enable_signal = True, FIXME + doc = f""" +Compute the final solution using the SOCC protocol +================================================== + +{docstring_params_header_main} + +{docstring_params_table_header} +| refreq_mesh | MeshReFreq or | -- | Grid of energy points :math:`\epsilon_k` used in CC regularization procedure. | +| | real 1D | | Either a TRIQS real-frequency mesh object or a strictly ordered list of points is allowed. | +| | array_like | | | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| verbosity | int | 1 on MPI rank 0, 0 otherwise. | Verbosity level (max level - 2). | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| good_chi_rel | float | 2.0 | Maximal ratio :math:`\chi/\chi_\mathrm{{min}}` for a particular solution to be selected. | +| | | | This criterion must be fulfilled together with the one set by `good_chi_abs`. | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| good_chi_abs | float | infinity | Maximal value of :math:`\chi` for a particular solution to be selected. | +| | | | This criterion must be fulfilled together with the one set by `good_chi_rel`. | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| default_model | Real 1D array | [] | Optional default model of the spectral function evaluated at energy points of `refreq_mesh` | +| | | | (:math:`A_T(\epsilon_k)`). | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| default_model_weights | Real 1D array | [] | Weights determining how much deviations from `default_model` are penalized at each energy point | +| | | | (:math:`T_k`). | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| max_iter | int | 20 | Maximum allowed number of parameter adjustment iterations. | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ + +{docstring_params_header_fine} + +{docstring_params_table_header} +| unity_sum_coeff | float | 1e6 | Coefficient of the term that enforces the unity sum constraint (:math:`\mathcal{{U}}`). | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| amp_penalty_max | float | 1e6 | Maximum value of the regularization parameter that penalizes negative values of | +| | | | the spectral function (:math:`\mathcal{{Q}}`). | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| amp_penalty_divisor | float | 100 | Divisor used to reduce the regularization parameter that penalizes negative values of | +| | | | the spectral function. | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| der_penalty_init | float | 1.0 | Initial value of the regularization parameters that penalize large derivatives of the solution | +| | | | (:math:`\mathcal{{D}}`). | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +| der_penalty_coeff | float | 2.0 | Coefficient used to increase the regularization parameters that penalize large derivatives of | +| | | | the solution (:math:`f`). | ++-----------------------+---------------+-------------------------------+-------------------------------------------------------------------------------------------------------+ +""") + # # SomCore.run() # diff --git a/test/c++/CMakeLists.txt b/test/c++/CMakeLists.txt index e9ccff964..192af7a57 100644 --- a/test/c++/CMakeLists.txt +++ b/test/c++/CMakeLists.txt @@ -93,7 +93,7 @@ foreach(test ${SIMPLE_TESTS}) endforeach() # MPI tests -set(MPI_TESTS spectral_stats.cpp) +set(MPI_TESTS global_index_map.cpp spectral_stats.cpp) foreach(test ${MPI_TESTS}) add_som_cpp_test(${test} 1 2 4) endforeach()