diff --git a/src/integrals/ao_integrals/uq_driver.cpp b/src/integrals/ao_integrals/uq_driver.cpp index 43ee1f87..081b73c1 100644 --- a/src/integrals/ao_integrals/uq_driver.cpp +++ b/src/integrals/ao_integrals/uq_driver.cpp @@ -25,9 +25,54 @@ using namespace tensorwrapper; namespace integrals::ao_integrals { namespace { +enum mean_type { none, max, geometric, harmonic }; + +mean_type get_mean_type(bool use_max, bool use_geom, bool use_harmonic) { + mean_type mean = none; + auto throw_error = []() { + throw std::runtime_error( + "Cannot use multiple mean types at once. Please select only one."); + }; + if(use_max) { mean = max; } + if(use_geom) { + if(mean != none) throw_error(); + mean = geometric; + } + if(use_harmonic) { + if(mean != none) throw_error(); + mean = harmonic; + } + return mean; +} + +template +FloatType geometric_mean(const std::span values) { + std::decay_t 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 +FloatType harmonic_mean(const std::span values) { + std::decay_t 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; +} + struct Kernel { using shape_type = buffer::Contiguous::shape_type; - Kernel(shape_type shape) : m_shape(std::move(shape)) {} + Kernel(shape_type shape, mean_type mean) : + m_shape(std::move(shape)), m_mean(mean) {} template Tensor operator()(const std::span t, @@ -49,10 +94,24 @@ struct Kernel { using uq_type = sigma::Uncertain; auto rv_buffer = buffer::make_contiguous(m_shape); auto rv_data = buffer::get_raw_data(rv_buffer); + + float_type max_error = 0.0; + if(m_mean == max) { + max_error = *std::max_element(error.begin(), error.end()); + } else if(m_mean == geometric) { + max_error = geometric_mean(error); + } else if(m_mean == harmonic) { + max_error = harmonic_mean(error); + } + uq_type max_uq{0.0, max_error}; + bool use_mean = (m_mean != none); for(std::size_t i = 0; i < t.size(); ++i) { - const auto elem = t[i]; - const auto elem_error = error[i]; - rv_data[i] = uq_type(elem, elem_error); + const auto elem = t[i]; + if(!use_mean) { + rv_data[i] = uq_type{elem, error[i]}; + } else { + rv_data[i] = elem + max_uq; + } } rv = tensorwrapper::Tensor(m_shape, std::move(rv_buffer)); #else @@ -63,6 +122,7 @@ struct Kernel { return rv; } shape_type m_shape; + mean_type m_mean; }; const auto desc = R"( @@ -81,10 +141,17 @@ MODULE_CTOR(UQDriver) { description(desc); add_submodule("ERIs"); add_submodule("ERI Error"); + add_input("Max Error").set_default(false); + add_input("Geometric Mean").set_default(false); + add_input("Harmonic Mean").set_default(false); } MODULE_RUN(UQDriver) { const auto& [braket] = eri_pt::unwrap_inputs(inputs); + bool use_max = inputs.at("Max Error").value(); + bool use_geom = inputs.at("Geometric Mean").value(); + bool use_harmonic = inputs.at("Harmonic Mean").value(); + auto mean = get_mean_type(use_max, use_geom, use_harmonic); auto& eri_mod = submods.at("ERIs").value(); auto tol = eri_mod.inputs().at("Threshold").value(); @@ -94,7 +161,7 @@ MODULE_RUN(UQDriver) { using buffer::visit_contiguous_buffer; shape::Smooth shape = t.buffer().layout().shape().as_smooth().make_smooth(); - Kernel k(shape); + Kernel k(shape, mean); auto t_buffer = make_contiguous(t.buffer()); auto error_buffer = make_contiguous(error.buffer()); auto t_w_error = visit_contiguous_buffer(k, t_buffer, error_buffer); diff --git a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp index 22ec0ff9..d19c678a 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp @@ -20,6 +20,13 @@ 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) { @@ -59,11 +66,97 @@ auto corr_answer(const simde::type::tensor& T) { } } +template +auto corr_max_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()); + double max_error = 0.0000170000000000; + corr_buffer.set_elem({0, 0, 0, 0}, FloatType{0.774606, max_error}); + corr_buffer.set_elem({0, 0, 0, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({0, 0, 1, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({0, 0, 1, 1}, FloatType{0.446701, max_error}); + corr_buffer.set_elem({0, 1, 0, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({0, 1, 0, 1}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({0, 1, 1, 0}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({0, 1, 1, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 0, 0, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 0, 0, 1}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({1, 0, 1, 0}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({1, 0, 1, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 1, 0, 0}, FloatType{0.446701, max_error}); + corr_buffer.set_elem({1, 1, 0, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 1, 1, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 1, 1, 1}, FloatType{0.774606, max_error}); + return T_corr; + } +} + +template +auto corr_geometric_mean_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()); + double max_error = 0.0000025712815907; + corr_buffer.set_elem({0, 0, 0, 0}, FloatType{0.774606, max_error}); + corr_buffer.set_elem({0, 0, 0, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({0, 0, 1, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({0, 0, 1, 1}, FloatType{0.446701, max_error}); + corr_buffer.set_elem({0, 1, 0, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({0, 1, 0, 1}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({0, 1, 1, 0}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({0, 1, 1, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 0, 0, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 0, 0, 1}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({1, 0, 1, 0}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({1, 0, 1, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 1, 0, 0}, FloatType{0.446701, max_error}); + corr_buffer.set_elem({1, 1, 0, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 1, 1, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 1, 1, 1}, FloatType{0.774606, max_error}); + return T_corr; + } +} + +template +auto corr_harmonic_mean_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()); + double max_error = 0.0000014571428571; + corr_buffer.set_elem({0, 0, 0, 0}, FloatType{0.774606, max_error}); + corr_buffer.set_elem({0, 0, 0, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({0, 0, 1, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({0, 0, 1, 1}, FloatType{0.446701, max_error}); + corr_buffer.set_elem({0, 1, 0, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({0, 1, 0, 1}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({0, 1, 1, 0}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({0, 1, 1, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 0, 0, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 0, 0, 1}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({1, 0, 1, 0}, FloatType{0.120666, max_error}); + corr_buffer.set_elem({1, 0, 1, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 1, 0, 0}, FloatType{0.446701, max_error}); + corr_buffer.set_elem({1, 1, 0, 1}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 1, 1, 0}, FloatType{0.265558, max_error}); + corr_buffer.set_elem({1, 1, 1, 1}, FloatType{0.774606, max_error}); + return T_corr; + } +} + +} // namespace + TEST_CASE("UQ Driver") { using float_type = tensorwrapper::types::udouble; - if constexpr(tensorwrapper::types::is_uncertain_v) { - using test_pt = simde::ERI4; + 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); @@ -85,11 +178,49 @@ TEST_CASE("UQ Driver") { // Make BraKet Input chemist::braket::BraKet braket(aos_squared, op, aos_squared); - // Call modules - auto T = mm.at("UQ Driver").run_as(braket); - - auto T_corr = corr_answer(T); using tensorwrapper::operations::approximately_equal; - REQUIRE(approximately_equal(T_corr, T, 1E-6)); + SECTION("No Mean") { + // Call modules + auto T = mm.at("UQ Driver").run_as(braket); + + auto T_corr = corr_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + SECTION("Max Error") { + mm.change_input("UQ Driver", "Max Error", true); + auto T = mm.at("UQ Driver").run_as(braket); + + auto T_corr = corr_max_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + SECTION("Geometric Mean") { + mm.change_input("UQ Driver", "Geometric Mean", true); + auto T = mm.at("UQ Driver").run_as(braket); + + auto T_corr = corr_geometric_mean_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + SECTION("Harmonic Mean") { + mm.change_input("UQ Driver", "Harmonic Mean", true); + auto T = mm.at("UQ Driver").run_as(braket); + + auto T_corr = corr_harmonic_mean_answer(T); + REQUIRE(approximately_equal(T_corr, T, 1E-6)); + } + SECTION("Invalid (T, T, F)") { + mm.change_input("UQ Driver", "Max Error", true); + mm.change_input("UQ Driver", "Geometric Mean", true); + REQUIRE_THROWS(mm.at("UQ Driver").run_as(braket)); + } + SECTION("Invalid (T, F, T)") { + mm.change_input("UQ Driver", "Max Error", true); + mm.change_input("UQ Driver", "Harmonic Mean", true); + REQUIRE_THROWS(mm.at("UQ Driver").run_as(braket)); + } + SECTION("Invalid (F, T, T)") { + mm.change_input("UQ Driver", "Geometric Mean", true); + mm.change_input("UQ Driver", "Harmonic Mean", true); + REQUIRE_THROWS(mm.at("UQ Driver").run_as(braket)); + } } }