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 2384b1d3..6f4640c4 100644 --- a/src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp +++ b/src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp @@ -29,17 +29,37 @@ namespace { template auto average_error(T&& strides, T&& nbf, T&& ao_i, Tensor&& error, utils::mean_type mean) { - std::vector buffer; +#ifdef ENABLE_SIGMA + using uq_type = sigma::Uncertain; + auto n_elements = nbf[0] * nbf[1] * nbf[2] * nbf[3]; + + if(mean == utils::mean_type::none) { + std::vector result; + result.reserve(n_elements); + 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]; + result.push_back(uq_type{0.0, error[loffset]}); + } + } + } + } + return result; + } + std::vector buffer; + buffer.reserve(n_elements); 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]); @@ -47,15 +67,20 @@ auto average_error(T&& strides, T&& nbf, T&& ao_i, Tensor&& error, } } } - - return utils::compute_mean(mean, buffer); + auto mean_value = utils::compute_mean(mean, buffer); + return std::vector(n_elements, uq_type{0.0, mean_value}); +#else + throw std::runtime_error("Sigma support not enabled!"); + return std::vector{}; +#endif } +#ifdef ENABLE_SIGMA template auto compute_block(T&& strides, T&& nbf, T&& ao_i, Tensor&& value, - FloatType error) { + const std::vector>& errors) { auto n_elements = nbf[0] * nbf[1] * nbf[2] * nbf[3]; - std::vector> buffer(n_elements, error); + std::vector> buffer(n_elements); for(std::size_t i = 0; i < nbf[0]; ++i) { auto ilocal = i * nbf[1] * nbf[2] * nbf[3]; @@ -70,15 +95,16 @@ auto compute_block(T&& strides, T&& nbf, T&& ao_i, Tensor&& value, 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]; + auto llocal = klocal + l; + auto loffset = koffset + (ao_i[3] + l) * strides[3]; + buffer[llocal] = errors[llocal] + value[loffset]; } } } } return buffer; } +#endif template void set_block(T&& strides, T&& nbf, @@ -191,13 +217,12 @@ struct Kernel { bool pair_gt = (c2eqc0 && centers[3] > centers[1]); if(pair_gt && all_same) break; - auto block_error = average_error( + auto block_errors = 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); + t, block_errors); // Set all symmetry equivalent blocks to `block` auto perms = get_permutations_with_sigma( 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 index e4c777b6..4c8bbafe 100644 --- 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 @@ -66,11 +66,36 @@ auto corr_answer(const simde::type::tensor& T) { } } +template +auto corr_answer_no_mean(const simde::type::tensor& T, + const simde::type::tensor& T_eri, + const simde::type::tensor& T_err) { + if constexpr(std::is_same_v) { + return true; + } else { + auto t_uq = eigen_tensor<4, FloatType>(T.buffer()); + auto t_eri = eigen_tensor<4, double>(T_eri.buffer()); + auto t_err = eigen_tensor<4, double>(T_err.buffer()); + + for(Eigen::Index i = 0; i < 2; ++i) + for(Eigen::Index j = 0; j < 2; ++j) + for(Eigen::Index k = 0; k < 2; ++k) + for(Eigen::Index l = 0; l < 2; ++l) { + REQUIRE(t_uq(i, j, k, l).mean() == + Catch::Approx(t_eri(i, j, k, l)).margin(1E-6)); + REQUIRE(t_uq(i, j, k, l).sd() == + Catch::Approx(t_err(i, j, k, l)).margin(1E-6)); + } + return true; + } +} + } // namespace TEST_CASE("UQ Atom Symm Blocked Driver") { using float_type = tensorwrapper::types::udouble; using test_pt = simde::ERI4; + using error_pt = integrals::property_types::Uncertainty; using tensorwrapper::operations::approximately_equal; if constexpr(tensorwrapper::types::is_uncertain_v) { @@ -80,7 +105,7 @@ TEST_CASE("UQ Atom Symm Blocked Driver") { 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"); + auto mod = mm.at("UQ Atom Symm Blocked Driver"); // Make Operator simde::type::v_ee_type op{}; @@ -96,7 +121,14 @@ TEST_CASE("UQ Atom Symm Blocked Driver") { // Make BraKet Input chemist::braket::BraKet braket(aos_squared, op, aos_squared); - SECTION("No Mean") { REQUIRE_THROWS(mod.run_as(braket)); } + SECTION("No Mean") { + auto T = mod.run_as(braket); + auto tol = 1.0e-6; + auto T_eri = mm.at("ERI4").run_as(braket); + auto T_err = + mm.at("Primitive Error Model").run_as(braket, tol); + REQUIRE(corr_answer_no_mean(T, T_eri, T_err)); + } SECTION("Max Error") { mod.change_input("Mean Type", "max"); auto T = mod.run_as(braket);