Skip to content

Add different means for QCUP #157

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ryanmrichard merged 2 commits into from
Mar 4, 2026
+210 −12
Merged

Add different means for QCUP #157

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8bae897
adds test for using max error
ryanmrichard Mar 3, 2026
1d681aa
adds error reduction methods
ryanmrichard Mar 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 72 additions & 5 deletions src/integrals/ao_integrals/uq_driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<typename FloatType>
FloatType geometric_mean(const std::span<FloatType> values) {
std::decay_t<FloatType> 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<typename FloatType>
FloatType harmonic_mean(const std::span<FloatType> values) {
std::decay_t<FloatType> 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<typename FloatType0, typename FloatType1>
Tensor operator()(const std::span<FloatType0> t,
Expand All @@ -49,10 +94,24 @@ struct Kernel {
using uq_type = sigma::Uncertain<float_type>;
auto rv_buffer = buffer::make_contiguous<uq_type>(m_shape);
auto rv_data = buffer::get_raw_data<uq_type>(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
Expand All @@ -63,6 +122,7 @@ struct Kernel {
return rv;
}
shape_type m_shape;
mean_type m_mean;
};

const auto desc = R"(
Expand All @@ -81,10 +141,17 @@ MODULE_CTOR(UQDriver) {
description(desc);
add_submodule<eri_pt>("ERIs");
add_submodule<error_pt>("ERI Error");
add_input<bool>("Max Error").set_default(false);
add_input<bool>("Geometric Mean").set_default(false);
add_input<bool>("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>();
bool use_geom = inputs.at("Geometric Mean").value<bool>();
bool use_harmonic = inputs.at("Harmonic Mean").value<bool>();
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<double>();
Expand All @@ -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);
Expand Down
145 changes: 138 additions & 7 deletions tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<typename FloatType>
auto corr_answer(const simde::type::tensor& T) {
if constexpr(std::is_same_v<FloatType, double>) {
Expand Down Expand Up @@ -59,11 +66,97 @@ auto corr_answer(const simde::type::tensor& T) {
}
}

template<typename FloatType>
auto corr_max_answer(const simde::type::tensor& T) {
if constexpr(std::is_same_v<FloatType, double>) {
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<typename FloatType>
auto corr_geometric_mean_answer(const simde::type::tensor& T) {
if constexpr(std::is_same_v<FloatType, double>) {
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<typename FloatType>
auto corr_harmonic_mean_answer(const simde::type::tensor& T) {
if constexpr(std::is_same_v<FloatType, double>) {
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<float_type>) {
using test_pt = simde::ERI4;
using test_pt = simde::ERI4;

if constexpr(tensorwrapper::types::is_uncertain_v<float_type>) {
auto rt = std::make_unique<parallelzone::runtime::RuntimeView>();
pluginplay::ModuleManager mm(std::move(rt), nullptr);
integrals::load_modules(mm);
Expand All @@ -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<test_pt>(braket);

auto T_corr = corr_answer<float_type>(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<test_pt>(braket);

auto T_corr = corr_answer<float_type>(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<test_pt>(braket);

auto T_corr = corr_max_answer<float_type>(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<test_pt>(braket);

auto T_corr = corr_geometric_mean_answer<float_type>(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<test_pt>(braket);

auto T_corr = corr_harmonic_mean_answer<float_type>(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<test_pt>(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<test_pt>(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<test_pt>(braket));
}
}
}