From ecfc347c3ddef640a115d92a375062f2d3c7b89b Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Thu, 12 Mar 2026 11:35:17 -0500 Subject: [PATCH 1/2] non-symmetric case works --- src/integrals/ao_integrals/ao_integrals.hpp | 5 + .../ao_integrals/uq_atom_blocked_driver.cpp | 213 +++++++++++++++ .../uq_atom_symm_blocked_driver.cpp | 258 ++++++++++++++++++ .../utils/uncertainty_reductions.hpp | 80 ++++++ .../test_uq_atom_blocked_driver.cpp | 121 ++++++++ .../utils/uncertainty_reductions.cpp | 90 ++++++ 6 files changed, 767 insertions(+) create mode 100644 src/integrals/ao_integrals/uq_atom_blocked_driver.cpp create mode 100644 src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp create mode 100644 src/integrals/utils/uncertainty_reductions.hpp create mode 100644 tests/cxx/unit/integrals/ao_integrals/test_uq_atom_blocked_driver.cpp create mode 100644 tests/cxx/unit/integrals/utils/uncertainty_reductions.cpp diff --git a/src/integrals/ao_integrals/ao_integrals.hpp b/src/integrals/ao_integrals/ao_integrals.hpp index 6fe9fb10..6ec507e8 100644 --- a/src/integrals/ao_integrals/ao_integrals.hpp +++ b/src/integrals/ao_integrals/ao_integrals.hpp @@ -27,6 +27,7 @@ DECLARE_MODULE(KDensityFitted); DECLARE_MODULE(DFIntegral); DECLARE_MODULE(CoulombMetric); DECLARE_MODULE(UQDriver); +DECLARE_MODULE(UQAtomBlockedDriver); inline void set_defaults(pluginplay::ModuleManager& mm) { mm.change_submod("AO integral driver", "Coulomb matrix", @@ -41,6 +42,9 @@ inline void set_defaults(pluginplay::ModuleManager& mm) { "Coulomb Metric"); mm.change_submod("UQ Driver", "ERIs", "ERI4"); mm.change_submod("UQ Driver", "ERI Error", "Primitive Error Model"); + mm.change_submod("UQ Atom Blocked Driver", "ERIs", "ERI4"); + mm.change_submod("UQ Atom Blocked Driver", "ERI Error", + "Primitive Error Model"); } inline void load_modules(pluginplay::ModuleManager& mm) { @@ -52,6 +56,7 @@ inline void load_modules(pluginplay::ModuleManager& mm) { mm.add_module("Density Fitting Integral"); mm.add_module("Coulomb Metric"); mm.add_module("UQ Driver"); + mm.add_module("UQ Atom Blocked Driver"); } } // namespace integrals::ao_integrals diff --git a/src/integrals/ao_integrals/uq_atom_blocked_driver.cpp b/src/integrals/ao_integrals/uq_atom_blocked_driver.cpp new file mode 100644 index 00000000..496f12a2 --- /dev/null +++ b/src/integrals/ao_integrals/uq_atom_blocked_driver.cpp @@ -0,0 +1,213 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "../utils/uncertainty_reductions.hpp" +#include "ao_integrals.hpp" +#include +#ifdef ENABLE_SIGMA +#include +#endif + +using namespace tensorwrapper; + +namespace integrals::ao_integrals { +namespace { + +template +auto average_error(T&& strides, T&& nbf, T&& ao_i, Tensor&& error, + utils::mean_type mean) { + std::vector buffer; + + for(std::size_t i = 0; i < nbf[0]; ++i) { + auto ioffset = (ao_i[0] + i) * strides[0]; + + for(std::size_t j = 0; j < nbf[1]; ++j) { + auto joffset = ioffset + (ao_i[1] + j) * strides[1]; + + for(std::size_t k = 0; k < nbf[2]; ++k) { + auto koffset = joffset + (ao_i[2] + k) * strides[2]; + + for(std::size_t l = 0; l < nbf[3]; ++l) { + auto loffset = koffset + (ao_i[3] + l) * strides[3]; + buffer.push_back(error[loffset]); + } + } + } + } + + return utils::compute_mean(mean, buffer); +} + +template +void update_block(T&& strides, T&& nbf, T&& ao_i, std::vector& out, + Tensor&& value, FloatType error) { + for(std::size_t i = 0; i < nbf[0]; ++i) { + auto ioffset = (ao_i[0] + i) * strides[0]; + + for(std::size_t j = 0; j < nbf[1]; ++j) { + auto joffset = ioffset + (ao_i[1] + j) * strides[1]; + + for(std::size_t k = 0; k < nbf[2]; ++k) { + auto koffset = joffset + (ao_i[2] + k) * strides[2]; + + for(std::size_t l = 0; l < nbf[3]; ++l) { + auto loffset = koffset + (ao_i[3] + l) * strides[3]; + out[loffset] = value[loffset] + error; + } + } + } + } +} + +struct Kernel { + using shape_type = buffer::Contiguous::shape_type; + Kernel(shape_type shape, std::array aos, + utils::mean_type mean) : + m_shape(std::move(shape)), m_aos(aos), m_mean(mean) {} + + template + Tensor operator()(const std::span t, + const std::span error) { + throw std::runtime_error( + "UQ Integrals Driver kernel only supports same float types"); + } + + template + auto operator()(const std::span t, + const std::span error) { + Tensor rv; + + using float_type = std::decay_t; + if constexpr(types::is_uncertain_v) { + throw std::runtime_error("Did not expect an uncertain type"); + } else { +#ifdef ENABLE_SIGMA + using tensorwrapper::buffer::make_contiguous; + + std::array n_centers{m_aos[0].size(), m_aos[1].size(), + m_aos[2].size(), m_aos[3].size()}; + + std::array centers{0, 0, 0, 0}; + std::array ao_i{0, 0, 0, 0}; + std::array nbf{0, 0, 0, 0}; + + using uq_type = sigma::Uncertain; + std::vector rv_data(m_shape.size()); + std::array strides{0, 0, 0, 1}; + strides[2] = strides[3] * m_aos[3].n_aos(); + strides[1] = strides[2] * m_aos[2].n_aos(); + strides[0] = strides[1] * m_aos[1].n_aos(); + + for(centers[0] = 0; centers[0] < n_centers[0]; ++centers[0]) { + nbf[0] = m_aos[0][centers[0]].n_aos(); + + ao_i[1] = 0; + for(centers[1] = 0; centers[1] < n_centers[1]; ++centers[1]) { + nbf[1] = m_aos[1][centers[1]].n_aos(); + + ao_i[2] = 0; + for(centers[2] = 0; centers[2] < n_centers[2]; + ++centers[2]) { + nbf[2] = m_aos[2][centers[2]].n_aos(); + + ao_i[3] = 0; + for(centers[3] = 0; centers[3] < n_centers[3]; + ++centers[3]) { + nbf[3] = m_aos[3][centers[3]].n_aos(); + + auto block_error = average_error( + strides, nbf, ao_i, error, m_mean); + uq_type max_uq{0.0, block_error}; + + update_block(strides, nbf, ao_i, rv_data, t, + max_uq); + + ao_i[3] += nbf[3]; + } + ao_i[2] += nbf[2]; + } + ao_i[1] += nbf[1]; + } + ao_i[0] += nbf[0]; + } + tensorwrapper::buffer::Contiguous t_w_contig(std::move(rv_data), + m_shape); + rv = tensorwrapper::Tensor(m_shape, std::move(t_w_contig)); +#else + throw std::runtime_error("Sigma support not enabled!"); +#endif + } + + return rv; + } + shape_type m_shape; + std::array m_aos; + utils::mean_type m_mean; +}; + +const auto desc = R"( +UQ Integrals Driver +------------------- + +)"; + +} // namespace + +using eri_pt = simde::ERI4; +using error_pt = integrals::property_types::Uncertainty; + +MODULE_CTOR(UQAtomBlockedDriver) { + satisfies_property_type(); + description(desc); + add_submodule("ERIs"); + add_submodule("ERI Error"); + add_input("Mean Type").set_default("none"); +} + +MODULE_RUN(UQAtomBlockedDriver) { + const auto& [braket] = eri_pt::unwrap_inputs(inputs); + auto mean_str = inputs.at("Mean Type").value(); + auto mean = utils::mean_from_string(mean_str); + + auto& eri_mod = submods.at("ERIs").value(); + auto tol = eri_mod.inputs().at("Threshold").value(); + + const auto& t = eri_mod.run_as(braket); + const auto& error = submods.at("ERI Error").run_as(braket, tol); + + using tensorwrapper::buffer::make_contiguous; + const auto& t_buffer = make_contiguous(t.buffer()); + const auto& e_buffer = make_contiguous(error.buffer()); + + const auto& bra = braket.bra(); + const auto& ket = braket.ket(); + const auto& mu = bra.first.ao_basis_set(); + const auto& nu = bra.second.ao_basis_set(); + const auto& lam = ket.first.ao_basis_set(); + const auto& sig = ket.second.ao_basis_set(); + + std::array aos{mu, nu, lam, sig}; + + using buffer::visit_contiguous_buffer; + shape::Smooth shape = t.buffer().layout().shape().as_smooth().make_smooth(); + + Kernel k(shape, aos, mean); + auto t_w_error = visit_contiguous_buffer(k, t_buffer, e_buffer); + + auto rv = results(); + return eri_pt::wrap_results(rv, t_w_error); +} +} // namespace integrals::ao_integrals diff --git a/src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp b/src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp new file mode 100644 index 00000000..cd741973 --- /dev/null +++ b/src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp @@ -0,0 +1,258 @@ +// /* +// * Copyright 2025 NWChemEx-Project +// * +// * 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 +// * +// * 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. +// */ + +// #include "../utils/uncertainty_reductions.hpp" +// #include "ao_integrals.hpp" +// #include +// #ifdef ENABLE_SIGMA +// #include +// #endif + +// using namespace tensorwrapper; + +// namespace integrals::ao_integrals { +// namespace { + +// template +// auto average_error(T&& strides, T&& nbf, T&& ao_i, Tensor&& error, +// utils::mean_type mean) { +// std::vector buffer; + +// for(std::size_t i = 0; i < nbf[0]; ++i) { +// auto ioffset = (ao_i[0] + i) * strides[0]; + +// for(std::size_t j = 0; j < nbf[1]; ++j) { +// auto joffset = ioffset + (ao_i[1] + j) * strides[1]; + +// for(std::size_t k = 0; k < nbf[2]; ++k) { +// auto koffset = joffset + (ao_i[2] + k) * strides[2]; + +// for(std::size_t l = 0; l < nbf[3]; ++l) { +// auto loffset = koffset + (ao_i[3] + l) * strides[3]; +// buffer.push_back(error[loffset]); +// } +// } +// } +// } + +// return utils::compute_mean(mean, buffer); +// } + +// template +// void update_block(T&& strides, T&& nbf, T&& ao_i, std::vector& +// out, +// Tensor&& value, FloatType error) { +// for(std::size_t i = 0; i < nbf[0]; ++i) { +// auto ioffset = (ao_i[0] + i) * strides[0]; + +// for(std::size_t j = 0; j < nbf[1]; ++j) { +// auto joffset = ioffset + (ao_i[1] + j) * strides[1]; + +// for(std::size_t k = 0; k < nbf[2]; ++k) { +// auto koffset = joffset + (ao_i[2] + k) * strides[2]; + +// for(std::size_t l = 0; l < nbf[3]; ++l) { +// auto loffset = koffset + (ao_i[3] + l) * strides[3]; +// out[loffset] = value[loffset] + error; +// } +// } +// } +// } +// } + +// template +// void copy_block(T&& strides, T&& nbf, T&& ao_i, T&& new_nbf, T&& new_ao_i, +// std::vector& out) { +// for(std::size_t i = 0; i < nbf[0]; ++i) { +// auto ioffset = (ao_i[0] + i) * strides[0]; + +// for(std::size_t j = 0; j < nbf[1]; ++j) { +// auto joffset = ioffset + (ao_i[1] + j) * strides[1]; + +// for(std::size_t k = 0; k < nbf[2]; ++k) { +// auto koffset = joffset + (ao_i[2] + k) * strides[2]; + +// for(std::size_t l = 0; l < nbf[3]; ++l) { +// auto loffset = koffset + (ao_i[3] + l) * strides[3]; +// out[loffset] = value[loffset] + error; +// } +// } +// } +// } +// } + +// struct Kernel { +// using shape_type = buffer::Contiguous::shape_type; +// Kernel(shape_type shape, std::array aos, +// utils::mean_type mean) : +// m_shape(std::move(shape)), m_aos(aos), m_mean(mean) {} + +// template +// Tensor operator()(const std::span t, +// const std::span error) { +// throw std::runtime_error( +// "UQ Integrals Driver kernel only supports same float types"); +// } + +// template +// auto operator()(const std::span t, +// const std::span error) { +// Tensor rv; + +// using float_type = std::decay_t; +// if constexpr(types::is_uncertain_v) { +// throw std::runtime_error("Did not expect an uncertain type"); +// } else { +// #ifdef ENABLE_SIGMA +// using tensorwrapper::buffer::make_contiguous; + +// std::array n_centers{m_aos[0].size(), m_aos[1].size(), +// m_aos[2].size(), m_aos[3].size()}; + +// std::array centers{0, 0, 0, 0}; +// std::array ao_i{0, 0, 0, 0}; +// std::array nbf{0, 0, 0, 0}; + +// using uq_type = sigma::Uncertain; +// std::vector rv_data(m_shape.size()); +// std::array strides{0, 0, 0, 1}; +// strides[2] = strides[3] * m_aos[3].n_aos(); +// strides[1] = strides[2] * m_aos[2].n_aos(); +// strides[0] = strides[1] * m_aos[1].n_aos(); + +// // If true, we can skip centers[1] > centers[0] +// bool mu_is_nu = (m_aos[0] == m_aos[1]); +// // If true, we can skip centers[3] > centers[2] +// bool lam_is_sig = (m_aos[2] == m_aos[3]); +// // If true, can skip centers[2] > centers[0] and centers[3] > +// // centers[1] +// bool all_same = (m_aos[0] == m_aos[2]) && mu_is_nu && lam_is_sig; + +// for(centers[0] = 0; centers[0] < n_centers[0]; ++centers[0]) { +// nbf[0] = m_aos[0][centers[0]].n_aos(); + +// ao_i[1] = 0; +// for(centers[1] = 0; centers[1] < n_centers[1]; ++centers[1]) +// { +// if(centers[1] > centers[0] && mu_is_nu) break; +// nbf[1] = m_aos[1][centers[1]].n_aos(); + +// ao_i[2] = 0; +// for(centers[2] = 0; centers[2] < n_centers[2]; +// ++centers[2]) { +// nbf[2] = m_aos[2][centers[2]].n_aos(); + +// ao_i[3] = 0; +// for(centers[3] = 0; centers[3] < n_centers[3]; +// ++centers[3]) { +// if(centers[3] > centers[2] && lam_is_sig) break; + +// nbf[3] = m_aos[3][centers[3]].n_aos(); + +// auto block_error = average_error( +// strides, nbf, ao_i, error, m_mean); +// uq_type max_uq{0.0, block_error}; + +// update_block(strides, nbf, ao_i, rv_data, t, +// max_uq); + +// if(mu_is_nu) { +// std::array new_nbf{ +// nbf[1], nbf[0], nbf[2], nbf[3]}; +// std::array new_ao_i{ +// ao_i[1], ao_i[0], ao_i[2], ao_i[3]}; +// copy_block(strides, nbf, ao_i, new_nbf, +// new_ao_i, rv_data); +// } + +// ao_i[3] += nbf[3]; +// } +// ao_i[2] += nbf[2]; +// } +// ao_i[1] += nbf[1]; +// } +// ao_i[0] += nbf[0]; +// } +// tensorwrapper::buffer::Contiguous t_w_contig(std::move(rv_data), +// m_shape); +// rv = tensorwrapper::Tensor(m_shape, std::move(t_w_contig)); +// #else +// throw std::runtime_error("Sigma support not enabled!"); +// #endif +// } + +// return rv; +// } +// shape_type m_shape; +// std::array m_aos; +// utils::mean_type m_mean; +// }; + +// const auto desc = R"( +// UQ Integrals Driver +// ------------------- + +// )"; + +// } // namespace + +// using eri_pt = simde::ERI4; +// using error_pt = integrals::property_types::Uncertainty; + +// MODULE_CTOR(UQAtomBlockedDriver) { +// satisfies_property_type(); +// description(desc); +// add_submodule("ERIs"); +// add_submodule("ERI Error"); +// add_input("Mean Type").set_default("none"); +// } + +// MODULE_RUN(UQAtomBlockedDriver) { +// const auto& [braket] = eri_pt::unwrap_inputs(inputs); +// auto mean_str = inputs.at("Mean Type").value(); +// auto mean = utils::mean_from_string(mean_str); + +// auto& eri_mod = submods.at("ERIs").value(); +// auto tol = eri_mod.inputs().at("Threshold").value(); + +// const auto& t = eri_mod.run_as(braket); +// const auto& error = submods.at("ERI Error").run_as(braket, +// tol); + +// using tensorwrapper::buffer::make_contiguous; +// const auto& t_buffer = make_contiguous(t.buffer()); +// const auto& e_buffer = make_contiguous(error.buffer()); + +// const auto& bra = braket.bra(); +// const auto& ket = braket.ket(); +// const auto& mu = bra.first.ao_basis_set(); +// const auto& nu = bra.second.ao_basis_set(); +// const auto& lam = ket.first.ao_basis_set(); +// const auto& sig = ket.second.ao_basis_set(); + +// std::array aos{mu, nu, lam, sig}; + +// using buffer::visit_contiguous_buffer; +// shape::Smooth shape = +// t.buffer().layout().shape().as_smooth().make_smooth(); + +// Kernel k(shape, aos, mean); +// auto t_w_error = visit_contiguous_buffer(k, t_buffer, e_buffer); + +// auto rv = results(); +// return eri_pt::wrap_results(rv, t_w_error); +// } +// } // namespace integrals::ao_integrals diff --git a/src/integrals/utils/uncertainty_reductions.hpp b/src/integrals/utils/uncertainty_reductions.hpp new file mode 100644 index 00000000..2ddfaf97 --- /dev/null +++ b/src/integrals/utils/uncertainty_reductions.hpp @@ -0,0 +1,80 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include +#include +#include + +namespace integrals::utils { + +enum class mean_type { none, max, geometric, harmonic }; + +template +auto geometric_mean(const ContainerType& values) { + using float_type = std::decay_t; + float_type log_sum = 0.0; + std::size_t non_zero_count = 0; + for(const auto& val : values) { + if(val == 0) continue; + log_sum += std::log(val); + ++non_zero_count; + } + return std::exp(log_sum / non_zero_count); +} + +template +auto harmonic_mean(const ContainerType& values) { + using float_type = std::decay_t; + float_type reciprocal_sum = 0.0; + std::size_t non_zero_count = 0; + for(const auto& val : values) { + if(val == 0) continue; + reciprocal_sum += 1 / val; + ++non_zero_count; + } + return non_zero_count / reciprocal_sum; +} + +template +auto compute_mean(mean_type mean, const ContainerType& span) { + if(mean == mean_type::max) { + return *std::max_element(span.begin(), span.end()); + } else if(mean == mean_type::harmonic) { + return harmonic_mean(span); + } else if(mean == mean_type::geometric) { + return geometric_mean(span); + } else { + throw std::runtime_error("Mean type not supported"); + } +} + +inline auto mean_from_string(const std::string& mean_str) { + if(mean_str == std::string("max")) { + return mean_type::max; + } else if(mean_str == std::string("geometric")) { + return mean_type::geometric; + } else if(mean_str == std::string("harmonic")) { + return mean_type::harmonic; + } else if(mean_str == std::string("none")) { + return mean_type::none; + } else { + throw std::runtime_error("Mean type not supported: " + mean_str); + } +} + +} // namespace integrals::utils diff --git a/tests/cxx/unit/integrals/ao_integrals/test_uq_atom_blocked_driver.cpp b/tests/cxx/unit/integrals/ao_integrals/test_uq_atom_blocked_driver.cpp new file mode 100644 index 00000000..48d26a8f --- /dev/null +++ b/tests/cxx/unit/integrals/ao_integrals/test_uq_atom_blocked_driver.cpp @@ -0,0 +1,121 @@ +/* + * Copyright 2022 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "../testing/testing.hpp" + +using namespace integrals::testing; + +using namespace tensorwrapper; + +namespace { + +// N.b. The "means" of the correct values are validated by comparing to Libint's +// results. With the exception of the "No Mean" values, the +// "standard deviations" are not validated, but seem reasonable (the "No Mean" +// values come from the unit tests of the underlying error module). + +template +auto corr_answer(const simde::type::tensor& T) { + if constexpr(std::is_same_v) { + return T; + } else { + simde::type::tensor T_corr(T); + auto& corr_buffer = buffer::make_contiguous(T_corr.buffer()); + corr_buffer.set_elem({0, 0, 0, 0}, FloatType{0.774606, 0}); + corr_buffer.set_elem({0, 0, 0, 1}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({0, 0, 1, 0}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({0, 0, 1, 1}, FloatType{0.446701, 0}); + corr_buffer.set_elem({0, 1, 0, 0}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({0, 1, 0, 1}, + FloatType{0.120666, 0.0000170000000000}); + corr_buffer.set_elem({0, 1, 1, 0}, + FloatType{0.120666, 0.0000170000000000}); + corr_buffer.set_elem({0, 1, 1, 1}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 0, 0, 0}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 0, 0, 1}, + FloatType{0.120666, 0.0000170000000000}); + corr_buffer.set_elem({1, 0, 1, 0}, + FloatType{0.120666, 0.0000170000000000}); + corr_buffer.set_elem({1, 0, 1, 1}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 1, 0, 0}, FloatType{0.446701, 0}); + corr_buffer.set_elem({1, 1, 0, 1}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 1, 1, 0}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 1, 1, 1}, FloatType{0.774606, 0}); + return T_corr; + } +} + +} // namespace + +TEST_CASE("UQ Atom Blocked Driver") { + using float_type = tensorwrapper::types::udouble; + using test_pt = simde::ERI4; + + if constexpr(tensorwrapper::types::is_uncertain_v) { + auto rt = std::make_unique(); + pluginplay::ModuleManager mm(std::move(rt), nullptr); + integrals::load_modules(mm); + integrals::set_defaults(mm); + REQUIRE(mm.count("UQ Atom Blocked Driver")); + mm.change_input("ERI4", "Threshold", 1.0e-6); + auto mod = mm.at("UQ Atom Blocked Driver"); + + // Get basis set + auto aobs = h2_sto3g_basis_set(); + + // Make AOS object + simde::type::aos aos(aobs); + simde::type::aos_squared aos_squared(aos, aos); + + // Make Operator + simde::type::v_ee_type op{}; + + // Make BraKet Input + chemist::braket::BraKet braket(aos_squared, op, aos_squared); + + using tensorwrapper::operations::approximately_equal; + SECTION("No Mean") { REQUIRE_THROWS(mod.run_as(braket)); } + SECTION("Max Error") { + mod.change_input("Mean Type", "max"); + auto T = mod.run_as(braket); + + auto T_corr = corr_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + SECTION("Geometric Mean") { + mod.change_input("Mean Type", "geometric"); + auto T = mod.run_as(braket); + + auto T_corr = corr_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + SECTION("Harmonic Mean") { + mod.change_input("Mean Type", "harmonic"); + auto T = mm.at("UQ Driver").run_as(braket); + + auto T_corr = corr_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + } +} diff --git a/tests/cxx/unit/integrals/utils/uncertainty_reductions.cpp b/tests/cxx/unit/integrals/utils/uncertainty_reductions.cpp new file mode 100644 index 00000000..fa49b577 --- /dev/null +++ b/tests/cxx/unit/integrals/utils/uncertainty_reductions.cpp @@ -0,0 +1,90 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "../testing/testing.hpp" +#include +using namespace integrals::utils; + +TEST_CASE("geometric_mean") { + SECTION("Basic Test") { + std::vector values{1.0, 10.0, 100.0}; + auto result = geometric_mean(values); + REQUIRE(result == Catch::Approx(10.0)); + } + + SECTION("Skips zero") { + std::vector values{1.0, 10.0, 100.0, 0.0}; + auto result = geometric_mean(values); + REQUIRE(result == Catch::Approx(10.0)); + } +} + +TEST_CASE("harmonic_mean") { + SECTION("Basic Test") { + std::vector values{1.0, 2.0, 3.0}; + auto result = harmonic_mean(values); + REQUIRE(result == Catch::Approx(1.63636363636)); + } + + SECTION("Skips zero") { + std::vector values{1.0, 2.0, 3.0, 0.0}; + auto result = harmonic_mean(values); + REQUIRE(result == Catch::Approx(1.63636363636)); + } +} + +TEST_CASE("compute_mean") { + std::vector values{1.0, 2.0, 3.0}; + SECTION("Max") { + auto result = compute_mean(mean_type::max, values); + REQUIRE(result == Catch::Approx(3.0)); + } + SECTION("Geometric") { + auto result = compute_mean(mean_type::geometric, values); + REQUIRE(result == Catch::Approx(1.81712059283)); + } + SECTION("Harmonic") { + auto result = compute_mean(mean_type::harmonic, values); + REQUIRE(result == Catch::Approx(1.63636363636)); + } + SECTION("None") { + using except_t = std::runtime_error; + REQUIRE_THROWS_AS(compute_mean(mean_type::none, values), except_t); + } +} + +TEST_CASE("mean_from_string") { + SECTION("Max") { + auto result = mean_from_string("max"); + REQUIRE(result == mean_type::max); + } + SECTION("Geometric") { + auto result = mean_from_string("geometric"); + REQUIRE(result == mean_type::geometric); + } + SECTION("Harmonic") { + auto result = mean_from_string("harmonic"); + REQUIRE(result == mean_type::harmonic); + } + SECTION("None") { + auto result = mean_from_string("none"); + REQUIRE(result == mean_type::none); + } + SECTION("Invalid") { + using except_t = std::runtime_error; + REQUIRE_THROWS_AS(mean_from_string("invalid"), except_t); + } +} From 31faf72b466be84516ac1ed20d78cefd193b4aee Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Sat, 14 Mar 2026 13:38:55 -0500 Subject: [PATCH 2/2] adds symmetric atom blocking --- .gitignore | 1 + src/integrals/ao_integrals/ao_integrals.hpp | 5 + .../uq_atom_symm_blocked_driver.cpp | 543 +++++++++--------- src/integrals/utils/get_permutations.hpp | 93 +++ .../utils/uncertainty_reductions.hpp | 16 +- .../test_uq_atom_symm_blocked_driver.cpp | 147 +++++ tests/cxx/unit/integrals/testing/ao_bases.hpp | 34 +- .../unit/integrals/utils/get_permutations.cpp | 364 ++++++++++++ .../utils/uncertainty_reductions.cpp | 45 ++ 9 files changed, 971 insertions(+), 277 deletions(-) create mode 100644 src/integrals/utils/get_permutations.hpp create mode 100644 tests/cxx/unit/integrals/ao_integrals/test_uq_atom_symm_blocked_driver.cpp create mode 100644 tests/cxx/unit/integrals/utils/get_permutations.cpp diff --git a/.gitignore b/.gitignore index 21a2bbcc..b57d034a 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,7 @@ # These are directories used by IDEs for storing settings .idea/ .vscode/ +.cache/ # These are common Python virtual enviornment directory names venv/ diff --git a/src/integrals/ao_integrals/ao_integrals.hpp b/src/integrals/ao_integrals/ao_integrals.hpp index 6ec507e8..683d9fd3 100644 --- a/src/integrals/ao_integrals/ao_integrals.hpp +++ b/src/integrals/ao_integrals/ao_integrals.hpp @@ -28,6 +28,7 @@ DECLARE_MODULE(DFIntegral); DECLARE_MODULE(CoulombMetric); DECLARE_MODULE(UQDriver); DECLARE_MODULE(UQAtomBlockedDriver); +DECLARE_MODULE(UQAtomSymmBlockedDriver); inline void set_defaults(pluginplay::ModuleManager& mm) { mm.change_submod("AO integral driver", "Coulomb matrix", @@ -45,6 +46,9 @@ inline void set_defaults(pluginplay::ModuleManager& mm) { mm.change_submod("UQ Atom Blocked Driver", "ERIs", "ERI4"); mm.change_submod("UQ Atom Blocked Driver", "ERI Error", "Primitive Error Model"); + mm.change_submod("UQ Atom Symm Blocked Driver", "ERIs", "ERI4"); + mm.change_submod("UQ Atom Symm Blocked Driver", "ERI Error", + "Primitive Error Model"); } inline void load_modules(pluginplay::ModuleManager& mm) { @@ -57,6 +61,7 @@ inline void load_modules(pluginplay::ModuleManager& mm) { mm.add_module("Coulomb Metric"); mm.add_module("UQ Driver"); mm.add_module("UQ Atom Blocked Driver"); + mm.add_module("UQ Atom Symm Blocked Driver"); } } // namespace integrals::ao_integrals diff --git a/src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp b/src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp index cd741973..2384b1d3 100644 --- a/src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp +++ b/src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp @@ -1,258 +1,285 @@ -// /* -// * Copyright 2025 NWChemEx-Project -// * -// * 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 -// * -// * 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. -// */ - -// #include "../utils/uncertainty_reductions.hpp" -// #include "ao_integrals.hpp" -// #include -// #ifdef ENABLE_SIGMA -// #include -// #endif - -// using namespace tensorwrapper; - -// namespace integrals::ao_integrals { -// namespace { - -// template -// auto average_error(T&& strides, T&& nbf, T&& ao_i, Tensor&& error, -// utils::mean_type mean) { -// std::vector buffer; - -// for(std::size_t i = 0; i < nbf[0]; ++i) { -// auto ioffset = (ao_i[0] + i) * strides[0]; - -// for(std::size_t j = 0; j < nbf[1]; ++j) { -// auto joffset = ioffset + (ao_i[1] + j) * strides[1]; - -// for(std::size_t k = 0; k < nbf[2]; ++k) { -// auto koffset = joffset + (ao_i[2] + k) * strides[2]; - -// for(std::size_t l = 0; l < nbf[3]; ++l) { -// auto loffset = koffset + (ao_i[3] + l) * strides[3]; -// buffer.push_back(error[loffset]); -// } -// } -// } -// } - -// return utils::compute_mean(mean, buffer); -// } - -// template -// void update_block(T&& strides, T&& nbf, T&& ao_i, std::vector& -// out, -// Tensor&& value, FloatType error) { -// for(std::size_t i = 0; i < nbf[0]; ++i) { -// auto ioffset = (ao_i[0] + i) * strides[0]; - -// for(std::size_t j = 0; j < nbf[1]; ++j) { -// auto joffset = ioffset + (ao_i[1] + j) * strides[1]; - -// for(std::size_t k = 0; k < nbf[2]; ++k) { -// auto koffset = joffset + (ao_i[2] + k) * strides[2]; - -// for(std::size_t l = 0; l < nbf[3]; ++l) { -// auto loffset = koffset + (ao_i[3] + l) * strides[3]; -// out[loffset] = value[loffset] + error; -// } -// } -// } -// } -// } - -// template -// void copy_block(T&& strides, T&& nbf, T&& ao_i, T&& new_nbf, T&& new_ao_i, -// std::vector& out) { -// for(std::size_t i = 0; i < nbf[0]; ++i) { -// auto ioffset = (ao_i[0] + i) * strides[0]; - -// for(std::size_t j = 0; j < nbf[1]; ++j) { -// auto joffset = ioffset + (ao_i[1] + j) * strides[1]; - -// for(std::size_t k = 0; k < nbf[2]; ++k) { -// auto koffset = joffset + (ao_i[2] + k) * strides[2]; - -// for(std::size_t l = 0; l < nbf[3]; ++l) { -// auto loffset = koffset + (ao_i[3] + l) * strides[3]; -// out[loffset] = value[loffset] + error; -// } -// } -// } -// } -// } - -// struct Kernel { -// using shape_type = buffer::Contiguous::shape_type; -// Kernel(shape_type shape, std::array aos, -// utils::mean_type mean) : -// m_shape(std::move(shape)), m_aos(aos), m_mean(mean) {} - -// template -// Tensor operator()(const std::span t, -// const std::span error) { -// throw std::runtime_error( -// "UQ Integrals Driver kernel only supports same float types"); -// } - -// template -// auto operator()(const std::span t, -// const std::span error) { -// Tensor rv; - -// using float_type = std::decay_t; -// if constexpr(types::is_uncertain_v) { -// throw std::runtime_error("Did not expect an uncertain type"); -// } else { -// #ifdef ENABLE_SIGMA -// using tensorwrapper::buffer::make_contiguous; - -// std::array n_centers{m_aos[0].size(), m_aos[1].size(), -// m_aos[2].size(), m_aos[3].size()}; - -// std::array centers{0, 0, 0, 0}; -// std::array ao_i{0, 0, 0, 0}; -// std::array nbf{0, 0, 0, 0}; - -// using uq_type = sigma::Uncertain; -// std::vector rv_data(m_shape.size()); -// std::array strides{0, 0, 0, 1}; -// strides[2] = strides[3] * m_aos[3].n_aos(); -// strides[1] = strides[2] * m_aos[2].n_aos(); -// strides[0] = strides[1] * m_aos[1].n_aos(); - -// // If true, we can skip centers[1] > centers[0] -// bool mu_is_nu = (m_aos[0] == m_aos[1]); -// // If true, we can skip centers[3] > centers[2] -// bool lam_is_sig = (m_aos[2] == m_aos[3]); -// // If true, can skip centers[2] > centers[0] and centers[3] > -// // centers[1] -// bool all_same = (m_aos[0] == m_aos[2]) && mu_is_nu && lam_is_sig; - -// for(centers[0] = 0; centers[0] < n_centers[0]; ++centers[0]) { -// nbf[0] = m_aos[0][centers[0]].n_aos(); - -// ao_i[1] = 0; -// for(centers[1] = 0; centers[1] < n_centers[1]; ++centers[1]) -// { -// if(centers[1] > centers[0] && mu_is_nu) break; -// nbf[1] = m_aos[1][centers[1]].n_aos(); - -// ao_i[2] = 0; -// for(centers[2] = 0; centers[2] < n_centers[2]; -// ++centers[2]) { -// nbf[2] = m_aos[2][centers[2]].n_aos(); - -// ao_i[3] = 0; -// for(centers[3] = 0; centers[3] < n_centers[3]; -// ++centers[3]) { -// if(centers[3] > centers[2] && lam_is_sig) break; - -// nbf[3] = m_aos[3][centers[3]].n_aos(); - -// auto block_error = average_error( -// strides, nbf, ao_i, error, m_mean); -// uq_type max_uq{0.0, block_error}; - -// update_block(strides, nbf, ao_i, rv_data, t, -// max_uq); - -// if(mu_is_nu) { -// std::array new_nbf{ -// nbf[1], nbf[0], nbf[2], nbf[3]}; -// std::array new_ao_i{ -// ao_i[1], ao_i[0], ao_i[2], ao_i[3]}; -// copy_block(strides, nbf, ao_i, new_nbf, -// new_ao_i, rv_data); -// } - -// ao_i[3] += nbf[3]; -// } -// ao_i[2] += nbf[2]; -// } -// ao_i[1] += nbf[1]; -// } -// ao_i[0] += nbf[0]; -// } -// tensorwrapper::buffer::Contiguous t_w_contig(std::move(rv_data), -// m_shape); -// rv = tensorwrapper::Tensor(m_shape, std::move(t_w_contig)); -// #else -// throw std::runtime_error("Sigma support not enabled!"); -// #endif -// } - -// return rv; -// } -// shape_type m_shape; -// std::array m_aos; -// utils::mean_type m_mean; -// }; - -// const auto desc = R"( -// UQ Integrals Driver -// ------------------- - -// )"; - -// } // namespace - -// using eri_pt = simde::ERI4; -// using error_pt = integrals::property_types::Uncertainty; - -// MODULE_CTOR(UQAtomBlockedDriver) { -// satisfies_property_type(); -// description(desc); -// add_submodule("ERIs"); -// add_submodule("ERI Error"); -// add_input("Mean Type").set_default("none"); -// } - -// MODULE_RUN(UQAtomBlockedDriver) { -// const auto& [braket] = eri_pt::unwrap_inputs(inputs); -// auto mean_str = inputs.at("Mean Type").value(); -// auto mean = utils::mean_from_string(mean_str); - -// auto& eri_mod = submods.at("ERIs").value(); -// auto tol = eri_mod.inputs().at("Threshold").value(); - -// const auto& t = eri_mod.run_as(braket); -// const auto& error = submods.at("ERI Error").run_as(braket, -// tol); - -// using tensorwrapper::buffer::make_contiguous; -// const auto& t_buffer = make_contiguous(t.buffer()); -// const auto& e_buffer = make_contiguous(error.buffer()); - -// const auto& bra = braket.bra(); -// const auto& ket = braket.ket(); -// const auto& mu = bra.first.ao_basis_set(); -// const auto& nu = bra.second.ao_basis_set(); -// const auto& lam = ket.first.ao_basis_set(); -// const auto& sig = ket.second.ao_basis_set(); - -// std::array aos{mu, nu, lam, sig}; - -// using buffer::visit_contiguous_buffer; -// shape::Smooth shape = -// t.buffer().layout().shape().as_smooth().make_smooth(); - -// Kernel k(shape, aos, mean); -// auto t_w_error = visit_contiguous_buffer(k, t_buffer, e_buffer); - -// auto rv = results(); -// return eri_pt::wrap_results(rv, t_w_error); -// } -// } // namespace integrals::ao_integrals +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "../utils/get_permutations.hpp" +#include "../utils/uncertainty_reductions.hpp" +#include "ao_integrals.hpp" +#include +#ifdef ENABLE_SIGMA +#include +#endif +using namespace tensorwrapper; + +namespace integrals::ao_integrals { +namespace { + +template +auto average_error(T&& strides, T&& nbf, T&& ao_i, Tensor&& error, + utils::mean_type mean) { + std::vector buffer; + + for(std::size_t i = 0; i < nbf[0]; ++i) { + auto ioffset = (ao_i[0] + i) * strides[0]; + + for(std::size_t j = 0; j < nbf[1]; ++j) { + auto joffset = ioffset + (ao_i[1] + j) * strides[1]; + + for(std::size_t k = 0; k < nbf[2]; ++k) { + auto koffset = joffset + (ao_i[2] + k) * strides[2]; + + for(std::size_t l = 0; l < nbf[3]; ++l) { + auto loffset = koffset + (ao_i[3] + l) * strides[3]; + buffer.push_back(error[loffset]); + } + } + } + } + + return utils::compute_mean(mean, buffer); +} + +template +auto compute_block(T&& strides, T&& nbf, T&& ao_i, Tensor&& value, + FloatType error) { + auto n_elements = nbf[0] * nbf[1] * nbf[2] * nbf[3]; + std::vector> buffer(n_elements, error); + + for(std::size_t i = 0; i < nbf[0]; ++i) { + auto ilocal = i * nbf[1] * nbf[2] * nbf[3]; + auto ioffset = (ao_i[0] + i) * strides[0]; + + for(std::size_t j = 0; j < nbf[1]; ++j) { + auto jlocal = ilocal + j * nbf[2] * nbf[3]; + auto joffset = ioffset + (ao_i[1] + j) * strides[1]; + + for(std::size_t k = 0; k < nbf[2]; ++k) { + auto klocal = jlocal + k * nbf[3]; + auto koffset = joffset + (ao_i[2] + k) * strides[2]; + + for(std::size_t l = 0; l < nbf[3]; ++l) { + auto llocal = klocal + l; + auto loffset = koffset + (ao_i[3] + l) * strides[3]; + buffer[llocal] += value[loffset]; + } + } + } + } + return buffer; +} + +template +void set_block(T&& strides, T&& nbf, + const std::array& permuted_ao_offsets, + const std::array& sigma, + const std::vector& block, + std::vector& out) { + // sigma[d] = the original mode that for what is now mode d. Therefore, + // sigma[d] maps us back to the original mode, e.g., if the permutation + // took 0, 1, 2, 3 to 3, 2, 1, 0 then sigma[0] = 3, sigma[1] = 2, + // sigma[2] = 1, sigma[3] = 0. + + // If cidx is a 4-tuple of indices using the original modes, then for + // output dimension d the new mode is cidx[sigma[d]]. + + // Here we iterate in canonical (i,j,k,l) order — the same order the block + // was filled — and then scatter to its permuted position in out. + std::size_t block_idx = 0; + for(std::size_t i = 0; i < nbf[0]; ++i) { + for(std::size_t j = 0; j < nbf[1]; ++j) { + for(std::size_t k = 0; k < nbf[2]; ++k) { + for(std::size_t l = 0; l < nbf[3]; ++l) { + // This is the index using the original modes + std::array cidx{i, j, k, l}; + const auto out_idx = + (permuted_ao_offsets[0] + cidx[sigma[0]]) * strides[0] + + (permuted_ao_offsets[1] + cidx[sigma[1]]) * strides[1] + + (permuted_ao_offsets[2] + cidx[sigma[2]]) * strides[2] + + (permuted_ao_offsets[3] + cidx[sigma[3]]) * strides[3]; + out[out_idx] = block[block_idx++]; + } + } + } + } +} + +struct Kernel { + using shape_type = buffer::Contiguous::shape_type; + Kernel(shape_type shape, std::array aos, + utils::mean_type mean) : + m_shape(std::move(shape)), m_aos(aos), m_mean(mean) {} + + template + Tensor operator()(const std::span t, + const std::span error) { + throw std::runtime_error( + "UQ Integrals Driver kernel only supports same float types"); + } + + template + auto operator()(const std::span t, + const std::span error) { + Tensor rv; + + using float_type = std::decay_t; + if constexpr(types::is_uncertain_v) { + throw std::runtime_error("Did not expect an uncertain type"); + } else { +#ifdef ENABLE_SIGMA + using tensorwrapper::buffer::make_contiguous; + using utils::get_permutations_with_sigma; + + std::array n_centers{m_aos[0].size(), m_aos[1].size(), + m_aos[2].size(), m_aos[3].size()}; + + std::array centers{0, 0, 0, 0}; + std::array ao_offsets{0, 0, 0, 0}; + std::array nbf{0, 0, 0, 0}; + + using uq_type = sigma::Uncertain; + std::vector rv_data(m_shape.size()); + std::array strides{0, 0, 0, 1}; + strides[2] = strides[3] * m_aos[3].n_aos(); + strides[1] = strides[2] * m_aos[2].n_aos(); + strides[0] = strides[1] * m_aos[1].n_aos(); + + // If true, we can skip centers[1] > centers[0] + bool mu_is_nu = (m_aos[0] == m_aos[1]); + // If true, we can skip centers[3] > centers[2] + bool lam_is_sig = (m_aos[2] == m_aos[3]); + // If true, can skip centers[2] > centers[0] and centers[3] > + // centers[1] + bool all_same = (m_aos[0] == m_aos[2]) && mu_is_nu && lam_is_sig; + + for(centers[0] = 0; centers[0] < n_centers[0]; ++centers[0]) { + nbf[0] = m_aos[0][centers[0]].n_aos(); + + ao_offsets[1] = 0; + for(centers[1] = 0; centers[1] < n_centers[1]; ++centers[1]) { + // We restrict our bra pairs to centers[0] <= centers[1] + if(centers[1] > centers[0] && mu_is_nu) break; + nbf[1] = m_aos[1][centers[1]].n_aos(); + + ao_offsets[2] = 0; + for(centers[2] = 0; centers[2] < n_centers[2]; + ++centers[2]) { + // (c2, c3) <= (c0, c1) is impossible if c2 > c0 + if(centers[2] > centers[0] && all_same) break; + bool c2eqc0 = centers[2] == centers[0]; + nbf[2] = m_aos[2][centers[2]].n_aos(); + + ao_offsets[3] = 0; + for(centers[3] = 0; centers[3] < n_centers[3]; + ++centers[3]) { + // Restrict ket pairs to centers[2] <= centers[3] + if(centers[3] > centers[2] && lam_is_sig) break; + + nbf[3] = m_aos[3][centers[3]].n_aos(); + // Skip (c2,c3) > (c0,c1) lexicographically + bool pair_gt = (c2eqc0 && centers[3] > centers[1]); + if(pair_gt && all_same) break; + + auto block_error = average_error( + strides, nbf, ao_offsets, error, m_mean); + uq_type max_uq{0.0, block_error}; + + // Compute (ab|cd) + auto block = compute_block(strides, nbf, ao_offsets, + t, max_uq); + + // Set all symmetry equivalent blocks to `block` + auto perms = get_permutations_with_sigma( + ao_offsets, mu_is_nu, lam_is_sig, all_same); + for(auto& [perm, sigma] : perms) { + set_block(strides, nbf, perm, sigma, block, + rv_data); + } + + ao_offsets[3] += nbf[3]; + } + ao_offsets[2] += nbf[2]; + } + ao_offsets[1] += nbf[1]; + } + ao_offsets[0] += nbf[0]; + } + tensorwrapper::buffer::Contiguous t_w_contig(std::move(rv_data), + m_shape); + rv = tensorwrapper::Tensor(m_shape, std::move(t_w_contig)); +#else + throw std::runtime_error("Sigma support not enabled!"); +#endif + } + + return rv; + } + shape_type m_shape; + std::array m_aos; + utils::mean_type m_mean; +}; + +const auto desc = R"( +UQ Integrals Driver +------------------- + +)"; + +} // namespace + +using eri_pt = simde::ERI4; +using error_pt = integrals::property_types::Uncertainty; + +MODULE_CTOR(UQAtomSymmBlockedDriver) { + satisfies_property_type(); + description(desc); + add_submodule("ERIs"); + add_submodule("ERI Error"); + add_input("Mean Type").set_default("none"); +} + +MODULE_RUN(UQAtomSymmBlockedDriver) { + const auto& [braket] = eri_pt::unwrap_inputs(inputs); + auto mean_str = inputs.at("Mean Type").value(); + auto mean = utils::mean_from_string(mean_str); + + auto& eri_mod = submods.at("ERIs").value(); + auto tol = eri_mod.inputs().at("Threshold").value(); + + const auto& t = eri_mod.run_as(braket); + const auto& error = submods.at("ERI Error").run_as(braket, tol); + + using tensorwrapper::buffer::make_contiguous; + const auto& t_buffer = make_contiguous(t.buffer()); + const auto& e_buffer = make_contiguous(error.buffer()); + + const auto& bra = braket.bra(); + const auto& ket = braket.ket(); + const auto& mu = bra.first.ao_basis_set(); + const auto& nu = bra.second.ao_basis_set(); + const auto& lam = ket.first.ao_basis_set(); + const auto& sig = ket.second.ao_basis_set(); + + std::array aos{mu, nu, lam, sig}; + + using buffer::visit_contiguous_buffer; + shape::Smooth shape = t.buffer().layout().shape().as_smooth().make_smooth(); + + Kernel k(shape, aos, mean); + auto t_w_error = visit_contiguous_buffer(k, t_buffer, e_buffer); + + auto rv = results(); + return eri_pt::wrap_results(rv, t_w_error); +} +} // namespace integrals::ao_integrals diff --git a/src/integrals/utils/get_permutations.hpp b/src/integrals/utils/get_permutations.hpp new file mode 100644 index 00000000..36dce63f --- /dev/null +++ b/src/integrals/utils/get_permutations.hpp @@ -0,0 +1,93 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include +#include +#include +#include + +namespace integrals::utils { + +/** @brief Works out all the permutations of index @p abcd that are symmetry + * equivalent. + * + * @return std::set of unique permuted arrays (original API, for backward + * compatibility with existing callers and unit tests). + */ +template +auto get_permutations(T&& abcd, bool p0_1, bool p2_3, bool p01_23) { + using element_type = std::decay_t; + using index_type = std::array; + + std::set permutations{abcd}; + + if(p0_1) permutations.insert({abcd[1], abcd[0], abcd[2], abcd[3]}); + if(p2_3) permutations.insert({abcd[0], abcd[1], abcd[3], abcd[2]}); + if(p01_23) permutations.insert({abcd[2], abcd[3], abcd[0], abcd[1]}); + if(p0_1 && p01_23) + permutations.insert({abcd[2], abcd[3], abcd[1], abcd[0]}); + if(p2_3 && p01_23) + permutations.insert({abcd[3], abcd[2], abcd[0], abcd[1]}); + if(p0_1 && p2_3) permutations.insert({abcd[1], abcd[0], abcd[3], abcd[2]}); + if(p0_1 && p2_3 && p01_23) + permutations.insert({abcd[3], abcd[2], abcd[1], abcd[0]}); + + return permutations; +} + +/** @brief Like get_permutations but also returns the index permutation sigma + * for each result, where sigma[d] is the canonical dimension that + * permuted dimension d came from: perm[d] == abcd[sigma[d]]. + * + * Deduplicates by permuted value (first sigma wins), so repeated values in + * @p abcd are handled correctly. + * + * @return vector of (permuted_abcd, sigma) pairs. + */ +template +auto get_permutations_with_sigma(T&& abcd, bool p0_1, bool p2_3, bool p01_23) { + using element_type = std::decay_t; + using index_type = std::array; + using sigma_type = std::array; + using pair_type = std::pair; + + std::map seen; + auto try_insert = [&](index_type perm, sigma_type sigma) { + seen.emplace(perm, sigma); + }; + + try_insert({abcd[0], abcd[1], abcd[2], abcd[3]}, {0, 1, 2, 3}); + if(p0_1) try_insert({abcd[1], abcd[0], abcd[2], abcd[3]}, {1, 0, 2, 3}); + if(p2_3) try_insert({abcd[0], abcd[1], abcd[3], abcd[2]}, {0, 1, 3, 2}); + if(p01_23) try_insert({abcd[2], abcd[3], abcd[0], abcd[1]}, {2, 3, 0, 1}); + if(p0_1 && p01_23) + try_insert({abcd[2], abcd[3], abcd[1], abcd[0]}, {2, 3, 1, 0}); + if(p2_3 && p01_23) + try_insert({abcd[3], abcd[2], abcd[0], abcd[1]}, {3, 2, 0, 1}); + if(p0_1 && p2_3) + try_insert({abcd[1], abcd[0], abcd[3], abcd[2]}, {1, 0, 3, 2}); + if(p0_1 && p2_3 && p01_23) + try_insert({abcd[3], abcd[2], abcd[1], abcd[0]}, {3, 2, 1, 0}); + + std::vector result; + result.reserve(seen.size()); + for(auto& [perm, sigma] : seen) result.emplace_back(perm, sigma); + return result; +} + +} // namespace integrals::utils diff --git a/src/integrals/utils/uncertainty_reductions.hpp b/src/integrals/utils/uncertainty_reductions.hpp index 2ddfaf97..e44c016a 100644 --- a/src/integrals/utils/uncertainty_reductions.hpp +++ b/src/integrals/utils/uncertainty_reductions.hpp @@ -24,6 +24,13 @@ namespace integrals::utils { enum class mean_type { none, max, geometric, harmonic }; +template +auto max_mean(const ContainerType& values) { + using float_type = std::decay_t; + if(values.empty()) return float_type{0}; + return *std::max_element(values.begin(), values.end()); +} + template auto geometric_mean(const ContainerType& values) { using float_type = std::decay_t; @@ -34,26 +41,31 @@ auto geometric_mean(const ContainerType& values) { log_sum += std::log(val); ++non_zero_count; } + if(non_zero_count == 0 || log_sum == 0) return float_type{0}; return std::exp(log_sum / non_zero_count); } template auto harmonic_mean(const ContainerType& values) { + // N.b. If val is very small then 1 / val can be very large and be infinity. + // Then non_zero_count / reciprocal_sum can be NaN. using float_type = std::decay_t; float_type reciprocal_sum = 0.0; std::size_t non_zero_count = 0; for(const auto& val : values) { - if(val == 0) continue; + if(val == 0 || std::isnan(1 / val)) continue; reciprocal_sum += 1 / val; ++non_zero_count; } + if(non_zero_count == 0 || reciprocal_sum == 0) return float_type{0}; + if(std::isnan(non_zero_count / reciprocal_sum)) return float_type{0}; return non_zero_count / reciprocal_sum; } template auto compute_mean(mean_type mean, const ContainerType& span) { if(mean == mean_type::max) { - return *std::max_element(span.begin(), span.end()); + return max_mean(span); } else if(mean == mean_type::harmonic) { return harmonic_mean(span); } else if(mean == mean_type::geometric) { diff --git a/tests/cxx/unit/integrals/ao_integrals/test_uq_atom_symm_blocked_driver.cpp b/tests/cxx/unit/integrals/ao_integrals/test_uq_atom_symm_blocked_driver.cpp new file mode 100644 index 00000000..e4c777b6 --- /dev/null +++ b/tests/cxx/unit/integrals/ao_integrals/test_uq_atom_symm_blocked_driver.cpp @@ -0,0 +1,147 @@ +/* + * Copyright 2022 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "../testing/testing.hpp" + +using namespace integrals::testing; + +using namespace tensorwrapper; + +namespace { + +// N.b. The "means" of the correct values are validated by comparing to Libint's +// results. With the exception of the "No Mean" values, the +// "standard deviations" are not validated, but seem reasonable (the "No Mean" +// values come from the unit tests of the underlying error module). + +template +auto corr_answer(const simde::type::tensor& T) { + if constexpr(std::is_same_v) { + return T; + } else { + simde::type::tensor T_corr(T); + auto& corr_buffer = buffer::make_contiguous(T_corr.buffer()); + corr_buffer.set_elem({0, 0, 0, 0}, FloatType{0.774606, 0}); + corr_buffer.set_elem({0, 0, 0, 1}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({0, 0, 1, 0}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({0, 0, 1, 1}, FloatType{0.446701, 0}); + corr_buffer.set_elem({0, 1, 0, 0}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({0, 1, 0, 1}, + FloatType{0.120666, 0.0000170000000000}); + corr_buffer.set_elem({0, 1, 1, 0}, + FloatType{0.120666, 0.0000170000000000}); + corr_buffer.set_elem({0, 1, 1, 1}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 0, 0, 0}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 0, 0, 1}, + FloatType{0.120666, 0.0000170000000000}); + corr_buffer.set_elem({1, 0, 1, 0}, + FloatType{0.120666, 0.0000170000000000}); + corr_buffer.set_elem({1, 0, 1, 1}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 1, 0, 0}, FloatType{0.446701, 0}); + corr_buffer.set_elem({1, 1, 0, 1}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 1, 1, 0}, + FloatType{0.265558, 0.0000010000000000}); + corr_buffer.set_elem({1, 1, 1, 1}, FloatType{0.774606, 0}); + return T_corr; + } +} + +} // namespace + +TEST_CASE("UQ Atom Symm Blocked Driver") { + using float_type = tensorwrapper::types::udouble; + using test_pt = simde::ERI4; + using tensorwrapper::operations::approximately_equal; + + if constexpr(tensorwrapper::types::is_uncertain_v) { + auto rt = std::make_unique(); + pluginplay::ModuleManager mm(std::move(rt)); + integrals::load_modules(mm); + integrals::set_defaults(mm); + REQUIRE(mm.count("UQ Atom Symm Blocked Driver")); + mm.change_input("ERI4", "Threshold", 1.0e-6); + auto mod = mm.at("UQ Atom Blocked Driver"); + + // Make Operator + simde::type::v_ee_type op{}; + + SECTION("H2") { + // Get basis set + auto aobs = h2_sto3g_basis_set(); + + // Make AOS object + simde::type::aos aos(aobs); + simde::type::aos_squared aos_squared(aos, aos); + + // Make BraKet Input + chemist::braket::BraKet braket(aos_squared, op, aos_squared); + + SECTION("No Mean") { REQUIRE_THROWS(mod.run_as(braket)); } + SECTION("Max Error") { + mod.change_input("Mean Type", "max"); + auto T = mod.run_as(braket); + + auto T_corr = corr_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + SECTION("Geometric Mean") { + mod.change_input("Mean Type", "geometric"); + auto T = mod.run_as(braket); + + auto T_corr = corr_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + SECTION("Harmonic Mean") { + mod.change_input("Mean Type", "harmonic"); + auto T = mod.run_as(braket); + + auto T_corr = corr_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + } + + SECTION("H2 Dimer") { + simde::type::nucleus h0("H", 1ul, 1836.15, 0.0, 0.0, 0.0); + simde::type::nucleus h1("H", 1ul, 1836.15, 0.0, 0.0, 1.39839); + simde::type::nucleus h2("H", 1ul, 1836.15, 0.0, 0.0, 4.39839); + simde::type::nucleus h3("H", 1ul, 1836.15, 0.0, 0.0, 5.79678); + auto aobs = apply_h2_sto3g_basis_set(std::vector{h0, h1, h2, h3}); + + // Make AOS object + simde::type::aos aos(aobs); + simde::type::aos_squared aos_squared(aos, aos); + + // Make BraKet Input + chemist::braket::BraKet braket(aos_squared, op, aos_squared); + const auto corr_key = "UQ Atom Blocked Driver"; + auto corr_mod = mm.at(corr_key); + SECTION("Max Error") { + mod.change_input("Mean Type", "max"); + corr_mod.change_input("Mean Type", "max"); + auto eris = mod.run_as(braket); + auto eris_corr = corr_mod.run_as(braket); + REQUIRE(approximately_equal(eris_corr, eris, 1E-6)); + } + } + } +} diff --git a/tests/cxx/unit/integrals/testing/ao_bases.hpp b/tests/cxx/unit/integrals/testing/ao_bases.hpp index f423fbb0..09284038 100644 --- a/tests/cxx/unit/integrals/testing/ao_bases.hpp +++ b/tests/cxx/unit/integrals/testing/ao_bases.hpp @@ -20,32 +20,32 @@ namespace integrals::testing { -inline simde::type::ao_basis_set h2_sto3g_basis_set() { +template +inline simde::type::ao_basis_set apply_h2_sto3g_basis_set( + const MoleculeType& mol) { using ao_basis_t = simde::type::ao_basis_set; using atomic_basis_t = simde::type::atomic_basis_set; using cg_t = simde::type::contracted_gaussian; - using point_t = simde::type::point; using doubles_t = std::vector; - auto mol = water_molecule(); - point_t r0 = mol[0].as_nucleus(); - point_t r1 = mol[1].as_nucleus(); - - doubles_t cs{0.1543289673, 0.5353281423, 0.4446345422}; - doubles_t es{3.425250914, 0.6239137298, 0.1688554040}; - cg_t cg0(cs.begin(), cs.end(), es.begin(), es.end(), r0); - cg_t cg1(cs.begin(), cs.end(), es.begin(), es.end(), r1); - atomic_basis_t h0("sto-3g", 1, r0); - atomic_basis_t h1("sto-3g", 1, r1); - h0.add_shell(chemist::ShellType::cartesian, 0, cg0); - h1.add_shell(chemist::ShellType::cartesian, 0, cg1); - ao_basis_t bs; - bs.add_center(h0); - bs.add_center(h1); + for(const auto& atom : mol) { + doubles_t cs{0.1543289673, 0.5353281423, 0.4446345422}; + doubles_t es{3.425250914, 0.6239137298, 0.1688554040}; + cg_t cg0(cs.begin(), cs.end(), es.begin(), es.end(), atom); + atomic_basis_t h0("sto-3g", 1, atom); + h0.add_shell(chemist::ShellType::cartesian, 0, cg0); + bs.add_center(h0); + } return bs; } +inline simde::type::ao_basis_set h2_sto3g_basis_set() { + auto mol = water_molecule(); + std::vector h2{mol[0].as_nucleus(), mol[1].as_nucleus()}; + return apply_h2_sto3g_basis_set(h2); +} + inline simde::type::ao_basis_set water_sto3g_basis_set() { using ao_basis_t = simde::type::ao_basis_set; using atomic_basis_t = simde::type::atomic_basis_set; diff --git a/tests/cxx/unit/integrals/utils/get_permutations.cpp b/tests/cxx/unit/integrals/utils/get_permutations.cpp new file mode 100644 index 00000000..90e10be2 --- /dev/null +++ b/tests/cxx/unit/integrals/utils/get_permutations.cpp @@ -0,0 +1,364 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "../testing/testing.hpp" +#include + +using namespace integrals::utils; + +/* Testing Notes: + * + * We assume that std::set works and thus won't allow us to insert duplicates. + * Thus a spot check for an input with all unique values, and one with a + * repeated index should suffice. This leaves 2^3 = 8 combinations of + * symmetry to check. + */ + +TEST_CASE("get_permutations") { + using type = std::array; + SECTION("No symmetry") { + SECTION("All unique") { + auto perms = + get_permutations(type{0, 1, 2, 3}, false, false, false); + REQUIRE(perms.size() == 1); + REQUIRE(perms.count({0, 1, 2, 3}) == 1); + } + SECTION("Three repeats") { + auto perms = + get_permutations(type{0, 0, 0, 1}, false, false, false); + REQUIRE(perms.size() == 1); + REQUIRE(perms.count({0, 0, 0, 1}) == 1); + } + } + SECTION("p0_1 symmetry") { + SECTION("All unique") { + auto perms = get_permutations(type{0, 1, 2, 3}, true, false, false); + REQUIRE(perms.size() == 2); + REQUIRE(perms.count({0, 1, 2, 3}) == 1); + REQUIRE(perms.count({1, 0, 2, 3}) == 1); + } + SECTION("One repeat (in symmetry)") { + auto perms = get_permutations(type{0, 0, 1, 2}, true, false, false); + REQUIRE(perms.size() == 1); + REQUIRE(perms.count({0, 0, 1, 2}) == 1); + } + SECTION("One repeat (not in symmetry)") { + auto perms = get_permutations(type{0, 1, 1, 2}, true, false, false); + REQUIRE(perms.size() == 2); + REQUIRE(perms.count({0, 1, 1, 2}) == 1); + REQUIRE(perms.count({1, 0, 1, 2}) == 1); + } + } + + SECTION("p2_3 symmetry") { + // Conceptually very similar to p0_1 symmetry so just spot check + auto perms = get_permutations(type{0, 1, 2, 3}, false, true, false); + REQUIRE(perms.size() == 2); + REQUIRE(perms.count({0, 1, 2, 3}) == 1); + REQUIRE(perms.count({0, 1, 3, 2}) == 1); + } + + SECTION("p01_23") { + // Conceptually very similar to p0_1 symmetry so just spot check + auto perms = get_permutations(type{0, 1, 2, 3}, false, false, true); + REQUIRE(perms.size() == 2); + REQUIRE(perms.count({0, 1, 2, 3}) == 1); + REQUIRE(perms.count({2, 3, 0, 1}) == 1); + } + + SECTION("P0_1 and p2_3") { + SECTION("All unique") { + auto perms = get_permutations(type{0, 1, 2, 3}, true, true, false); + REQUIRE(perms.size() == 4); + REQUIRE(perms.count({0, 1, 2, 3}) == 1); + REQUIRE(perms.count({1, 0, 2, 3}) == 1); + REQUIRE(perms.count({0, 1, 3, 2}) == 1); + REQUIRE(perms.count({1, 0, 3, 2}) == 1); + } + SECTION("One repeat") { + auto perms = get_permutations(type{0, 0, 1, 2}, true, true, false); + REQUIRE(perms.size() == 2); + REQUIRE(perms.count({0, 0, 1, 2}) == 1); + REQUIRE(perms.count({0, 0, 2, 1}) == 1); + } + SECTION("Two repeats (one in each symmetry group)") { + auto perms = get_permutations(type{0, 0, 1, 1}, true, true, false); + REQUIRE(perms.size() == 1); + REQUIRE(perms.count({0, 0, 1, 1}) == 1); + } + } + SECTION("P0_1 and p01_23") { + // Conceptually very similar to p0_1 && p2_3 symmetry so just spot check + auto perms = get_permutations(type{0, 1, 2, 3}, true, false, true); + REQUIRE(perms.size() == 4); + REQUIRE(perms.count({0, 1, 2, 3}) == 1); + REQUIRE(perms.count({1, 0, 2, 3}) == 1); + REQUIRE(perms.count({2, 3, 0, 1}) == 1); + REQUIRE(perms.count({2, 3, 1, 0}) == 1); + } + SECTION("P2_3 and p01_23") { + // Conceptually very similar to p0_1 && p2_3 symmetry so just spot check + auto perms = get_permutations(type{0, 1, 2, 3}, false, true, true); + REQUIRE(perms.size() == 4); + REQUIRE(perms.count({0, 1, 2, 3}) == 1); + REQUIRE(perms.count({0, 1, 3, 2}) == 1); + REQUIRE(perms.count({2, 3, 0, 1}) == 1); + REQUIRE(perms.count({3, 2, 0, 1}) == 1); + } + SECTION("P0_1 and p2_3 and p01_23") { + SECTION("All unique") { + auto perms = get_permutations(type{0, 1, 2, 3}, true, true, true); + REQUIRE(perms.size() == 8); + REQUIRE(perms.count({0, 1, 2, 3}) == 1); + REQUIRE(perms.count({1, 0, 2, 3}) == 1); + REQUIRE(perms.count({0, 1, 3, 2}) == 1); + REQUIRE(perms.count({1, 0, 3, 2}) == 1); + REQUIRE(perms.count({2, 3, 0, 1}) == 1); + REQUIRE(perms.count({3, 2, 0, 1}) == 1); + REQUIRE(perms.count({2, 3, 1, 0}) == 1); + REQUIRE(perms.count({3, 2, 1, 0}) == 1); + } + SECTION("One repeat") { + auto perms = get_permutations(type{0, 0, 1, 2}, true, true, true); + REQUIRE(perms.size() == 4); + REQUIRE(perms.count({0, 0, 1, 2}) == 1); + REQUIRE(perms.count({0, 0, 2, 1}) == 1); + REQUIRE(perms.count({1, 2, 0, 0}) == 1); + REQUIRE(perms.count({2, 1, 0, 0}) == 1); + } + } +} + +/* Testing Notes for get_permutations_with_sigma: + * + * For every returned (perm, sigma) pair the invariant perm[d] == abcd[sigma[d]] + * must hold for all d in {0,1,2,3}. We verify this invariant explicitly for + * every pair in every test case, in addition to checking the expected set of + * permuted arrays and their counts. + * + * Deduplication behaviour: when two symmetry operations produce the same + * permuted array (because abcd has repeated values) only one entry is kept and + * its sigma is the one that was inserted first (identity wins over later + * operations). + */ + +namespace { +// Helper: verify perm[d] == abcd[sigma[d]] for all d +template +bool sigma_consistent(const Index& abcd, const Index& perm, + const Sigma& sigma) { + for(std::size_t d = 0; d < 4; ++d) + if(perm[d] != abcd[sigma[d]]) return false; + return true; +} + +// Helper: find a (perm, sigma) pair by its perm value +template +auto find_perm(const Pairs& pairs, const Index& target) { + for(auto& [p, s] : pairs) + if(p == target) return s; + // Return identity sigma as sentinel — test will fail on sigma check + return std::array{0, 1, 2, 3}; +} + +// Helper: count how many pairs have a given perm value +template +std::size_t count_perm(const Pairs& pairs, const Index& target) { + std::size_t n = 0; + for(auto& [p, s] : pairs) + if(p == target) ++n; + return n; +} +} // namespace + +TEST_CASE("get_permutations_with_sigma") { + using type = std::array; + using sigma_type = std::array; + + SECTION("No symmetry") { + SECTION("All unique") { + type abcd{0, 1, 2, 3}; + auto pairs = get_permutations_with_sigma(abcd, false, false, false); + REQUIRE(pairs.size() == 1); + REQUIRE(count_perm(pairs, type{0, 1, 2, 3}) == 1); + auto s = find_perm(pairs, type{0, 1, 2, 3}); + REQUIRE(s == sigma_type{0, 1, 2, 3}); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + SECTION("Three repeats") { + type abcd{0, 0, 0, 1}; + auto pairs = get_permutations_with_sigma(abcd, false, false, false); + REQUIRE(pairs.size() == 1); + REQUIRE(count_perm(pairs, type{0, 0, 0, 1}) == 1); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + } + + SECTION("p0_1 symmetry") { + SECTION("All unique") { + type abcd{0, 1, 2, 3}; + auto pairs = get_permutations_with_sigma(abcd, true, false, false); + REQUIRE(pairs.size() == 2); + REQUIRE(count_perm(pairs, type{0, 1, 2, 3}) == 1); + REQUIRE(count_perm(pairs, type{1, 0, 2, 3}) == 1); + REQUIRE(find_perm(pairs, type{0, 1, 2, 3}) == + sigma_type{0, 1, 2, 3}); + REQUIRE(find_perm(pairs, type{1, 0, 2, 3}) == + sigma_type{1, 0, 2, 3}); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + SECTION("One repeat (in symmetry) — deduplicates to 1 entry") { + type abcd{0, 0, 1, 2}; + auto pairs = get_permutations_with_sigma(abcd, true, false, false); + REQUIRE(pairs.size() == 1); + REQUIRE(count_perm(pairs, type{0, 0, 1, 2}) == 1); + // Identity sigma wins + REQUIRE(find_perm(pairs, type{0, 0, 1, 2}) == + sigma_type{0, 1, 2, 3}); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + SECTION("One repeat (not in symmetry)") { + type abcd{0, 1, 1, 2}; + auto pairs = get_permutations_with_sigma(abcd, true, false, false); + REQUIRE(pairs.size() == 2); + REQUIRE(count_perm(pairs, type{0, 1, 1, 2}) == 1); + REQUIRE(count_perm(pairs, type{1, 0, 1, 2}) == 1); + REQUIRE(find_perm(pairs, type{0, 1, 1, 2}) == + sigma_type{0, 1, 2, 3}); + REQUIRE(find_perm(pairs, type{1, 0, 1, 2}) == + sigma_type{1, 0, 2, 3}); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + } + + SECTION("p2_3 symmetry") { + type abcd{0, 1, 2, 3}; + auto pairs = get_permutations_with_sigma(abcd, false, true, false); + REQUIRE(pairs.size() == 2); + REQUIRE(count_perm(pairs, type{0, 1, 2, 3}) == 1); + REQUIRE(count_perm(pairs, type{0, 1, 3, 2}) == 1); + REQUIRE(find_perm(pairs, type{0, 1, 2, 3}) == sigma_type{0, 1, 2, 3}); + REQUIRE(find_perm(pairs, type{0, 1, 3, 2}) == sigma_type{0, 1, 3, 2}); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + + SECTION("p01_23 symmetry") { + type abcd{0, 1, 2, 3}; + auto pairs = get_permutations_with_sigma(abcd, false, false, true); + REQUIRE(pairs.size() == 2); + REQUIRE(count_perm(pairs, type{0, 1, 2, 3}) == 1); + REQUIRE(count_perm(pairs, type{2, 3, 0, 1}) == 1); + REQUIRE(find_perm(pairs, type{0, 1, 2, 3}) == sigma_type{0, 1, 2, 3}); + REQUIRE(find_perm(pairs, type{2, 3, 0, 1}) == sigma_type{2, 3, 0, 1}); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + + SECTION("p0_1 and p2_3") { + SECTION("All unique") { + type abcd{0, 1, 2, 3}; + auto pairs = get_permutations_with_sigma(abcd, true, true, false); + REQUIRE(pairs.size() == 4); + REQUIRE(count_perm(pairs, type{0, 1, 2, 3}) == 1); + REQUIRE(count_perm(pairs, type{1, 0, 2, 3}) == 1); + REQUIRE(count_perm(pairs, type{0, 1, 3, 2}) == 1); + REQUIRE(count_perm(pairs, type{1, 0, 3, 2}) == 1); + REQUIRE(find_perm(pairs, type{0, 1, 2, 3}) == + sigma_type{0, 1, 2, 3}); + REQUIRE(find_perm(pairs, type{1, 0, 2, 3}) == + sigma_type{1, 0, 2, 3}); + REQUIRE(find_perm(pairs, type{0, 1, 3, 2}) == + sigma_type{0, 1, 3, 2}); + REQUIRE(find_perm(pairs, type{1, 0, 3, 2}) == + sigma_type{1, 0, 3, 2}); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + SECTION("One repeat — deduplicates") { + type abcd{0, 0, 1, 2}; + auto pairs = get_permutations_with_sigma(abcd, true, true, false); + REQUIRE(pairs.size() == 2); + REQUIRE(count_perm(pairs, type{0, 0, 1, 2}) == 1); + REQUIRE(count_perm(pairs, type{0, 0, 2, 1}) == 1); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + SECTION( + "Two repeats (one in each symmetry group) — deduplicates to 1") { + type abcd{0, 0, 1, 1}; + auto pairs = get_permutations_with_sigma(abcd, true, true, false); + REQUIRE(pairs.size() == 1); + REQUIRE(count_perm(pairs, type{0, 0, 1, 1}) == 1); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + } + + SECTION("p0_1 and p01_23") { + type abcd{0, 1, 2, 3}; + auto pairs = get_permutations_with_sigma(abcd, true, false, true); + REQUIRE(pairs.size() == 4); + REQUIRE(count_perm(pairs, type{0, 1, 2, 3}) == 1); + REQUIRE(count_perm(pairs, type{1, 0, 2, 3}) == 1); + REQUIRE(count_perm(pairs, type{2, 3, 0, 1}) == 1); + REQUIRE(count_perm(pairs, type{2, 3, 1, 0}) == 1); + REQUIRE(find_perm(pairs, type{0, 1, 2, 3}) == sigma_type{0, 1, 2, 3}); + REQUIRE(find_perm(pairs, type{1, 0, 2, 3}) == sigma_type{1, 0, 2, 3}); + REQUIRE(find_perm(pairs, type{2, 3, 0, 1}) == sigma_type{2, 3, 0, 1}); + REQUIRE(find_perm(pairs, type{2, 3, 1, 0}) == sigma_type{2, 3, 1, 0}); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + + SECTION("p2_3 and p01_23") { + type abcd{0, 1, 2, 3}; + auto pairs = get_permutations_with_sigma(abcd, false, true, true); + REQUIRE(pairs.size() == 4); + REQUIRE(count_perm(pairs, type{0, 1, 2, 3}) == 1); + REQUIRE(count_perm(pairs, type{0, 1, 3, 2}) == 1); + REQUIRE(count_perm(pairs, type{2, 3, 0, 1}) == 1); + REQUIRE(count_perm(pairs, type{3, 2, 0, 1}) == 1); + REQUIRE(find_perm(pairs, type{0, 1, 2, 3}) == sigma_type{0, 1, 2, 3}); + REQUIRE(find_perm(pairs, type{0, 1, 3, 2}) == sigma_type{0, 1, 3, 2}); + REQUIRE(find_perm(pairs, type{2, 3, 0, 1}) == sigma_type{2, 3, 0, 1}); + REQUIRE(find_perm(pairs, type{3, 2, 0, 1}) == sigma_type{3, 2, 0, 1}); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + + SECTION("p0_1 and p2_3 and p01_23") { + SECTION("All unique — all 8 permutations") { + type abcd{0, 1, 2, 3}; + auto pairs = get_permutations_with_sigma(abcd, true, true, true); + REQUIRE(pairs.size() == 8); + const std::vector> expected{ + {{0, 1, 2, 3}, {0, 1, 2, 3}}, {{1, 0, 2, 3}, {1, 0, 2, 3}}, + {{0, 1, 3, 2}, {0, 1, 3, 2}}, {{1, 0, 3, 2}, {1, 0, 3, 2}}, + {{2, 3, 0, 1}, {2, 3, 0, 1}}, {{2, 3, 1, 0}, {2, 3, 1, 0}}, + {{3, 2, 0, 1}, {3, 2, 0, 1}}, {{3, 2, 1, 0}, {3, 2, 1, 0}}, + }; + for(auto& [exp_perm, exp_sigma] : expected) { + REQUIRE(count_perm(pairs, exp_perm) == 1); + REQUIRE(find_perm(pairs, exp_perm) == exp_sigma); + } + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + SECTION("One repeat — deduplicates to 4 entries") { + type abcd{0, 0, 1, 2}; + auto pairs = get_permutations_with_sigma(abcd, true, true, true); + REQUIRE(pairs.size() == 4); + REQUIRE(count_perm(pairs, type{0, 0, 1, 2}) == 1); + REQUIRE(count_perm(pairs, type{0, 0, 2, 1}) == 1); + REQUIRE(count_perm(pairs, type{1, 2, 0, 0}) == 1); + REQUIRE(count_perm(pairs, type{2, 1, 0, 0}) == 1); + for(auto& [p, sig] : pairs) REQUIRE(sigma_consistent(abcd, p, sig)); + } + } +} diff --git a/tests/cxx/unit/integrals/utils/uncertainty_reductions.cpp b/tests/cxx/unit/integrals/utils/uncertainty_reductions.cpp index fa49b577..44bf60fc 100644 --- a/tests/cxx/unit/integrals/utils/uncertainty_reductions.cpp +++ b/tests/cxx/unit/integrals/utils/uncertainty_reductions.cpp @@ -18,6 +18,29 @@ #include using namespace integrals::utils; +TEST_CASE("max_mean") { + SECTION("Basic Test") { + std::vector values{1.0, 10.0, 100.0}; + auto result = max_mean(values); + REQUIRE(result == Catch::Approx(100.0)); + } + SECTION("Contains zero") { + std::vector values{1.0, 10.0, 100.0, 0.0}; + auto result = max_mean(values); + REQUIRE(result == Catch::Approx(100.0)); + } + SECTION("All zero") { + std::vector values{0.0, 0.0, 0.0}; + auto result = max_mean(values); + REQUIRE(result == Catch::Approx(0.0)); + } + SECTION("Empty") { + std::vector values{}; + auto result = max_mean(values); + REQUIRE(result == Catch::Approx(0.0)); + } +} + TEST_CASE("geometric_mean") { SECTION("Basic Test") { std::vector values{1.0, 10.0, 100.0}; @@ -30,6 +53,16 @@ TEST_CASE("geometric_mean") { auto result = geometric_mean(values); REQUIRE(result == Catch::Approx(10.0)); } + SECTION("All zero") { + std::vector values{0.0, 0.0, 0.0}; + auto result = geometric_mean(values); + REQUIRE(result == Catch::Approx(0.0)); + } + SECTION("No elements") { + std::vector values{}; + auto result = geometric_mean(values); + REQUIRE(result == Catch::Approx(0.0)); + } } TEST_CASE("harmonic_mean") { @@ -44,6 +77,18 @@ TEST_CASE("harmonic_mean") { auto result = harmonic_mean(values); REQUIRE(result == Catch::Approx(1.63636363636)); } + + SECTION("All zero") { + std::vector values{0.0, 0.0, 0.0}; + auto result = harmonic_mean(values); + REQUIRE(result == Catch::Approx(0.0)); + } + + SECTION("No elements") { + std::vector values{}; + auto result = harmonic_mean(values); + REQUIRE(result == Catch::Approx(0.0)); + } } TEST_CASE("compute_mean") {