Skip to content

Affine #168

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 7 commits into from
Jun 16, 2026
Merged

Affine #168

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
909dd65
works with affine forms
ryanmrichard Jun 10, 2026
db4e08a
integrals works!
ryanmrichard Jun 12, 2026
879116c
changes for affine
ryanmrichard Jun 16, 2026
b187ee5
clean up unit tests
ryanmrichard Jun 16, 2026
3142894
remove unused typedef
ryanmrichard Jun 16, 2026
0d20ef5
change static asserts to exceptions to satisfy compilers
ryanmrichard Jun 16, 2026
6a940a9
Fix incorrect logic for when Sigma is not enabled.
ryanmrichard Jun 16, 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
3 changes: 1 addition & 2 deletions src/integrals/ao_integrals/coulomb_metric.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,7 @@ struct Kernel {
template<typename FloatType>
simde::type::tensor operator()(const std::span<FloatType> I) {
using clean_type = std::decay_t<FloatType>;
if constexpr(types::is_uncertain_v<clean_type> ||
types::is_interval_v<clean_type>) {
if constexpr(types::is_uq_type_v<clean_type>) {
throw std::runtime_error(
"integrals::ao_integrals::CoulombMetric: UQ types not supported");
} else {
Expand Down
37 changes: 14 additions & 23 deletions src/integrals/ao_integrals/uq_atom_symm_blocked_driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,11 +33,11 @@ auto average_error(T&& strides, T&& nbf, T&& ao_i, Tensor&& error,
"integrals::ao_integrals::UQAtomSymmBlockedDriver: ";

#ifdef ENABLE_SIGMA
using float_type = typename UQType::value_t;
using demangler_type = ::utilities::printing::Demangler;
using float_type = typename UQType::value_t;

auto n_elements = nbf[0] * nbf[1] * nbf[2] * nbf[3];

using tensorwrapper::types::construct_uq_type;
if(mean == utils::mean_type::none) {
std::vector<UQType> result;
result.reserve(n_elements);
Expand All @@ -50,16 +50,8 @@ auto average_error(T&& strides, T&& nbf, T&& ao_i, Tensor&& error,
for(std::size_t l = 0; l < nbf[3]; ++l) {
auto loffset = koffset + (ao_i[3] + l) * strides[3];
auto ei = std::fabs(error[loffset]);
if constexpr(types::is_interval_v<UQType>) {
result.push_back(UQType{-ei, ei});
} else if constexpr(types::is_uncertain_v<UQType>) {
result.push_back(UQType{0.0, ei});
} else {
auto type0 = demangler_type::demangle<UQType>();
throw std::runtime_error(
error_base + "Invalid UQ type " + type0 +
" for non-mean reduction");
}
auto elem = construct_uq_type<UQType>(0.0, ei);
result.push_back(elem);
}
}
}
Expand All @@ -83,15 +75,8 @@ auto average_error(T&& strides, T&& nbf, T&& ao_i, Tensor&& error,
}
}
auto mean_value = utils::compute_mean(mean, buffer);
if constexpr(types::is_interval_v<UQType>) {
return std::vector<UQType>(n_elements, UQType{-mean_value, mean_value});
} else if constexpr(types::is_uncertain_v<UQType>) {
return std::vector<UQType>(n_elements, UQType{0.0, mean_value});
} else {
auto type0 = demangler_type::demangle<UQType>();
throw std::runtime_error(error_base + "Invalid UQ type " + type0 +
" for mean reduction");
}
auto value = construct_uq_type<UQType>(0.0, mean_value);
return std::vector<UQType>(n_elements, value);
#else
throw std::runtime_error(error_base + "Sigma support not enabled!");
return std::vector<UQType>{};
Expand Down Expand Up @@ -191,8 +176,7 @@ struct Kernel {
Tensor rv;

using float_type = std::decay_t<FloatType>;
if constexpr(types::is_uncertain_v<float_type> ||
types::is_interval_v<float_type>) {
if constexpr(types::is_uq_type_v<float_type>) {
auto type0 = demangler_type::demangle<FloatType>();
throw std::runtime_error(
m_error_base + "Expected non-UQ floating point type, but got " +
Expand Down Expand Up @@ -348,6 +332,13 @@ MODULE_RUN(UQAtomSymmBlockedDriver) {
} else if(uq_type == "interval") {
Kernel<tensorwrapper::types::interval_type> k(shape, aos, mean);
t_w_error = visit_contiguous_buffer(k, t_buffer, e_buffer);
} else if(uq_type == "affine") {
Kernel<tensorwrapper::types::affine_type> k(shape, aos, mean);
t_w_error = visit_contiguous_buffer(k, t_buffer, e_buffer);
} else if(uq_type == "thresholded affine") {
Kernel<tensorwrapper::types::thresholded_affine_type> k(shape, aos,
mean);
t_w_error = visit_contiguous_buffer(k, t_buffer, e_buffer);
} else {
throw std::runtime_error(
"integrals::ao_integrals::UQAtomSymmBlockedDriver: Invalid UQ type "
Expand Down
38 changes: 14 additions & 24 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 @@ -52,46 +52,30 @@ struct Kernel {
Tensor rv;

using float_type = std::decay_t<FloatType>;
if constexpr(types::is_uncertain_v<float_type> ||
types::is_interval_v<float_type>) {
if constexpr(types::is_uq_type_v<float_type>) {
auto type0 = demangler_type::demangle<FloatType>();
throw std::runtime_error(
m_error_base + "Expected non-UQ floating point type, but got " +
type0);
} else {
#ifdef ENABLE_SIGMA
using uq_type = UQType<float_type>;
using uq_type = UQType<float_type>;
using tensorwrapper::types::construct_uq_type;
using tensorwrapper::types::fabs;
auto rv_buffer = buffer::make_contiguous<uq_type>(m_shape);
auto rv_data = buffer::get_raw_data<uq_type>(rv_buffer);

bool use_mean = (m_mean != utils::mean_type::none);
if(!use_mean) {
for(std::size_t i = 0; i < t.size(); ++i) {
const auto elem = t[i];
if constexpr(types::is_interval_v<uq_type>) {
auto ei = std::fabs(error[i]);
rv_data[i] = uq_type{elem - ei, elem + ei};
} else if constexpr(types::is_uncertain_v<uq_type>) {
rv_data[i] = uq_type{elem, error[i]};
} else {
auto type0 = demangler_type::demangle<uq_type>();
throw std::runtime_error(m_error_base +
"Invalid UQ type " + type0 +
" for non-mean reduction");
}
const auto ei = std::fabs(error[i]);
rv_data[i] = construct_uq_type<uq_type>(elem, ei);
}
} else {
float_type max_error = utils::compute_mean(m_mean, error);
uq_type max_uq;
if constexpr(types::is_interval_v<uq_type>) {
max_uq = uq_type(-max_error, max_error);
} else if constexpr(types::is_uncertain_v<uq_type>) {
max_uq = uq_type(0.0, max_error);
} else {
auto type0 = demangler_type::demangle<uq_type>();
throw std::runtime_error(m_error_base + "Invalid UQ type " +
type0 + " for mean reduction");
}
auto max_uq = construct_uq_type<uq_type>(0.0, max_error);

for(std::size_t i = 0; i < t.size(); ++i) {
const auto elem = t[i];
rv_data[i] = uq_type(elem) + max_uq;
Expand Down Expand Up @@ -155,6 +139,12 @@ MODULE_RUN(UQDriver) {
} else if(uq_type == "interval") {
Kernel<tensorwrapper::types::interval_type> k(shape, mean);
t_w_error = visit_contiguous_buffer(k, t_buffer, error_buffer);
} else if(uq_type == "affine") {
Kernel<tensorwrapper::types::affine_type> k(shape, mean);
t_w_error = visit_contiguous_buffer(k, t_buffer, error_buffer);
} else if(uq_type == "thresholded affine") {
Kernel<tensorwrapper::types::thresholded_affine_type> k(shape, mean);
t_w_error = visit_contiguous_buffer(k, t_buffer, error_buffer);
} else {
throw std::runtime_error(
"integrals::ao_integrals::UQDriver: Invalid UQ type name " + uq_type);
Expand Down
8 changes: 8 additions & 0 deletions src/integrals/libint/libint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,14 @@ TEMPLATED_MODULE_RUN(Libint, BraKetType) {
} else if(uq_type == "interval") {
using interval_type = tensorwrapper::types::interval_type<float_type>;
t = detail_::fill_tensor<N, interval_type>(basis_sets, op, rv, thresh);
} else if(uq_type == "affine") {
using affine_type = tensorwrapper::types::affine_type<float_type>;
t = detail_::fill_tensor<N, affine_type>(basis_sets, op, rv, thresh);
} else if(uq_type == "thresholded affine") {
using thresholded_affine_type =
tensorwrapper::types::thresholded_affine_type<float_type>;
t = detail_::fill_tensor<N, thresholded_affine_type>(basis_sets, op, rv,
thresh);
} else {
throw std::runtime_error(
"integrals::libint::Libint: Invalid UQ type name " + uq_type);
Expand Down
4 changes: 2 additions & 2 deletions src/integrals/utils/uncertainty_reductions.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@ auto max_mean(const ContainerType& values) {

template<typename ContainerType>
auto geometric_mean(const ContainerType& values) {
using float_type = std::decay_t<decltype(values[0])>;
float_type log_sum = 0.0;
using float_type = std::decay_t<decltype(values[0])>;
float_type log_sum(0.0);
std::size_t non_zero_count = 0;
for(const auto& val : values) {
if(val == float_type{0}) continue;
Expand Down
143 changes: 62 additions & 81 deletions tests/cxx/unit/integrals/ao_integrals/test_uq_atom_symm_blocked_driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,20 +22,6 @@ using namespace tensorwrapper;

namespace {

template<typename UQType, typename FloatType>
auto elem(FloatType value, FloatType error) {
if constexpr(tensorwrapper::types::is_uncertain_v<UQType>) {
return UQType{value, error};
} else if constexpr(tensorwrapper::types::is_interval_v<UQType>) {
return UQType{value - error, value + error};
} else {
::utilities::printing::Demangler demangler;
auto type0 = demangler.template demangle<UQType>();
throw std::runtime_error(
"UQ Atom Symm Blocked Driver Test: Invalid UQ type " + type0);
}
}

// 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"
Expand All @@ -49,34 +35,38 @@ auto corr_answer(const simde::type::tensor& T) {
using value_t = typename UQType::value_t;
simde::type::tensor T_corr(T);
auto& corr_buffer = buffer::make_contiguous(T_corr.buffer());
corr_buffer.set_elem({0, 0, 0, 0}, elem<UQType>(0.774606, value_t(0)));
corr_buffer.set_elem({0, 0, 0, 1},
elem<UQType>(0.265558, 0.0000010000000000));
corr_buffer.set_elem({0, 0, 1, 0},
elem<UQType>(0.265558, 0.0000010000000000));
corr_buffer.set_elem({0, 0, 1, 1}, elem<UQType>(0.446701, value_t(0)));
corr_buffer.set_elem({0, 1, 0, 0},
elem<UQType>(0.265558, 0.0000010000000000));
corr_buffer.set_elem({0, 1, 0, 1},
elem<UQType>(0.120666, 0.0000170000000000));
corr_buffer.set_elem({0, 1, 1, 0},
elem<UQType>(0.120666, 0.0000170000000000));
corr_buffer.set_elem({0, 1, 1, 1},
elem<UQType>(0.265558, 0.0000010000000000));
corr_buffer.set_elem({1, 0, 0, 0},
elem<UQType>(0.265558, 0.0000010000000000));
corr_buffer.set_elem({1, 0, 0, 1},
elem<UQType>(0.120666, 0.0000170000000000));
corr_buffer.set_elem({1, 0, 1, 0},
elem<UQType>(0.120666, 0.0000170000000000));
corr_buffer.set_elem({1, 0, 1, 1},
elem<UQType>(0.265558, 0.0000010000000000));
corr_buffer.set_elem({1, 1, 0, 0}, elem<UQType>(0.446701, value_t(0)));
corr_buffer.set_elem({1, 1, 0, 1},
elem<UQType>(0.265558, 0.0000010000000000));
corr_buffer.set_elem({1, 1, 1, 0},
elem<UQType>(0.265558, 0.0000010000000000));
corr_buffer.set_elem({1, 1, 1, 1}, elem<UQType>(0.774606, value_t(0)));
corr_buffer.set_elem(
{0, 0, 0, 0}, types::construct_uq_type<UQType>(0.774606, value_t(0)));
corr_buffer.set_elem({0, 0, 0, 1}, types::construct_uq_type<UQType>(
0.265558, 0.0000010000000000));
corr_buffer.set_elem({0, 0, 1, 0}, types::construct_uq_type<UQType>(
0.265558, 0.0000010000000000));
corr_buffer.set_elem(
{0, 0, 1, 1}, types::construct_uq_type<UQType>(0.446701, value_t(0)));
corr_buffer.set_elem({0, 1, 0, 0}, types::construct_uq_type<UQType>(
0.265558, 0.0000010000000000));
corr_buffer.set_elem({0, 1, 0, 1}, types::construct_uq_type<UQType>(
0.120666, 0.0000170000000000));
corr_buffer.set_elem({0, 1, 1, 0}, types::construct_uq_type<UQType>(
0.120666, 0.0000170000000000));
corr_buffer.set_elem({0, 1, 1, 1}, types::construct_uq_type<UQType>(
0.265558, 0.0000010000000000));
corr_buffer.set_elem({1, 0, 0, 0}, types::construct_uq_type<UQType>(
0.265558, 0.0000010000000000));
corr_buffer.set_elem({1, 0, 0, 1}, types::construct_uq_type<UQType>(
0.120666, 0.0000170000000000));
corr_buffer.set_elem({1, 0, 1, 0}, types::construct_uq_type<UQType>(
0.120666, 0.0000170000000000));
corr_buffer.set_elem({1, 0, 1, 1}, types::construct_uq_type<UQType>(
0.265558, 0.0000010000000000));
corr_buffer.set_elem(
{1, 1, 0, 0}, types::construct_uq_type<UQType>(0.446701, value_t(0)));
corr_buffer.set_elem({1, 1, 0, 1}, types::construct_uq_type<UQType>(
0.265558, 0.0000010000000000));
corr_buffer.set_elem({1, 1, 1, 0}, types::construct_uq_type<UQType>(
0.265558, 0.0000010000000000));
corr_buffer.set_elem(
{1, 1, 1, 1}, types::construct_uq_type<UQType>(0.774606, value_t(0)));
return T_corr;
}
}
Expand Down Expand Up @@ -112,11 +102,14 @@ auto corr_answer_no_mean(const simde::type::tensor& T,

} // namespace

TEST_CASE("UQ Atom Symm Blocked Driver") {
using uncertain_type = tensorwrapper::types::udouble;
using interval_type = tensorwrapper::types::interval_type<double>;
using test_pt = simde::ERI4;
using error_pt = integrals::property_types::Uncertainty<test_pt>;
using uq_types =
std::tuple<tensorwrapper::types::interval_type<double>,
tensorwrapper::types::affine_type<double>,
tensorwrapper::types::thresholded_affine_type<double>>;
TEMPLATE_LIST_TEST_CASE("UQ Atom Symm Blocked Driver", "", uq_types) {
using float_type = TestType;
using test_pt = simde::ERI4;
using error_pt = integrals::property_types::Uncertainty<test_pt>;
using tensorwrapper::operations::approximately_equal;

auto rt = std::make_unique<parallelzone::runtime::RuntimeView>();
Expand All @@ -141,60 +134,48 @@ TEST_CASE("UQ Atom Symm Blocked Driver") {
// Make BraKet Input
chemist::braket::BraKet braket(aos_squared, op, aos_squared);

if constexpr(tensorwrapper::types::is_uncertain_v<uncertain_type>) {
SECTION("No Mean") {
auto T = mod.run_as<test_pt>(braket);
auto tol = 1.0e-6;
auto T_eri = mm.at("ERI4").run_as<test_pt>(braket);
auto T_err =
mm.at("Primitive Error Model").run_as<error_pt>(braket, tol);
REQUIRE(corr_answer_no_mean<uncertain_type>(T, T_eri, T_err));
}
SECTION("Max Error") {
mod.change_input("Mean Type", "max");
auto T = mod.run_as<test_pt>(braket);

auto T_corr = corr_answer<uncertain_type>(T);
REQUIRE(approximately_equal(T_corr, T, 1E-6));
std::string uq_type = []() -> std::string {
if constexpr(tensorwrapper::types::is_uncertain_v<float_type>) {
return "uncertain";
} else if constexpr(tensorwrapper::types::is_interval_v<
float_type>) {
return "interval";
} else if constexpr(tensorwrapper::types::is_affine_v<float_type>) {
return "affine";
} else if constexpr(tensorwrapper::types::is_thresholded_affine_v<
float_type>) {
return "thresholded affine";
} else {
return ""; // Won't be used because of static assert below.
}
SECTION("Geometric Mean") {
mod.change_input("Mean Type", "geometric");
auto T = mod.run_as<test_pt>(braket);
}();

auto T_corr = corr_answer<uncertain_type>(T);
REQUIRE(approximately_equal(T_corr, T, 1E-6));
}
}
if constexpr(tensorwrapper::types::is_interval_v<interval_type>) {
if constexpr(tensorwrapper::types::is_uq_type_v<float_type>) {
// The errors for the integrals are on the order of 1e-5.
// Subtracting
// the intervals in approximately_equal treats the intervals as
// independent and doubles them instead of cancelling them. This in
// turn means center + radius on the difference inside
// approximately_equal will be between 1e-5 and 1e-4.
// Subtracting results in differences between 1e-5 and 1e-4.
SECTION("No Mean") {
mod.change_input("UQ Type", "interval");
mod.change_input("UQ Type", uq_type);
auto T = mod.run_as<test_pt>(braket);
auto tol = 1.0e-6;
auto T_eri = mm.at("ERI4").run_as<test_pt>(braket);
auto T_err =
mm.at("Primitive Error Model").run_as<error_pt>(braket, tol);
REQUIRE(corr_answer_no_mean<interval_type>(T, T_eri, T_err));
REQUIRE(corr_answer_no_mean<float_type>(T, T_eri, T_err));
}
SECTION("Max Error") {
mod.change_input("UQ Type", "interval");
mod.change_input("UQ Type", uq_type);
mod.change_input("Mean Type", "max");
auto T = mod.run_as<test_pt>(braket);

auto T_corr = corr_answer<interval_type>(T);
auto T_corr = corr_answer<float_type>(T);
REQUIRE(approximately_equal(T_corr, T, 1E-4));
}
SECTION("Geometric Mean") {
mod.change_input("UQ Type", "interval");
mod.change_input("UQ Type", uq_type);
mod.change_input("Mean Type", "geometric");
auto T = mod.run_as<test_pt>(braket);

auto T_corr = corr_answer<interval_type>(T);
auto T_corr = corr_answer<float_type>(T);
REQUIRE(approximately_equal(T_corr, T, 1E-4));
}
}
Expand Down
Loading
Loading