numpy compat: full ufunc coverage, sorting, safe reductions, tolerance comparisons for mp dtypes - #306
Merged
Merged
Conversation
…/argmax slots np.abs, np.conj, the transcendental family (exp/log/trig/hyperbolic + inverses), power, sign, reciprocal, minimum/maximum/fmin/fmax, the rounding family (floor/ceil/trunc/rint half-to-even via mpfr_rint), remainder/fmod/floor_divide (numpy sign semantics), arctan2/hypot/ copysign, and the isnan/isinf/isfinite/signbit predicates now have guarded loops on the real_mp/complex_mp numpy dtypes -- previously "ufunc not supported" TypeErrors. All loops read through value_or_zero per the uninitialized-slot doctrine (ADR-0006) and call the same boost::multiprecision free functions the multiprec scalar functions bind. absolute on complex outputs real_mp. New HardenCompare/HardenArgMinMax fill the compare/argmax/argmin PyArray_ArrFuncs slots for real_mp: np.sort/argsort/searchsorted/ median/argmax/argmin now work (complex stays unordered by design). Also fixes a copy-paste bug where the scalar free function mp.imag() returned the REAL part, and adds array overloads of multiprec real/imag/arg: numpy's ndarray .real/.imag attributes silently return wrong values for legacy user dtypes (numpy cannot know complex_mp is complex-like), so these are the sanctioned array component accessors and the np.angle replacement. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…he buffer Stock eigenpy's getitem (and our heal-on-read specializations, which copied it) returns boost::ref(slot) -- a Python scalar ALIASING numpy array storage. Scalars extracted from temporary arrays dangled once the array was freed: np.sum(v) kept past the statement read freed memory (zeros, garbage precision, MPFR assertion SIGABRT on str()), and np.mean returned a silently-wrong 0 through the same mechanism. This is the root of the ADR-0031 / #259 hazard class, previously worked around consumer-by-consumer with copy-at-extraction. Returning an owned copy kills the class at the source: an indexed element is a durable value, as numpy users expect, and the identity-seeded reductions (np.sum/np.prod/np.mean) are now genuinely safe rather than accidentally readable. Adds python/test/classes/numpy_ufuncs_test.py: element-wise agreement of every new ufunc with the multiprec scalar functions, numpy semantics corners (rint half-to-even, remainder sign-of-divisor, nan propagation in minimum/maximum vs fmin/fmax), sort/argmax/searchsorted/ median, reductions regression, precision preservation through every loop shape, unwritten-slot safety, the mp.imag regression, the array real/imag/arg accessors, and named regressions for the dangling-scalar fix. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The gotchas page now documents what IS supported (the full ufunc set,
sorting, reductions -- verified on numpy 2.3/2.4 and regression-tested)
and shrinks the gotchas to the real, permanent edges: the float64
boundary (by design -- double->mp casts stay unsafe, reaffirmed
2026-07-08: the user has to think about float-literal intent; the
sanctioned closeness idiom is np.all(np.abs(a-b) <= real_mp('1e-10'))),
complex component access (.real/.imag/np.angle silently wrong for
legacy user dtypes -- use the new mp.real/imag/arg array overloads),
complex ordering (deliberately unimplemented), and no mp->int casts.
The reductions section becomes a version-qualified historical note with
the initial= idiom kept as the old-numpy fallback.
ADR-0051 records the whole decision set, including the owned-copy
getitem that root-causes the ADR-0031/#259 aliasing class. Also
updates the stale endgame_test comment, the multiprec docstring
(and its pre-rename type names), and a stale dtype name in the
precision-models tutorial. All 180 sphinx doctests pass.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…conversion Comparisons against a double tolerance -- np.abs(a - b) < 1e-10 -- are safe: an ordering yields a bool, so no float value ever flows into a multiprecision computation, and the comparison is exact (boost compares the number against the double directly; 1e-22 is not lost against a 1e-30 tolerance). This matches the C++ solvers' double ToleranceT and the line the scalar bindings already drew (GreatLessVisitor<T,double> bound, double equality deliberately not). Register mixed mp-vs-float64 ordering loops (< <= > >=, both operand orders, real type only). Mixed EQUALITY with a float stays unregistered (exact equality against a float literal is the 0.1-intent trap the unsafe double->mp cast exists to block), as does mixed arithmetic; np.isclose/np.allclose still raise (they COMPUTE with float64 tolerances internally). The float boundary protects polynomial-system construction, not tolerance checks. Also register mp -> complex128 casts as unsafe, alongside the pre-existing mp -> double: arr.astype(complex) / astype(float) are the explicit, conscious truncations (astype(complex) had never actually been registered). Registration-order fix: casts now register BEFORE the ufunc loops -- registering the mixed loops makes numpy query the mp<->double casts, and a cast first registered after being queried is permanently ignored (numpy RuntimeWarning "registered/modified ... after the cast had been used"). Docs: the gotchas page float64-boundary section now states the values-vs-comparisons line and the plain-float closeness idiom; ADR-0051 decision 5 amended. Tests cover both operand orders, exactness of the mixed compare, the still-blocked equality/arithmetic/ isclose, and astype. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…f-to-even) numpy defines rint for complex dtypes (rounds real and imaginary parts independently) and np.round dispatches through it; complex_mp had no loop, so np.round(complex_mp array) raised. Register the component-wise rint for complex_mp, matching complex128 semantics exactly (half-to-even via mpfr_rint in MPFR_RNDN). Rounds out the np.abs/np.round ask of #301. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Contributor
Author
|
Contributor
Author
|
Fable make many hours of work, very simple. I love using this system. Let me focus on the math!!! |
…docs page reframed Reconciles #305's interim numpy helpers with the binding-level ufunc coverage from this branch: - bertini.real/imag/abs/conj/sum/norm/is_real keep their API but now take the native C++ ufunc loops when handed an mp-dtype array (returning proper mp-dtype arrays that keep working with sort and reductions), falling back to element-wise work for lists and mixed input. bertini.round keeps its decimal-digit Decimal semantics. - The "Known gotchas" page is retired. Its replacement, docs/source/numpy.rst ("Multiprecision numbers and NumPy"), documents the interop as the feature it now is: what works (everything), the float64 boundary as deliberate digit protection with the tolerance- comparison carve-out, the complex component accessors, and the boundaries by design. The one remaining hazard -- ndarray .real/.imag lie on complex_mp arrays, a numpy hardwiring no binding can reach -- is a warning admonition, not a page of dragons. - All cross-references updated (module docstrings, test comments). 769 tests pass (both suites merged), 169 doctests, doclint clean. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
numpy's ndarray.real/.imag are C getsets gated on PyArray_ISCOMPLEX -- a hardwired builtin-type-number check (verified against numpy 2.4.x getset.c). No user dtype, legacy or NEP-42-style, can make the attributes correct: on complex_mp arrays .real returns the complex values and .imag returns zeros, silently. Decision: a crash is better than incorrect values. Defend every spelling within reach: - Solution overrides .real/.imag at the subclass level (a python property shadows the C getset), so solve results are simply correct -- for the mp dtype and for double solves alike. - bertini._numpy_guard, installed at import: np.real/np.imag raise a TypeError naming the right tool on plain mp-complex arrays (and on lists that would convert to them inside numpy); np.angle raises for EVERY mp-complex input, since it branches on dtype and dies in arctan2 where not even a subclass property can reach; all other inputs pass through to numpy untouched. - The raw .real/.imag attributes on a plain self-built ndarray remain the one spelling nothing can reach -- documented with a warning and a pinning test. Docs page and ADR-0051 updated accordingly. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…bulary `from bertini.operators import *` now gives functions that work on EVERYTHING, dispatched per argument: a symbolic Variable/expression builds a function-tree node; multiprecision scalars, mp-dtype numpy arrays, lists, and plain python numbers take the numeric path through the native precision-preserving loops. No more remembering that arg lives in bertini.multiprec while imag is at the top level. - sin/cos/tan/asin/acos/atan/exp/log/sqrt: fully polymorphic (symbolic twin exists). - sinh/cosh/tanh/asinh/acosh/atanh and abs/arg/real/imag/conj/round/ sum/norm/is_real: numeric, with a clear TypeError on symbolic input. - E/Pi/I ride along. abs/round/sum shadow the builtins only inside this opt-in star-import (they fall back to builtin behavior on plain python input); `from bertini import *` still never shadows. - the top-level elementary functions are rebound to the polymorphic versions (a strict superset: bertini.sin(x) now also accepts numbers and arrays), and bertini.arg + the hyperbolics join the top level. - new _numpy_helpers.arg rides mp.arg on complex arrays (real arrays go through the exact real->complex cast). Also fixes a latent #305 bug the new tests caught: _numpy_helpers' scalar fallbacks called bare abs()/round(), which resolve to the module's own shadowing functions at module scope -> infinite recursion on plain python input. Now explicitly builtins. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Contributor
Author
|
and a while later, i understand that numpy has a significant BUG for complex user types: imag and real will always return incorrect values. what?!?!?!? #307 is the lasting memory of it, and I should probably do the work to help numpy fix this. for now, users will get an exception if they try to, because i would rather have an error, than an incorrect value/result. |
Contributor
Author
|
here's the final state of the PR's work, in picture form (assuming i get green on CI and don't have to do more work): 
|
….5 reduce UAF) CI segfaulted on Linux py3.12-3.14 (numpy 2.5.1; py3.10/3.11 got numpy 2.2/2.4 and passed). Root cause, pinned locally with valgrind after reproducing in a numpy-2.5.1 venv: numpy 2.5 initializes the accumulator of an IDENTITYLESS reduce (np.min/np.max) by memcpy of element 0, so the accumulator slot and v[0] share one mpfr allocation. Our loops' `res = Op::apply(...)` move-assignment freed the slot's old limbs -- freeing v[0]'s storage out from under it (np.max(v) also silently rewrote v[0] through the shared limbs: v[0] went from 3 to 3.5). The freed block was re-read by the next reduce: use-after-free, double-free, corrupted mimalloc metadata, and a SIGSEGV three tests later inside np.median -- the classic action-at-a-distance crash. Fix: slot_write is now the only sanctioned store into an mp output slot. It computes the value first (the slot may also be an input), memsets the slot to BMP's uninitialized sentinel, and move-assigns the fresh value in -- no existing allocation is ever freed or written through. Applied to every loop that writes mp slots (binary, unary, cross-type unary, matmul, dotfunc); bool outputs stay plain stores. Same crash-into-bounded-leak doctrine as HardenSetitem (ADR-0006); recorded as ADR-0051 decision 7. Verified: numpy 2.5.1 -- min/max/median correct, input arrays intact, valgrind 0 errors (previously UAF pair between the maximum and minimum loops + fatal invalid read inside mi_malloc); numpy 2.4.6 -- full suite green. New regression test runs the reduces repeatedly and asserts the input survives. Also fixes the all-Windows wheel failure (separate, infrastructural): clang rejects a PCH whose input mtime changed even with identical content, which happens when the wheel build re-configures into a shared bld/. Add -Xclang -fno-pch-timestamp to the two PCH targets for clang builds. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
ofloveandhate
added a commit
that referenced
this pull request
Jul 9, 2026
…ows PCH -fno-pch-timestamp fix)
ofloveandhate
added a commit
that referenced
this pull request
Jul 9, 2026
Bump `VERSION` from `3.1.0.dev1` to **`3.1.0rc1`** — release candidate for 3.1.0. `VERSION` is the single source of truth (scikit-build-core reads it dynamically; `publish.yml`'s `check_version` asserts a release tag matches it). This one-line change is docs-adjacent: `VERSION` is in the workflow `paths-ignore`, so this PR runs **no CI** and needs no approvals — it can merge straight in. 3.1.0.dev1 → 3.1.0rc1 gathers the 3.1.0 line so far: the UI quality-of-life batch (#305, issues #293–#304 + solution group projection), the numpy-compat overhaul (#306), and `ZeroDimConfig.recall` (#308).
ofloveandhate
added a commit
that referenced
this pull request
Jul 9, 2026
The changelog had drifted badly: the last entry was **2.0.1**, while **2.0.2**, the entire **3.0.0** modernization, and the **3.1.0** line had all shipped — with their notes living only in commit messages, PRs, and GitHub Releases. This consolidates them back into `CHANGELOG.md`, newest-first, in the existing *Keep a Changelog* format. ### New entries - **[3.1.0] – 2026-07-09** — NumPy interop for the mp dtypes (#306), the Python UI quality-of-life batch (#293–#304, #305), `ZeroDimConfig.recall` (#308), prebuilt CI deps (ADR-0049, #282), and docs-store Pages (ADR-0050, #291, #292). - **[3.0.0] – 2026-07-07** — reworked from the hand-written v3.0.0 release notes (~70 PRs; full themed index in #238) into Added / Changed / Fixed sections. - **[2.0.2] – 2026-05-22** — the packaging/CI maintenance entry that was never recorded. Older 1.0.x / 2.0.1 entries and the commented template are untouched. ### Why now (load-bearing) `publish.yml`'s `github-release` job builds the release body from the **top** `CHANGELOG.md` block. That block was the stale **[2.0.1]** — so a final `v3.1.0` tag would have published 2.0.1's notes as the 3.1.0 release. With this merged, the extraction yields exactly the **[3.1.0]** block (verified locally against the workflow's extraction logic). Docs-only (`**/*.md` → `paths-ignore`), so this runs no CI. **Merge before tagging `v3.1.0`.**
ofloveandhate
added a commit
that referenced
this pull request
Jul 9, 2026
…ns (v3.1.0 docs deploy) (#313) The **v3.1.0** versioned-docs deploy failed. `build_docs.yml` runs `sphinx-build -b html -W` (warnings-as-errors), and three docstrings wrote absolute-value / norm notation with **bare pipes**: ``` ERROR: Undefined substitution referenced: "p_i - q_i" — bertini.multiprec.is_distinct_up_to ERROR: Undefined substitution referenced: "p_i - q_i" — bertini.is_distinct_up_to (top-level re-export) ERROR: Undefined substitution referenced: "imag" — bertini.operators.is_real ``` In reStructuredText `|word|` is a **substitution reference**, so Sphinx tried to resolve undefined substitutions `p_i - q_i` and `imag` and errored out. `is_distinct_up_to` (#305) and `is_real` (#306) both landed in the 3.1.0 line. ### Fix Reword the two source strings to `abs(...)` — `max_i abs(p_i - q_i)` and `abs(imag) < tol` — which read cleanly in Sphinx HTML *and* plain `help()` and carry no RST metacharacters. The other `|…|` docstrings in the tree are already safe (wrapped in ``inline literals`` or inside code blocks) and are left untouched. ### Not a package problem Docstrings only — the released wheels are unaffected. **3.1.0 is already on PyPI and the GitHub Release is published**; only the docs-site deploy failed. ### Verified Reproduced the exact failing command locally (`sphinx-build -b html -W --keep-going`): now **`build succeeded`**, zero substitution errors. ### Root-cause note (for the post-3.1.0 CI pass) The `-b html -W` build runs **only at docs-deploy time**, not in PR CI (PR CI runs `-b doctest`, which passed). So this whole class of RST/docstring error is invisible until release. Running the `-W` html build in PR CI would have caught it — worth adding alongside #312.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What this PR delivers
The numpy gotchas for the
real_mp/complex_mpdtypes are gone. The scary "Known gotchas" page is retired; its replacement,docs/source/numpy.rst("Multiprecision numbers and NumPy"), documents the interop as the feature it now is. Full decision record in ADR-0051; includes the post-#305 reconciliation.Full ufunc coverage
np.abs,np.conj,np.exp/log/log2/log10/expm1/log1p, all trig/hyperbolic + inverses,power,sign,reciprocal,minimum/maximum/fmin/fmax,floor/ceil/trunc/rint,isnan/isinf/isfinite/signbit,arctan2,hypot,copysign,mod/fmod/floor_divide— previously allTypeError: ufunc not supported. Every loop reads throughvalue_or_zero(ADR-0006 doctrine) and calls the same boost::multiprecision function themultiprecscalar functions bind, sonp.f(a)[i] == mp.f(a[i])exactly. numpy-semantics corners done deliberately:rintis half-to-even via directmpfr_rint(boost's is half-away; registered for complex too, component-wise, sonp.roundworks on complex arrays),modtakes the divisor's sign,minimum/maximumpropagate nan whilefmin/fmaxignore it.Sorting & reductions — and the bug hiding under them
HardenCompare/HardenArgMinMaxfill the dtype slots forreal_mp:np.sort,np.argsort,np.searchsorted,np.median,np.argmax/argmin,np.max/min(complex stays unordered by design).np.sum/np.prod/np.meanwork (verified numpy 2.3.2 & 2.4.6, regression-tested) — but investigating them exposed the real hazard:getitemreturnedboost::refinto the numpy buffer (stock eigenpy behavior), so a scalar extracted from a temporary array dangled —np.meanreturned a silently wrong0, ands = np.sum(v)kept past the statement SIGABRT'd instr(). This was the root of the ADR-0031 / #259 aliasing class.getitemnow returns an owned copy, killing the class at the source.The float64 boundary: closed for values, open for tolerance comparisons
double→mpcasts stay unsafe — floats must not pollute polynomial-system construction. But tolerance orderings against a double are registered (np.abs(a - b) < 1e-10, both operand orders): a comparison yields a bool, no float value enters an mp computation, and the compare is exact (1e-22 is not lost against a 1e-30 tolerance). Mirrors the C++ solvers' doubleToleranceTand the scalarGreatLessVisitor<T,double>precedent. Mixed equality and arithmetic stay blocked;np.isclose/allclosestill raise (they compute with float64 tolerances).arr.astype(float)/astype(complex)are the explicit conscious truncations (mp→complex128cast registered unsafe; it never had been). Registration-order note: casts must register before the ufunc loops, or numpy permanently ignores them.Component access on complex arrays — "a crash is better than incorrect values"
numpy's
ndarray.real/.imagare C getsets gated onPyArray_ISCOMPLEX, a hardwired builtin-type check with no user-dtype hook (verified against numpy 2.4.xgetset.c; neither the legacy API nor NEP 42 offers one) — on mp-complex arrays they return silently wrong values. Defended everywhere reachable:Solutionoverrides.real/.imagat the subclass level — solve results are simply correct;bertini._numpy_guard, installed at import, wrapsnp.real/np.imag/np.angleto raise aTypeErrornaming the right tool on plain mp-complex input (pass-through otherwise;angleraises for every mp-complex input since it branches on dtype);bertini.real/imag/argandmp.real/imag/arg(en route: fixed the scalarmp.imag, which had returned the real part since it was written — copy-paste bug);One namespace for the whole vocabulary
bertini.operatorsis now polymorphic:from bertini.operators import *givessin/cos/exp/... that dispatch per argument (symbolic node for a Variable/expression; numeric through the native precision-preserving loops for mp scalars, numpy containers, lists, and plain numbers), plusabs/arg/real/imag/conj/round/sum/norm/is_realandE/Pi/I. The top-level elementary functions are rebound to the polymorphic versions (a strict superset), andbertini.arg+ the hyperbolics join the top level. No more remembering thatarglives inbertini.multiprecwhileimagis top-level.#305 reconciliation
Develop (the UI QoL batch) is merged in. Its interim
bertini.real/imag/abs/conj/sum/norm/is_realhelpers keep their API but now ride the native ufunc loops on mp-dtype arrays, falling back to element-wise work for lists. Two latent #305 bugs fixed: the helpers' scalar fallbacks called bareabs()/round(), which resolve to the module's own shadowing names → infinite recursion on plain python input; and the docs claim that the numpy boundary "cannot be patched in the bindings" (this PR is the patch).Tests / verification
numpy_ufuncs_test.py(~120) +operators_test.py(17): exact agreement with the scalar functions, semantics corners, precision preservation through every loop shape (including the mixed real/complex division path — loops re-tag viaat_precision_of), unwritten-slot safety per ADR-0006, sorting incl. nan-wins argmax, reductions, tolerance comparisons, astype, guard behavior, polymorphic dispatch, and named regressions for the dangling-scalar,mp.imag, and shadowed-builtin bugs.Linked issues
.realin numpy does NOT get the real part #298 —.realin numpy doesn't get the real part: Solution overrides + guarded numpy functions +bertini.real/imag/mp.real/imag/argaccessors; the scalarmp.imagbug found under it is fixed. The unreachable plain-attribute residue is numpy ndarray .real/.imag return silently wrong values on complex_mp arrays — no user-dtype hook exists in numpy #307.np.abs(complex_mp)/np.round(complex_mp)and the rest of the ufunc family now just work.🤖 Generated with Claude Code