From f7f336f47bc95ba6bc2bbeca0e5808765b9c5472 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 7 Feb 2025 14:52:57 -0600 Subject: [PATCH 1/5] adds sigma support --- CMakeLists.txt | 2 + src/integrals/ao_integrals/ao_integrals.cpp | 145 +++++++++--------- .../ao_integrals/lib_int_visitor.hpp | 57 +++++++ src/integrals/uncertain_types.hpp | 16 ++ .../ao_integrals/test_ao_integrals.hpp | 17 +- .../ao_integrals/test_arbitrary_operator.cpp | 113 ++++++++++---- 6 files changed, 241 insertions(+), 109 deletions(-) create mode 100644 src/integrals/ao_integrals/lib_int_visitor.hpp create mode 100644 src/integrals/uncertain_types.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index c55816cf..805b326f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,7 @@ nwx_cxx_api_docs("${INTEGRALS_SOURCE_DIR}" "${INTEGRALS_INCLUDE_DIR}") cmaize_option_list( BUILD_TESTING OFF "Should we build the tests?" BUILD_PYBIND11_PYBINDINGS ON "Build pybind11 python3 bindings?" + ENABLE_SIGMA OFF "Should we enable Sigma for uncertainty tracking?" ) # Can't build from github due to Libint setup. @@ -56,6 +57,7 @@ cmaize_find_or_build_dependency( FIND_TARGET nwx::simde CMAKE_ARGS BUILD_TESTING=OFF BUILD_PYBIND11_PYBINDINGS=${BUILD_PYBIND11_PYBINDINGS} + BUILD_SIGMA=${BUILD_SIGMA} ) cmaize_add_library( diff --git a/src/integrals/ao_integrals/ao_integrals.cpp b/src/integrals/ao_integrals/ao_integrals.cpp index b477759c..03b8ecb3 100644 --- a/src/integrals/ao_integrals/ao_integrals.cpp +++ b/src/integrals/ao_integrals/ao_integrals.cpp @@ -13,99 +13,52 @@ * See the License for the specific language governing permissions and * limitations under the License. */ - +#include "../uncertain_types.hpp" #include "ao_integrals.hpp" #include "detail_/get_basis_sets.hpp" #include "detail_/libint_op.hpp" #include "detail_/make_engine.hpp" #include "detail_/make_libint_basis_set.hpp" #include "detail_/shells2ord.hpp" +#include "lib_int_visitor.hpp" #include namespace integrals::ao_integrals { - -class LibIntVisitor : public chemist::qm_operator::OperatorVisitor { -public: - using s_e_type = simde::type::s_e_type; - using t_e_type = simde::type::t_e_type; - using v_ee_type = simde::type::v_ee_type; - using v_en_type = simde::type::v_en_type; - - LibIntVisitor(const std::vector& bases, double thresh, - std::size_t deriv = 0) : - m_bases(bases), m_thresh(thresh), m_deriv(deriv){}; - - void run(const s_e_type& S_e) { - m_engine = detail_::make_engine(m_bases, S_e, m_thresh, m_deriv); - } - - void run(const t_e_type& T_e) { - m_engine = detail_::make_engine(m_bases, T_e, m_thresh, m_deriv); - } - - void run(const v_en_type& V_en) { - m_engine = detail_::make_engine(m_bases, V_en, m_thresh, m_deriv); - } - - void run(const v_ee_type& V_ee) { - m_engine = detail_::make_engine(m_bases, V_ee, m_thresh, m_deriv); +namespace { + +template +auto build_eigen_buffer(const std::vector& basis_sets, + double thresh) { + FloatType initial_value; + if constexpr(std::is_same_v) { + initial_value = 0.0; + } else { // Presumably sigma::UDouble + initial_value = FloatType(0.0, thresh); } - - libint2::Engine& engine() { return m_engine; } - -private: - const std::vector& m_bases; - double m_thresh; - std::size_t m_deriv; - libint2::Engine m_engine; -}; - -template -TEMPLATED_MODULE_CTOR(AOIntegral, BraKetType) { - using my_pt = simde::EvaluateBraKet; - satisfies_property_type(); - description("Computes integrals with Libint"); - - add_input("Threshold") - .set_default(1.0E-16) - .set_description( - "The target precision with which the integrals will be computed"); -} - -template -TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) { - using my_pt = simde::EvaluateBraKet; - - const auto& [braket] = my_pt::unwrap_inputs(inputs); - auto thresh = inputs.at("Threshold").value(); - auto bra = braket.bra(); - auto ket = braket.ket(); - auto& op = braket.op(); - - // Gather information from Bra, Ket, and Op - auto basis_sets = detail_::get_basis_sets(bra, ket); - constexpr int N = detail_::get_n(bra, ket); - - // Dimensional information - std::vector dims_shells(N); Eigen::array dims_bfs; - for(auto i = 0; i < N; ++i) { - dims_shells[i] = basis_sets[i].size(); - dims_bfs[i] = basis_sets[i].nbf(); - } + for(auto i = 0; i < N; ++i) dims_bfs[i] = basis_sets[i].nbf(); - // Build tensor inputs - using tensor_t = simde::type::tensor; using shape_t = tensorwrapper::shape::Smooth; using layout_t = tensorwrapper::layout::Physical; - using buffer_t = tensorwrapper::buffer::Eigen; + using buffer_t = tensorwrapper::buffer::Eigen; using data_t = typename buffer_t::data_type; shape_t s{dims_bfs.begin(), dims_bfs.end()}; layout_t l(s); data_t d(dims_bfs); buffer_t b{d, l}; - b.value().setZero(); + b.value().setConstant(initial_value); + return b; +} + +template +auto fill_tensor(const std::vector& basis_sets, + const chemist::qm_operator::OperatorBase& op, double thresh) { + // Dimensional information + std::vector dims_shells(N); + for(auto i = 0; i < N; ++i) dims_shells[i] = basis_sets[i].size(); + + auto b = build_eigen_buffer(basis_sets, thresh); // Make libint engine LibIntVisitor visitor(basis_sets, thresh); @@ -140,7 +93,51 @@ TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) { } } - tensor_t t(s, b); + return simde::type::tensor(b.layout().shape().clone(), b); +} + +} // namespace + +template +TEMPLATED_MODULE_CTOR(AOIntegral, BraKetType) { + using my_pt = simde::EvaluateBraKet; + satisfies_property_type(); + description("Driver for computing integrals with Libint"); + + add_input("Threshold") + .set_default(1.0E-16) + .set_description( + "The target precision with which the integrals will be computed"); + + add_input("With UQ?").set_default(false); +} + +template +TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) { + using my_pt = simde::EvaluateBraKet; + + const auto& [braket] = my_pt::unwrap_inputs(inputs); + auto thresh = inputs.at("Threshold").value(); + auto with_uq = inputs.at("With UQ?").value(); + auto bra = braket.bra(); + auto ket = braket.ket(); + auto& op = braket.op(); + + // Gather information from Bra, Ket, and Op + auto basis_sets = detail_::get_basis_sets(bra, ket); + constexpr int N = detail_::get_n(bra, ket); + + simde::type::tensor t; + if(with_uq) { + if constexpr(!std::is_same_v) { + t = fill_tensor(basis_sets, op, thresh); + } else { + throw std::runtime_error("Sigma support not enabled!"); + } + } else { + t = fill_tensor(basis_sets, op, thresh); + } + auto rv = results(); return my_pt::wrap_results(rv, t); } diff --git a/src/integrals/ao_integrals/lib_int_visitor.hpp b/src/integrals/ao_integrals/lib_int_visitor.hpp new file mode 100644 index 00000000..d03778e5 --- /dev/null +++ b/src/integrals/ao_integrals/lib_int_visitor.hpp @@ -0,0 +1,57 @@ +/* + * Copyright 2024 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 + +namespace integrals::ao_integrals { + +class LibIntVisitor : public chemist::qm_operator::OperatorVisitor { +public: + using s_e_type = simde::type::s_e_type; + using t_e_type = simde::type::t_e_type; + using v_ee_type = simde::type::v_ee_type; + using v_en_type = simde::type::v_en_type; + + LibIntVisitor(const std::vector& bases, double thresh, + std::size_t deriv = 0) : + m_bases(bases), m_thresh(thresh), m_deriv(deriv) {}; + + void run(const s_e_type& S_e) { + m_engine = detail_::make_engine(m_bases, S_e, m_thresh, m_deriv); + } + + void run(const t_e_type& T_e) { + m_engine = detail_::make_engine(m_bases, T_e, m_thresh, m_deriv); + } + + void run(const v_en_type& V_en) { + m_engine = detail_::make_engine(m_bases, V_en, m_thresh, m_deriv); + } + + void run(const v_ee_type& V_ee) { + m_engine = detail_::make_engine(m_bases, V_ee, m_thresh, m_deriv); + } + + libint2::Engine& engine() { return m_engine; } + +private: + const std::vector& m_bases; + double m_thresh; + std::size_t m_deriv; + libint2::Engine m_engine; +}; + +} // namespace integrals::ao_integrals \ No newline at end of file diff --git a/src/integrals/uncertain_types.hpp b/src/integrals/uncertain_types.hpp new file mode 100644 index 00000000..ca6b4d8c --- /dev/null +++ b/src/integrals/uncertain_types.hpp @@ -0,0 +1,16 @@ +#pragma once + +#ifdef ENABLE_SIGMA +#include +#endif + +namespace integrals::type { +#ifdef ENABLE_SIGMA +using uncertain_float = sigma::UFloat; +using uncertain_double = sigma::UDouble; +#else +using uncertain_float = void; +using uncertain_double = void; +#endif + +} // namespace integrals::type \ No newline at end of file diff --git a/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp b/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp index d11b1a15..902445aa 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp @@ -21,20 +21,21 @@ namespace test { -template +template auto eigen_buffer(const tensorwrapper::buffer::BufferBase& buffer) { - return static_cast&>(buffer); + return static_cast&>( + buffer); } -template -auto trace(const TensorType& t) { - Eigen::Tensor trace = t.value().trace(); +template +auto trace(const tensorwrapper::buffer::Eigen& t) { + Eigen::Tensor trace = t.value().trace(); return trace.coeff(); } -template -auto norm(const TensorType& t) { - Eigen::Tensor norm = +template +auto norm(const tensorwrapper::buffer::Eigen& t) { + Eigen::Tensor norm = t.value().square().sum().sqrt(); return norm.coeff(); } diff --git a/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp b/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp index 1b63467f..4fe28bf4 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp @@ -14,8 +14,12 @@ * limitations under the License. */ +#include "integrals/uncertain_types.hpp" #include "test_ao_integrals.hpp" +using udouble = integrals::type::uncertain_double; +constexpr bool has_sigma = !std::is_same_v; + TEST_CASE("OperatorBase") { using aos_t = simde::type::aos; using aos_squared_t = simde::type::aos_squared; @@ -48,16 +52,33 @@ TEST_CASE("OperatorBase") { braket_t braket(aos, op_base, aos); // Call module - auto T = mm.at("Evaluate 2-Index BraKet").run_as(braket); - - // Check output - auto t = test::eigen_buffer<2>(T.buffer()); - REQUIRE(test::trace(t) == - Catch::Approx(124.7011973877891364).margin(1.0e-16)); - REQUIRE(test::norm(t) == - Catch::Approx(90.2562579028763707).margin(1.0e-16)); + auto& mod = mm.at("Evaluate 2-Index BraKet"); + + SECTION("No UQ") { + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<2>(T.buffer()); + REQUIRE(test::trace(t) == + Catch::Approx(124.7011973877891364).margin(1.0e-16)); + REQUIRE(test::norm(t) == + Catch::Approx(90.2562579028763707).margin(1.0e-16)); + } + + SECTION("With UQ") { + if constexpr(has_sigma) { + mod.change_input("With UQ?", true); + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<2, udouble>(T.buffer()); + REQUIRE(test::trace(t).mean() == + Catch::Approx(124.7011973877891364).margin(1.0e-16)); + REQUIRE(test::norm(t).mean() == + Catch::Approx(90.2562579028763707).margin(1.0e-16)); + } + } } - SECTION("3-Index") { using braket_t = simde::type::braket; using test_pt = simde::EvaluateBraKet; @@ -65,15 +86,34 @@ TEST_CASE("OperatorBase") { // Make BraKet Input braket_t braket(aos, op_base, aos_squared); - // Call module - auto T = mm.at("Evaluate 3-Index BraKet").run_as(braket); - - // Check output - auto t = test::eigen_buffer<3>(T.buffer()); - REQUIRE(test::trace(t) == - Catch::Approx(16.8245948391706577).margin(1.0e-16)); - REQUIRE(test::norm(t) == - Catch::Approx(20.6560572032543597).margin(1.0e-16)); + auto& mod = mm.at("Evaluate 3-Index BraKet"); + + SECTION("No UQ") { + // Call module + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<3>(T.buffer()); + REQUIRE(test::trace(t) == + Catch::Approx(16.8245948391706577).margin(1.0e-16)); + REQUIRE(test::norm(t) == + Catch::Approx(20.6560572032543597).margin(1.0e-16)); + } + + SECTION("With UQ") { + if constexpr(has_sigma) { + mod.change_input("With UQ?", true); + // Call module + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<3, udouble>(T.buffer()); + REQUIRE(test::trace(t).mean() == + Catch::Approx(16.8245948391706577).margin(1.0e-16)); + REQUIRE(test::norm(t).mean() == + Catch::Approx(20.6560572032543597).margin(1.0e-16)); + } + } } SECTION("4-Index") { @@ -84,14 +124,33 @@ TEST_CASE("OperatorBase") { // Make BraKet Input braket_t braket(aos_squared, op_base, aos_squared); - // Call module - auto T = mm.at("Evaluate 4-Index BraKet").run_as(braket); - - // Check output - auto t = test::eigen_buffer<4>(T.buffer()); - REQUIRE(test::trace(t) == - Catch::Approx(9.7919608941952063).margin(1.0e-16)); - REQUIRE(test::norm(t) == - Catch::Approx(7.7796143419802553).margin(1.0e-16)); + auto& mod = mm.at("Evaluate 4-Index BraKet"); + + SECTION("No UQ") { + // Call module + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<4>(T.buffer()); + REQUIRE(test::trace(t) == + Catch::Approx(9.7919608941952063).margin(1.0e-16)); + REQUIRE(test::norm(t) == + Catch::Approx(7.7796143419802553).margin(1.0e-16)); + } + + SECTION("With UQ") { + if constexpr(has_sigma) { + // Call module + mod.change_input("With UQ?", true); + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<4, udouble>(T.buffer()); + REQUIRE(test::trace(t).mean() == + Catch::Approx(9.7919608941952063).margin(1.0e-16)); + REQUIRE(test::norm(t).mean() == + Catch::Approx(7.7796143419802553).margin(1.0e-16)); + } + } } } From 910b1e4240386d040f50c845c2dd56eae8f8309e Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 7 Feb 2025 14:54:59 -0600 Subject: [PATCH 2/5] add sigma to CI --- .github/enable_sigma.cmake | 1 + .github/workflows/pull_request.yaml | 10 ++++++++++ 2 files changed, 11 insertions(+) create mode 100644 .github/enable_sigma.cmake diff --git a/.github/enable_sigma.cmake b/.github/enable_sigma.cmake new file mode 100644 index 00000000..2fbe8f6a --- /dev/null +++ b/.github/enable_sigma.cmake @@ -0,0 +1 @@ +set(ENABLE_SIGMA ON) \ No newline at end of file diff --git a/.github/workflows/pull_request.yaml b/.github/workflows/pull_request.yaml index f08bc7f1..f5da17e7 100644 --- a/.github/workflows/pull_request.yaml +++ b/.github/workflows/pull_request.yaml @@ -29,3 +29,13 @@ jobs: compilers: '["gcc-11", "clang-14"]' doc_target: 'integrals_cxx_api' secrets: inherit + Test-Enable-Sigma: + uses: NWChemEx/.github/.github/workflows/common_pull_request.yaml@master + with: + config_file: '' + source_dir: '' + format_python: false + compilers: '["gcc-11", "clang-14"]' + repo_toolchain: ".github/enable_sigma.cmake" + doc_target: '' + secrets: inherit From c0750b288f8e49234c229c54d155a0fbb2ebd751 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Fri, 7 Feb 2025 20:58:41 +0000 Subject: [PATCH 3/5] Committing clang-format changes --- src/integrals/ao_integrals/lib_int_visitor.hpp | 2 +- src/integrals/uncertain_types.hpp | 16 ++++++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/integrals/ao_integrals/lib_int_visitor.hpp b/src/integrals/ao_integrals/lib_int_visitor.hpp index d03778e5..1ddb2564 100644 --- a/src/integrals/ao_integrals/lib_int_visitor.hpp +++ b/src/integrals/ao_integrals/lib_int_visitor.hpp @@ -27,7 +27,7 @@ class LibIntVisitor : public chemist::qm_operator::OperatorVisitor { LibIntVisitor(const std::vector& bases, double thresh, std::size_t deriv = 0) : - m_bases(bases), m_thresh(thresh), m_deriv(deriv) {}; + m_bases(bases), m_thresh(thresh), m_deriv(deriv){}; void run(const s_e_type& S_e) { m_engine = detail_::make_engine(m_bases, S_e, m_thresh, m_deriv); diff --git a/src/integrals/uncertain_types.hpp b/src/integrals/uncertain_types.hpp index ca6b4d8c..aadca69b 100644 --- a/src/integrals/uncertain_types.hpp +++ b/src/integrals/uncertain_types.hpp @@ -1,3 +1,19 @@ +/* + * 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. + */ + #pragma once #ifdef ENABLE_SIGMA From 73686c2e06716915740fd50e4f58c765eb9a4d2c Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Sun, 9 Feb 2025 22:12:38 -0600 Subject: [PATCH 4/5] fix when sigma is disabled --- src/integrals/ao_integrals/ao_integrals.cpp | 8 +++---- src/integrals/uncertain_types.hpp | 12 +++++++--- .../ao_integrals/test_ao_integrals.hpp | 1 - .../ao_integrals/test_arbitrary_operator.cpp | 23 +++++++++++++------ 4 files changed, 29 insertions(+), 15 deletions(-) diff --git a/src/integrals/ao_integrals/ao_integrals.cpp b/src/integrals/ao_integrals/ao_integrals.cpp index 03b8ecb3..4e04e335 100644 --- a/src/integrals/ao_integrals/ao_integrals.cpp +++ b/src/integrals/ao_integrals/ao_integrals.cpp @@ -36,7 +36,7 @@ auto build_eigen_buffer(const std::vector& basis_sets, initial_value = FloatType(0.0, thresh); } Eigen::array dims_bfs; - for(auto i = 0; i < N; ++i) dims_bfs[i] = basis_sets[i].nbf(); + for(decltype(N) i = 0; i < N; ++i) dims_bfs[i] = basis_sets[i].nbf(); using shape_t = tensorwrapper::shape::Smooth; using layout_t = tensorwrapper::layout::Physical; @@ -56,7 +56,7 @@ auto fill_tensor(const std::vector& basis_sets, const chemist::qm_operator::OperatorBase& op, double thresh) { // Dimensional information std::vector dims_shells(N); - for(auto i = 0; i < N; ++i) dims_shells[i] = basis_sets[i].size(); + for(decltype(N) i = 0; i < N; ++i) dims_shells[i] = basis_sets[i].size(); auto b = build_eigen_buffer(basis_sets, thresh); @@ -83,7 +83,7 @@ auto fill_tensor(const std::vector& basis_sets, // Increment index shells[N - 1] += 1; - for(auto i = 1; i < N; ++i) { + for(decltype(N) i = 1; i < N; ++i) { if(shells[N - i] >= dims_shells[N - i]) { // Reset this dimension and increment the next one // shells[0] accumulates until we reach the end @@ -129,7 +129,7 @@ TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) { simde::type::tensor t; if(with_uq) { - if constexpr(!std::is_same_v) { + if constexpr(integrals::type::has_sigma()) { t = fill_tensor(basis_sets, op, thresh); } else { throw std::runtime_error("Sigma support not enabled!"); diff --git a/src/integrals/uncertain_types.hpp b/src/integrals/uncertain_types.hpp index aadca69b..f7160372 100644 --- a/src/integrals/uncertain_types.hpp +++ b/src/integrals/uncertain_types.hpp @@ -15,18 +15,24 @@ */ #pragma once - +#include #ifdef ENABLE_SIGMA #include #endif namespace integrals::type { #ifdef ENABLE_SIGMA + +static constexpr bool has_sigma() { return true; } + using uncertain_float = sigma::UFloat; using uncertain_double = sigma::UDouble; #else -using uncertain_float = void; -using uncertain_double = void; + +static constexpr bool has_sigma() { return false; } + +using uncertain_float = float; +using uncertain_double = double; #endif } // namespace integrals::type \ No newline at end of file diff --git a/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp b/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp index 902445aa..7a51fccc 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp @@ -39,5 +39,4 @@ auto norm(const tensorwrapper::buffer::Eigen& t) { t.value().square().sum().sqrt(); return norm.coeff(); } - } // namespace test \ No newline at end of file diff --git a/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp b/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp index 4fe28bf4..e692c7d6 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp @@ -18,7 +18,16 @@ #include "test_ao_integrals.hpp" using udouble = integrals::type::uncertain_double; -constexpr bool has_sigma = !std::is_same_v; +constexpr bool has_sigma = integrals::type::has_sigma(); + +template +auto unwrap_uq(InputType&& uq) { + if constexpr(has_sigma) { + return uq.mean(); + } else { + return uq; + } +} TEST_CASE("OperatorBase") { using aos_t = simde::type::aos; @@ -72,9 +81,9 @@ TEST_CASE("OperatorBase") { // Check output auto t = test::eigen_buffer<2, udouble>(T.buffer()); - REQUIRE(test::trace(t).mean() == + REQUIRE(unwrap_uq(test::trace(t)) == Catch::Approx(124.7011973877891364).margin(1.0e-16)); - REQUIRE(test::norm(t).mean() == + REQUIRE(unwrap_uq(test::norm(t)) == Catch::Approx(90.2562579028763707).margin(1.0e-16)); } } @@ -108,9 +117,9 @@ TEST_CASE("OperatorBase") { // Check output auto t = test::eigen_buffer<3, udouble>(T.buffer()); - REQUIRE(test::trace(t).mean() == + REQUIRE(unwrap_uq(test::trace(t)) == Catch::Approx(16.8245948391706577).margin(1.0e-16)); - REQUIRE(test::norm(t).mean() == + REQUIRE(unwrap_uq(test::norm(t)) == Catch::Approx(20.6560572032543597).margin(1.0e-16)); } } @@ -146,9 +155,9 @@ TEST_CASE("OperatorBase") { // Check output auto t = test::eigen_buffer<4, udouble>(T.buffer()); - REQUIRE(test::trace(t).mean() == + REQUIRE(unwrap_uq(test::trace(t)) == Catch::Approx(9.7919608941952063).margin(1.0e-16)); - REQUIRE(test::norm(t).mean() == + REQUIRE(unwrap_uq(test::norm(t)) == Catch::Approx(7.7796143419802553).margin(1.0e-16)); } } From 477dde7875d8d6af57239095cb7ac8c941aeebc1 Mon Sep 17 00:00:00 2001 From: "Jonathan M. Waldrop" Date: Mon, 10 Feb 2025 11:40:45 -0600 Subject: [PATCH 5/5] add standard deviation testing for uncertain values --- src/integrals/ao_integrals/ao_integrals.cpp | 2 +- .../ao_integrals/test_arbitrary_operator.cpp | 35 +++++++++++++++---- 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/src/integrals/ao_integrals/ao_integrals.cpp b/src/integrals/ao_integrals/ao_integrals.cpp index 4e04e335..8d66aaeb 100644 --- a/src/integrals/ao_integrals/ao_integrals.cpp +++ b/src/integrals/ao_integrals/ao_integrals.cpp @@ -77,7 +77,7 @@ auto fill_tensor(const std::vector& basis_sets, auto ord = detail_::shells2ord(basis_sets, shells); auto n_ord = ord.size(); for(decltype(n_ord) i_ord = 0; i_ord < n_ord; ++i_ord) { - b.value().data()[ord[i_ord]] = vals[i_ord]; + b.value().data()[ord[i_ord]] += vals[i_ord]; } } diff --git a/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp b/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp index e692c7d6..f20d312d 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp @@ -21,7 +21,7 @@ using udouble = integrals::type::uncertain_double; constexpr bool has_sigma = integrals::type::has_sigma(); template -auto unwrap_uq(InputType&& uq) { +auto unwrap_mean(InputType&& uq) { if constexpr(has_sigma) { return uq.mean(); } else { @@ -29,6 +29,15 @@ auto unwrap_uq(InputType&& uq) { } } +template +auto unwrap_sd(InputType&& uq) { + if constexpr(has_sigma) { + return uq.sd(); + } else { + return 0.0; + } +} + TEST_CASE("OperatorBase") { using aos_t = simde::type::aos; using aos_squared_t = simde::type::aos_squared; @@ -81,10 +90,14 @@ TEST_CASE("OperatorBase") { // Check output auto t = test::eigen_buffer<2, udouble>(T.buffer()); - REQUIRE(unwrap_uq(test::trace(t)) == + REQUIRE(unwrap_mean(test::trace(t)) == Catch::Approx(124.7011973877891364).margin(1.0e-16)); - REQUIRE(unwrap_uq(test::norm(t)) == + REQUIRE(unwrap_sd(test::trace(t)) == + Catch::Approx(7e-16).margin(1.0e-16)); + REQUIRE(unwrap_mean(test::norm(t)) == Catch::Approx(90.2562579028763707).margin(1.0e-16)); + REQUIRE(unwrap_sd(test::norm(t)) == + Catch::Approx(3e-16).margin(1.0e-16)); } } } @@ -117,10 +130,14 @@ TEST_CASE("OperatorBase") { // Check output auto t = test::eigen_buffer<3, udouble>(T.buffer()); - REQUIRE(unwrap_uq(test::trace(t)) == + REQUIRE(unwrap_mean(test::trace(t)) == Catch::Approx(16.8245948391706577).margin(1.0e-16)); - REQUIRE(unwrap_uq(test::norm(t)) == + REQUIRE(unwrap_sd(test::trace(t)) == + Catch::Approx(7e-16).margin(1.0e-16)); + REQUIRE(unwrap_mean(test::norm(t)) == Catch::Approx(20.6560572032543597).margin(1.0e-16)); + REQUIRE(unwrap_sd(test::norm(t)) == + Catch::Approx(7e-16).margin(1.0e-16)); } } } @@ -155,10 +172,14 @@ TEST_CASE("OperatorBase") { // Check output auto t = test::eigen_buffer<4, udouble>(T.buffer()); - REQUIRE(unwrap_uq(test::trace(t)) == + REQUIRE(unwrap_mean(test::trace(t)) == Catch::Approx(9.7919608941952063).margin(1.0e-16)); - REQUIRE(unwrap_uq(test::norm(t)) == + REQUIRE(unwrap_sd(test::trace(t)) == + Catch::Approx(7e-16).margin(1.0e-16)); + REQUIRE(unwrap_mean(test::norm(t)) == Catch::Approx(7.7796143419802553).margin(1.0e-16)); + REQUIRE(unwrap_sd(test::norm(t)) == + Catch::Approx(11e-16).margin(1.0e-16)); } } }