Skip to content

Fix container type and empty-vector defects in the composite quadratures #105

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
wavefunction91 merged 2 commits into from
Oct 8, 2026
Merged

Fix container type and empty-vector defects in the composite quadratures #105

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
dbbf796
Fix container type and empty-vector defects in the composite quadratures
susilehtola Aug 29, 2026
7ab221b
Close the remaining empty-partition UB and cover both fixes
Oct 6, 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
17 changes: 14 additions & 3 deletions include/integratorxx/composite_quadratures/pruned_spherical_quadrature.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include <integratorxx/types.hpp>
#include <integratorxx/type_traits.hpp>

#include <stdexcept>
#include <utility>
#include <cassert>
#include <memory>
Expand Down Expand Up @@ -48,6 +49,9 @@ class RadialGridPartition {

template <typename RadialQuad>
void finalize( const RadialQuad& rq ) {
if( partition_idx_.empty() )

@wavefunction91 wavefunction91 Oct 5, 2026

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The same default-constructed partition still hits UB through the iterators: end(), cend() and end() const (lines 141/148/155) compute partition_idx_.end()-1, which on an empty vector steps before begin(). With -D_GLIBCXX_DEBUG on this head, RadialGridPartition<LebedevLaikov<double>>{}.end() aborts with "attempt to retreat a past-the-end iterator". Since this PR is about the empty-vector cases, could those be made safe here too? A REQUIRE_THROWS next to the default-constructed rgp in test/composite_quadratures.cxx:202 would also pin the new finalize() behaviour.

@susilehtola susilehtola Oct 6, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right, and the iterators were the live path. Fixed in 7ab221b: end(), cend() and end() const present an empty range when partition_idx_ is empty instead of computing end()-1. I reproduced your abort first — -D_GLIBCXX_DEBUG on master gives "attempt to retreat a past-the-end iterator" from RadialGridPartition<LebedevLaikov<double>>{}.end(), and the patched version reports an empty range.

The new "Unfinalized Partition" section pins REQUIRE_THROWS( rgp.finalize(r) ) together with begin() == end() for all three overloads.

throw std::runtime_error(
"RadialGridPartition: cannot finalize without any angular quadrature");
if( partition_idx_.back() != rq.npts() )
partition_idx_.emplace_back( rq.npts() );
}
Expand Down Expand Up @@ -134,20 +138,27 @@ class RadialGridPartition {
return iterator( partition_idx_.begin(), quads_.begin() );
}
iterator end() {
// partition_idx_ carries one entry past the last partition (the npts
// sentinel added by finalize), so the end index iterator is one before its
// end. An unfinalized partition holds no entries at all, and end()-1 would
// step before begin(); the range is simply empty in that case.
if( partition_idx_.empty() ) return begin();
return iterator( partition_idx_.end()-1, quads_.end() );
}

const_iterator cbegin() const {
return const_iterator( partition_idx_.cbegin(), quads_.cbegin() );
}
const_iterator cend() const {
if( partition_idx_.empty() ) return cbegin();
return const_iterator( partition_idx_.cend()-1, quads_.cend() );
}

const_iterator begin() const {
return const_iterator( partition_idx_.cbegin(), quads_.cbegin() );
}
const_iterator end() const {
if( partition_idx_.empty() ) return begin();
return const_iterator( partition_idx_.cend()-1, quads_.cend() );
}
};
Expand Down Expand Up @@ -183,13 +194,13 @@ class PrunedSphericalQuadrature :
const point_container& sph_points_adaptor() const override {
return quad_base_type::points();
}
point_container& sph_points_adaptor() override {;
point_container& sph_points_adaptor() override {
return quad_base_type::points();
}
const weight_container& sph_weights_adaptor() const override {;
const weight_container& sph_weights_adaptor() const override {
return quad_base_type::weights();
}
weight_container& sph_weights_adaptor() override {;
weight_container& sph_weights_adaptor() override {
return quad_base_type::weights();
}

Expand Down
6 changes: 3 additions & 3 deletions include/integratorxx/composite_quadratures/spherical_quadrature.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,13 +133,13 @@ class SphericalQuadrature :
const point_container& sph_points_adaptor() const override {
return quad_base_type::points();
}
point_container& sph_points_adaptor() override {;
point_container& sph_points_adaptor() override {
return quad_base_type::points();
}
const weight_container& sph_weights_adaptor() const override {;
const weight_container& sph_weights_adaptor() const override {
return quad_base_type::weights();
}
weight_container& sph_weights_adaptor() override {;
weight_container& sph_weights_adaptor() override {
return quad_base_type::weights();
}

Expand Down
17 changes: 7 additions & 10 deletions include/integratorxx/composite_quadratures/sub_quadrature.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,16 +49,13 @@ struct quadrature_traits<
inline static std::tuple<point_container,weight_container>
generate( range_idx_type idx_range, const Quadrature<DerivedQuad>& quad ) {
const auto [begin_idx, end_idx] = idx_range;
assert( begin_idx >= 0 );
assert( begin_idx <= quad.npts() );
assert( end_idx >= 0 );
assert( end_idx <= quad.npts() );
assert( begin_idx < end_idx );

point_container points( quad.points().begin() + begin_idx,
quad.points().begin() + end_idx );
point_container weights( quad.weights().begin() + begin_idx,
quad.weights().begin() + end_idx );
assert( begin_idx < end_idx );
assert( end_idx <= quad.npts() );

point_container points ( quad.points().begin() + begin_idx,
quad.points().begin() + end_idx );
weight_container weights( quad.weights().begin() + begin_idx,

@wavefunction91 wavefunction91 Oct 5, 2026

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nothing in the suite instantiates SubQuadrature, so the mixed-precision case this fixes isn't covered, and it would regress silently. Could we add a small case, e.g. SubQuadrature<Becke<float,double>> over a slice of a Becke<float,double>, checking that the points and (double) weights match the parent's slice?

@susilehtola susilehtola Oct 6, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added as a "Sub Quadratures" case over Becke<float,double>, checking the slice against the parent exactly, plus static_asserts on both container types.

One correction to the premise, in the good direction: it would not have regressed silently. With master's sub_quadrature.hpp the instantiation does not compile at all — make_tuple yields tuple<vector<float>,vector<float>>, which will not convert to tuple<vector<float>,vector<double>>. So the case was simply absent rather than passing for the wrong reason.

quad.weights().begin() + end_idx );

return std::make_tuple( points, weights );
}
Expand Down
45 changes: 45 additions & 0 deletions test/composite_quadratures.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include <integratorxx/quadratures/s2.hpp>
#include <integratorxx/composite_quadratures/spherical_quadrature.hpp>
#include <integratorxx/composite_quadratures/pruned_spherical_quadrature.hpp>
#include <integratorxx/composite_quadratures/sub_quadrature.hpp>
#include <integratorxx/batch/spherical_micro_batcher.hpp>
#include <numeric>

Expand Down Expand Up @@ -239,6 +240,29 @@ TEST_CASE( "Pruned Spherical Quadratures", "[sph-quad]" ) {
}


SECTION("Unfinalized Partition") {

size_t nrad = 4;
IntegratorXX::MuraKnowles<double,double> r(nrad);
IntegratorXX::LebedevLaikov<double> q(6);
IntegratorXX::RadialGridPartition<decltype(q)> rgp;

// Nothing has been added, so there is no last index to compare against
// rq.npts(): finalize must say so rather than read back() off an empty
// vector.
REQUIRE_THROWS( rgp.finalize(r) );

// The iterators must present an empty range over it. partition_idx_ is
// empty, so the end()-1 that skips the npts sentinel would step before
// begin().
CHECK( rgp.begin() == rgp.end() );
CHECK( rgp.cbegin() == rgp.cend() );
const auto& crgp = rgp;
CHECK( crgp.begin() == crgp.end() );

}


SECTION("Validity") {

size_t nrad = 12;
Expand Down Expand Up @@ -300,9 +324,30 @@ TEST_CASE( "Pruned Spherical Quadratures", "[sph-quad]" ) {
}


TEST_CASE( "Sub Quadratures", "[sph-quad]" ) {

// The slice's weights must be built in the parent's weight_container, not in
// its point_container. The two coincide for the usual <double,double> case,
// which is why nothing caught it; under mixed precision the weights would be
// narrowed to the point type.
IntegratorXX::Becke<float,double> r(20);

const std::pair<size_t,size_t> range{ 5ul, 12ul };
IntegratorXX::SubQuadrature<decltype(r)> sq( range, r );

static_assert( std::is_same_v< decltype(sq)::weight_container,
decltype(r)::weight_container >,
"sub-quadrature weights must keep the parent's weight type" );
static_assert( std::is_same_v< decltype(sq)::point_container,
decltype(r)::point_container >,
"sub-quadrature points must keep the parent's point type" );

REQUIRE( sq.npts() == range.second - range.first );

// The slice is a copy, so it must agree with the parent exactly
for( size_t i = 0; i < sq.npts(); ++i ) {
CHECK( sq.points()[i] == r.points()[i + range.first] );
CHECK( sq.weights()[i] == r.weights()[i + range.first] );
}

}
Loading