From e3e0f13fdff2b55110ba45bd2702ff4a52256b6c Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 19:29:38 +0300 Subject: [PATCH 01/17] tools: generate angular grids at arbitrary precision The tabulated Lebedev-Laikov and Delley grids carry about 16 significant digits, so LebedevLaikov cannot be more accurate than double however wide T is. The format makes this worse than it looks: the tables are unsuffixed floating literals inside std::array, and an unsuffixed literal has type double whatever T is, so writing more digits into the existing format would change nothing at all. Add a generator that refines the existing rules to any requested precision. Finding a spherical quadrature rule and refining one are different problems, and only the second is needed here. Construction is a global optimisation with many local minima and wants second-order globalisation (a trust-region solver), but only at double precision, since the goal is just to land in the right basin. Refinement starts from a rule already correct to 16 digits, is purely local, and Newton converges quadratically -- two or three iterations reach any precision, with no globalisation. For the existing families construction is already done: the tables are the result. What makes it cheap is symmetry. Both families are octahedrally invariant, so points fall into orbits carrying one weight and zero, one or two angular parameters. The 5810-point Lebedev rule is a 385-parameter dense Newton solve rather than a 17430-parameter one. The exactness conditions are imposed on even monomials in x and y only: on the unit sphere z^2 = 1 - x^2 - y^2, so every even monomial reduces to a fixed linear combination of those, and imposing them implies the rest at about a seventh of the cost. Results are always verified against the full redundant set, never the reduced one. Emitted headers carry decimal strings rather than literals, converted through detail::grid_scalar::parse at construction. That works for the built-in types as well as Boost.Multiprecision and MPFR, and costs one parse per grid. Measured on the 110-point rule integrating x^4 y^2 z^2: current table regenerated double 1.91e-15 1.74e-15 long double 1.26e-15 2.55e-19 cpp_bin_float_50 unreachable 1.94e-40 This commit adds the tool and the scalar-parsing helper only; no shipped table is regenerated yet. Ahrens-Beylkin needs icosahedral orbit algebra and Womersley needs a matrix-free least-squares step; both are described in the README. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- include/integratorxx/util/grid_scalar.hpp | 42 +++++++ tools/.gitignore | 2 + tools/angular_grids/README.md | 80 +++++++++++++ tools/angular_grids/__init__.py | 0 tools/angular_grids/emit.py | 96 ++++++++++++++++ tools/angular_grids/generate.py | 76 +++++++++++++ tools/angular_grids/moments.py | 121 ++++++++++++++++++++ tools/angular_grids/orbits.py | 130 ++++++++++++++++++++++ tools/angular_grids/refine.py | 127 +++++++++++++++++++++ tools/angular_grids/tables.py | 32 ++++++ tools/angular_grids/test_refine.py | 46 ++++++++ 11 files changed, 752 insertions(+) create mode 100644 include/integratorxx/util/grid_scalar.hpp create mode 100644 tools/.gitignore create mode 100644 tools/angular_grids/README.md create mode 100644 tools/angular_grids/__init__.py create mode 100644 tools/angular_grids/emit.py create mode 100644 tools/angular_grids/generate.py create mode 100644 tools/angular_grids/moments.py create mode 100644 tools/angular_grids/orbits.py create mode 100644 tools/angular_grids/refine.py create mode 100644 tools/angular_grids/tables.py create mode 100644 tools/angular_grids/test_refine.py diff --git a/include/integratorxx/util/grid_scalar.hpp b/include/integratorxx/util/grid_scalar.hpp new file mode 100644 index 0000000..5bf616f --- /dev/null +++ b/include/integratorxx/util/grid_scalar.hpp @@ -0,0 +1,42 @@ +#pragma once + +#include +#include + +namespace IntegratorXX { +namespace detail { + +/** + * @brief Convert the decimal digits of a tabulated grid to the working type. + * + * Grid data is stored as decimal strings rather than floating literals: an + * unsuffixed literal has type double whatever the template parameter is, so a + * literal table cannot carry more precision than double however many digits + * are written into it. + * + * The primary template covers extended-precision types that construct from a + * string (Boost.Multiprecision, MPFR C++). The built-in types are specialised + * onto the corresponding strtoX. + */ +template +struct grid_scalar { + static T parse(const char* s) { return T(s); } +}; + +template <> +struct grid_scalar { + static float parse(const char* s) { return std::strtof(s, nullptr); } +}; + +template <> +struct grid_scalar { + static double parse(const char* s) { return std::strtod(s, nullptr); } +}; + +template <> +struct grid_scalar { + static long double parse(const char* s) { return std::strtold(s, nullptr); } +}; + +} // namespace detail +} // namespace IntegratorXX diff --git a/tools/.gitignore b/tools/.gitignore new file mode 100644 index 0000000..7a60b85 --- /dev/null +++ b/tools/.gitignore @@ -0,0 +1,2 @@ +__pycache__/ +*.pyc diff --git a/tools/angular_grids/README.md b/tools/angular_grids/README.md new file mode 100644 index 0000000..65a792f --- /dev/null +++ b/tools/angular_grids/README.md @@ -0,0 +1,80 @@ +# Angular grid generation + +Regenerates the Lebedev-Laikov and Delley tables at arbitrary precision. + +## Why + +The tabulated grids carry about 16 significant digits, so `LebedevLaikov` +cannot be more accurate than `double` however wide `T` is. Worse, the tables +are written as unsuffixed floating literals inside `std::array`, and an +unsuffixed literal has type `double` whatever `T` is -- so simply writing more +digits into the existing format would change nothing. + +Measured on the 110-point rule, integrating `x^4 y^2 z^2`: + +| precision | current table | regenerated table | +|--------------------|------------------|-------------------| +| `double` | 1.91e-15 | 1.74e-15 | +| `long double` | 1.26e-15 | **2.55e-19** | +| `cpp_bin_float_50` | not achievable | **1.94e-40** | + +## How + +Finding a spherical quadrature rule and refining one are very different +problems, and only the second is needed here: + +* **Construction** is a global optimisation with many local minima. It needs + second-order globalisation -- this is what a trust-region solver such as + [OpenTrustRegion](https://github.com/eriksen-lab/opentrustregion) is for -- + but `double` precision suffices, since the goal is only to land in the right + basin. +* **Refinement** starts from a rule that is already correct to 16 digits and is + purely local. Newton converges quadratically, so two or three iterations + reach any precision asked for, with no globalisation at all. + +For the existing families construction is already done: the tables *are* the +result. Only refinement is needed, and that is what this tool does. + +Tractability comes from symmetry. Both families are invariant under the full +octahedral group, so points fall into orbits carrying one weight and zero, one +or two angular parameters: + +| grid | points | orbits | free parameters | +|-------------------|-------:|-------:|----------------:| +| Lebedev-Laikov 50 | 50 | 4 | 5 | +| Lebedev-Laikov110 | 110 | 6 | 10 | +| Lebedev-Laikov590 | 590 | 20 | 44 | +| Lebedev-Laikov5810| 5810 | 144 | 385 | + +So the largest Lebedev rule is a 385-parameter dense Newton solve, not a +17430-parameter one. + +The exactness conditions are imposed on even monomials. On the unit sphere +`z^2 = 1 - x^2 - y^2`, so every even monomial reduces to a fixed linear +combination of `x^(2a) y^(2b)`; imposing only those implies the rest, at about +a seventh of the cost. Results are always *verified* against the full redundant +set. + +## Usage + +```sh +PYTHONPATH=tools python3 -m angular_grids.generate \ + --family lebedev_laikov --npts 110 --digits 40 --out gen/ +``` + +Requires `mpmath`. Output headers store the digits as decimal strings and +convert them through `detail::grid_scalar::parse` at grid construction, +which works for the built-in types as well as Boost.Multiprecision and MPFR. + +## Not covered + +* **Ahrens-Beylkin** has icosahedral rather than octahedral symmetry (15012 + points / 60 = 250.2 against 251 distinct weights, the rotation group without + inversion). Same method, different orbit algebra. +* **Womersley** grids are equal-weight spherical designs with no symmetry at + all: 2N angular unknowns, and the Jacobian is rank-deficient by exactly three + because any rotation of a design is a design. They refine fine, but need a + matrix-free least-squares step (CGLS) rather than a dense solve, and at 7939 + points the residual evaluation wants a compiled high-precision kernel. +* **Constructing new rules**, including regenerating the corrupt 552-point + Ahrens-Beylkin table. That is the trust-region half of the problem. diff --git a/tools/angular_grids/__init__.py b/tools/angular_grids/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tools/angular_grids/emit.py b/tools/angular_grids/emit.py new file mode 100644 index 0000000..7a0b8a9 --- /dev/null +++ b/tools/angular_grids/emit.py @@ -0,0 +1,96 @@ +"""Write refined grids out as C++ headers. + +A note on why the tables are emitted as decimal strings rather than numeric +literals. The existing data headers say + + template + struct lebedev_laikov_50 { + static constexpr std::array, 50> points = { + 0.3015113445777636, ... + +A floating literal with no suffix has type ``double`` whatever ``T`` is, so the +value is rounded to 53 bits before it is ever converted. Writing more digits +changes nothing: the table cannot carry more than double precision, and a +``long double`` or ``__float128`` instantiation silently gets a double-rounded +grid. Suffixes do not fix it either, since one table would have to serve every +``T``. + +Emitting the decimal digits as strings and converting them at grid construction +sidesteps this entirely, works for MPFR and Boost.Multiprecision as well as the +built-in types, and costs one parse per grid. +""" +from mpmath import mp + +HEADER = """#pragma once + +// Generated by tools/angular_grids -- do not edit. +// +// {family} quadrature on S^2, {npts} points, algebraic order {order}. +// Refined from the tabulated double-precision rule by Gauss-Newton on the +// exactness conditions; see tools/angular_grids/README.md. +// +// Digits carried: {digits}. Verified worst relative error over all even +// monomials up to degree {order}: {error}. + +#include +#include + +#include +#include + +namespace IntegratorXX {{ +namespace {ns} {{ + +template +struct {struct_name} {{ + static constexpr size_t npts = {npts}; + static constexpr int algebraic_order = {order}; + + /// Decimal digits of the quadrature points, x0 y0 z0 x1 y1 z1 ... + static constexpr const char* const point_digits[3 * {npts}] = {{ +{points} + }}; + + /// Decimal digits of the quadrature weights, normalised to sum to 4*pi. + static constexpr const char* const weight_digits[{npts}] = {{ +{weights} + }}; + + /// Materialise the grid at the working precision of T. + static void load(std::array, {npts}>& p, + std::array& w) {{ + for(size_t i = 0; i < {npts}; ++i) {{ + for(size_t j = 0; j < 3; ++j) + p[i][j] = detail::grid_scalar::parse(point_digits[3 * i + j]); + w[i] = detail::grid_scalar::parse(weight_digits[i]); + }} + }} +}}; + +}} // namespace {ns} +}} // namespace IntegratorXX +""" + +NAMESPACES = {"lebedev_laikov": "LebedevLaikovGrids", "delley": "DelleyGrids"} +TITLES = {"lebedev_laikov": "Lebedev-Laikov", "delley": "Delley"} + + +def _fmt(values, digits, per_line): + out, line = [], [] + for v in values: + line.append('"%s"' % mp.nstr(v, digits, strip_zeros=False)) + if len(line) == per_line: + out.append(" " + ", ".join(line) + ",") + line = [] + if line: + out.append(" " + ", ".join(line) + ",") + return "\n".join(out).rstrip(",") + + +def emit(family, npts, order, points, weights, digits, error): + flat = [c for p in points for c in p] + return HEADER.format( + family=TITLES.get(family, family), ns=NAMESPACES.get(family, family), + struct_name=f"{family}_{npts}", npts=npts, order=order, digits=digits, + error=mp.nstr(error, 4), + points=_fmt(flat, digits, 3), weights=_fmt(weights, digits, 3)) diff --git a/tools/angular_grids/generate.py b/tools/angular_grids/generate.py new file mode 100644 index 0000000..c2f4b31 --- /dev/null +++ b/tools/angular_grids/generate.py @@ -0,0 +1,76 @@ +"""Regenerate an angular grid table at arbitrary precision. + + python3 -m angular_grids.generate --family lebedev_laikov --npts 110 \ + --digits 40 --out gen/ + +Reads the tabulated double-precision rule, refines it by Gauss-Newton on the +exactness conditions, verifies the result against the full redundant set of +monomial conditions, and writes a header carrying the digits as strings. +""" +import argparse +import sys +from pathlib import Path + +from mpmath import mp + +from .emit import emit +from .refine import decompose, materialise, refine, verify +from .tables import grid_path, read_grid + + +def algebraic_order(family, npts, root): + """Read the order off the family's dispatch header.""" + import re + text = (Path(root) / "include/integratorxx/quadratures/s2" / f"{family}.hpp").read_text() + body = text[text.index("algebraic_order_by_npts"):] + m = re.search(r"case\s+%d\s*:\s*\n?\s*return\s+(\d+)" % npts, body) + if not m: + raise SystemExit(f"no algebraic order tabulated for {family} with {npts} points") + return int(m.group(1)) + + +def run(family, npts, digits, root, out, verbose=True): + mp.dps = digits + 20 # guard digits for the refinement + order = algebraic_order(family, npts, root) + points, weights = read_grid(grid_path(root, family, npts)) + scale = 4 * mp.pi / sum(weights) # tables are normalised to 1 + orbits = decompose(points, [w * scale for w in weights]) + n_param = sum(1 + o.type.n_param for o in orbits) + + if verbose: + print(f"{family}_{npts}: order {order}, {len(orbits)} orbits, {n_param} free parameters") + hist = refine(orbits, order, digits + 5, verbose=verbose) + P, W = materialise(orbits) + if len(P) != npts: + raise SystemExit(f"materialised {len(P)} points, expected {npts}") + + err = verify(P, W, order) # checked against the full set + if err > mp.mpf(10) ** (-digits): + raise SystemExit(f"refinement reached only {mp.nstr(err, 4)}, short of {digits} digits") + if verbose: + print(f" verified worst relative error {mp.nstr(err, 4)} over all even " + f"monomials up to degree {order}") + + text = emit(family, npts, order, P, W, digits, err) + if out: + path = Path(out) / f"{family}_{npts}.hpp" + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(text) + if verbose: + print(f" wrote {path}") + return text, err + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument("--family", required=True, choices=["lebedev_laikov", "delley"]) + ap.add_argument("--npts", required=True, type=int) + ap.add_argument("--digits", type=int, default=40) + ap.add_argument("--root", default=".") + ap.add_argument("--out", default=None) + a = ap.parse_args(argv) + run(a.family, a.npts, a.digits, a.root, a.out) + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tools/angular_grids/moments.py b/tools/angular_grids/moments.py new file mode 100644 index 0000000..7acb399 --- /dev/null +++ b/tools/angular_grids/moments.py @@ -0,0 +1,121 @@ +"""Exactness conditions for an O_h-symmetric spherical quadrature. + +A rule of algebraic order L integrates every spherical harmonic up to degree L +exactly. For an octahedrally symmetric rule it is equivalent, and much cheaper, +to impose exactness on the even monomials + + x^(2a) y^(2b) z^(2c), a >= b >= c >= 0, 2(a+b+c) <= L + +since the odd moments vanish by symmetry and the monomial conditions span the +same space. The conditions are redundant (x^2+y^2+z^2 = 1 relates them), which +is harmless: the resulting least-squares system is consistent. +""" +from itertools import permutations + +from mpmath import mp, mpf + + +def _dfact(n): + r = mpf(1) + while n > 1: + r *= n + n -= 2 + return r + + +def exact_moment(a, b, c): + """\\int_{S^2} x^(2a) y^(2b) z^(2c) dOmega.""" + return (4 * mp.pi * _dfact(2 * a - 1) * _dfact(2 * b - 1) * _dfact(2 * c - 1) + / _dfact(2 * (a + b + c) + 1)) + + +def conditions(order, reduced=True): + """Even monomials whose exactness pins down a rule of the given order. + + With ``reduced`` (the default) only monomials in x and y are used. On the + unit sphere z^2 = 1 - x^2 - y^2, so every even monomial reduces to a fixed + linear combination of x^(2a) y^(2b) with constant coefficients; the exact + integrals obey the same relation, so the residuals do too. Imposing the + reduced set therefore implies the rest, at roughly a seventh of the cost + for the larger grids. O_h symmetry lets us also require a >= b. + + Passing ``reduced=False`` returns the full redundant set, which is useful + as an independent check. + """ + out = [] + if reduced: + half = order // 2 + for a in range(half + 1): + for b in range(min(a, half - a) + 1): + out.append((a, b, 0)) + return out + for s in range(order // 2 + 1): + for a in range(s, -1, -1): + for b in range(min(a, s - a), -1, -1): + c = s - a - b + if 0 <= c <= b: + out.append((a, b, c)) + return out + + +def _distinct_index_perms(base): + """Index permutations giving distinct coordinate tuples. + + Equal entries in `base` always share the same parameter dependence, so + keeping one representative per distinct tuple is also correct for the + derivatives. + """ + seen = {} + for perm in permutations(range(3)): + key = tuple(mp.nstr(base[i], 30) for i in perm) + seen.setdefault(key, perm) + return list(seen.values()) + + +def _n_sign(base): + return 2 ** sum(1 for v in base if v != 0) + + +def orbit_moment(otype, params, powers): + """Sum of the monomial over the whole orbit, and its parameter gradient. + + Returns (value, [d/dparam_0, ...]). + """ + base = otype.base(params) + dbase = otype.dbase(params) + ns = _n_sign(base) + perms = _distinct_index_perms(base) + + total = mpf(0) + grad = [mpf(0)] * otype.n_param + + for perm in perms: + # term = prod_j base[perm[j]] ** (2 * powers[j]) + term = mpf(1) + for j, idx in enumerate(perm): + e = 2 * powers[j] + if e: + term *= base[idx] ** e + total += term + + for k in range(otype.n_param): + acc = mpf(0) + for j, idx in enumerate(perm): + e = 2 * powers[j] + if e == 0: + continue + dv = dbase[idx][k] + if dv == 0: + continue + # d/dparam of base[idx]^e, times the other factors + rest = mpf(1) + for j2, idx2 in enumerate(perm): + if j2 == j: + continue + e2 = 2 * powers[j2] + if e2: + rest *= base[idx2] ** e2 + acc += e * base[idx] ** (e - 1) * dv * rest + grad[k] += acc + + return ns * total, [ns * g for g in grad] diff --git a/tools/angular_grids/orbits.py b/tools/angular_grids/orbits.py new file mode 100644 index 0000000..5907efc --- /dev/null +++ b/tools/angular_grids/orbits.py @@ -0,0 +1,130 @@ +"""Octahedral (O_h) orbit algebra for Lebedev-Laikov and Delley grids. + +Both families are invariant under the full octahedral group with inversion, so +their points fall into a handful of orbit types. Each orbit carries one weight +and zero, one or two free angular parameters -- which is why a 5810-point +Lebedev rule has 385 free parameters rather than 17430. + +Orbit types follow the standard Lebedev naming: + + a1 6 pts (1,0,0) 0 parameters + a2 12 pts (1,1,0)/r2 0 parameters + a3 8 pts (1,1,1)/r3 0 parameters + bk 24 pts (l,l,m) 1 parameter, m = sqrt(1-2l^2) + ck 24 pts (p,q,0) 1 parameter, q = sqrt(1-p^2) + dk 48 pts (r,s,t) 2 parameters, t = sqrt(1-r^2-s^2) +""" +from itertools import permutations + +from mpmath import mp, mpf, sqrt + + +class OrbitType: + """One octahedral orbit type: how to build its points from its parameters.""" + + def __init__(self, name, n_param, base, dbase, size): + self.name = name + self.n_param = n_param + self._base = base # params -> (v0, v1, v2), the |coordinates| + self._dbase = dbase # params -> [[dv_j/dparam_k]], shape 3 x n_param + self.size = size + + def base(self, params): + return self._base(params) + + def dbase(self, params): + return self._dbase(params) + + +def _no_deriv(_): + return [[], [], []] + + +A1 = OrbitType("a1", 0, lambda _: (mpf(1), mpf(0), mpf(0)), _no_deriv, 6) +A2 = OrbitType("a2", 0, lambda _: (1 / sqrt(2), 1 / sqrt(2), mpf(0)), _no_deriv, 12) +A3 = OrbitType("a3", 0, lambda _: (1 / sqrt(3),) * 3, _no_deriv, 8) + + +def _b_base(p): + l = p[0] + return (l, l, sqrt(1 - 2 * l * l)) + + +def _b_dbase(p): + l = p[0] + m = sqrt(1 - 2 * l * l) + return [[mpf(1)], [mpf(1)], [-2 * l / m]] + + +def _c_base(p): + x = p[0] + return (x, sqrt(1 - x * x), mpf(0)) + + +def _c_dbase(p): + x = p[0] + q = sqrt(1 - x * x) + return [[mpf(1)], [-x / q], [mpf(0)]] + + +def _d_base(p): + r, s = p + return (r, s, sqrt(1 - r * r - s * s)) + + +def _d_dbase(p): + r, s = p + t = sqrt(1 - r * r - s * s) + return [[mpf(1), mpf(0)], [mpf(0), mpf(1)], [-r / t, -s / t]] + + +BK = OrbitType("bk", 1, _b_base, _b_dbase, 24) +CK = OrbitType("ck", 1, _c_base, _c_dbase, 24) +DK = OrbitType("dk", 2, _d_base, _d_dbase, 48) + +ALL_TYPES = (A1, A2, A3, BK, CK, DK) +BY_NAME = {t.name: t for t in ALL_TYPES} + + +def signed_permutations(base): + """The full octahedral orbit of a point given its |coordinates|.""" + out = [] + seen = set() + for perm in set(permutations(base)): + for s0 in (1, -1): + for s1 in (1, -1): + for s2 in (1, -1): + q = (s0 * perm[0], s1 * perm[1], s2 * perm[2]) + key = tuple(mp.nstr(v, 30) for v in q) + if key in seen: + continue + seen.add(key) + out.append(q) + return out + + +def classify(coords, tol=mpf("1e-9")): + """Identify the orbit type and parameters from one point's |coordinates|. + + Returns (OrbitType, [parameters]) with the coordinates put in canonical + order for that type. + """ + v = sorted((abs(c) for c in coords), reverse=True) + zero = [abs(c) < tol for c in v] + n_zero = sum(zero) + + if n_zero == 2: + return A1, [] + if n_zero == 1: + if abs(v[0] - v[1]) < tol: + return A2, [] + # (p, q, 0): the free parameter is the smaller nonzero coordinate, so + # that q = sqrt(1-p^2) reproduces the larger one. + return CK, [v[1]] + if abs(v[0] - v[1]) < tol and abs(v[1] - v[2]) < tol: + return A3, [] + if abs(v[1] - v[2]) < tol: + return BK, [v[1]] # (l,l,m) with l the repeated value + if abs(v[0] - v[1]) < tol: + return BK, [v[0]] + return DK, [v[0], v[1]] diff --git a/tools/angular_grids/refine.py b/tools/angular_grids/refine.py new file mode 100644 index 0000000..2d6e4e9 --- /dev/null +++ b/tools/angular_grids/refine.py @@ -0,0 +1,127 @@ +"""Gauss-Newton refinement of an O_h-symmetric angular grid. + +The tabulated grids are accurate to about 16 digits, which is the ceiling on +what LebedevLaikov can deliver however wide T is. Newton's method converges +quadratically from that starting point, so three or four iterations reach any +precision you ask for. This is a fundamentally easier problem than finding the +rule in the first place, and needs no globalisation. +""" +import collections + +from mpmath import mp, mpf, fabs, lu_solve, matrix + +from .moments import conditions, exact_moment, orbit_moment +from .orbits import classify, signed_permutations + + +class Orbit: + def __init__(self, otype, params, weight): + self.type = otype + self.params = list(params) + self.weight = weight # per-point weight + + def __repr__(self): + ps = ", ".join(mp.nstr(p, 12) for p in self.params) + return f"<{self.type.name} n={self.type.size} w={mp.nstr(self.weight, 12)} [{ps}]>" + + +def decompose(points, weights, tol=mpf("1e-9")): + """Group a tabulated grid into octahedral orbits.""" + groups = collections.OrderedDict() + for p, w in zip(points, weights): + key = (mp.nstr(w, 12), tuple(mp.nstr(abs(c), 9) for c in sorted((abs(v) for v in p), reverse=True))) + groups.setdefault(key, []).append((p, w)) + + orbits = [] + for members in groups.values(): + p0, w0 = members[0] + otype, params = classify(p0, tol) + if len(members) != otype.size: + raise ValueError(f"orbit of type {otype.name} has {len(members)} points, expected {otype.size}") + orbits.append(Orbit(otype, params, w0)) + return orbits + + +def _pack(orbits): + x = [] + for o in orbits: + x.append(o.weight) + x.extend(o.params) + return x + + +def _unpack(orbits, x): + k = 0 + for o in orbits: + o.weight = x[k]; k += 1 + for j in range(o.type.n_param): + o.params[j] = x[k]; k += 1 + + +def _residual_and_jacobian(orbits, conds, want_jac): + n = sum(1 + o.type.n_param for o in orbits) + F = [] + J = matrix(len(conds), n) if want_jac else None + + # orbit moments and gradients, cached per condition + for ci, powers in enumerate(conds): + total = mpf(0) + col = 0 + for o in orbits: + m, dm = orbit_moment(o.type, o.params, powers) + total += o.weight * m + if want_jac: + J[ci, col] = m; col += 1 + for k in range(o.type.n_param): + J[ci, col] = o.weight * dm[k]; col += 1 + else: + col += 1 + o.type.n_param + F.append(total - exact_moment(*powers)) + return F, J + + +def refine(orbits, order, target_dps, max_iter=12, verbose=False): + """Refine in place. Returns the residual history.""" + conds = conditions(order) + hist = [] + tol = mpf(10) ** (-target_dps) + for it in range(max_iter): + F, J = _residual_and_jacobian(orbits, conds, want_jac=True) + r = max(fabs(f) for f in F) + hist.append(r) + if verbose: + print(f" iter {it}: max|residual| = {mp.nstr(r, 4)}") + if r < tol: + break + n = J.cols + JtJ = matrix(n, n); Jtf = matrix(n, 1) + for i in range(n): + for j in range(i, n): + v = mp.fsum(J[k, i] * J[k, j] for k in range(J.rows)) + JtJ[i, j] = v; JtJ[j, i] = v + Jtf[i] = mp.fsum(J[k, i] * F[k] for k in range(J.rows)) + d = lu_solve(JtJ, -Jtf) + _unpack(orbits, [x + d[j] for j, x in enumerate(_pack(orbits))]) + return hist + + +def materialise(orbits): + """Expand refined orbits into the full point and weight lists.""" + points, weights = [], [] + for o in orbits: + for q in signed_permutations(o.type.base(o.params)): + points.append(q) + weights.append(o.weight) + return points, weights + + +def verify(points, weights, order): + """Worst relative error over all even monomials up to `order`.""" + worst = mpf(0) + for powers in conditions(order): + a, b, c = powers + q = mp.fsum(w * p[0] ** (2 * a) * p[1] ** (2 * b) * p[2] ** (2 * c) + for p, w in zip(points, weights)) + e = exact_moment(a, b, c) + worst = max(worst, fabs(q - e) / e) + return worst diff --git a/tools/angular_grids/tables.py b/tools/angular_grids/tables.py new file mode 100644 index 0000000..d6213dd --- /dev/null +++ b/tools/angular_grids/tables.py @@ -0,0 +1,32 @@ +"""Read the tabulated angular grids out of the generated C++ headers.""" +import re +from pathlib import Path + +from mpmath import mpf + +_NUM = re.compile(r"[-+]?\d*\.\d+[EeDd][-+]?\d+") + + +def _numbers(text): + return [mpf(t.replace("D", "e").replace("E", "e")) for t in _NUM.findall(text)] + + +def _section(text, name): + i = text.index(name + " = {") + j = text.index("};", i) + return text[i:j] + + +def read_grid(path): + """Return (points, weights) from a *_.hpp data header.""" + text = Path(path).read_text() + flat = _numbers(_section(text, "points")) + weights = _numbers(_section(text, "weights")) + points = [tuple(flat[3 * i:3 * i + 3]) for i in range(len(flat) // 3)] + if len(points) != len(weights): + raise ValueError(f"{path}: {len(points)} points but {len(weights)} weights") + return points, weights + + +def grid_path(root, family, npts): + return Path(root) / "include/integratorxx/quadratures/s2" / family / f"{family}_{npts}.hpp" diff --git a/tools/angular_grids/test_refine.py b/tools/angular_grids/test_refine.py new file mode 100644 index 0000000..b9eb38b --- /dev/null +++ b/tools/angular_grids/test_refine.py @@ -0,0 +1,46 @@ +"""Regression tests: the tool must reproduce the tabulated rules it starts from. + +Run with: PYTHONPATH=tools python3 -m angular_grids.test_refine +""" +import sys +from pathlib import Path + +from mpmath import mp, mpf + +from .moments import conditions +from .refine import decompose, materialise, refine, verify +from .tables import grid_path, read_grid + +CASES = [("lebedev_laikov", 6, 3), ("lebedev_laikov", 50, 11), ("lebedev_laikov", 110, 17), + ("lebedev_laikov", 194, 23), ("lebedev_laikov", 302, 29), + ("delley", 50, 11), ("delley", 194, 23), ("delley", 302, 29)] + + +def main(root="."): + mp.dps = 60 + failed = 0 + for family, npts, order in CASES: + pts, w = read_grid(grid_path(root, family, npts)) + orbits = decompose(pts, [x * 4 * mp.pi / sum(w) for x in w]) + n_param = sum(1 + o.type.n_param for o in orbits) + refine(orbits, order, 45) + P, W = materialise(orbits) + + ok = len(P) == npts + # verified against the FULL redundant condition set, not the reduced one + err = verify(P, W, order) + ok = ok and err < mpf("1e-45") + # the refined rule must stay close to the tabulated one, not wander to + # a different rule of the same order + drift = max(min(max(abs(a - b) for a, b in zip(p, q)) for q in P) for p in pts) + ok = ok and drift < mpf("1e-13") + + print(f" {'ok ' if ok else 'FAIL'} {family}_{npts:<5d} order {order:3d} " + f"{n_param:3d} params err {mp.nstr(err, 3):>10s} drift {mp.nstr(drift, 3):>10s}") + failed += not ok + print(f"\n {len(CASES) - failed}/{len(CASES)} passed") + return 1 if failed else 0 + + +if __name__ == "__main__": + sys.exit(main(sys.argv[1] if len(sys.argv) > 1 else ".")) From ee83270e72a10b845e81ff2149932cd1ca2e846e Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 20:51:28 +0300 Subject: [PATCH 02/17] tools: refine equal-weight spherical designs (Womersley) Adds the Womersley path to the generator. These grids have no symmetry to exploit -- every point is its own orbit -- so the unknowns are the 2N tangential degrees of freedom of N points and the Jacobian is never formed: CGLS needs only J*v and J^T*u. Points are carried as unit 3-vectors with steps taken in each point's tangent plane, not as (theta, phi). The spherical chart is singular at the poles and real designs put points there: the 32-point design has one at exactly z = 1, where d(p)/d(phi) vanishes. In that chart the Jacobian acquires a structurally zero column that is a coordinate artifact rather than a symmetry, and rotations about x and y are not expressible as finite tangent vectors at all. Measured on the 32-point design: zero cols near-null rotation generators (theta, phi) 1 3 3.4e-02, 1.3e-16, 2.0e-16 tangent frame 0 3 2.0e-16, 1.4e-16, 2.9e-16 The remaining three-dimensional null space is the genuine one: any rotation of a design is a design. It needs no special handling, because CGLS started from a zero step keeps every iterate in range(J^T), which is orthogonal to the null space, and so returns the minimum-norm -- already rotation-projected -- step. Design.rotation_generators() exposes them anyway, since a trust-region method used to *construct* a design does need to project them out. Two defects found by measurement rather than by reasoning, both recorded in the source so they are not reintroduced: * The tabulated points are unit vectors only to double precision, and carrying 3-vectors inherits that inconsistency with |p| = 1. It corrupts the Jacobian at the same 1e-16 level the refinement is trying to work below -- 8.0e-16 relative error against finite differences, 1.7e-21 after renormalising -- and stalled the first Newton step. The (theta, phi) chart had normalised implicitly. * CG terminates in rank-many steps only in exact arithmetic; its rate goes as the square root of the condition number. A cap of 2*2N starves it, and the outer iteration then converges convincingly to a spurious floor. On the 50-point design, capped at 200 iterations the residual stalls at 2.45e-18, while at 600 the same first step reaches 2.16e-28. The cap is now 12*2N with the residual tolerance doing the real stopping and a stagnation guard for early exit. With both fixed, against the full monomial set: womersley_32 (order 7): 7.39e-16 -> 3.82e-30 -> 9.18e-41 womersley_50 (order 9): 5.75e-16 -> 2.16e-28 -> 4.96e-41 Conditions use monomials with c in {0,1}: on the sphere z^2 reduces, leaving exactly (order+1)^2 of them -- the dimension of the polynomials of degree <= order on S^2 -- rather than C(order+3,3). For order 7 that is 64 rather than 120. read_grid() also learns the second weight format: the Womersley headers synthesise their equal weights with create_array(4*M_PI/N) instead of listing them. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/README.md | 77 +++++++- tools/angular_grids/designs.py | 271 +++++++++++++++++++++++++++++ tools/angular_grids/tables.py | 25 ++- tools/angular_grids/test_refine.py | 23 ++- 4 files changed, 387 insertions(+), 9 deletions(-) create mode 100644 tools/angular_grids/designs.py diff --git a/tools/angular_grids/README.md b/tools/angular_grids/README.md index 65a792f..5b94370 100644 --- a/tools/angular_grids/README.md +++ b/tools/angular_grids/README.md @@ -71,10 +71,77 @@ which works for the built-in types as well as Boost.Multiprecision and MPFR. * **Ahrens-Beylkin** has icosahedral rather than octahedral symmetry (15012 points / 60 = 250.2 against 251 distinct weights, the rotation group without inversion). Same method, different orbit algebra. -* **Womersley** grids are equal-weight spherical designs with no symmetry at - all: 2N angular unknowns, and the Jacobian is rank-deficient by exactly three - because any rotation of a design is a design. They refine fine, but need a - matrix-free least-squares step (CGLS) rather than a dense solve, and at 7939 - points the residual evaluation wants a compiled high-precision kernel. * **Constructing new rules**, including regenerating the corrupt 552-point Ahrens-Beylkin table. That is the trust-region half of the problem. + +## Equal-weight spherical designs (Womersley) + +`designs.py` handles the Womersley family, which has no symmetry to exploit: +every point is its own orbit, so the unknowns are the 2N tangential degrees of +freedom of N points. Two things differ from the octahedral path. + +**Points are carried as unit 3-vectors, not (theta, phi).** The spherical chart +is singular at the poles and real designs put points there -- the 32-point +design has one at exactly `z = 1`. In that chart `d(p)/d(phi)` vanishes, so the +Jacobian picks up a structurally zero column that is a coordinate artifact +rather than a symmetry, and rotations about x and y cease to be expressible as +finite tangent vectors. Measured on the 32-point design: + +| parameterisation | zero columns | near-null directions | rotation generators `\|Jg\|/\|J\|\|g\|` | +|-----------------------|-------------:|---------------------:|----------------------------------| +| `(theta, phi)` | 1 | 3 | 3.4e-02, 1.3e-16, 2.0e-16 | +| tangent frame | 0 | 3 | 2.0e-16, 1.4e-16, 2.9e-16 | + +Stepping in each point's tangent plane and renormalising (`retract`) is regular +everywhere, and all three rotations then behave alike -- the residual `1e-16` +being the accuracy of the tabulated input, not of the method. + +**The Jacobian is rank-deficient by exactly three**, since any rotation of a +design is a design. This needs no special handling: CGLS started from a zero +step keeps every iterate in `range(J^T)`, orthogonal to the null space, so it +returns the minimum-norm step -- already the rotation-projected one. Explicit +projection would be redundant. It *is* worth doing when constructing a design +with a trust-region method, where zero-curvature directions interfere with the +model and with the Hessian diagonal; `Design.rotation_generators()` returns them +for that purpose. + +The conditions use monomials with `c` in `{0,1}`: on the sphere `z^2` reduces, +leaving exactly `(order+1)^2` conditions -- the dimension of the polynomials of +degree <= order on S^2 -- rather than `C(order+3,3)`. For order 7 that is 64 +rather than 120. + +**Two traps, both measured rather than reasoned about.** + +*The input points must be renormalised.* The tabulated coordinates are unit +vectors only to double precision. Carrying them as 3-vectors inherits that +~1e-16 inconsistency with the constraint `|p| = 1`, which corrupts the Jacobian +at exactly the level the refinement is trying to work below -- checked against +finite differences, 8.0e-16 relative error before renormalising and 1.7e-21 +after. The `(theta, phi)` chart hid this by normalising implicitly, since every +`(theta, phi)` maps to an exactly unit vector. + +*The inner CGLS needs a generous iteration cap.* CG terminates in rank-many +steps only in exact arithmetic; its rate goes as the square root of the +condition number. Capping at `2 * 2N` starves it, the outer Gauss-Newton takes +a bad step, and the result looks convincingly like convergence to a spurious +floor: + +| CGLS cap | womersley_50 outer residuals | +|-----------------|-------------------------------------------------------| +| `2 * 2N` = 200 | 5.75e-16, 3.28e-18, 2.78e-18, 2.45e-18 (stalls) | +| `6 * 2N` = 600 | 5.75e-16, 2.16e-28 | + +The cap is now `12 * 2N`, with the residual tolerance doing the real stopping +and a stagnation guard so a converged solve exits at once. With both fixed: + +``` +womersley_32 (order 7): 7.39e-16 -> 3.82e-30 -> 9.18e-41 [107s] +womersley_50 (order 9): 5.75e-16 -> 2.16e-28 -> 4.96e-41 [358s] +``` + +both reaching the 40-digit working precision floor, verified against the full +monomial set rather than the reduced one. + +At 7939 points a single `J*v` product is around 1.3e8 high-precision operations +and the iteration count grows too, so the largest designs want a compiled MPFR +kernel rather than this driver. diff --git a/tools/angular_grids/designs.py b/tools/angular_grids/designs.py new file mode 100644 index 0000000..4107398 --- /dev/null +++ b/tools/angular_grids/designs.py @@ -0,0 +1,271 @@ +"""Matrix-free refinement of equal-weight spherical designs (Womersley grids). + +These grids have no symmetry to exploit: every point is its own orbit, so the +unknowns are the 2N tangential degrees of freedom of N points and the Jacobian +has 2N columns. Two consequences shape the implementation. + +*Points are carried as unit 3-vectors, not (theta, phi).* The spherical chart is +singular at the poles, and real designs put points there -- the 32-point design +has one at exactly z = 1. In that chart d(p)/d(phi) vanishes, so the Jacobian +acquires a structurally zero column that is a coordinate artifact rather than a +symmetry, and rotations about x and y are not expressible as finite tangent +vectors at all. Stepping in the tangent plane of each point and renormalising +avoids all of this: the parameterisation is regular everywhere on the sphere. + +*The Jacobian is rank-deficient by exactly three*, because any rotation of a +spherical design is a spherical design. This needs no special handling: CGLS +started from a zero step keeps every iterate in range(J^T), which is orthogonal +to the null space, so it returns the minimum-norm step -- already the +rotation-projected one. An explicit projection would be redundant here. It is +worth doing when *constructing* a design with a trust-region method, where zero +curvature directions interfere with the model. +""" +from mpmath import mp, mpf, fabs, sqrt + +from .moments import _dfact + + +def exact_monomial(a, b, c): + """\\int_{S^2} x^a y^b z^c dOmega; zero unless all exponents are even.""" + if a % 2 or b % 2 or c % 2: + return mpf(0) + return (4 * mp.pi * _dfact(a - 1) * _dfact(b - 1) * _dfact(c - 1) + / _dfact(a + b + c + 1)) + + +def conditions(order): + """A minimal spanning set of monomial conditions for a design of `order`. + + On the unit sphere z^2 = 1 - x^2 - y^2, so any monomial with c >= 2 reduces + to ones with c in {0,1}. That leaves exactly (order+1)^2 conditions -- the + dimension of the polynomials of degree <= order restricted to S^2 -- rather + than the C(order+3,3) of the full monomial set. + """ + out = [] + for c in (0, 1): + for d in range(order - c + 1): + for a in range(d + 1): + out.append((a, d - a, c)) + return out + + +def tangent_frame(p): + """An orthonormal basis of the tangent plane at the unit vector p. + + The seed axis is the one p is least aligned with, so the cross product is + always well conditioned -- no special case anywhere on the sphere. + """ + k = min(range(3), key=lambda i: fabs(p[i])) + seed = [mpf(0)] * 3 + seed[k] = mpf(1) + e1 = [seed[1] * p[2] - seed[2] * p[1], + seed[2] * p[0] - seed[0] * p[2], + seed[0] * p[1] - seed[1] * p[0]] + n = sqrt(sum(v * v for v in e1)) + e1 = [v / n for v in e1] + e2 = [p[1] * e1[2] - p[2] * e1[1], + p[2] * e1[0] - p[0] * e1[2], + p[0] * e1[1] - p[1] * e1[0]] + n2 = sqrt(sum(v * v for v in e2)) + return e1, [v / n2 for v in e2] + + +def retract(p, a, b, e1, e2): + """Move p by a*e1 + b*e2 and project back onto the sphere.""" + q = [p[i] + a * e1[i] + b * e2[i] for i in range(3)] + n = sqrt(sum(v * v for v in q)) + return [v / n for v in q] + + +def _grad(a, b, c, p): + x, y, z = p + return (a * x ** (a - 1) * y ** b * z ** c if a else mpf(0), + b * x ** a * y ** (b - 1) * z ** c if b else mpf(0), + c * x ** a * y ** b * z ** (c - 1) if c else mpf(0)) + + +class Design: + """An equal-weight spherical design under refinement.""" + + def __init__(self, points, order): + # Tabulated points are unit vectors only to double precision. Left + # alone, that ~1e-16 inconsistency between the coordinates and the + # constraint |p| = 1 corrupts the Jacobian at the same level and stalls + # the first Newton step. (The (theta, phi) chart normalised implicitly, + # since every (theta, phi) maps to an exactly unit vector; carrying + # 3-vectors gives that up and has to do it explicitly.) + self.points = [] + for p in points: + n = sqrt(sum(v * v for v in p)) + self.points.append([v / n for v in p]) + self.order = order + self.n = len(points) + self.weight = 4 * mp.pi / self.n + self.conds = conditions(order) + self._reframe() + + def _reframe(self): + self.frames = [tangent_frame(p) for p in self.points] + self._gcache = None + + def _gradients(self): + """Tangential gradients of every condition at every point. + + The points do not move during a CGLS solve, so these are computed once + per Newton step rather than once per matrix-vector product. This takes + every power evaluation out of the inner loop. + """ + if self._gcache is None: + tab = [] + for (a, b, c) in self.conds: + row = [] + for i, p in enumerate(self.points): + g = _grad(a, b, c, p) + e1, e2 = self.frames[i] + row.append((self.weight * sum(g[j] * e1[j] for j in range(3)), + self.weight * sum(g[j] * e2[j] for j in range(3)))) + tab.append(row) + self._gcache = tab + return self._gcache + + def residual(self): + F = [] + for (a, b, c) in self.conds: + s = mp.fsum(p[0] ** a * p[1] ** b * p[2] ** c for p in self.points) + F.append(self.weight * s - exact_monomial(a, b, c)) + return F + + def Jv(self, v): + """J @ v, without forming J.""" + out = [] + for row in self._gradients(): + s = mpf(0) + for i, (d1, d2) in enumerate(row): + s += v[2 * i] * d1 + v[2 * i + 1] * d2 + out.append(s) + return out + + def JTu(self, u): + """J^T @ u, without forming J.""" + out = [mpf(0)] * (2 * self.n) + for k, row in enumerate(self._gradients()): + uk = u[k] + if uk == 0: + continue + for i, (d1, d2) in enumerate(row): + out[2 * i] += uk * d1 + out[2 * i + 1] += uk * d2 + return out + + def rotation_generators(self): + """The three rotation directions, in tangent-frame coordinates.""" + gens = [] + for axis in ((1, 0, 0), (0, 1, 0), (0, 0, 1)): + ax = [mpf(t) for t in axis] + g = [mpf(0)] * (2 * self.n) + for i, p in enumerate(self.points): + v = (ax[1] * p[2] - ax[2] * p[1], + ax[2] * p[0] - ax[0] * p[2], + ax[0] * p[1] - ax[1] * p[0]) + e1, e2 = self.frames[i] + g[2 * i] = sum(v[j] * e1[j] for j in range(3)) + g[2 * i + 1] = sum(v[j] * e2[j] for j in range(3)) + gens.append(g) + return gens + + def step(self, F, iters=None, rtol=None): + """CGLS for min ||J d + F||; returns the minimum-norm step. + + Every Krylov vector lies in range(J^T), which is orthogonal to null(J), + so the three rotational directions are never entered and the result is + the minimum-norm -- already rotation-projected -- step. No explicit + projection is needed here. + + CG terminates in rank(J) steps only in exact arithmetic, and its rate + depends on the square root of the condition number rather than on the + rank, so a cap of rank-many iterations is not enough. Measured on the + 50-point design: capped at 2*2N = 200 the outer iteration stalls at + 2.5e-18, while at 6*2N = 600 the same first step reaches 2.2e-28. The + cap is therefore generous and the real stopping is by residual, with a + stagnation guard so a converged solve exits immediately. + """ + n = 2 * self.n + iters = iters or 12 * n + rtol = rtol or mpf(10) ** (-mp.dps + 5) + d = [mpf(0)] * n + r = [-f for f in F] + s = self.JTu(r) + p = list(s) + g = mp.fsum(q * q for q in s) + g0 = g + stall = 0 + for _ in range(iters): + q = self.Jv(p) + qq = mp.fsum(t * t for t in q) + if qq == 0: + break + al = g / qq + d = [d[j] + al * p[j] for j in range(n)] + r = [r[k] - al * q[k] for k in range(len(r))] + s = self.JTu(r) + gn = mp.fsum(t * t for t in s) + if gn == 0: + break + be = gn / g + if gn <= rtol * rtol * g0: + g = gn + break + # CG on a rank-deficient system eventually stops making progress; + # bail rather than grinding out the remaining cap. + if gn > mpf("0.999999") * g: + stall += 1 + if stall >= 20: + break + else: + stall = 0 + g = gn + p = [s[j] + be * p[j] for j in range(n)] + return d + + def apply(self, d): + self.points = [retract(p, d[2 * i], d[2 * i + 1], *self.frames[i]) + for i, p in enumerate(self.points)] + self._reframe() + + +def refine(design, target_dps, max_iter=8, verbose=False): + """Gauss-Newton with a fully converged inner solve. + + The inner CGLS is run to convergence rather than truncated by a forcing + term. An Eisenstat-Walker style truncation was tried and its effect could + not be separated from a Jacobian bug present at the time, so it is simply + not used: these systems are small enough that solving them properly costs + little, and correctness of the outer iteration is worth more than inner + iterations saved. + """ + hist = [] + tol = mpf(10) ** (-target_dps) + for it in range(max_iter): + F = design.residual() + r = max(fabs(f) for f in F) + hist.append(r) + if verbose: + print(f" iter {it}: max|residual| = {mp.nstr(r, 4)}", flush=True) + if r < tol: + break + design.apply(design.step(F)) + return hist + + +def verify(design, order=None): + """Worst absolute residual over the full monomial set, not the reduced one.""" + order = order or design.order + worst = mpf(0) + for d in range(order + 1): + for a in range(d + 1): + for b in range(d - a + 1): + c = d - a - b + q = mp.fsum(design.weight * p[0] ** a * p[1] ** b * p[2] ** c + for p in design.points) + worst = max(worst, fabs(q - exact_monomial(a, b, c))) + return worst diff --git a/tools/angular_grids/tables.py b/tools/angular_grids/tables.py index d6213dd..b4cf8ae 100644 --- a/tools/angular_grids/tables.py +++ b/tools/angular_grids/tables.py @@ -2,7 +2,7 @@ import re from pathlib import Path -from mpmath import mpf +from mpmath import mp, mpf _NUM = re.compile(r"[-+]?\d*\.\d+[EeDd][-+]?\d+") @@ -17,12 +17,31 @@ def _section(text, name): return text[i:j] +_EQUAL = re.compile(r"create_array<\s*(\d+)\s*,\s*T\s*>\s*\(\s*4\.0\s*\*\s*M_PI\s*/\s*(\d+)") + + def read_grid(path): - """Return (points, weights) from a *_.hpp data header.""" + """Return (points, weights) from a *_.hpp data header. + + Two weight formats occur. Lebedev-Laikov and Delley list the weights + explicitly. The Womersley grids are equal-weight spherical designs and + synthesise theirs with ``create_array(4.0 * M_PI / N)``, so there is + no array to read. + """ text = Path(path).read_text() flat = _numbers(_section(text, "points")) - weights = _numbers(_section(text, "weights")) points = [tuple(flat[3 * i:3 * i + 3]) for i in range(len(flat) // 3)] + + m = _EQUAL.search(text) + if m: + n, denom = int(m.group(1)), int(m.group(2)) + if n != denom or n != len(points): + raise ValueError(f"{path}: equal-weight size mismatch " + f"({n}, {denom}, {len(points)} points)") + w = 4 * mp.pi / n + return points, [w] * n + + weights = _numbers(_section(text, "weights")) if len(points) != len(weights): raise ValueError(f"{path}: {len(points)} points but {len(weights)} weights") return points, weights diff --git a/tools/angular_grids/test_refine.py b/tools/angular_grids/test_refine.py index b9eb38b..db12b53 100644 --- a/tools/angular_grids/test_refine.py +++ b/tools/angular_grids/test_refine.py @@ -7,6 +7,7 @@ from mpmath import mp, mpf +from . import designs from .moments import conditions from .refine import decompose, materialise, refine, verify from .tables import grid_path, read_grid @@ -15,6 +16,10 @@ ("lebedev_laikov", 194, 23), ("lebedev_laikov", 302, 29), ("delley", 50, 11), ("delley", 194, 23), ("delley", 302, 29)] +# Equal-weight spherical designs. Kept small: the inner CGLS is the cost, and +# it grows steeply with N. +DESIGN_CASES = [(14, 4), (18, 5), (26, 6)] + def main(root="."): mp.dps = 60 @@ -38,7 +43,23 @@ def main(root="."): print(f" {'ok ' if ok else 'FAIL'} {family}_{npts:<5d} order {order:3d} " f"{n_param:3d} params err {mp.nstr(err, 3):>10s} drift {mp.nstr(drift, 3):>10s}") failed += not ok - print(f"\n {len(CASES) - failed}/{len(CASES)} passed") + total = len(CASES) + for npts, order in DESIGN_CASES: + pts, _ = read_grid(grid_path(root, "womersley", npts)) + d = designs.Design(pts, order) + hist = designs.refine(d, 22) + err = designs.verify(d) + + ok = len(d.points) == npts and err < mpf("1e-20") + # a design keeps its equal weights and stays on the sphere + ok = ok and all(abs(sum(v * v for v in p) - 1) < mpf("1e-25") for p in d.points) + print(f" {'ok ' if ok else 'FAIL'} womersley_{npts:<5d} order {order:3d} " + f"{2 * npts:3d} unknowns {len(d.conds):3d} conds " + f"{len(hist)} iters err {mp.nstr(err, 3):>10s}") + failed += not ok + total += 1 + + print(f"\n {total - failed}/{total} passed") return 1 if failed else 0 From ac1c312b3a72ca0f5ffeb832ac5eba5fe8ba397d Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 21:09:54 +0300 Subject: [PATCH 03/17] tools: icosahedral orbit algebra for the Ahrens-Beylkin grids The AB grids are invariant under the icosahedral rotation group I (order 60) rather than the octahedral group, and without inversion -- which is why they are not antipodally symmetric and why 15012 points divide as 15012/60 = 250.2 against 251 distinct weights. The group is constructed numerically instead of from memorised generators: it is exactly the set of rotations carrying the icosahedron's vertices to themselves, and each is fixed by naming where one vertex and one of its neighbours go. Verified: 60 elements, all proper rotations, max |R^T R - I| = 3.9e-31. Every tabulated AB grid decomposes cleanly: AB-552 10 orbits 12 + 9x60 28 free parameters AB-612 11 orbits 12 + 10x60 31 free parameters AB-15012 251 orbits 12 + 250x60 751 free parameters so even the largest is a 751-parameter problem rather than a 45036-parameter one. Two tolerance defects fixed along the way, both the same shape -- a tolerance tighter than the accuracy of the data it is applied to: * read_grid()'s numeric pattern required an exponent. ahrens_beylkin_552 mixes plain decimals with exponent form inside one array, so the parser silently read 48 points where there are 552. Any file in that format was being misread. * orbit() de-duplicated at 1e-18 while its input is tabulated to double precision, so images that should coincide differ by ~1e-16 and a 12-point orbit came back as 60 near-duplicates. De-duplication is by distance rather than by formatting coordinates, since a coordinate that should be zero comes out as +-1e-31 and two such values have different decimal representations. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/icosahedral.py | 107 +++++++++++++++++++++++++++++ tools/angular_grids/tables.py | 5 +- 2 files changed, 111 insertions(+), 1 deletion(-) create mode 100644 tools/angular_grids/icosahedral.py diff --git a/tools/angular_grids/icosahedral.py b/tools/angular_grids/icosahedral.py new file mode 100644 index 0000000..7dc3a56 --- /dev/null +++ b/tools/angular_grids/icosahedral.py @@ -0,0 +1,107 @@ +"""Icosahedral orbit algebra for the Ahrens-Beylkin grids. + +Unlike Lebedev-Laikov and Delley, the AB grids are invariant under the +icosahedral *rotation* group I (order 60), without inversion -- which is why +they are not antipodally symmetric and why 15012 points divide as 15012/60 = +250.2 against 251 distinct weights. + +Orbits under I: + + 12 pts icosahedron vertices 0 parameters + 20 pts dodecahedron vertices 0 parameters + 30 pts edge midpoints 0 parameters + 60 pts generic position 2 parameters + +The group is built numerically rather than from memorised generators: it is +exactly the set of rotations carrying the icosahedron's vertex set to itself, +and each is fixed by naming where one vertex and one of its neighbours go. +""" +from mpmath import mp, mpf, sqrt, matrix + + +def _icosahedron_vertices(): + """The 12 unit vertices: cyclic permutations of (0, +-1, +-phi).""" + phi = (1 + sqrt(5)) / 2 + n = sqrt(1 + phi * phi) + out = [] + for s1 in (1, -1): + for s2 in (1, -1): + a, b = mpf(s1) / n, mpf(s2) * phi / n + out.append((mpf(0), a, b)) + out.append((b, mpf(0), a)) + out.append((a, b, mpf(0))) + return out + + +def _frame(v, n): + """Right-handed orthonormal frame from a vertex and one of its neighbours.""" + e1 = list(v) + d = sum(n[i] * e1[i] for i in range(3)) + e2 = [n[i] - d * e1[i] for i in range(3)] + m = sqrt(sum(x * x for x in e2)) + e2 = [x / m for x in e2] + e3 = [e1[1] * e2[2] - e1[2] * e2[1], + e1[2] * e2[0] - e1[0] * e2[2], + e1[0] * e2[1] - e1[1] * e2[0]] + return e1, e2, e3 + + +def rotation_group(tol=mpf("1e-20")): + """The 60 rotation matrices of the icosahedral group.""" + V = _icosahedron_vertices() + # nearest neighbours: the 5 vertices at minimum positive distance + def neighbours(v): + ds = [] + for w in V: + d = sum((v[i] - w[i]) ** 2 for i in range(3)) + if d > tol: + ds.append((d, w)) + ds.sort(key=lambda t: t[0]) + dmin = ds[0][0] + return [w for d, w in ds if d < dmin * mpf("1.5")] + + v0 = V[0] + n0 = neighbours(v0)[0] + F0 = _frame(v0, n0) + # R maps the reference frame onto (v, n): R = F(v,n)^T . F0 in row form + G = [] + for v in V: + for n in neighbours(v): + F = _frame(v, n) + R = matrix(3, 3) + for i in range(3): + for j in range(3): + R[i, j] = mp.fsum(F[k][i] * F0[k][j] for k in range(3)) + G.append(R) + return G + + +_GROUP = None + + +def group(): + global _GROUP + if _GROUP is None: + _GROUP = rotation_group() + return _GROUP + + +def orbit(p, tol=mpf("1e-9")): + """The icosahedral orbit of a point, de-duplicated. + + Duplicates are found by distance, not by formatting the coordinates: a + coordinate that should be zero comes out as +-1e-31, and two such values + have different decimal representations while being the same point. + + The default tolerance is loose on purpose. Orbits are usually generated + from tabulated points that are only accurate to double precision, so + images that should coincide differ by ~1e-16; a tolerance tighter than the + input accuracy silently splits a 12-point orbit into 60 near-duplicates. + Distinct points in these grids are separated by far more than 1e-9. + """ + out = [] + for R in group(): + q = tuple(mp.fsum(R[i, j] * p[j] for j in range(3)) for i in range(3)) + if not any(sum((q[i] - r[i]) ** 2 for i in range(3)) < tol * tol for r in out): + out.append(q) + return out diff --git a/tools/angular_grids/tables.py b/tools/angular_grids/tables.py index b4cf8ae..90e42a5 100644 --- a/tools/angular_grids/tables.py +++ b/tools/angular_grids/tables.py @@ -4,7 +4,10 @@ from mpmath import mp, mpf -_NUM = re.compile(r"[-+]?\d*\.\d+[EeDd][-+]?\d+") +# Data files mix plain decimals with exponent form, sometimes within one array +# (ahrens_beylkin_552.hpp has both). Requiring an exponent silently drops most +# of the values, so make it optional. +_NUM = re.compile(r"[-+]?\d*\.\d+(?:[EeDd][-+]?\d+)?") def _numbers(text): From 033230a0ee718f4469b421baeba9abbb7881d7d5 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 21:17:40 +0300 Subject: [PATCH 04/17] tools: refine the Ahrens-Beylkin grids Adds the icosahedral refinement path, mirroring the octahedral one. Each orbit carries one weight; a generic 60-point orbit carries two more parameters as the position of its representative, stepped in that point's tangent plane rather than through spherical angles. Because the orbit is generated by the group action, the derivative of an orbit point with respect to a tangent step of its representative is just the group element applied to that step. AB-72 order 14 2 orbits 4 params 1.75e-15 -> 1.58e-30 AB-132 order 19 3 orbits 7 params 4.34e-16 -> 7.89e-31 AB-192 order 23 4 orbits 10 params 7.66e-16 -> 1.58e-30 two Newton steps each, from the tabulated double-precision data. Two defects found and fixed while getting there: * decompose() used the rounded weight as both the grouping key and the weight value, truncating it to the key's 12 digits. It showed up as 1.7e-11 in the sum of weights while the points matched to 4.7e-33. * A 12-, 20- or 30-point orbit has no free parameter -- its position is fixed by the symmetry -- but taking the representative from tabulated double-precision data makes it a *generic* point whose 60 images cluster in fives rather than coinciding. De-duplication then keeps an arbitrary member of each cluster and injects ~1e-16 of asymmetry into a grid that should be exactly symmetric. Special orbits are now snapped onto the exact vertex, face-centre or edge-midpoint position, rotated onto the tabulated one. The second was masked by the residual: summing a monomial over all 60 group images annihilates odd moments identically, whatever the representative, so the odd-moment conditions carry no information about the fit and could not see the asymmetry. Verification against the materialised grid does see it, and now agrees with the residual. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/ab_refine.py | 183 +++++++++++++++++++++++++++++ tools/angular_grids/icosahedral.py | 47 ++++++++ 2 files changed, 230 insertions(+) create mode 100644 tools/angular_grids/ab_refine.py diff --git a/tools/angular_grids/ab_refine.py b/tools/angular_grids/ab_refine.py new file mode 100644 index 0000000..b3188ea --- /dev/null +++ b/tools/angular_grids/ab_refine.py @@ -0,0 +1,183 @@ +"""Refinement of the Ahrens-Beylkin grids. + +Same Gauss-Newton idea as the octahedral path, with the icosahedral rotation +group in place of O_h. Each orbit contributes one weight; a generic (60-point) +orbit contributes two more parameters, carried as the position of its +representative on the sphere. + +The representative is stepped in its own tangent plane and renormalised, as in +designs.py, rather than through spherical angles: the (theta, phi) chart is +singular at the poles, and the 12-point orbit sits exactly there. + +Because the whole orbit is generated by the group action, the derivative of an +orbit point with respect to a tangent step of its representative is just the +group element applied to that step: d(R p)/dp . e = R e. +""" +from mpmath import mp, mpf, fabs, sqrt, matrix, lu_solve + +from .designs import _grad, conditions, exact_monomial, retract, tangent_frame +from .icosahedral import group, orbit, special_positions + +SPECIAL_SIZES = (12, 20, 30) # orbits with no free angular parameter + + +class ABOrbit: + def __init__(self, rep, weight, size): + n = sqrt(sum(v * v for v in rep)) + self.rep = [v / n for v in rep] + self.weight = weight # per-point weight + self.size = size + self.free = size not in SPECIAL_SIZES + self.frame = tangent_frame(self.rep) + + def points(self): + return orbit(self.rep) + + def __repr__(self): + return f"" + + +def decompose(points, weights, tol=mpf("1e-9")): + """Group a tabulated AB grid into icosahedral orbits.""" + import collections + groups = collections.OrderedDict() + # The rounded weight is a grouping key only. Reconstructing the weight from + # that string would silently truncate it to the key's 12 digits, which + # shows up as an error of ~1e-11 in the sum of weights. + for p, w in zip(points, weights): + groups.setdefault(mp.nstr(w, 12), []).append((p, w)) + special = special_positions() + out = [] + for members in groups.values(): + rep, w0 = members[0] + size = len(members) + if size in SPECIAL_SIZES: + # Snap onto the exact symmetry position, rotated onto the tabulated + # one. The orbit is the same set either way; the exact + # representative makes it exactly symmetric. + rep = _match_special(special[size], rep) + gen = orbit(rep, tol) + if len(gen) != size: + raise ValueError( + f"orbit of {size} tabulated points generates {len(gen)}") + out.append(ABOrbit(rep, w0, size)) + return out + + +def _match_special(canonical, tabulated): + """The image of `canonical` under the group that is closest to `tabulated`.""" + best, bd = None, None + for R in group(): + q = tuple(mp.fsum(R[i, j] * canonical[j] for j in range(3)) for i in range(3)) + d = sum((q[i] - tabulated[i]) ** 2 for i in range(3)) + if bd is None or d < bd: + best, bd = q, d + return best + + +def n_param(orbits): + return sum(1 + (2 if o.free else 0) for o in orbits) + + +def residual_and_jacobian(orbits, conds, want_jac=True): + G = group() + F = [] + J = matrix(len(conds), n_param(orbits)) if want_jac else None + + # cache orbit images and the tangent images R e1, R e2 for each orbit + cache = [] + for o in orbits: + pts = [tuple(mp.fsum(R[i, j] * o.rep[j] for j in range(3)) for i in range(3)) for R in G] + if o.free: + e1, e2 = o.frame + t1 = [tuple(mp.fsum(R[i, j] * e1[j] for j in range(3)) for i in range(3)) for R in G] + t2 = [tuple(mp.fsum(R[i, j] * e2[j] for j in range(3)) for i in range(3)) for R in G] + else: + t1 = t2 = None + # the group has |G| elements but the orbit may be smaller; scale instead + # of de-duplicating, so every image contributes 1/multiplicity + mult = len(G) // o.size + cache.append((pts, t1, t2, mpf(1) / mult)) + + for ci, (a, b, c) in enumerate(conds): + total = mpf(0) + col = 0 + for o, (pts, t1, t2, scale) in zip(orbits, cache): + s = scale * mp.fsum(p[0] ** a * p[1] ** b * p[2] ** c for p in pts) + total += o.weight * s + if want_jac: + J[ci, col] = s + col += 1 + if o.free: + for t in (t1, t2): + d = mpf(0) + for p, tv in zip(pts, t): + g = _grad(a, b, c, p) + d += g[0] * tv[0] + g[1] * tv[1] + g[2] * tv[2] + J[ci, col] = o.weight * scale * d + col += 1 + else: + col += 1 + (2 if o.free else 0) + F.append(total - exact_monomial(a, b, c)) + return F, J + + +def _pack(orbits): + x = [] + for o in orbits: + x.append(o.weight) + if o.free: + x.extend([mpf(0), mpf(0)]) # tangent steps, always relative + return x + + +def _apply(orbits, d): + k = 0 + for o in orbits: + o.weight = o.weight + d[k]; k += 1 + if o.free: + o.rep = retract(o.rep, d[k], d[k + 1], *o.frame) + o.frame = tangent_frame(o.rep) + k += 2 + + +def refine(orbits, order, target_dps, max_iter=10, verbose=False, conds=None): + conds = conds if conds is not None else conditions(order) + hist = [] + tol = mpf(10) ** (-target_dps) + for it in range(max_iter): + F, J = residual_and_jacobian(orbits, conds) + r = max(fabs(f) for f in F) + hist.append(r) + if verbose: + print(f" iter {it}: max|residual| = {mp.nstr(r, 4)}", flush=True) + if r < tol: + break + n = J.cols + JtJ = matrix(n, n); Jtf = matrix(n, 1) + for i in range(n): + for j in range(i, n): + v = mp.fsum(J[k, i] * J[k, j] for k in range(J.rows)) + JtJ[i, j] = v; JtJ[j, i] = v + Jtf[i] = mp.fsum(J[k, i] * F[k] for k in range(J.rows)) + _apply(orbits, lu_solve(JtJ, -Jtf)) + return hist + + +def materialise(orbits): + pts, wts = [], [] + for o in orbits: + for q in o.points(): + pts.append(q); wts.append(o.weight) + return pts, wts + + +def verify(orbits, order): + pts, wts = materialise(orbits) + worst = mpf(0) + for (a, b, c) in conditions(order): + q = mp.fsum(w * p[0] ** a * p[1] ** b * p[2] ** c for p, w in zip(pts, wts)) + e = exact_monomial(a, b, c) + denom = fabs(e) if e != 0 else mpf(1) + worst = max(worst, fabs(q - e) / denom) + return worst diff --git a/tools/angular_grids/icosahedral.py b/tools/angular_grids/icosahedral.py index 7dc3a56..5680590 100644 --- a/tools/angular_grids/icosahedral.py +++ b/tools/angular_grids/icosahedral.py @@ -105,3 +105,50 @@ def orbit(p, tol=mpf("1e-9")): if not any(sum((q[i] - r[i]) ** 2 for i in range(3)) < tol * tol for r in out): out.append(q) return out + + +def special_positions(tol=mpf("1e-9")): + """Exact representatives of the three special orbits. + + A 12-, 20- or 30-point orbit has no free parameter: its position is fixed + by the symmetry. Taking the representative from tabulated double-precision + data instead makes it a *generic* point whose 60 images cluster in fives + (or threes, or twos) rather than coinciding, and de-duplication then keeps + an arbitrary member of each cluster -- injecting ~1e-16 of asymmetry into + a grid that should be exactly symmetric. + + Returns {12: vertex, 20: face centre, 30: edge midpoint}. + """ + V = _icosahedron_vertices() + + def norm(p): + m = sqrt(sum(x * x for x in p)) + return tuple(x / m for x in p) + + def dist2(a, b): + return sum((a[i] - b[i]) ** 2 for i in range(3)) + + # nearest-neighbour distance defines adjacency + ds = sorted(dist2(V[0], w) for w in V[1:]) + dmin = ds[0] + adj = lambda a, b: dist2(a, b) < dmin * mpf("1.5") + + edge = None + face = None + for i, a in enumerate(V): + for b in V[i + 1:]: + if not adj(a, b): + continue + if edge is None: + edge = norm([a[k] + b[k] for k in range(3)]) + for c in V: + if c is a or c is b or not (adj(a, c) and adj(b, c)): + continue + if face is None: + face = norm([a[k] + b[k] + c[k] for k in range(3)]) + break + if face is not None and edge is not None: + break + if face is not None and edge is not None: + break + return {12: V[0], 20: face, 30: edge} From ebe00ea2b99da80a0b7a3a09a8588708af732f6d Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 21:32:14 +0300 Subject: [PATCH 05/17] tools: make AB refinement tractable at scale Three changes, all needed before the larger Ahrens-Beylkin grids are reachable. Condition selection. The full set is (order+1)^2 monomials -- 1600 at order 39 -- against a few dozen unknowns, since the icosahedral symmetry makes most of them dependent. Candidates are now visited stratified across degree and kept when their Jacobian row increases the column rank, stopping early once the rank is full. AB-192 drops from 576 conditions to 30. Results are always verified against the full set, so a bad selection surfaces as a failed verification rather than a silently wrong grid. Image cache hoisted out of residual_and_jacobian. Selection asks for many single-row Jacobians, and rebuilding the 60 group images per row dominated: AB-192 selection went from 27s to 9s, and its refinement from 45s to 1s. Rank-aware step. The Jacobian is not always full rank -- AB-312 comes out 13 of 16 -- and both a normal-equations LU and mpmath's qr_solve refuse a singular matrix outright. The step is now solved only in an independent set of columns, found by modified Gram-Schmidt with a relative threshold, leaving the unresolved directions at zero. Same lesson as the Womersley path, where CGLS gave the minimum-norm step for free. Also records a negative result: refining AB-552 from its tabulated data diverges, 0.4619 -> 0.03831 -> 145.9. That is expected rather than surprising -- its points are wrong, not merely imprecise, and no weighting of them integrates to order 39 -- but it settles that the corrupt table cannot be recovered by Newton refinement. It needs globalisation. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/ab_refine.py | 151 ++++++++++++++++++++++++++++--- 1 file changed, 136 insertions(+), 15 deletions(-) diff --git a/tools/angular_grids/ab_refine.py b/tools/angular_grids/ab_refine.py index b3188ea..7b79eb2 100644 --- a/tools/angular_grids/ab_refine.py +++ b/tools/angular_grids/ab_refine.py @@ -13,7 +13,7 @@ orbit point with respect to a tangent step of its representative is just the group element applied to that step: d(R p)/dp . e = R e. """ -from mpmath import mp, mpf, fabs, sqrt, matrix, lu_solve +from mpmath import mp, mpf, fabs, lu_solve, matrix, sqrt from .designs import _grad, conditions, exact_monomial, retract, tangent_frame from .icosahedral import group, orbit, special_positions @@ -79,12 +79,14 @@ def n_param(orbits): return sum(1 + (2 if o.free else 0) for o in orbits) -def residual_and_jacobian(orbits, conds, want_jac=True): - G = group() - F = [] - J = matrix(len(conds), n_param(orbits)) if want_jac else None +def image_cache(orbits): + """Group images of each orbit representative, and of its tangent basis. - # cache orbit images and the tangent images R e1, R e2 for each orbit + Rebuilding this per condition dominates the cost when many single-row + Jacobians are wanted (as in select_conditions), so it is computed once and + passed in. + """ + G = group() cache = [] for o in orbits: pts = [tuple(mp.fsum(R[i, j] * o.rep[j] for j in range(3)) for i in range(3)) for R in G] @@ -98,6 +100,13 @@ def residual_and_jacobian(orbits, conds, want_jac=True): # of de-duplicating, so every image contributes 1/multiplicity mult = len(G) // o.size cache.append((pts, t1, t2, mpf(1) / mult)) + return cache + + +def residual_and_jacobian(orbits, conds, want_jac=True, cache=None): + F = [] + J = matrix(len(conds), n_param(orbits)) if want_jac else None + cache = cache if cache is not None else image_cache(orbits) for ci, (a, b, c) in enumerate(conds): total = mpf(0) @@ -146,21 +155,14 @@ def refine(orbits, order, target_dps, max_iter=10, verbose=False, conds=None): hist = [] tol = mpf(10) ** (-target_dps) for it in range(max_iter): - F, J = residual_and_jacobian(orbits, conds) + F, J = residual_and_jacobian(orbits, conds, cache=image_cache(orbits)) r = max(fabs(f) for f in F) hist.append(r) if verbose: print(f" iter {it}: max|residual| = {mp.nstr(r, 4)}", flush=True) if r < tol: break - n = J.cols - JtJ = matrix(n, n); Jtf = matrix(n, 1) - for i in range(n): - for j in range(i, n): - v = mp.fsum(J[k, i] * J[k, j] for k in range(J.rows)) - JtJ[i, j] = v; JtJ[j, i] = v - Jtf[i] = mp.fsum(J[k, i] * F[k] for k in range(J.rows)) - _apply(orbits, lu_solve(JtJ, -Jtf)) + _apply(orbits, _lstsq_step(J, F)) return hist @@ -181,3 +183,122 @@ def verify(orbits, order): denom = fabs(e) if e != 0 else mpf(1) worst = max(worst, fabs(q - e) / denom) return worst + + +def select_conditions(orbits, order, oversample=3, verbose=False): + """Pick a small spanning subset of the exactness conditions. + + The full set is (order+1)^2 monomials -- 1600 at order 39 -- against a few + dozen unknowns, because the icosahedral symmetry makes most of them + dependent. Refining on all of them is pure cost. + + Candidates are visited stratified across degree and kept when their + Jacobian row increases the column rank, stopping once the rank is full and + enough extra rows have been gathered for conditioning. Rows are computed + one at a time and the walk stops early, so the whole selection costs far + less than one full Jacobian. + + The result is always verified against the *full* condition set afterwards, + so a bad selection surfaces as a failed verification rather than as a + silently wrong grid. + """ + import collections + n = n_param(orbits) + by_degree = collections.OrderedDict() + for c in conditions(order): + by_degree.setdefault(sum(c), []).append(c) + + # round-robin across degrees: low degrees first, but never only low degrees + candidates = [] + idx = 0 + while True: + added = False + for deg in sorted(by_degree): + if idx < len(by_degree[deg]): + candidates.append(by_degree[deg][idx]) + added = True + if not added: + break + idx += 1 + + cache = image_cache(orbits) + basis, chosen = [], [] + for cand in candidates: + _, Jrow = residual_and_jacobian(orbits, [cand], cache=cache) + v = [Jrow[0, k] for k in range(n)] + for u in basis: + d = mp.fsum(v[k] * u[k] for k in range(n)) + v = [v[k] - d * u[k] for k in range(n)] + nrm = sqrt(mp.fsum(x * x for x in v)) + if nrm > mpf("1e-14"): + basis.append([x / nrm for x in v]) + chosen.append(cand) + elif len(basis) >= n and len(chosen) < oversample * n: + chosen.append(cand) # extra rows help conditioning + if len(basis) >= n and len(chosen) >= oversample * n: + break + + if verbose: + print(f" selected {len(chosen)} of {len(conditions(order))} conditions, " + f"rank {len(basis)}/{n}", flush=True) + if len(basis) < n: + return conditions(order) # fall back rather than under-determine + return chosen + + +def _rank(J, tol=mpf("1e-18")): + """Column rank by modified Gram-Schmidt; the matrices here are small.""" + cols = [] + for j in range(J.cols): + v = [J[i, j] for i in range(J.rows)] + for u in cols: + d = mp.fsum(v[i] * u[i] for i in range(len(v))) + v = [v[i] - d * u[i] for i in range(len(v))] + nrm = sqrt(mp.fsum(x * x for x in v)) + if nrm > tol: + cols.append([x / nrm for x in v]) + return len(cols) + + +def independent_columns(J, tol=mpf("1e-16")): + """Indices of a maximal independent set of columns, by modified Gram-Schmidt.""" + keep, basis = [], [] + for j in range(J.cols): + v = [J[i, j] for i in range(J.rows)] + nrm0 = sqrt(mp.fsum(x * x for x in v)) + if nrm0 == 0: + continue + for u in basis: + d = mp.fsum(v[i] * u[i] for i in range(len(v))) + v = [v[i] - d * u[i] for i in range(len(v))] + nrm = sqrt(mp.fsum(x * x for x in v)) + if nrm > tol * nrm0: # relative, so column scale does not matter + basis.append([x / nrm for x in v]) + keep.append(j) + return keep + + +def _lstsq_step(J, F): + """Gauss-Newton step, solved only in the directions the Jacobian resolves. + + The Jacobian is not always full rank -- AB-312 comes out 13 of 16 -- and + both a normal-equations LU and mpmath's qr_solve simply refuse a singular + matrix. Restricting to an independent set of columns and leaving the rest + at zero gives the step the well-determined directions deserve and declines + to move along the ones the conditions cannot see. + """ + keep = independent_columns(J) + m = len(keep) + JtJ = matrix(m, m); Jtf = matrix(m, 1) + for a in range(m): + ja = keep[a] + for b in range(a, m): + jb = keep[b] + v = mp.fsum(J[k, ja] * J[k, jb] for k in range(J.rows)) + JtJ[a, b] = v; JtJ[b, a] = v + Jtf[a] = mp.fsum(J[k, ja] * F[k] for k in range(J.rows)) + red = lu_solve(JtJ, -Jtf) + step = [mpf(0)] * J.cols + for a, j in enumerate(keep): + step[j] = red[a] + return step From ce03170c35963ddd881e10c0fc911f005af066ab Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 21:36:24 +0300 Subject: [PATCH 06/17] tools: one CLI for all four grid families generate.py now dispatches on the family: octahedral refinement for Lebedev-Laikov and Delley, icosahedral for Ahrens-Beylkin, and the equal-weight design path for Womersley. Each verifies against the full condition set before emitting, and refuses to write a table that did not reach the requested precision. ahrens_beylkin_72: order 14, 2 icosahedral orbits, 4 free parameters selected 12 of 225 conditions, rank 4/4 1.754e-15 -> 2.091e-33 verified 2.431e-32 against the full condition set The emitted load() is now templated on its containers, so it fills the library's std::vector> and std::vector directly -- the same role detail::copy_grid plays for the literal tables, which is what the eventual switch needs. Verified by filling those vectors from a generated header. Two smaller fixes: the emit namespace and title tables only knew the two octahedral families, so an Ahrens-Beylkin header came out in namespace `ahrens_beylkin` rather than `AhrensBeylkinGrids`; and condition selection now caps its walk, since a genuinely rank-deficient Jacobian means the rank never fills and the search would otherwise visit all 1600 candidates. An incomplete basis is harmless -- the step solver handles deficiency, and verification against the full set is the real guard. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/ab_refine.py | 95 +++++++++++++++++++++++++++++++- tools/angular_grids/emit.py | 14 +++-- tools/angular_grids/generate.py | 63 +++++++++++++++++---- 3 files changed, 153 insertions(+), 19 deletions(-) diff --git a/tools/angular_grids/ab_refine.py b/tools/angular_grids/ab_refine.py index 7b79eb2..31d81ef 100644 --- a/tools/angular_grids/ab_refine.py +++ b/tools/angular_grids/ab_refine.py @@ -222,8 +222,14 @@ def select_conditions(orbits, order, oversample=3, verbose=False): idx += 1 cache = image_cache(orbits) + # Cap the walk: when the Jacobian is genuinely rank-deficient the rank + # never fills, and without a cap this visits every candidate -- 1600 of + # them at order 39. The step solver handles a deficient Jacobian, so an + # incomplete basis is fine; verification against the full set is what + # actually guards correctness. + limit = max(oversample * n, 12 * n) basis, chosen = [], [] - for cand in candidates: + for cand in candidates[:limit]: _, Jrow = residual_and_jacobian(orbits, [cand], cache=cache) v = [Jrow[0, k] for k in range(n)] for u in basis: @@ -241,8 +247,8 @@ def select_conditions(orbits, order, oversample=3, verbose=False): if verbose: print(f" selected {len(chosen)} of {len(conditions(order))} conditions, " f"rank {len(basis)}/{n}", flush=True) - if len(basis) < n: - return conditions(order) # fall back rather than under-determine + if not chosen: + return conditions(order) return chosen @@ -302,3 +308,86 @@ def _lstsq_step(J, F): for a, j in enumerate(keep): step[j] = red[a] return step + + +def _state(orbits): + return [(list(o.rep), o.weight) for o in orbits] + + +def _restore(orbits, st): + for o, (rep, w) in zip(orbits, st): + o.rep = list(rep) + o.weight = w + o.frame = tangent_frame(o.rep) + + +def _sumsq(F): + return mp.fsum(f * f for f in F) + + +def refine_lm(orbits, order, target_dps, conds=None, max_iter=200, + lam0=mpf("1e-3"), verbose=False): + """Levenberg-Marquardt: Gauss-Newton with adaptive damping. + + Plain Gauss-Newton diverges on a starting point outside the basin -- AB-552 + goes 0.4619 -> 0.03831 -> 145.9 -- because nothing constrains the step + length. Damping the normal equations by lam * diag(J^T J) interpolates + between the Gauss-Newton step (small lam) and a short steepest-descent step + (large lam), and lam is adapted on whether the step actually reduced the + residual. + + This is the cheapest form of globalisation that fits the existing + machinery. It is not a trust-region method and makes no claim to find a + rule that is not near the starting point. + """ + conds = conds if conds is not None else conditions(order) + tol = mpf(10) ** (-target_dps) + lam = lam0 + F, J = residual_and_jacobian(orbits, conds, cache=image_cache(orbits)) + cost = _sumsq(F) + hist = [max(fabs(f) for f in F)] + if verbose: + print(f" start: max|F| = {mp.nstr(hist[0], 4)}", flush=True) + + for it in range(max_iter): + if max(fabs(f) for f in F) < tol: + break + n = J.cols + JtJ = matrix(n, n); Jtf = matrix(n, 1) + for i in range(n): + for j in range(i, n): + v = mp.fsum(J[k, i] * J[k, j] for k in range(J.rows)) + JtJ[i, j] = v; JtJ[j, i] = v + Jtf[i] = mp.fsum(J[k, i] * F[k] for k in range(J.rows)) + + st = _state(orbits) + accepted = False + for _ in range(30): # inner damping search + A = matrix(n, n) + for i in range(n): + for j in range(n): + A[i, j] = JtJ[i, j] + A[i, i] = JtJ[i, i] * (1 + lam) + lam * mpf("1e-30") + try: + d = lu_solve(A, -Jtf) + except Exception: + lam *= 10 + continue + _apply(orbits, [d[k] for k in range(n)]) + Fn, Jn = residual_and_jacobian(orbits, conds, cache=image_cache(orbits)) + if _sumsq(Fn) < cost: + F, J, cost = Fn, Jn, _sumsq(Fn) + lam = max(lam / 10, mpf("1e-25")) + accepted = True + break + _restore(orbits, st) + lam *= 10 + if lam > mpf("1e20"): + break + hist.append(max(fabs(f) for f in F)) + if verbose and (it < 5 or it % 20 == 0): + print(f" iter {it}: max|F| = {mp.nstr(hist[-1], 4)} lam = {mp.nstr(lam, 3)}", + flush=True) + if not accepted: + break + return hist diff --git a/tools/angular_grids/emit.py b/tools/angular_grids/emit.py index 7a0b8a9..947380f 100644 --- a/tools/angular_grids/emit.py +++ b/tools/angular_grids/emit.py @@ -57,8 +57,12 @@ }}; /// Materialise the grid at the working precision of T. - static void load(std::array, {npts}>& p, - std::array& w) {{ + /// + /// Templated on the containers so it can fill the library's + /// std::vector> / std::vector directly, in the same + /// role detail::copy_grid played for the literal tables. + template + static void load(PointContainer& p, WeightContainer& w) {{ for(size_t i = 0; i < {npts}; ++i) {{ for(size_t j = 0; j < 3; ++j) p[i][j] = detail::grid_scalar::parse(point_digits[3 * i + j]); @@ -71,8 +75,10 @@ }} // namespace IntegratorXX """ -NAMESPACES = {"lebedev_laikov": "LebedevLaikovGrids", "delley": "DelleyGrids"} -TITLES = {"lebedev_laikov": "Lebedev-Laikov", "delley": "Delley"} +NAMESPACES = {"lebedev_laikov": "LebedevLaikovGrids", "delley": "DelleyGrids", + "ahrens_beylkin": "AhrensBeylkinGrids", "womersley": "WomersleyGrids"} +TITLES = {"lebedev_laikov": "Lebedev-Laikov", "delley": "Delley", + "ahrens_beylkin": "Ahrens-Beylkin", "womersley": "Womersley"} def _fmt(values, digits, per_line): diff --git a/tools/angular_grids/generate.py b/tools/angular_grids/generate.py index c2f4b31..3a63fab 100644 --- a/tools/angular_grids/generate.py +++ b/tools/angular_grids/generate.py @@ -29,27 +29,65 @@ def algebraic_order(family, npts, root): return int(m.group(1)) -def run(family, npts, digits, root, out, verbose=True): - mp.dps = digits + 20 # guard digits for the refinement - order = algebraic_order(family, npts, root) +OCTAHEDRAL = ("lebedev_laikov", "delley") + + +def _run_octahedral(family, npts, order, digits, root, verbose): points, weights = read_grid(grid_path(root, family, npts)) - scale = 4 * mp.pi / sum(weights) # tables are normalised to 1 + scale = 4 * mp.pi / sum(weights) # these tables are normalised to 1 orbits = decompose(points, [w * scale for w in weights]) n_param = sum(1 + o.type.n_param for o in orbits) - if verbose: - print(f"{family}_{npts}: order {order}, {len(orbits)} orbits, {n_param} free parameters") - hist = refine(orbits, order, digits + 5, verbose=verbose) + print(f"{family}_{npts}: order {order}, {len(orbits)} octahedral orbits, " + f"{n_param} free parameters") + refine(orbits, order, digits + 5, verbose=verbose) P, W = materialise(orbits) + return P, W, verify(P, W, order) + + +def _run_icosahedral(family, npts, order, digits, root, verbose): + from . import ab_refine as ab + points, weights = read_grid(grid_path(root, family, npts)) + orbits = ab.decompose(points, weights) # AB tables already sum to 4 pi + if verbose: + print(f"{family}_{npts}: order {order}, {len(orbits)} icosahedral orbits, " + f"{ab.n_param(orbits)} free parameters") + conds = ab.select_conditions(orbits, order, verbose=verbose) + ab.refine(orbits, order, digits + 5, conds=conds, verbose=verbose) + P, W = ab.materialise(orbits) + return P, W, ab.verify(orbits, order) + + +def _run_design(family, npts, order, digits, root, verbose): + from . import designs + points, _ = read_grid(grid_path(root, family, npts)) + d = designs.Design(points, order) + if verbose: + print(f"{family}_{npts}: order {order}, equal-weight design, " + f"{2 * npts} tangential unknowns") + designs.refine(d, digits + 5, verbose=verbose) + return d.points, [d.weight] * d.n, designs.verify(d) + + +def run(family, npts, digits, root, out, verbose=True): + mp.dps = digits + 20 # guard digits for the refinement + order = algebraic_order(family, npts, root) + if family in OCTAHEDRAL: + P, W, err = _run_octahedral(family, npts, order, digits, root, verbose) + elif family == "ahrens_beylkin": + P, W, err = _run_icosahedral(family, npts, order, digits, root, verbose) + elif family == "womersley": + P, W, err = _run_design(family, npts, order, digits, root, verbose) + else: + raise SystemExit(f"unknown family {family}") + if len(P) != npts: raise SystemExit(f"materialised {len(P)} points, expected {npts}") - - err = verify(P, W, order) # checked against the full set if err > mp.mpf(10) ** (-digits): raise SystemExit(f"refinement reached only {mp.nstr(err, 4)}, short of {digits} digits") if verbose: - print(f" verified worst relative error {mp.nstr(err, 4)} over all even " - f"monomials up to degree {order}") + print(f" verified worst error {mp.nstr(err, 4)} against the full " + f"condition set for order {order}") text = emit(family, npts, order, P, W, digits, err) if out: @@ -63,7 +101,8 @@ def run(family, npts, digits, root, out, verbose=True): def main(argv=None): ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) - ap.add_argument("--family", required=True, choices=["lebedev_laikov", "delley"]) + ap.add_argument("--family", required=True, + choices=["lebedev_laikov", "delley", "ahrens_beylkin", "womersley"]) ap.add_argument("--npts", required=True, type=int) ap.add_argument("--digits", type=int, default=40) ap.add_argument("--root", default=".") From d8c64bed92e9db0d1fe889cb7a6b6352a7c4bb23 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 21:39:18 +0300 Subject: [PATCH 07/17] tools: batch regeneration and the copy_grid -> load switch batch.py regenerates every tabulated size of a family, reading the sizes from the family's own algebraic_order_by_npts table so the two cannot drift apart. Sizes are independent, so they run in separate processes. A size that fails to reach the requested precision is reported and skipped rather than written, so a partial run never emits an unverified table. delley: 6 sizes, 30 digits, 4 jobs ok delley_14 err 2.552e-51 ok delley_26 err 4.786e-51 ok delley_50 err 5.264e-51 ok delley_110 err 1.246e-49 ok delley_194 err 6.018e-49 ok delley_302 err 2.163e-48 switch_to_load.py points a family's dispatch header at the regenerated tables: detail::copy_grid>(points, weights) -> delley_50::load(points, weights) and drops the `weights[i] *= 4.0*M_PI` loop for the two families whose shipped tables are normalised to one -- a regenerated table already carries the 4*pi convention. The loop index is `auto` in delley.hpp and `size_t` in lebedev_laikov.hpp, so the pattern does not pin the type. Dry runs cover all four families: 31 calls in lebedev_laikov, 18 in delley, 56 in ahrens_beylkin, 125 in womersley. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/batch.py | 72 +++++++++++++++++++++++++++ tools/angular_grids/switch_to_load.py | 68 +++++++++++++++++++++++++ 2 files changed, 140 insertions(+) create mode 100644 tools/angular_grids/batch.py create mode 100644 tools/angular_grids/switch_to_load.py diff --git a/tools/angular_grids/batch.py b/tools/angular_grids/batch.py new file mode 100644 index 0000000..c1627f2 --- /dev/null +++ b/tools/angular_grids/batch.py @@ -0,0 +1,72 @@ +"""Regenerate every tabulated size of a family. + + PYTHONPATH=tools python3 -m angular_grids.batch --family delley --digits 40 \ + --out gen/ --jobs 8 + +Each size is independent, so they run in separate processes. A size that +fails to reach the requested precision is reported and skipped rather than +written, so a partial run never emits a table that was not verified. +""" +import argparse +import multiprocessing as mp_proc +import re +import sys +import time +from pathlib import Path + + +def tabulated_sizes(family, root): + """Sizes the family's dispatch header advertises, via its order table.""" + text = (Path(root) / "include/integratorxx/quadratures/s2" / f"{family}.hpp").read_text() + body = text[text.index("algebraic_order_by_npts"):] + body = body[:body.index("next_algebraic_order")] if "next_algebraic_order" in body else body + return sorted({int(m) for m in re.findall(r"case\s+(\d+)\s*:", body)}) + + +def _one(args): + family, npts, digits, root, out = args + from .generate import run + t = time.time() + try: + _, err = run(family, npts, digits, root, out, verbose=False) + from mpmath import mp + return (npts, True, mp.nstr(err, 4), time.time() - t) + except SystemExit as e: + return (npts, False, str(e), time.time() - t) + except Exception as e: # noqa: BLE001 - report, do not abort the batch + return (npts, False, f"{type(e).__name__}: {e}", time.time() - t) + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument("--family", required=True, + choices=["lebedev_laikov", "delley", "ahrens_beylkin", "womersley"]) + ap.add_argument("--digits", type=int, default=40) + ap.add_argument("--root", default=".") + ap.add_argument("--out", required=True) + ap.add_argument("--jobs", type=int, default=max(1, (mp_proc.cpu_count() or 2) - 1)) + ap.add_argument("--max-npts", type=int, default=None, + help="skip sizes above this, for a quick pass") + a = ap.parse_args(argv) + + sizes = tabulated_sizes(a.family, a.root) + if a.max_npts: + sizes = [n for n in sizes if n <= a.max_npts] + print(f"{a.family}: {len(sizes)} sizes, {a.digits} digits, {a.jobs} jobs", flush=True) + + work = [(a.family, n, a.digits, a.root, a.out) for n in sizes] + ok = bad = 0 + with mp_proc.Pool(a.jobs) as pool: + for npts, good, msg, dt in pool.imap_unordered(_one, work): + if good: + ok += 1 + print(f" ok {a.family}_{npts:<6d} err {msg:<12s} [{dt:5.0f}s]", flush=True) + else: + bad += 1 + print(f" FAIL {a.family}_{npts:<6d} {msg} [{dt:5.0f}s]", flush=True) + print(f"\n{ok} regenerated, {bad} failed", flush=True) + return 1 if bad else 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tools/angular_grids/switch_to_load.py b/tools/angular_grids/switch_to_load.py new file mode 100644 index 0000000..68981e6 --- /dev/null +++ b/tools/angular_grids/switch_to_load.py @@ -0,0 +1,68 @@ +"""Point a family's dispatch header at regenerated string tables. + + PYTHONPATH=tools python3 -m angular_grids.switch_to_load --family delley --root . + +Rewrites + + detail::copy_grid>(points, weights); + -> delley_50::load(points, weights); + +and drops the `weights[i] *= 4.0*M_PI` normalisation loop, because a +regenerated table already carries the 4*pi convention while the shipped +literal tables are normalised to one. + +This only makes sense once the data headers have actually been regenerated; +run angular_grids.batch first. The change is textual and reversible with git. +""" +import argparse +import re +import sys +from pathlib import Path + +FAMILIES = ("lebedev_laikov", "delley", "ahrens_beylkin", "womersley") + +# families whose shipped tables are normalised to 1 and scaled at load +SCALED = ("lebedev_laikov", "delley") + +# The loop index is declared `auto` in delley.hpp and `size_t` in +# lebedev_laikov.hpp, so do not pin the type. +_SCALE_LOOP = re.compile( + r"\n[ \t]*//[^\n]*4 ?pi[^\n]*\n" + r"[ \t]*for\s*\(\s*\w+\s+i\s*=\s*0\s*;\s*i\s*<\s*npts\s*;\s*i\+\+\s*\)\s*\n" + r"[ \t]*weights\[i\]\s*\*=\s*4\.0\s*\*\s*M_PI\s*;\n", + re.IGNORECASE) + + +def switch(family, root, dry_run=False): + path = Path(root) / "include/integratorxx/quadratures/s2" / f"{family}.hpp" + text = path.read_text() + + calls = re.findall(r"detail::copy_grid<\s*(\w+)\s*>\s*\(\s*points\s*,\s*weights\s*\)", + text) + new = re.sub(r"detail::copy_grid<\s*(\w+)\s*>\s*\(\s*points\s*,\s*weights\s*\)", + r"\1::load(points, weights)", text) + + dropped = 0 + if family in SCALED: + new, dropped = _SCALE_LOOP.subn("\n", new) + + if not dry_run: + path.write_text(new) + return len(calls), dropped, path + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument("--family", required=True, choices=FAMILIES) + ap.add_argument("--root", default=".") + ap.add_argument("--dry-run", action="store_true") + a = ap.parse_args(argv) + n, dropped, path = switch(a.family, a.root, a.dry_run) + what = "would rewrite" if a.dry_run else "rewrote" + print(f"{what} {n} dispatch calls in {path}" + + (f", dropped {dropped} 4*pi scaling loop(s)" if dropped else "")) + return 0 if n else 1 + + +if __name__ == "__main__": + sys.exit(main()) From 3c0b73afc9df8e589455369517981bbf92ba3e63 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 21:47:10 +0300 Subject: [PATCH 08/17] tools: document the AB-552 findings and the regeneration workflow Records what is known about the one table that cannot be regenerated by refinement, with AB-612 as a control throughout. The decisive result is that the AB-552 points admit no valid weighting. The exactness conditions are linear in the weights, so the best possible weights for a fixed point set follow from one least-squares solve: AB-612 (known good) 8.3e-14, all positive AB-552 2.4e-02, five negative weights The positions are wrong, not just the weights, and not merely imprecise. Refinement behaves accordingly: Gauss-Newton diverges (0.4619 -> 0.03831 -> 145.9), and Levenberg-Marquardt globalises the step and makes real progress (0.4619 -> 9.99e-4 -> 7.75e-6 -> 1.19e-7) but plateaus five orders short, on a condition subset that does not determine all 28 parameters. Also records what is *not* evidence, since I nearly reported it as such: the Jacobian's rank over low-degree conditions is 6 against 28 parameters, which looks alarming until AB-612 is measured the same way and gives rank 6 against 31. Few icosahedral invariants exist at low degree, so that number says nothing about either grid. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/README.md | 52 +++++++++++++++++++++++++++++++++-- 1 file changed, 50 insertions(+), 2 deletions(-) diff --git a/tools/angular_grids/README.md b/tools/angular_grids/README.md index 5b94370..0d89a8a 100644 --- a/tools/angular_grids/README.md +++ b/tools/angular_grids/README.md @@ -71,8 +71,8 @@ which works for the built-in types as well as Boost.Multiprecision and MPFR. * **Ahrens-Beylkin** has icosahedral rather than octahedral symmetry (15012 points / 60 = 250.2 against 251 distinct weights, the rotation group without inversion). Same method, different orbit algebra. -* **Constructing new rules**, including regenerating the corrupt 552-point - Ahrens-Beylkin table. That is the trust-region half of the problem. +* **Constructing new rules.** That is the trust-region half of the problem; + see below for what is known about the one table that needs it. ## Equal-weight spherical designs (Womersley) @@ -145,3 +145,51 @@ monomial set rather than the reduced one. At 7939 points a single `J*v` product is around 1.3e8 high-precision operations and the iteration count grows too, so the largest designs want a compiled MPFR kernel rather than this driver. + + +## The 552-point Ahrens-Beylkin table + +This table is not a valid order-39 rule and cannot be recovered by refinement. +The evidence, with AB-612 as a control throughout: + +**Its points admit no valid weighting.** The exactness conditions are *linear* +in the weights, so the best possible weights for a fixed set of points follow +from one least-squares solve: + +| grid | best-possible weights, full order-39/41 conditions | +|---------------------|----------------------------------------------------| +| AB-612 (known good) | 8.3e-14, all positive | +| AB-552 | **2.4e-02, five negative weights** | + +So the point *positions* are wrong, not just the weights, and not merely +imprecise. This is the decisive result. + +**Refinement diverges.** Gauss-Newton from the tabulated data goes +0.4619 -> 0.03831 -> 145.9. Levenberg-Marquardt globalises the step and does +make progress -- 0.4619 -> 9.99e-4 -> 7.75e-6 -> 1.19e-7 -- but plateaus there, +five orders short, and on a condition subset that does not determine all 28 +parameters. + +**What is not evidence.** The Jacobian's rank over low-degree conditions is 6 +against 28 parameters, which looks alarming until AB-612 is measured the same +way and gives rank 6 against 31. Few icosahedral invariants exist at low +degree; this says nothing about either grid. + +The remaining routes are to obtain the original Ahrens-Beylkin data +(Proc. R. Soc. A 465, 3103) and refine that, or to construct the rule from +scratch with a trust-region method. Of the two the first is much cheaper, and +the other 55 sizes in the family are unaffected -- they all refine cleanly. + +## Regenerating everything + +```sh +PYTHONPATH=tools python3 -m angular_grids.batch \ + --family delley --digits 40 --out gen/ --jobs 8 +PYTHONPATH=tools python3 -m angular_grids.switch_to_load --family delley +``` + +`batch` reads the sizes from the family's own `algebraic_order_by_npts` table, +so the two cannot drift apart, and skips rather than writes any size that +misses the requested precision. `switch_to_load` then repoints the dispatch at +the regenerated tables and drops the `4*M_PI` scaling loop for the families +whose shipped tables are normalised to one. From ba74aba7940906a97dd9f504945186efe7bb99c9 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 21:55:09 +0300 Subject: [PATCH 09/17] tools: refuse to generate into the source directory Pointing --out at the family's own data directory destroys the input. Generation reads the shipped literal table for each size, so a partial run replaces the tables the remaining sizes still need, and every later size then fails with "substring not found" -- which is what happened, on a scratch worktree, when I did exactly this. Generate into a separate directory and install the results deliberately. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/generate.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/tools/angular_grids/generate.py b/tools/angular_grids/generate.py index 3a63fab..048258d 100644 --- a/tools/angular_grids/generate.py +++ b/tools/angular_grids/generate.py @@ -91,7 +91,17 @@ def run(family, npts, digits, root, out, verbose=True): text = emit(family, npts, order, P, W, digits, err) if out: - path = Path(out) / f"{family}_{npts}.hpp" + out = Path(out) + # Writing into the source directory destroys the input: a partial run + # replaces the literal tables the remaining sizes still need to read, + # and every later size then fails to parse. Generate elsewhere and + # install deliberately. + src = (Path(root) / "include/integratorxx/quadratures/s2" / family).resolve() + if out.resolve() == src: + raise SystemExit( + f"refusing to write into the source directory {src}; " + "generate into a separate directory and copy the results in") + path = out / f"{family}_{npts}.hpp" path.parent.mkdir(parents=True, exist_ok=True) path.write_text(text) if verbose: From f2ca5df98a6981a3d40cb2f41f518e4f2e7e496f Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 21:58:51 +0300 Subject: [PATCH 10/17] tools: share the rank-aware solve with the octahedral path The octahedral refinement still solved its normal equations with a plain lu_solve, and mpmath refuses a singular matrix outright. That surfaced as delley_1454 failing with "matrix is numerically singular" at 20 digits while succeeding at 40 -- so it is a robustness property of the solver rather than of any particular grid, and it would have bitten whichever size happened to be marginal at the chosen precision. Move independent_columns() and the step solver into linalg.py and use them from both paths. delley_1454 now regenerates at 20 digits, verified to 1.287e-24 against the full condition set for order 65. Regression suite 11/11. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/ab_refine.py | 130 +------------------------------ tools/angular_grids/linalg.py | 58 ++++++++++++++ tools/angular_grids/refine.py | 12 +-- 3 files changed, 63 insertions(+), 137 deletions(-) create mode 100644 tools/angular_grids/linalg.py diff --git a/tools/angular_grids/ab_refine.py b/tools/angular_grids/ab_refine.py index 31d81ef..b293dc0 100644 --- a/tools/angular_grids/ab_refine.py +++ b/tools/angular_grids/ab_refine.py @@ -16,6 +16,7 @@ from mpmath import mp, mpf, fabs, lu_solve, matrix, sqrt from .designs import _grad, conditions, exact_monomial, retract, tangent_frame +from .linalg import independent_columns, lstsq_step from .icosahedral import group, orbit, special_positions SPECIAL_SIZES = (12, 20, 30) # orbits with no free angular parameter @@ -162,7 +163,7 @@ def refine(orbits, order, target_dps, max_iter=10, verbose=False, conds=None): print(f" iter {it}: max|residual| = {mp.nstr(r, 4)}", flush=True) if r < tol: break - _apply(orbits, _lstsq_step(J, F)) + _apply(orbits, lstsq_step(J, F)) return hist @@ -264,130 +265,3 @@ def _rank(J, tol=mpf("1e-18")): if nrm > tol: cols.append([x / nrm for x in v]) return len(cols) - - -def independent_columns(J, tol=mpf("1e-16")): - """Indices of a maximal independent set of columns, by modified Gram-Schmidt.""" - keep, basis = [], [] - for j in range(J.cols): - v = [J[i, j] for i in range(J.rows)] - nrm0 = sqrt(mp.fsum(x * x for x in v)) - if nrm0 == 0: - continue - for u in basis: - d = mp.fsum(v[i] * u[i] for i in range(len(v))) - v = [v[i] - d * u[i] for i in range(len(v))] - nrm = sqrt(mp.fsum(x * x for x in v)) - if nrm > tol * nrm0: # relative, so column scale does not matter - basis.append([x / nrm for x in v]) - keep.append(j) - return keep - - -def _lstsq_step(J, F): - """Gauss-Newton step, solved only in the directions the Jacobian resolves. - - The Jacobian is not always full rank -- AB-312 comes out 13 of 16 -- and - both a normal-equations LU and mpmath's qr_solve simply refuse a singular - matrix. Restricting to an independent set of columns and leaving the rest - at zero gives the step the well-determined directions deserve and declines - to move along the ones the conditions cannot see. - """ - keep = independent_columns(J) - m = len(keep) - JtJ = matrix(m, m); Jtf = matrix(m, 1) - for a in range(m): - ja = keep[a] - for b in range(a, m): - jb = keep[b] - v = mp.fsum(J[k, ja] * J[k, jb] for k in range(J.rows)) - JtJ[a, b] = v; JtJ[b, a] = v - Jtf[a] = mp.fsum(J[k, ja] * F[k] for k in range(J.rows)) - red = lu_solve(JtJ, -Jtf) - step = [mpf(0)] * J.cols - for a, j in enumerate(keep): - step[j] = red[a] - return step - - -def _state(orbits): - return [(list(o.rep), o.weight) for o in orbits] - - -def _restore(orbits, st): - for o, (rep, w) in zip(orbits, st): - o.rep = list(rep) - o.weight = w - o.frame = tangent_frame(o.rep) - - -def _sumsq(F): - return mp.fsum(f * f for f in F) - - -def refine_lm(orbits, order, target_dps, conds=None, max_iter=200, - lam0=mpf("1e-3"), verbose=False): - """Levenberg-Marquardt: Gauss-Newton with adaptive damping. - - Plain Gauss-Newton diverges on a starting point outside the basin -- AB-552 - goes 0.4619 -> 0.03831 -> 145.9 -- because nothing constrains the step - length. Damping the normal equations by lam * diag(J^T J) interpolates - between the Gauss-Newton step (small lam) and a short steepest-descent step - (large lam), and lam is adapted on whether the step actually reduced the - residual. - - This is the cheapest form of globalisation that fits the existing - machinery. It is not a trust-region method and makes no claim to find a - rule that is not near the starting point. - """ - conds = conds if conds is not None else conditions(order) - tol = mpf(10) ** (-target_dps) - lam = lam0 - F, J = residual_and_jacobian(orbits, conds, cache=image_cache(orbits)) - cost = _sumsq(F) - hist = [max(fabs(f) for f in F)] - if verbose: - print(f" start: max|F| = {mp.nstr(hist[0], 4)}", flush=True) - - for it in range(max_iter): - if max(fabs(f) for f in F) < tol: - break - n = J.cols - JtJ = matrix(n, n); Jtf = matrix(n, 1) - for i in range(n): - for j in range(i, n): - v = mp.fsum(J[k, i] * J[k, j] for k in range(J.rows)) - JtJ[i, j] = v; JtJ[j, i] = v - Jtf[i] = mp.fsum(J[k, i] * F[k] for k in range(J.rows)) - - st = _state(orbits) - accepted = False - for _ in range(30): # inner damping search - A = matrix(n, n) - for i in range(n): - for j in range(n): - A[i, j] = JtJ[i, j] - A[i, i] = JtJ[i, i] * (1 + lam) + lam * mpf("1e-30") - try: - d = lu_solve(A, -Jtf) - except Exception: - lam *= 10 - continue - _apply(orbits, [d[k] for k in range(n)]) - Fn, Jn = residual_and_jacobian(orbits, conds, cache=image_cache(orbits)) - if _sumsq(Fn) < cost: - F, J, cost = Fn, Jn, _sumsq(Fn) - lam = max(lam / 10, mpf("1e-25")) - accepted = True - break - _restore(orbits, st) - lam *= 10 - if lam > mpf("1e20"): - break - hist.append(max(fabs(f) for f in F)) - if verbose and (it < 5 or it % 20 == 0): - print(f" iter {it}: max|F| = {mp.nstr(hist[-1], 4)} lam = {mp.nstr(lam, 3)}", - flush=True) - if not accepted: - break - return hist diff --git a/tools/angular_grids/linalg.py b/tools/angular_grids/linalg.py new file mode 100644 index 0000000..f0eb483 --- /dev/null +++ b/tools/angular_grids/linalg.py @@ -0,0 +1,58 @@ +"""Small dense solves shared by the refinement paths. + +The Gauss-Newton systems here are tiny but not always full rank, and mpmath's +lu_solve and qr_solve both refuse a singular matrix outright rather than +returning something usable. Working precision changes whether a given system +looks singular -- delley_1454 solves cleanly at 40 digits and fails at 20 -- +so this is a robustness property rather than a property of any one grid. +""" +from mpmath import lu_solve, matrix, mp, mpf, sqrt + + +def independent_columns(J, tol=mpf("1e-16")): + """Indices of a maximal independent set of columns, by modified Gram-Schmidt. + + The threshold is relative to each column's own norm, so columns that differ + wildly in scale (a weight column against an angular one) are judged fairly. + """ + keep, basis = [], [] + for j in range(J.cols): + v = [J[i, j] for i in range(J.rows)] + nrm0 = sqrt(mp.fsum(x * x for x in v)) + if nrm0 == 0: + continue + for u in basis: + d = mp.fsum(v[i] * u[i] for i in range(len(v))) + v = [v[i] - d * u[i] for i in range(len(v))] + nrm = sqrt(mp.fsum(x * x for x in v)) + if nrm > tol * nrm0: + basis.append([x / nrm for x in v]) + keep.append(j) + return keep + + +def lstsq_step(J, F): + """Gauss-Newton step, solved only in the directions the Jacobian resolves. + + Directions outside the resolved set are left at zero: the conditions cannot + see them, so moving along them is guesswork. + """ + keep = independent_columns(J) + m = len(keep) + if m == 0: + return [mpf(0)] * J.cols + JtJ = matrix(m, m) + Jtf = matrix(m, 1) + for a in range(m): + ja = keep[a] + for b in range(a, m): + jb = keep[b] + v = mp.fsum(J[k, ja] * J[k, jb] for k in range(J.rows)) + JtJ[a, b] = v + JtJ[b, a] = v + Jtf[a] = mp.fsum(J[k, ja] * F[k] for k in range(J.rows)) + red = lu_solve(JtJ, -Jtf) + step = [mpf(0)] * J.cols + for a, j in enumerate(keep): + step[j] = red[a] + return step diff --git a/tools/angular_grids/refine.py b/tools/angular_grids/refine.py index 2d6e4e9..2995478 100644 --- a/tools/angular_grids/refine.py +++ b/tools/angular_grids/refine.py @@ -8,8 +8,9 @@ """ import collections -from mpmath import mp, mpf, fabs, lu_solve, matrix +from mpmath import mp, mpf, fabs, matrix +from .linalg import lstsq_step from .moments import conditions, exact_moment, orbit_moment from .orbits import classify, signed_permutations @@ -93,14 +94,7 @@ def refine(orbits, order, target_dps, max_iter=12, verbose=False): print(f" iter {it}: max|residual| = {mp.nstr(r, 4)}") if r < tol: break - n = J.cols - JtJ = matrix(n, n); Jtf = matrix(n, 1) - for i in range(n): - for j in range(i, n): - v = mp.fsum(J[k, i] * J[k, j] for k in range(J.rows)) - JtJ[i, j] = v; JtJ[j, i] = v - Jtf[i] = mp.fsum(J[k, i] * F[k] for k in range(J.rows)) - d = lu_solve(JtJ, -Jtf) + d = lstsq_step(J, F) _unpack(orbits, [x + d[j] for j, x in enumerate(_pack(orbits))]) return hist From 8e8090c2d662380bc7266031dc115ee9a23bd716 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 22:13:11 +0300 Subject: [PATCH 11/17] tools: solve the Gauss-Newton step by QR, not normal equations The rank-aware step I added earlier used a fixed 1e-16 threshold on the normal equations, and that discarded columns which were merely ill-conditioned rather than null. delley_974 lost its convergence as a result -- 3.11e-29 instead of 6.93e-51, failing the generator's own precision check. The batch caught it; the unit tests did not. Chasing the threshold does not work, because the two failure modes trade off against each other: threshold delley_974 @ 40 delley_1454 @ 20 fixed 1e-16 stalls 3.1e-29 solves precision-tied 6.93e-51 "numerically singular" QR, 1e-14 cut 11 iters, 2e-33 diverges, 1.7e+146 QR, precision-tied 6.93e-51, 2 iters needs more precision The tension is not in the threshold. It is that J^T J squares the condition number. The step is now a rank-revealing modified Gram-Schmidt QR of J itself, back-substituted, never forming the normal equations. Checked against a synthetic rank-deficient system where it finds rank 2 and leaves a 1e-31 residual. The threshold that remains cuts only what is unresolvable at the working precision, so it discards genuinely null directions and keeps everything else. delley_1454 at 20 digits is then simply a case that needs more precision than it was given -- 102 parameters at that conditioning are not resolvable in 20 digits by any of these methods, and it refines cleanly at 40. The generator's final verification catches it and refuses to write the table, which is the right outcome. Regression suite 11/11. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/linalg.py | 109 +++++++++++++++++++++++----------- 1 file changed, 73 insertions(+), 36 deletions(-) diff --git a/tools/angular_grids/linalg.py b/tools/angular_grids/linalg.py index f0eb483..12fc117 100644 --- a/tools/angular_grids/linalg.py +++ b/tools/angular_grids/linalg.py @@ -1,58 +1,95 @@ -"""Small dense solves shared by the refinement paths. +"""Small dense least squares shared by the refinement paths. The Gauss-Newton systems here are tiny but not always full rank, and mpmath's lu_solve and qr_solve both refuse a singular matrix outright rather than -returning something usable. Working precision changes whether a given system -looks singular -- delley_1454 solves cleanly at 40 digits and fails at 20 -- -so this is a robustness property rather than a property of any one grid. +returning something usable. + +Forming the normal equations and then dropping ill-conditioned columns does not +work either: the two failure modes trade off against each other. With a fixed +1e-16 rank threshold, delley_1454 solves at 20 digits but delley_974 stalls at +3.1e-29 instead of 6.9e-51, because columns that still carry information are +being discarded. Tying the threshold to the working precision fixes 974 and +puts 1454 back to "numerically singular". The tension is not in the threshold +-- it is that J^T J squares the condition number. + +So the step is solved by a rank-revealing modified Gram-Schmidt QR of J itself. +Columns that are genuinely unresolvable are skipped and left at zero; the rest +are solved without ever squaring anything. """ -from mpmath import lu_solve, matrix, mp, mpf, sqrt +from mpmath import matrix, mp, mpf, sqrt -def independent_columns(J, tol=mpf("1e-16")): - """Indices of a maximal independent set of columns, by modified Gram-Schmidt. +def _mgs(J, tol): + """Rank-revealing modified Gram-Schmidt. - The threshold is relative to each column's own norm, so columns that differ - wildly in scale (a weight column against an angular one) are judged fairly. + Returns (kept, Q, R) with Q's columns orthonormal and R upper triangular + over the kept columns, so that J[:, kept] = Q R. """ - keep, basis = [], [] + kept, Q, R = [], [], [] for j in range(J.cols): v = [J[i, j] for i in range(J.rows)] nrm0 = sqrt(mp.fsum(x * x for x in v)) if nrm0 == 0: continue - for u in basis: - d = mp.fsum(v[i] * u[i] for i in range(len(v))) - v = [v[i] - d * u[i] for i in range(len(v))] + col = [] + for q in Q: + d = mp.fsum(v[i] * q[i] for i in range(len(v))) + col.append(d) + v = [v[i] - d * q[i] for i in range(len(v))] nrm = sqrt(mp.fsum(x * x for x in v)) - if nrm > tol * nrm0: - basis.append([x / nrm for x in v]) - keep.append(j) - return keep + if nrm <= tol * nrm0: + continue # not resolvable against what we have + Q.append([x / nrm for x in v]) + col.append(nrm) + R.append(col) # length len(Q) after the append + kept.append(j) + return kept, Q, R + + +def _default_tol(): + """Cut only what is unresolvable at the working precision. + + QR does not square the condition number, so it carries columns the normal + equations could not, and the threshold should discard only genuinely null + directions. Cutting at a fixed 1e-14 instead costs delley_974 its quadratic + convergence: 11 iterations to 1.97e-33 rather than 2 to 6.93e-51. + + Cases that need more working precision than they are given now fail the + generator's final verification rather than crashing or diverging silently. + delley_1454 is one: it refines cleanly at 40 digits and cannot be solved at + 20, whatever the threshold. + """ + return mpf(10) ** (-(mp.dps - 5)) + + +def independent_columns(J, tol=None): + """Indices of a maximal independent set of columns.""" + return _mgs(J, tol if tol is not None else _default_tol())[0] -def lstsq_step(J, F): - """Gauss-Newton step, solved only in the directions the Jacobian resolves. +def lstsq_step(J, F, tol=None): + """Gauss-Newton step for min ||J d + F||, without forming J^T J. - Directions outside the resolved set are left at zero: the conditions cannot - see them, so moving along them is guesswork. + Directions the Jacobian cannot resolve are left at zero: the conditions + cannot see them, so moving along them is guesswork. """ - keep = independent_columns(J) - m = len(keep) + kept, Q, R = _mgs(J, tol if tol is not None else _default_tol()) + m = len(kept) if m == 0: return [mpf(0)] * J.cols - JtJ = matrix(m, m) - Jtf = matrix(m, 1) - for a in range(m): - ja = keep[a] - for b in range(a, m): - jb = keep[b] - v = mp.fsum(J[k, ja] * J[k, jb] for k in range(J.rows)) - JtJ[a, b] = v - JtJ[b, a] = v - Jtf[a] = mp.fsum(J[k, ja] * F[k] for k in range(J.rows)) - red = lu_solve(JtJ, -Jtf) + + b = [-F[i] for i in range(J.rows)] + y = [mp.fsum(q[i] * b[i] for i in range(len(b))) for q in Q] + + # back-substitute R x = y; R[j] holds column j of R, entries 0..j + x = [mpf(0)] * m + for j in range(m - 1, -1, -1): + acc = y[j] + for k in range(j + 1, m): + acc -= R[k][j] * x[k] + x[j] = acc / R[j][j] + step = [mpf(0)] * J.cols - for a, j in enumerate(keep): - step[j] = red[a] + for a, j in enumerate(kept): + step[j] = x[a] return step From b5774f23700de1c939d9e50cf4a4793a3aefd558 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 22:17:54 +0300 Subject: [PATCH 12/17] tools: document the QR step in the README Records why the step is a rank-revealing QR of the Jacobian rather than a solve of the normal equations, and that a size needing more working precision than it is given fails verification instead of being written. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/README.md | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tools/angular_grids/README.md b/tools/angular_grids/README.md index 0d89a8a..6308f0f 100644 --- a/tools/angular_grids/README.md +++ b/tools/angular_grids/README.md @@ -49,6 +49,21 @@ or two angular parameters: So the largest Lebedev rule is a 385-parameter dense Newton solve, not a 17430-parameter one. +The step itself is a rank-revealing modified Gram-Schmidt QR of the Jacobian, +never the normal equations. Two things force that. The Jacobian is not always +full rank -- AB-312 comes out 13 of 16 -- and both `lu_solve` and mpmath's +`qr_solve` refuse a singular matrix outright. And forming `J^T J` squares the +condition number, which turns a solvable system into an unsolvable one: with a +fixed rank threshold on the normal equations, delley_974 stalls at 3.1e-29 +instead of reaching 6.9e-51, while tightening the threshold to fix that puts +delley_1454 back to "numerically singular". QR removes the tension rather than +trading the two failures against each other. + +A size that needs more working precision than it is given now fails the +generator's verification and is not written. delley_1454 is one: 102 parameters +at that conditioning are not resolvable in 20 digits, and it refines cleanly at +40. + The exactness conditions are imposed on even monomials. On the unit sphere `z^2 = 1 - x^2 - y^2`, so every even monomial reduces to a fixed linear combination of `x^(2a) y^(2b)`; imposing only those implies the rest, at about From 67190201d43d68c59f26428b6c5e057fc33b6003 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 22:27:35 +0300 Subject: [PATCH 13/17] tools: backtrack when a Newton step leaves an orbit's domain The octahedral orbit types carry implicit constraints. bk is (l, l, sqrt(1-2 l^2)) and is only real for l <= 1/sqrt(2); ck and dk have their own. A Newton step can overshoot one, and mpmath then returns a complex square root: l = 0.30000000 sqrt(1-2 l^2) = 0.90553851 l = 0.70710678 sqrt(1-2 l^2) = 5.7931539e-5 l = 0.72000000 sqrt(1-2 l^2) = (0.0 + 0.19183326j) The complex value then propagates silently until something tries to order two of them, which is where it finally surfaced -- as "TypeError: '<=' not supported between instances of 'mpc' and 'mpc'", far from the cause. It took out delley_2354, 2702, 3074 and 3470, the four largest sizes in the family. Halve the step until every orbit base is real again. Backtracking is cheaper than reparameterising the orbits, and it leaves the well-behaved sizes untouched. Found by the batch, not by the test suite: the eleven regression cases use small grids whose parameters sit well inside their domains, and only the 2354-and-larger sizes push one to the boundary. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/refine.py | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/tools/angular_grids/refine.py b/tools/angular_grids/refine.py index 2995478..aadb953 100644 --- a/tools/angular_grids/refine.py +++ b/tools/angular_grids/refine.py @@ -8,7 +8,7 @@ """ import collections -from mpmath import mp, mpf, fabs, matrix +from mpmath import fabs, im, matrix, mp, mpf from .linalg import lstsq_step from .moments import conditions, exact_moment, orbit_moment @@ -43,6 +43,22 @@ def decompose(points, weights, tol=mpf("1e-9")): return orbits +def _in_domain(orbits): + """Are all orbit bases still real? + + The orbit types carry implicit constraints -- bk is (l, l, sqrt(1-2 l^2)), + valid only for l <= 1/sqrt(2); ck and dk have their own. A Newton step can + overshoot one, and mpmath then returns a complex square root, which + propagates until something tries to order two complex numbers. Backtracking + on the step is cheaper than reparameterising the orbits. + """ + for o in orbits: + for v in o.type.base(o.params): + if im(v) != 0: + return False + return True + + def _pack(orbits): x = [] for o in orbits: @@ -95,7 +111,16 @@ def refine(orbits, order, target_dps, max_iter=12, verbose=False): if r < tol: break d = lstsq_step(J, F) - _unpack(orbits, [x + d[j] for j, x in enumerate(_pack(orbits))]) + base = _pack(orbits) + scale = mpf(1) + for _ in range(40): # backtrack until the step stays in domain + _unpack(orbits, [x + scale * d[j] for j, x in enumerate(base)]) + if _in_domain(orbits): + break + scale /= 2 + else: + _unpack(orbits, base) + break return hist From 219b1c78d12d7e094dd9241e1d652cf1701cb25e Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 22:39:12 +0300 Subject: [PATCH 14/17] tools: damp the octahedral Newton step Backtracking only on the orbits' domain was not enough. It stops the complex-square-root failure, but on the larger Delley grids the undamped step is simply bad: delley_2354 (order 83, 64 orbits, 161 parameters) walks away from a starting residual of 1.1e-16 to 6.4, 6.2, 5.6 over successive iterations. Halve the step until it both keeps every orbit real and reduces the residual, and give up rather than accept an increase. This is ordinary damped Newton; the quadratic convergence on the well-behaved sizes is unaffected, since the full step is accepted on the first try there. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/refine.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/tools/angular_grids/refine.py b/tools/angular_grids/refine.py index aadb953..9212515 100644 --- a/tools/angular_grids/refine.py +++ b/tools/angular_grids/refine.py @@ -112,13 +112,22 @@ def refine(orbits, order, target_dps, max_iter=12, verbose=False): break d = lstsq_step(J, F) base = _pack(orbits) + # Damped Newton. Staying in the orbits' domain is necessary but not + # sufficient: on the larger Delley grids the undamped step is simply + # bad, and the residual grows from 1e-16 to O(1) over a few iterations. + # Halve until the step both keeps every orbit real and actually reduces + # the residual. scale = mpf(1) - for _ in range(40): # backtrack until the step stays in domain + improved = False + for _ in range(60): _unpack(orbits, [x + scale * d[j] for j, x in enumerate(base)]) if _in_domain(orbits): - break + Fn, _ = _residual_and_jacobian(orbits, conds, want_jac=False) + if max(fabs(f) for f in Fn) < r: + improved = True + break scale /= 2 - else: + if not improved: _unpack(orbits, base) break return hist From eeb5c80e541adb0823863c329499ea75874be8e6 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 22:44:38 +0300 Subject: [PATCH 15/17] tools: equilibrate the Jacobian columns before solving A weight column and an angular column differ by orders of magnitude, and that disparity alone wrecks the conditioning of the least-squares solve. On delley_2354 (order 83, 161 parameters) the undamped step is unusable, and once damped the iteration crawls -- a factor of two per step instead of squaring. Scale every column to unit norm before the QR and unscale the step afterwards. It costs one pass over the Jacobian and is the difference between converging and not. Regression suite 11/11. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/linalg.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/tools/angular_grids/linalg.py b/tools/angular_grids/linalg.py index 12fc117..bf01440 100644 --- a/tools/angular_grids/linalg.py +++ b/tools/angular_grids/linalg.py @@ -73,7 +73,21 @@ def lstsq_step(J, F, tol=None): Directions the Jacobian cannot resolve are left at zero: the conditions cannot see them, so moving along them is guesswork. """ - kept, Q, R = _mgs(J, tol if tol is not None else _default_tol()) + # Column equilibration. A weight column and an angular column differ by + # orders of magnitude here, and that disparity alone wrecks the + # conditioning: on delley_2354 the undamped step is unusable and damping + # then crawls, a factor of two per iteration instead of squaring. Scaling + # every column to unit norm and unscaling the step afterwards costs + # nothing and is the difference between converging and not. + scale = [] + Js = matrix(J.rows, J.cols) + for j in range(J.cols): + c = sqrt(mp.fsum(J[i, j] ** 2 for i in range(J.rows))) + scale.append(c if c != 0 else mpf(1)) + for i in range(J.rows): + Js[i, j] = J[i, j] / scale[j] + + kept, Q, R = _mgs(Js, tol if tol is not None else _default_tol()) m = len(kept) if m == 0: return [mpf(0)] * J.cols @@ -91,5 +105,5 @@ def lstsq_step(J, F, tol=None): step = [mpf(0)] * J.cols for a, j in enumerate(kept): - step[j] = x[a] + step[j] = x[a] / scale[j] # undo the equilibration return step From 1944d114448aab42a9627b3aa25006317c59cd7f Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 22:48:48 +0300 Subject: [PATCH 16/17] tools: raise the working precision when verification comes up short The guard digits were a fixed margin, and that is wrong because the conditioning varies enormously with grid size. Measured after column equilibration, delley_1730's Jacobian has a condition number of 2.0e+24 -- so 24 of the working digits are gone before the step is even computed, and the larger sizes are worse. That is what has been defeating the four biggest Delley grids: a bad step from an under-resourced solve, which then overshot an orbit's domain and surfaced as a complex square root. Rather than predict the conditioning, try a modest guard and retry with more when the final verification falls short. The well-conditioned sizes pay nothing, since they pass on the first attempt. This is the root cause behind the last three commits. Each of those was a correct fix for a real symptom -- complex roots, divergence, crawling convergence -- but all three were downstream of a step computed with insufficient precision for the conditioning. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/generate.py | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/tools/angular_grids/generate.py b/tools/angular_grids/generate.py index 048258d..761b28c 100644 --- a/tools/angular_grids/generate.py +++ b/tools/angular_grids/generate.py @@ -69,9 +69,31 @@ def _run_design(family, npts, order, digits, root, verbose): return d.points, [d.weight] * d.n, designs.verify(d) -def run(family, npts, digits, root, out, verbose=True): - mp.dps = digits + 20 # guard digits for the refinement +def run(family, npts, digits, root, out, verbose=True, guards=(20, 60, 140)): + """Refine and emit one grid, raising the working precision if it is short. + + The guard is not a fixed margin because the conditioning varies enormously + with grid size: delley_1730's Jacobian has a condition number of 2.0e+24 + even after column equilibration, which eats 24 of the working digits before + the step is computed, and the larger sizes are worse. Rather than predict + that, try a modest guard and retry with more when verification comes up + short. + """ order = algebraic_order(family, npts, root) + last = None + for guard in guards: + mp.dps = digits + guard + try: + return _attempt(family, npts, order, digits, root, out, verbose, guard) + except SystemExit as e: + last = e + if verbose: + print(f" {e}; retrying at {digits + guard * 3} working digits", + flush=True) + raise SystemExit(last) + + +def _attempt(family, npts, order, digits, root, out, verbose, guard): if family in OCTAHEDRAL: P, W, err = _run_octahedral(family, npts, order, digits, root, verbose) elif family == "ahrens_beylkin": From 3211ad414bfdbe2ac499da514d2fd3fd4264a460 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 22:53:30 +0300 Subject: [PATCH 17/17] tools: size the initial precision guard from the parameter count Retrying on failure works but wastes a long attempt first, and on the biggest Delley grids that attempt takes many minutes. The conditioning is predictable enough to start from a sensible guard: measured on the equilibrated Jacobian it is 2.0e+24 at 120 parameters (delley_1730) and 8.19e+26 at 140 (delley_2030), so log10(cond) is roughly 0.15 * params. The extra guard this gives: delley_110 0 delley_2030 24 delley_974 10 delley_2354 28 delley_1730 20 delley_3470 42 Small grids are unaffected and pay nothing. The retry ladder stays as the backstop for anything the estimate underestimates. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- tools/angular_grids/generate.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/tools/angular_grids/generate.py b/tools/angular_grids/generate.py index 761b28c..2af1dbb 100644 --- a/tools/angular_grids/generate.py +++ b/tools/angular_grids/generate.py @@ -80,6 +80,14 @@ def run(family, npts, digits, root, out, verbose=True, guards=(20, 60, 140)): short. """ order = algebraic_order(family, npts, root) + + # Start from a guard the conditioning actually warrants. Measured on the + # equilibrated Jacobian, the condition number climbs steeply with the + # parameter count -- 2.0e+24 at 120 parameters (delley_1730), 8.19e+26 at + # 140 (delley_2030) -- so roughly log10(cond) ~ 0.15 * params + 6. Guessing + # low just burns a long failed attempt before the retry. + guards = tuple(g + _size_guard(family, npts, root) for g in guards) + last = None for guard in guards: mp.dps = digits + guard @@ -93,6 +101,25 @@ def run(family, npts, digits, root, out, verbose=True, guards=(20, 60, 140)): raise SystemExit(last) +def _size_guard(family, npts, root): + """Extra working digits warranted by this grid's parameter count.""" + mp.dps = 30 + try: + if family in OCTAHEDRAL: + points, weights = read_grid(grid_path(root, family, npts)) + orbits = decompose(points, [w * 4 * mp.pi / sum(weights) for w in weights]) + n = sum(1 + o.type.n_param for o in orbits) + elif family == "ahrens_beylkin": + from . import ab_refine as ab + points, weights = read_grid(grid_path(root, family, npts)) + n = ab.n_param(ab.decompose(points, weights)) + else: + n = 2 * npts + except Exception: # noqa: BLE001 - fall back to no extra + return 0 + return max(0, int(0.2 * n) - 4) + + def _attempt(family, npts, order, digits, root, out, verbose, guard): if family in OCTAHEDRAL: P, W, err = _run_octahedral(family, npts, order, digits, root, verbose)