diff --git a/SeQuant/core/eval/result.hpp b/SeQuant/core/eval/result.hpp index 6f3cc48d5e..74b839d6ac 100644 --- a/SeQuant/core/eval/result.hpp +++ b/SeQuant/core/eval/result.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -355,17 +356,13 @@ auto biorthogonal_nns_project_ta(TA::DistArray const& arr, SEQUANT_ASSERT(bra_rank <= rank); size_t const ket_rank = rank - bra_rank; - if (rank <= 4) { - return arr; - } + // Residuals of rank 4 or less have no redundancy and don't require NNS + // projection + if (rank <= 4) return arr; using numeric_type = typename TA::DistArray::numeric_type; - size_t factorial_ket = 1; - for (size_t i = 2; i <= ket_rank; ++i) { - factorial_ket *= i; - } - numeric_type norm_factor = numeric_type(1) / numeric_type(factorial_ket); + const auto& nns_p_coeffs = nns_projection_weights(ket_rank); TA::DistArray result; @@ -375,43 +372,28 @@ auto biorthogonal_nns_project_ta(TA::DistArray const& arr, const auto lannot = ords_to_annot(perm); - auto process_permutations = [&lannot](const TA::DistArray& input_arr, - size_t range_rank, perm_t range_perm, - const std::string& other_annot, - bool is_bra) -> TA::DistArray { - if (range_rank <= 1) return input_arr; - TA::DistArray result; + if (ket_rank > 2 && !nns_p_coeffs.empty()) { + const auto bra_annot = bra_rank == 0 ? "" : ords_to_annot(bra_perm); - auto callback = [&]([[maybe_unused]] int parity) { - const auto range_annot = ords_to_annot(range_perm); - const auto annot = other_annot.empty() - ? range_annot - : (is_bra ? range_annot + "," + other_annot - : other_annot + "," + range_annot); + size_t num_perms = nns_p_coeffs.size(); + for (size_t perm_rank = 0; perm_rank < num_perms; ++perm_rank) { + perm_t permuted_ket = + compute_permuted_indices(ket_perm, perm_rank, ket_rank); + + numeric_type coeff = nns_p_coeffs[perm_rank]; + + const auto ket_annot = ords_to_annot(permuted_ket); + const auto annot = + bra_annot.empty() ? ket_annot : bra_annot + "," + ket_annot; - // ignore parity, all permutations get same coefficient - numeric_type p_ = 1; if (result.is_initialized()) { - result(lannot) += p_ * input_arr(annot); + result(lannot) += coeff * arr(annot); } else { - result(lannot) = p_ * input_arr(annot); + result(lannot) = coeff * arr(annot); } - }; - antisymmetric_permutation(ParticleRange{range_perm.begin(), range_rank}, - callback); - return result; - }; - - // identity term with coefficient +1 - result(lannot) = arr(lannot); - - // process only ket permutations with coefficient norm_factor - if (ket_rank > 1) { - const auto bra_annot = bra_rank == 0 ? "" : ords_to_annot(bra_perm); - auto ket_result = - process_permutations(arr, ket_rank, ket_perm, bra_annot, false); - - result(lannot) -= norm_factor * ket_result(lannot); + } + } else { + result(lannot) = arr(lannot); } TA::DistArray::wait_for_lazy_cleanup(result.world()); @@ -428,64 +410,52 @@ auto biorthogonal_nns_project_ta(TA::DistArray const& arr, template auto biorthogonal_nns_project_btas(btas::Tensor const& arr, size_t bra_rank) { - using ranges::views::concat; using ranges::views::iota; size_t const rank = arr.rank(); SEQUANT_ASSERT(bra_rank <= rank); size_t const ket_rank = rank - bra_rank; - if (rank <= 4) { - return arr; - } + // Residuals of rank 4 or less have no redundancy and don't require NNS + // projection + if (rank <= 4) return arr; using numeric_type = typename btas::Tensor::numeric_type; - size_t factorial_ket = 1; - for (size_t i = 2; i <= ket_rank; ++i) { - factorial_ket *= i; - } - numeric_type norm_factor = numeric_type(1) / numeric_type(factorial_ket); + const auto& nns_p_coeffs = nns_projection_weights(ket_rank); + btas::Tensor result; + + perm_t perm = iota(size_t{0}, rank) | ranges::to; perm_t bra_perm = iota(size_t{0}, bra_rank) | ranges::to; perm_t ket_perm = iota(bra_rank, rank) | ranges::to; - const auto lannot = iota(size_t{0}, rank) | ranges::to; - auto process_permutations = [&lannot](const btas::Tensor& input_arr, - size_t range_rank, perm_t range_perm, - const perm_t& other_perm, bool is_bra) { - if (range_rank <= 1) return input_arr; - btas::Tensor result{input_arr.range()}; - result.fill(0); + if (ket_rank > 2 && !nns_p_coeffs.empty()) { + bool result_initialized = false; - auto callback = [&]([[maybe_unused]] int parity) { - const auto annot = - is_bra ? concat(range_perm, other_perm) | ranges::to() - : concat(other_perm, range_perm) | ranges::to(); + size_t num_perms = nns_p_coeffs.size(); + for (size_t perm_rank = 0; perm_rank < num_perms; ++perm_rank) { + perm_t permuted_ket = + compute_permuted_indices(ket_perm, perm_rank, ket_rank); - // ignore parity, all permutations get same coefficient - numeric_type p_ = 1; - btas::Tensor temp; - btas::permute(input_arr, lannot, temp, annot); - btas::scal(p_, temp); - result += temp; - }; + numeric_type coeff = nns_p_coeffs[perm_rank]; - antisymmetric_permutation(ParticleRange{range_perm.begin(), range_rank}, - callback); - return result; - }; + perm_t annot = bra_perm; + annot.insert(annot.end(), permuted_ket.begin(), permuted_ket.end()); - // identity term with coefficient +1 - auto result = arr; + btas::Tensor temp; + btas::permute(arr, annot, temp, perm); + btas::scal(coeff, temp); - // process only ket permutations with coefficient norm_factor - if (ket_rank > 1) { - const auto bra_annot = bra_rank == 0 ? perm_t{} : bra_perm; - auto ket_result = - process_permutations(arr, ket_rank, ket_perm, bra_annot, false); + if (result_initialized) { + result += temp; + } else { + result = temp; + result_initialized = true; + } + } - btas::scal(norm_factor, ket_result); - result -= ket_result; + } else { + result = arr; } return result; @@ -903,8 +873,14 @@ class ResultTensorTA final : public Result { [[nodiscard]] ResultPtr biorthogonal_nns_project( size_t bra_rank) const override { - return eval_result( - biorthogonal_nns_project_ta(get(), bra_rank)); + if constexpr (std::is_integral_v) { + SEQUANT_ABORT( + "biorthogonal_nns_project is not supported for integral numeric " + "types"); + } else { + return eval_result( + biorthogonal_nns_project_ta(get(), bra_rank)); + } } private: @@ -1166,8 +1142,14 @@ class ResultTensorBTAS final : public Result { [[nodiscard]] ResultPtr biorthogonal_nns_project( [[maybe_unused]] size_t bra_rank) const override { - return eval_result>( - biorthogonal_nns_project_btas(get(), bra_rank)); + if constexpr (std::is_integral_v) { + SEQUANT_ABORT( + "biorthogonal_nns_project is not supported for integral numeric " + "types"); + } else { + return eval_result>( + biorthogonal_nns_project_btas(get(), bra_rank)); + } } private: diff --git a/SeQuant/domain/mbpt/biorthogonalization.cpp b/SeQuant/domain/mbpt/biorthogonalization.cpp index 9c372b6f2c..a850259f5d 100644 --- a/SeQuant/domain/mbpt/biorthogonalization.cpp +++ b/SeQuant/domain/mbpt/biorthogonalization.cpp @@ -9,6 +9,7 @@ #include #include +#include #include #include @@ -29,6 +30,108 @@ struct compare_first_less { using IndexPair = std::pair; using ParticlePairings = container::svector; +std::vector hardcoded_biorth_coeffs_first_row( + std::size_t n_particles) { + switch (n_particles) { + case 1: + return std::vector{ratio(1, 2)}; + + case 2: + return std::vector{ratio(1, 3), ratio(1, 6)}; + + case 3: + return std::vector{ratio(17, 120), ratio(-7, 120), + ratio(-1, 120), ratio(-1, 120), + ratio(-1, 120), ratio(-7, 120)}; + + case 4: + return std::vector{ + ratio(43, 840), ratio(-19, 1680), ratio(-19, 1680), + ratio(-1, 105), ratio(-19, 1680), ratio(-19, 1680), + ratio(13, 840), ratio(1, 120), ratio(-1, 105), + ratio(1, 120), ratio(-1, 105), ratio(-19, 1680), + ratio(-1, 105), ratio(1, 120), ratio(1, 120), + ratio(13, 840), ratio(-1, 105), ratio(-1, 105), + ratio(1, 120), ratio(-19, 1680), ratio(-19, 1680), + ratio(13, 840), ratio(-19, 1680), ratio(1, 120)}; + + case 5: + return std::vector{ + ratio(59, 3780), ratio(-5, 3024), ratio(-5, 3024), + ratio(-5, 3024), ratio(-31, 7560), ratio(-5, 3024), + ratio(-5, 3024), ratio(-23, 30240), ratio(19, 7560), + ratio(37, 15120), ratio(-5, 3024), ratio(-23, 30240), + ratio(-5, 3024), ratio(19, 7560), ratio(37, 15120), + ratio(-31, 7560), ratio(37, 15120), ratio(37, 15120), + ratio(-31, 7560), ratio(-5, 3024), ratio(-5, 3024), + ratio(-23, 30240), ratio(-23, 30240), ratio(-23, 30240), + ratio(-13, 7560), ratio(-5, 3024), ratio(-5, 3024), + ratio(19, 7560), ratio(-23, 30240), ratio(37, 15120), + ratio(19, 7560), ratio(-23, 30240), ratio(19, 7560), + ratio(-23, 30240), ratio(-13, 7560), ratio(37, 15120), + ratio(-13, 7560), ratio(-13, 7560), ratio(37, 15120), + ratio(-23, 30240), ratio(-31, 7560), ratio(-13, 7560), + ratio(37, 15120), ratio(37, 15120), ratio(19, 7560), + ratio(37, 15120), ratio(37, 15120), ratio(-13, 7560), + ratio(-13, 7560), ratio(-23, 30240), ratio(-31, 7560), + ratio(37, 15120), ratio(-31, 7560), ratio(37, 15120), + ratio(-5, 3024), ratio(-5, 3024), ratio(-23, 30240), + ratio(19, 7560), ratio(-5, 3024), ratio(37, 15120), + ratio(-31, 7560), ratio(37, 15120), ratio(37, 15120), + ratio(-13, 7560), ratio(19, 7560), ratio(37, 15120), + ratio(37, 15120), ratio(-13, 7560), ratio(-13, 7560), + ratio(-23, 30240), ratio(37, 15120), ratio(-13, 7560), + ratio(37, 15120), ratio(-13, 7560), ratio(-23, 30240), + ratio(19, 7560), ratio(-23, 30240), ratio(-23, 30240), + ratio(19, 7560), ratio(-13, 7560), ratio(-31, 7560), + ratio(37, 15120), ratio(-13, 7560), ratio(37, 15120), + ratio(19, 7560), ratio(-31, 7560), ratio(-31, 7560), + ratio(37, 15120), ratio(37, 15120), ratio(-5, 3024), + ratio(37, 15120), ratio(-13, 7560), ratio(37, 15120), + ratio(-13, 7560), ratio(-23, 30240), ratio(-5, 3024), + ratio(19, 7560), ratio(-23, 30240), ratio(-5, 3024), + ratio(37, 15120), ratio(-5, 3024), ratio(-23, 30240), + ratio(-23, 30240), ratio(-23, 30240), ratio(-13, 7560), + ratio(19, 7560), ratio(19, 7560), ratio(-23, 30240), + ratio(-23, 30240), ratio(-13, 7560), ratio(-5, 3024), + ratio(19, 7560), ratio(-5, 3024), ratio(-23, 30240), + ratio(37, 15120), ratio(37, 15120), ratio(-13, 7560), + ratio(-13, 7560), ratio(37, 15120), ratio(-23, 30240)}; + + default: + throw std::runtime_error( + "hardcoded biorthogonal coefficients only available for ranks 1-5, " + "requested rank is : " + + std::to_string(n_particles)); + } +} + +Eigen::Matrix +make_hardcoded_biorth_coeffs_matrix( + const std::vector& first_row, std::size_t n_particles) { + const auto n = first_row.size(); + Eigen::Matrix M(n, n); + + for (std::size_t row = 0; row < n; ++row) { + for (std::size_t col = 0; col < n; ++col) { + perm::Permutation row_perm = perm::unrank(n - 1 - row, n_particles); + perm::Permutation col_perm = perm::unrank(col, n_particles); + + col_perm->preMultiply(row_perm); + + std::size_t source_idx = perm::rank(col_perm, n_particles); + M(row, col) = first_row[source_idx]; + } + } + return M; +} + +Eigen::Matrix +hardcoded_biorth_coeffs_matrix(std::size_t n_particles) { + auto first_row = hardcoded_biorth_coeffs_first_row(n_particles); + return make_hardcoded_biorth_coeffs_matrix(first_row, n_particles); +} + ResultExpr biorthogonal_transform_copy(const ResultExpr& expr, double threshold) { container::svector wrapper = {expr.clone()}; @@ -309,7 +412,7 @@ void biorthogonal_transform(container::svector& result_exprs, // like R^{IJ}_{AB} and the index pairing of the result is what determines // the required symmetrization. Hence, the symmetrization operator must not // be changed when transforming from one representation into the other. - assert(std::all_of( + SEQUANT_ASSERT(std::all_of( result_exprs.begin(), result_exprs.end(), [](const ResultExpr& res) { bool found = false; res.expression()->visit( @@ -336,12 +439,7 @@ void biorthogonal_transform(container::svector& result_exprs, ranges::to>(); const std::size_t n_particles = externals.front().size(); - - Eigen::MatrixXd coefficients = compute_biorth_coeffs(n_particles, threshold); - auto num_perms = factorial(n_particles); - SEQUANT_ASSERT(num_perms == coefficients.rows()); - SEQUANT_ASSERT(num_perms == coefficients.cols()); auto original_exprs = result_exprs | ranges::views::transform([](const ResultExpr& res) { @@ -349,6 +447,59 @@ void biorthogonal_transform(container::svector& result_exprs, }) | ranges::to>(); + auto memoize = [](container::map, + std::optional>& cache, + std::mutex& mutex, std::condition_variable& cv, + std::pair key, + auto compute_fn) -> const T& { + { + std::unique_lock lock(mutex); + auto [it, inserted] = cache.try_emplace(key, std::nullopt); + if (!inserted) { + cv.wait(lock, [&] { return it->second.has_value(); }); + return it->second.value(); + } + } + + T result = compute_fn(); + + { + std::lock_guard lock(mutex); + cache[key] = std::move(result); + cv.notify_all(); + return cache[key].value(); + } + }; + + using HardcodedMatrix = + Eigen::Matrix; + using ComputedMatrix = Eigen::MatrixXd; + using CacheKey = std::pair; + + static std::mutex cache_mutex; + static std::condition_variable cache_cv; + static container::map> + hardcoded_cache; + static container::map> computed_cache; + + constexpr std::size_t max_rank_hardcoded_biorth_coeffs = 5; + CacheKey key{n_particles, threshold}; + + const HardcodedMatrix* hardcoded_coefficients = nullptr; + const ComputedMatrix* computed_coefficients = nullptr; + + if (n_particles <= max_rank_hardcoded_biorth_coeffs) { + hardcoded_coefficients = + &memoize(hardcoded_cache, cache_mutex, cache_cv, key, + [&] { return hardcoded_biorth_coeffs_matrix(n_particles); }); + } else { + computed_coefficients = + &memoize(computed_cache, cache_mutex, cache_cv, key, + [&] { return compute_biorth_coeffs(n_particles, threshold); }); + SEQUANT_ASSERT(num_perms == computed_coefficients->rows()); + SEQUANT_ASSERT(num_perms == computed_coefficients->cols()); + } + for (std::size_t i = 0; i < result_exprs.size(); ++i) { result_exprs.at(i).expression() = ex(0); perm::Permutation reference = perm::unrank(ranks.at(i), n_particles); @@ -358,9 +509,14 @@ void biorthogonal_transform(container::svector& result_exprs, perm::Permutation perm = perm::unrank(rank, n_particles); perm->postMultiply(reference); + sequant::rational coeff = + (n_particles <= max_rank_hardcoded_biorth_coeffs) + ? (*hardcoded_coefficients)(ranks.at(i), rank) + : to_rational((*computed_coefficients)(ranks.at(i), rank), + threshold); + result_exprs.at(i).expression() += - ex( - to_rational(coefficients(ranks.at(i), rank), threshold)) * + ex(coeff) * create_expr_for(externals.at(i), perm, externals, original_exprs); } @@ -389,5 +545,39 @@ ExprPtr biorthogonal_transform( return res.expression(); } +Eigen::MatrixXd compute_nns_p_matrix(std::size_t n_particles, + double threshold) { + auto perm_ovlp_mat = permutational_overlap_matrix(n_particles); + auto normalized_pinv = compute_biorth_coeffs(n_particles, threshold); + + auto nns = perm_ovlp_mat * normalized_pinv; + + return nns; +} + +std::vector compute_nns_p_coeffs(std::size_t n_particles, + double threshold) { + Eigen::MatrixXd nns_matrix = compute_nns_p_matrix(n_particles, threshold); + std::size_t num_perms = nns_matrix.rows(); + + std::vector coeffs; + coeffs.reserve(num_perms); + for (std::size_t i = 0; i < num_perms; ++i) { + coeffs.push_back(nns_matrix(num_perms - 1, i)); + } + return coeffs; +} + +container::svector compute_permuted_indices( + const container::svector& indices, size_t perm_rank, + size_t n_particles) { + perm::Permutation perm_obj = perm::unrank(perm_rank, n_particles); + + container::svector permuted_indices(n_particles); + for (size_t i = 0; i < n_particles; ++i) { + permuted_indices[i] = indices[perm_obj[i]]; + } + return permuted_indices; +} } // namespace sequant diff --git a/SeQuant/domain/mbpt/biorthogonalization.hpp b/SeQuant/domain/mbpt/biorthogonalization.hpp index 62c415dd86..b448a884ee 100644 --- a/SeQuant/domain/mbpt/biorthogonalization.hpp +++ b/SeQuant/domain/mbpt/biorthogonalization.hpp @@ -5,11 +5,39 @@ #include #include +#include + namespace sequant { namespace { static constexpr double default_biorth_threshold = 1e-12; -} +} // namespace + +// clang-format off +/// \brief Provides the first row of the biorthogonal coefficients matrix, +/// hardcoded from Mathematica to avoid numerical precision loss. +/// +/// The Myrvold-Ruskey unrank1 algorithm (doi.org/10.1016/S0020-0190(01)00141-7) +/// is used to order permutations, then the permutational overlap matrix M is +/// constructed with elements (-2)^{c} × (-1)^{n_particles}, where c is the +/// number of cycles in the relative permutation. +/// +/// The biorthogonal coefficients are obtained from the normalized pseudoinverse +/// of M: first compute M_pinv (the pseudoinverse), then normalize it by the +/// factor ((n_particles)!/rank(M)). +/// Finally, biorthogonal coefficients = normalized_M_pinv · e_1, +/// where e_1 is the first unit vector. +/// See [DOI 10.48550/ARXIV.1805.00565](https://doi.org/10.48550/ARXIV.1805.00565) +/// for more details. +/// +/// \param n_particles The rank of external index pairs +/// +/// \return Vector of rational coefficients representing the first row +/// +/// \throw std::runtime_error if n_particles is not in the range [1,5] +// clang-format on +std::vector hardcoded_biorth_coeffs_first_row( + std::size_t n_particles); [[nodiscard]] ResultExpr biorthogonal_transform_copy( const ResultExpr& expr, double threshold = default_biorth_threshold); @@ -27,12 +55,168 @@ void biorthogonal_transform(container::svector& exprs, /// performs symbolic biorthogonal transform of CC-like equation using ///(for rank-3 and higher /// Wang-Knizia biorthogonalization (https://arxiv.org/abs/1805.00565) is used +/// uses hardcoded coefficients for ranks 1-6, computed coefficients for higher +/// ranks [[nodiscard]] ExprPtr biorthogonal_transform( const ExprPtr& expr, const container::svector>& ext_index_groups = {}, double threshold = default_biorth_threshold); +/// \brief Computes the non-null space (NNS) projection coefficients +/// +/// \param n_particles The rank of external index pairs +/// \param threshold The threshold to compute the pseudoinverse matrix +/// (set to default_biorth_threshold) +/// +/// \return Vector of computed NNS projection coefficients +[[nodiscard]] std::vector compute_nns_p_coeffs( + std::size_t n_particles, double threshold = default_biorth_threshold); + +/// \brief Provides permuted indices using libperm unrank function +/// +/// \param indices The indices to permute +/// \param perm_rank The rank of the permutation +/// \param n_particles The rank of external index pairs +/// +/// \return The permuted indices +container::svector compute_permuted_indices( + const container::svector& indices, size_t perm_rank, + size_t n_particles); + +/// \brief Provides one row of the NNS projector matrix, +/// hardcoded from Mathematica to avoid numerical precision loss. +/// +/// The NNS projector weights are obtained from the normalized pseudoinverse +/// of M: first compute M_pinv (the pseudoinverse), then normalize it by the +/// factor ((n_particles)!/rank(M)). +//// Finally, NNS projector = normalized_M_pinv · M. +/// +/// \param n_particles The rank of external index pairs +/// +/// \return Vector of NNS projector weights representing the last row +/// +/// \throw std::runtime_error if n_particles is not in the range [1,5] +template + requires(std::floating_point || meta::is_complex_v) +std::optional> hardcoded_nns_projector(std::size_t n_particles) { + switch (n_particles) { + case 1: + return std::vector{T(1) / T(1)}; + + case 2: + return std::vector{T(1) / T(1), T(1) / T(1)}; + + case 3: + return std::vector{T(-1) / T(5), T(-1) / T(5), T(-1) / T(5), + T(-1) / T(5), T(-1) / T(5), T(1) / T(1)}; + + case 4: + return std::vector{ + T(1) / T(7), T(1) / T(7), T(1) / T(7), T(-1) / T(14), + T(1) / T(7), T(1) / T(7), T(1) / T(7), T(-1) / T(14), + T(-1) / T(14), T(-1) / T(14), T(1) / T(7), T(-2) / T(7), + T(-1) / T(14), T(1) / T(7), T(-1) / T(14), T(-2) / T(7), + T(1) / T(7), T(-1) / T(14), T(-1) / T(14), T(-2) / T(7), + T(-2) / T(7), T(-2) / T(7), T(-2) / T(7), T(1) / T(1)}; + + case 5: + return std::vector{ + T(-1) / T(14), T(-1) / T(14), T(-1) / T(14), T(-1) / T(14), + T(2) / T(21), T(-1) / T(14), T(-1) / T(14), T(-1) / T(14), + T(-1) / T(14), T(2) / T(21), T(-1) / T(14), T(-1) / T(14), + T(-1) / T(14), T(-1) / T(14), T(2) / T(21), T(2) / T(21), + T(2) / T(21), T(2) / T(21), T(-1) / T(21), T(0) / T(1), + T(-1) / T(14), T(-1) / T(14), T(-1) / T(14), T(-1) / T(14), + T(2) / T(21), T(-1) / T(14), T(-1) / T(14), T(-1) / T(14), + T(-1) / T(14), T(2) / T(21), T(-1) / T(14), T(-1) / T(14), + T(-1) / T(14), T(-1) / T(14), T(2) / T(21), T(2) / T(21), + T(2) / T(21), T(2) / T(21), T(-1) / T(21), T(0) / T(1), + T(2) / T(21), T(2) / T(21), T(-1) / T(21), T(2) / T(21), + T(0) / T(1), T(2) / T(21), T(2) / T(21), T(-1) / T(21), + T(2) / T(21), T(0) / T(1), T(-1) / T(21), T(-1) / T(21), + T(-1) / T(21), T(-1) / T(21), T(1) / T(7), T(0) / T(1), + T(0) / T(1), T(1) / T(7), T(1) / T(7), T(-1) / T(3), + T(2) / T(21), T(-1) / T(21), T(2) / T(21), T(2) / T(21), + T(0) / T(1), T(-1) / T(21), T(-1) / T(21), T(-1) / T(21), + T(-1) / T(21), T(1) / T(7), T(2) / T(21), T(-1) / T(21), + T(2) / T(21), T(2) / T(21), T(0) / T(1), T(0) / T(1), + T(1) / T(7), T(0) / T(1), T(1) / T(7), T(-1) / T(3), + T(-1) / T(21), T(-1) / T(21), T(-1) / T(21), T(-1) / T(21), + T(1) / T(7), T(-1) / T(21), T(2) / T(21), T(2) / T(21), + T(2) / T(21), T(0) / T(1), T(-1) / T(21), T(2) / T(21), + T(2) / T(21), T(2) / T(21), T(0) / T(1), T(1) / T(7), + T(0) / T(1), T(0) / T(1), T(1) / T(7), T(-1) / T(3), + T(0) / T(1), T(1) / T(7), T(1) / T(7), T(0) / T(1), + T(-1) / T(3), T(1) / T(7), T(0) / T(1), T(1) / T(7), + T(0) / T(1), T(-1) / T(3), T(1) / T(7), T(1) / T(7), + T(0) / T(1), T(0) / T(1), T(-1) / T(3), T(-1) / T(3), + T(-1) / T(3), T(-1) / T(3), T(-1) / T(3), T(1) / T(1)}; + + default: + return std::nullopt; + } +} + +/// \brief Provides NNS projection weights for a given rank +/// +/// \tparam T The numeric type (must be floating point or complex) +/// \param n_particles The rank of external index pairs +/// \param threshold The threshold to compute the pseudoinverse matrix +/// (set to default_biorth_threshold) +/// +/// \return (memoized) Vector of hrdcoded/computed NNS projection weights +template + requires(std::floating_point || meta::is_complex_v) +[[nodiscard]] const std::vector& nns_projection_weights( + std::size_t n_particles, double threshold = default_biorth_threshold) { + static const std::vector empty_vec{}; + + if (n_particles < 3) { + return empty_vec; + } + + using CacheKey = std::pair; + using CacheValue = std::optional>; + + static std::mutex cache_mutex; + static std::condition_variable cache_cv; + static container::map cache; + + CacheKey key{n_particles, threshold}; + + { + std::unique_lock lock(cache_mutex); + auto [it, inserted] = cache.try_emplace(key, std::nullopt); + if (!inserted) { + cache_cv.wait(lock, [&] { return it->second.has_value(); }); + return it->second.value(); + } + } + + std::vector nns_p_coeffs; + + constexpr std::size_t max_rank_hardcoded_nns_projector = 5; + if (n_particles <= max_rank_hardcoded_nns_projector) { + auto hardcoded_coeffs = hardcoded_nns_projector(n_particles); + if (hardcoded_coeffs) { + nns_p_coeffs = std::move(hardcoded_coeffs.value()); + } + } else { + auto coeffs = compute_nns_p_coeffs(n_particles, threshold); + nns_p_coeffs.reserve(coeffs.size()); + for (const auto& c : coeffs) { + nns_p_coeffs.push_back(static_cast(c)); + } + } + + { + std::lock_guard lock(cache_mutex); + cache[key] = std::move(nns_p_coeffs); + cache_cv.notify_all(); + return cache[key].value(); + } +} } // namespace sequant #endif diff --git a/SeQuant/domain/mbpt/spin.cpp b/SeQuant/domain/mbpt/spin.cpp index ae5b3f0acf..90c5ec8a82 100644 --- a/SeQuant/domain/mbpt/spin.cpp +++ b/SeQuant/domain/mbpt/spin.cpp @@ -1169,28 +1169,20 @@ ExprPtr closed_shell_CC_spintrace_v2(ExprPtr const& expr, st_expr; } simplify(st_expr); - // expanding S after spintracing and biorthogonalization, to avoid dealing - // with large number of terms + st_expr = S_maps(st_expr); // canonicalizer must be called before hash-filter to combine terms canonicalize(st_expr); - // apply hash filter method to get unique set of terms st_expr = WK_biorthogonalization_filter(st_expr, ext_idxs); - // add S tensor again + st_expr = ex(Tensor{L"S", bra(std::move(kixs)), ket(std::move(bixs))}) * st_expr; - rational combined_factor; - if (ext_idxs.size() <= 2) { - combined_factor = rational(1, factorial(ext_idxs.size())); - } else { - auto fact_n = factorial(ext_idxs.size()); - combined_factor = - rational(1, fact_n - 1); // this is (1/fact_n) * (fact_n/(fact_n-1)) - } - st_expr = ex(combined_factor) * st_expr; + const auto nf = ex( + rational{1, factorial(ext_idxs.size())}); // normalization factor for S + st_expr = nf * st_expr; } simplify(st_expr); diff --git a/tests/unit/test_eval_btas.cpp b/tests/unit/test_eval_btas.cpp index 1d0db339ad..cbdae80816 100644 --- a/tests/unit/test_eval_btas.cpp +++ b/tests/unit/test_eval_btas.cpp @@ -401,14 +401,13 @@ TEST_CASE("eval_with_btas", "[eval_btas]") { BTensorD perm_sum{r2.range()}; perm_sum.fill(0); - perm_sum += r2; perm_sum += BTensorD{permute(r2, {0, 1, 2, 3, 5, 4})}; perm_sum += BTensorD{permute(r2, {0, 1, 2, 4, 3, 5})}; perm_sum += BTensorD{permute(r2, {0, 1, 2, 4, 5, 3})}; perm_sum += BTensorD{permute(r2, {0, 1, 2, 5, 3, 4})}; perm_sum += BTensorD{permute(r2, {0, 1, 2, 5, 4, 3})}; - btas::scal(1.0 / 6.0, perm_sum); + btas::scal(1.0 / 5.0, perm_sum); man2 -= perm_sum; REQUIRE(norm(eval2) == Catch::Approx(norm(man2))); diff --git a/tests/unit/test_eval_ta.cpp b/tests/unit/test_eval_ta.cpp index 721034a6a3..d3e6deca9a 100644 --- a/tests/unit/test_eval_ta.cpp +++ b/tests/unit/test_eval_ta.cpp @@ -510,7 +510,7 @@ TEST_CASE("eval_with_tiledarray", "[eval]") { } SECTION("Biorthogonal Cleanup") { - // low-rank residuals: skip cleanup + // low-rank residuals: skip nns auto expr1 = parse_antisymm(L"R_{a1, a2}^{i1, i2}"); auto eval1 = eval_biorthogonal_nns_project(expr1, "a_1,a_2,i_1,i_2"); auto const& arr1 = yield(L"R{a1,a2;i1,i2}"); @@ -524,8 +524,8 @@ TEST_CASE("eval_with_tiledarray", "[eval]") { REQUIRE(norm(zero1) == Catch::Approx(0).margin( 100 * std::numeric_limits::epsilon())); - // high-rank residuals: cleanup applies: - // result = identity - (1/ket_rank!) * sum_of_ket_permutations + // for rank 3 residual, nns applies: + // result = NNS_P * sum_of_ket_permutations auto expr2 = parse_antisymm(L"R_{a1, a2, a3}^{i1, i2, i3}"); auto eval2 = eval_biorthogonal_nns_project(expr2, "a_1,a_2,a_3,i_1,i_2,i_3"); @@ -534,15 +534,56 @@ TEST_CASE("eval_with_tiledarray", "[eval]") { auto man2 = TArrayD{}; man2("0,1,2,3,4,5") = arr2("0,1,2,3,4,5") - - (1.0 / 6.0) * - (arr2("0,1,2,3,4,5") + arr2("0,1,2,3,5,4") + arr2("0,1,2,4,3,5") + - arr2("0,1,2,4,5,3") + arr2("0,1,2,5,3,4") + arr2("0,1,2,5,4,3")); + (1.0 / 5.0) * + (arr2("0,1,2,3,5,4") + arr2("0,1,2,4,3,5") + arr2("0,1,2,4,5,3") + + arr2("0,1,2,5,3,4") + arr2("0,1,2,5,4,3")); REQUIRE(norm(man2) == Catch::Approx(norm(eval2))); TArrayD zero2; zero2("0,1,2,3,4,5") = man2("0,1,2,3,4,5") - eval2("0,1,2,3,4,5"); REQUIRE(norm(zero1) == Catch::Approx(0).margin( 100 * std::numeric_limits::epsilon())); + + // for rank 4 residual, nns applies: + // result = NNS_P * sum_of_ket_permutations + auto expr3 = parse_antisymm(L"R_{a1, a2, a3, a4}^{i1, i2, i3, i4}"); + auto eval3 = eval_biorthogonal_nns_project( + expr3, "a_1,a_2,a_3,a_4,i_1,i_2,i_3,i_4"); + auto const& arr3 = yield(L"R{a1,a2,a3,a4;i1,i2,i3,i4}"); + + auto man3 = TArrayD{}; + man3("0,1,2,3,4,5,6,7") = 1.0 * arr3("0,1,2,3,4,5,6,7") + + -4.0 / 14.0 * arr3("0,1,2,3,4,5,7,6") + + -4.0 / 14.0 * arr3("0,1,2,3,4,6,5,7") + + -1.0 / 14.0 * arr3("0,1,2,3,4,6,7,5") + + -1.0 / 14.0 * arr3("0,1,2,3,4,7,5,6") + + -4.0 / 14.0 * arr3("0,1,2,3,4,7,6,5") + + -4.0 / 14.0 * arr3("0,1,2,3,5,4,6,7") + + 2.0 / 14.0 * arr3("0,1,2,3,5,4,7,6") + + -1.0 / 14.0 * arr3("0,1,2,3,5,6,4,7") + + 2.0 / 14.0 * arr3("0,1,2,3,5,6,7,4") + + 2.0 / 14.0 * arr3("0,1,2,3,5,7,4,6") + + -1.0 / 14.0 * arr3("0,1,2,3,5,7,6,4") + + -1.0 / 14.0 * arr3("0,1,2,3,6,4,5,7") + + 2.0 / 14.0 * arr3("0,1,2,3,6,4,7,5") + + -4.0 / 14.0 * arr3("0,1,2,3,6,5,4,7") + + -1.0 / 14.0 * arr3("0,1,2,3,6,5,7,4") + + 2.0 / 14.0 * arr3("0,1,2,3,6,7,4,5") + + 2.0 / 14.0 * arr3("0,1,2,3,6,7,5,4") + + 2.0 / 14.0 * arr3("0,1,2,3,7,4,5,6") + + -1.0 / 14.0 * arr3("0,1,2,3,7,4,6,5") + + -1.0 / 14.0 * arr3("0,1,2,3,7,5,4,6") + + -4.0 / 14.0 * arr3("0,1,2,3,7,5,6,4") + + 2.0 / 14.0 * arr3("0,1,2,3,7,6,4,5") + + 2.0 / 14.0 * arr3("0,1,2,3,7,6,5,4"); + + REQUIRE(norm(man3) == Catch::Approx(norm(eval3))); + TArrayD zero3; + zero3("0,1,2,3,4,5,6,7") = + man3("0,1,2,3,4,5,6,7") - eval3("0,1,2,3,4,5,6,7"); + REQUIRE(norm(zero3) == + Catch::Approx(0).margin(1000 * + std::numeric_limits::epsilon())); } SECTION("Others") { diff --git a/tests/unit/test_spin.cpp b/tests/unit/test_spin.cpp index e3aeb180c3..da85c53262 100644 --- a/tests/unit/test_spin.cpp +++ b/tests/unit/test_spin.cpp @@ -1090,8 +1090,8 @@ SECTION("Closed-shell spintrace CCSDT terms") { "g{a_1,a_2;a_4,a_5}:N-C-S * t{a_3,a_4,a_5;i_1,i_2,i_3}:N-C-S + 2 " "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_1,i_3,i_2}:N-C-S")); - // the new efficient method, spintracing with partial expansion, then - // expanding by S_map ( this method is used in + // the new efficient method, does spintracing with partial expansion, then + // expanding by S_map (this method is used in // closed_shell_CC_spintrace_v2) auto result_2 = closed_shell_spintrace( input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}}); @@ -1124,19 +1124,13 @@ SECTION("Closed-shell spintrace CCSDT terms") { "g{a_1,a_3;a_4,a_5}:N-C-S * t{a_2,a_4,a_5;i_1,i_3,i_2}:N-C-S")); } - SECTION("ppl term in optimal") { // results in 1 term + SECTION("most expensive terms in CCSDT") { // results in 1 term const auto input = ex(ExprPtrList{ parse_expr(L"1/24 A{i_1,i_2,i_3;a_1,a_2,a_3} * " L"g{a_1,a_2;a_4,a_5} * t{a_3,a_4,a_5;i_1,i_2,i_3}", Symmetry::Antisymm)}); auto result = closed_shell_CC_spintrace_v2(input); - // multiply the result by 6/5 to revert the rescaling factor - result *= ex(rational{5, 6}); - - // There is a problem with casting a single term to Sum - // REQUIRE(result->size()== 1); // it needs to be checked - REQUIRE_THAT( result, EquivalentTo(