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..6308f0f --- /dev/null +++ b/tools/angular_grids/README.md @@ -0,0 +1,210 @@ +# 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 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 +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. +* **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) + +`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. + + +## 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. 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/ab_refine.py b/tools/angular_grids/ab_refine.py new file mode 100644 index 0000000..b293dc0 --- /dev/null +++ b/tools/angular_grids/ab_refine.py @@ -0,0 +1,267 @@ +"""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, 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 + + +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 image_cache(orbits): + """Group images of each orbit representative, and of its tangent basis. + + 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] + 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)) + 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) + 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, 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 + _apply(orbits, lstsq_step(J, F)) + 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 + + +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) + # 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[:limit]: + _, 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 not chosen: + return conditions(order) + 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) 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/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/emit.py b/tools/angular_grids/emit.py new file mode 100644 index 0000000..947380f --- /dev/null +++ b/tools/angular_grids/emit.py @@ -0,0 +1,102 @@ +"""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. + /// + /// 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]); + w[i] = detail::grid_scalar::parse(weight_digits[i]); + }} + }} +}}; + +}} // namespace {ns} +}} // namespace IntegratorXX +""" + +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): + 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..2af1dbb --- /dev/null +++ b/tools/angular_grids/generate.py @@ -0,0 +1,174 @@ +"""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)) + + +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) # 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)} 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, 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) + + # 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 + 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 _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) + 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}") + 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 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: + 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: + 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", "ahrens_beylkin", "womersley"]) + 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/icosahedral.py b/tools/angular_grids/icosahedral.py new file mode 100644 index 0000000..5680590 --- /dev/null +++ b/tools/angular_grids/icosahedral.py @@ -0,0 +1,154 @@ +"""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 + + +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} diff --git a/tools/angular_grids/linalg.py b/tools/angular_grids/linalg.py new file mode 100644 index 0000000..bf01440 --- /dev/null +++ b/tools/angular_grids/linalg.py @@ -0,0 +1,109 @@ +"""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. + +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 matrix, mp, mpf, sqrt + + +def _mgs(J, tol): + """Rank-revealing modified Gram-Schmidt. + + Returns (kept, Q, R) with Q's columns orthonormal and R upper triangular + over the kept columns, so that J[:, kept] = Q R. + """ + 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 + 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: + 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, tol=None): + """Gauss-Newton step for min ||J d + F||, without forming J^T J. + + Directions the Jacobian cannot resolve are left at zero: the conditions + cannot see them, so moving along them is guesswork. + """ + # 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 + + 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(kept): + step[j] = x[a] / scale[j] # undo the equilibration + return step 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..9212515 --- /dev/null +++ b/tools/angular_grids/refine.py @@ -0,0 +1,155 @@ +"""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 fabs, im, matrix, mp, mpf + +from .linalg import lstsq_step +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 _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: + 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 + 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) + improved = False + for _ in range(60): + _unpack(orbits, [x + scale * d[j] for j, x in enumerate(base)]) + if _in_domain(orbits): + Fn, _ = _residual_and_jacobian(orbits, conds, want_jac=False) + if max(fabs(f) for f in Fn) < r: + improved = True + break + scale /= 2 + if not improved: + _unpack(orbits, base) + break + 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/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()) diff --git a/tools/angular_grids/tables.py b/tools/angular_grids/tables.py new file mode 100644 index 0000000..90e42a5 --- /dev/null +++ b/tools/angular_grids/tables.py @@ -0,0 +1,54 @@ +"""Read the tabulated angular grids out of the generated C++ headers.""" +import re +from pathlib import Path + +from mpmath import mp, mpf + +# 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): + 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] + + +_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. + + 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")) + 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 + + +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..db12b53 --- /dev/null +++ b/tools/angular_grids/test_refine.py @@ -0,0 +1,67 @@ +"""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 . import designs +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)] + +# 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 + 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 + 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 + + +if __name__ == "__main__": + sys.exit(main(sys.argv[1] if len(sys.argv) > 1 else "."))