diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index e5b171634..0f558f1d2 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -570,6 +570,15 @@ target_precompile_headers(bertini2_exe PRIVATE ) +# clang validates a PCH by the mtime of its input; when the wheel build configures +# twice in one job (test build + wheel build sharing bld/), cmake regenerates an +# IDENTICAL cmake_pch.cxx with a fresh mtime and clang refuses the still-valid PCH +# ("has been modified since the precompiled header was built: mtime changed"). +# -fno-pch-timestamp exists precisely for this (see the ccache docs on PCH). +if(CMAKE_CXX_COMPILER_ID MATCHES "Clang") + target_compile_options(bertini2_exe PRIVATE "$<$:SHELL:-Xclang -fno-pch-timestamp>") +endif() + # todo: this should be made a devmode thing #target_compile_options(bertini2 PRIVATE -Wall -Wextra) diff --git a/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md new file mode 100644 index 000000000..c1dccf848 --- /dev/null +++ b/docs/adr/0051-numpy-ufunc-coverage-and-owned-getitem.md @@ -0,0 +1,161 @@ +# ADR-0051: Full numpy ufunc coverage for the mp dtypes; getitem returns owned copies + +**Status:** Accepted +**Date:** 2026-07-08 + +## Context + +The eigenpy-registered numpy dtypes for `real_mp` / `complex_mp` covered only the +arithmetic core (add/subtract/multiply/divide, equality + real orderings, +negative/square/sqrt, matmul, and the hardened `dotfunc`). Everything else — +`np.abs`, `np.conj`, the transcendental family, `power`, `sign`, `minimum/maximum`, +rounding, the `isnan` predicates — raised `ufunc ... not supported`; `np.sort` / +`np.argmax` failed on empty dtype slots ("type does not have compare function" / +"data type not ordered"); and the docs declared the identity-seeded reductions +(`np.sum`/`np.prod`/`np.mean`) permanently unsupported after a build-dependent +`SystemError` (documented 2026-06-29, "not something Bertini can patch"). + +Investigating the reductions on current numpy (2.3.2 and 2.4.6) showed the +`SystemError` no longer reproduces — but exposed something worse hiding behind them: + +**eigenpy's `getitem` returns `boost::ref(slot)` — a Python scalar aliasing numpy +array storage.** Our heal-on-read `getitem` specializations (ADR-0006) had copied +that behavior faithfully. Consequences: + +- A scalar extracted from a *temporary* array dangles once the array is freed. + `s = np.sum(v)` is exactly that: the reduce result is a temporary 0-d array, + `s` aliased its buffer, and reading `s` later gave zeros, garbage precision, or an + MPFR assertion SIGABRT inside `str()`. +- `np.mean` returned a **silently wrong `0`** through the same mechanism. +- This is the root of the ADR-0031 / #259 hazard class ("indexing an eigenpy Vec + returns an aliasing view"), previously mitigated consumer-by-consumer with a + copy-at-extraction rule in Python code. + +## Decision + +1. **Register guarded loops for the full ufunc set** on both dtypes + (`python_bindings/include/eigenpy_interaction.hpp`, `registerGuardedUfunct`): + + - both dtypes: `absolute` (complex → `real_mp` output), `conjugate`, `sign` + (numpy-2 semantics: complex sign is `z/|z|`), `positive`, `reciprocal`, + `power`, `exp`, `log`, `log10`, full trig/hyperbolic + inverses, + `isnan`/`isinf`/`isfinite` (→ bool); + - real only (ordering- or domain-dependent): `greater`/`less`/... , `fabs`, + `exp2`, `log2`, `expm1`, `log1p`, `cbrt`, `floor`, `ceil`, `trunc`, `rint`, + `signbit`, `arctan2`, `hypot`, `copysign`, `fmod`, `remainder`, + `floor_divide`, `minimum`, `maximum`, `fmin`, `fmax`. + + Every loop reads through `value_or_zero` (the ADR-0006 uninitialized-slot + doctrine) and calls the same boost::multiprecision free function the + `bertini.multiprec` scalar function binds, so `np.f(a)[i] == mp.f(a[i])` + exactly. Deliberate semantic choices, matching numpy's float64 behavior: + `rint` rounds half-to-even (direct `mpfr_rint` in `MPFR_RNDN`; boost's `rint` + rounds half away), `remainder`/`mod` takes the sign of the divisor (`fmod` + keeps C semantics), `minimum`/`maximum` propagate nan while `fmin`/`fmax` + ignore it. + +2. **Fill the `compare`/`argmax`/`argmin` dtype slots for `real_mp`** + (`HardenCompare`, `HardenArgMinMax`, same install pattern as `HardenDotfunc`). + Enables `np.sort`/`argsort`/`searchsorted`/`unique`/`median`/`argmax`/`argmin`. + Complex stays unordered on purpose — numpy's lexicographic complex ordering is + historical baggage we do not reproduce. + +3. **`getitem` returns an owned copy, never `boost::ref`.** An indexed element is + a durable value, as numpy users expect. This kills the ADR-0031/#259 hazard + class at the source (that ADR's copy-at-extraction rule in Python remains good + hygiene but is no longer load-bearing), and it is what makes the reductions + *actually* safe rather than accidentally readable. Cost: one mp copy per + element read. + +4. **Reductions are supported and regression-tested**, on numpy ≥ 2.3 (verified + 2.3.2 and 2.4.6; `python/test/classes/numpy_ufuncs_test.py::TestReductions` + pins them in CI on all three platforms). The docs note the historical + `SystemError` and keep the `initial=` idiom as the fallback for older numpy. + `pyproject.toml` keeps `numpy` unpinned. + +5. **The float64 boundary stays closed for VALUES, open for tolerance + comparisons** (both decided 2026-07-08). `double → mp` casts remain + registered *unsafe*, so float64 scalars/arrays do not silently promote into + mp arrays: the conversion itself is bit-exact, but a promoted float64 `0.1` + is not the decimal `0.1` the user typed — the user has to think. However, + mixed mp-vs-float64 **ordering** loops (`<`, `<=`, `>`, `>=`, both operand + orders, real only) ARE registered: `np.abs(a - b) < 1e-10` is safe — the + result is a bool, no float flows into an mp value, the comparison is exact + (boost compares the number against the double directly), and it matches the + C++ solvers' double `ToleranceT` and the scalar `GreatLessVisitor` + precedent. Mixed EQUALITY stays unregistered (exact equality against a + float literal is the 0.1-intent trap; not bound at scalar level either), as + does mixed arithmetic; `np.isclose`/`np.allclose` still raise (they *compute* + with float64 tolerances internally). `mp → complex128` casts are registered + unsafe alongside the pre-existing `mp → double`, so `arr.astype(complex)` / + `astype(float)` are the explicit, conscious truncations. Registration-order + note: casts must be registered 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). + +6. **Component access on complex arrays** goes through new array overloads of + `multiprec.real`/`imag`/`arg` (returning `real_mp` arrays). numpy cannot know + a legacy user dtype is complex-like — `ndarray.real`/`.imag` are C getsets + gated on `PyArray_ISCOMPLEX`, a hardwired builtin-type-number check (verified + against numpy 2.4.x `getset.c`; neither the legacy user-dtype API nor the + NEP 42 new DType API offers a complex-like hook) — so on mp-complex arrays + `.real` returns the complex values and `.imag` returns zeros, silently. Not + fixable at the attribute; defended everywhere reachable (decided 2026-07-08, + "a crash would be better than incorrect values"): + + - `bertini.records.Solution` overrides `.real`/`.imag` at the subclass level — + solve results are simply correct; + - importing bertini wraps the module functions `np.real`/`np.imag`/`np.angle` + (`bertini._numpy_guard`): plain mp-complex arrays (and lists that would + convert to them) raise a `TypeError` naming the right tool; `np.angle` + raises for every mp-complex input since it branches on dtype and no subclass + property can reach it; everything else passes through untouched; + - the raw `.real`/`.imag` attributes on a plain self-built ndarray remain the + single lying spelling — documented with a warning and a pinning test. + + (En route, fixed a copy-paste bug: the scalar `mp.imag` was bound to + `boost::multiprecision::real` and returned the real part.) + +7. **Loops write output slots only through `slot_write`** (added 2026-07-08, after + CI segfaults on numpy 2.5.1): numpy may hand a loop an output slot that + **bitwise-aliases** another slot's mpfr allocation — numpy 2.5 initializes the + accumulator of an *identityless* reduce (`np.min`/`np.max`) by `memcpy` of + element 0, so the accumulator and `v[0]` share one set of limbs. A plain + BMP assignment move-frees the slot's old limbs (freeing `v[0]`'s storage: + use-after-free → double-free → corrupted allocator → SIGSEGV several calls + later — the CI crash landed in `np.median`, three tests after the damage) + and writes through shared storage (`np.max(v)` silently rewrote `v[0]`). + `slot_write` computes the value first, memsets the slot to BMP's + uninitialized sentinel, and move-assigns the fresh value in — no existing + allocation is ever freed or written through. Same crash-into-bounded-leak + trade as `HardenSetitem`. Diagnosed with valgrind (UAF pair between the + `maximum` and `minimum` reduce loops); regression test + `test_identityless_reduce_does_not_corrupt_input` runs the reduces + repeatedly and asserts the input array survives. numpy < 2.5 never + aliased, which is why every pre-2.5 environment was green. + +## Consequences + +- The "Known gotchas" docs page shrinks to the real, permanent edges: the float64 + boundary (by design), complex component access (numpy limitation), complex + ordering (by design), no mp→int casts. The reductions section becomes a + historical note. +- `python/test/classes/numpy_ufuncs_test.py` pins: element-wise agreement with the + scalar functions, the numpy-semantics corners above, precision preservation + through every loop shape (including `sign`/`reciprocal`, which cross the mixed + real/complex division path with the known boost precision-mis-tagging hazard — + loops re-tag via `at_precision_of`), unwritten-slot safety per loop shape, + sorting/arg-extrema (nan-wins semantics), reductions, and the dangling-scalar + regressions. +- Not upstreamed to eigenpy (the guarded loops already diverge; see ADR-0006). + Upstreaming the owned-copy `getitem` would fix the aliasing class for all + eigenpy user types and may be worth an issue later. + +## Relation to prior ADRs + +- **ADR-0006** — the slot-guard doctrine these loops follow; its "known gaps" list + shrinks (compare/argmax slots now filled and guarded). +- **ADR-0031 / #259** — root-caused and fixed at the binding level by the owned-copy + `getitem`; the Python-side copy rule is now belt-and-suspenders. +- **ADR-0001/0008** — unrelated eigenpy hazards, unchanged. diff --git a/docs/adr/README.md b/docs/adr/README.md index c50fe2370..878340bc5 100644 --- a/docs/adr/README.md +++ b/docs/adr/README.md @@ -67,3 +67,4 @@ Each ADR follows the template: | [0047](0047-casual-records-surface.md) | The casual records surface: bertini.solve/save/load, Solution = points that remember, CLI records-on-by-default beside b1 files (no flag) | Python + CLI / records | | [0048](0048-cauchy-endgame-security-and-operating-zone.md) | Cauchy divergence handling: security check watches the ENDPOINT, truncates only in the operating zone, no pole-growth truncation (the acceptance gate alone cures junk-success); restores cyclic-6's 156 solutions (refines #70) | Core / endgames | | [0050](0050-docs-deploy-from-branch-not-deploy-pages.md) | Docs publish by serving the docs-store branch directly (Pages branch-source), NOT actions/deploy-pages — which failed structurally with BlobNotFound on the versioned store; custom domain via a /CNAME file; /stable/ is a redirect, .doctrees dropped | CI / docs | +| [0051](0051-numpy-ufunc-coverage-and-owned-getitem.md) | Full numpy ufunc coverage for the mp dtypes (abs/conj/transcendentals/min-max/rounding/predicates, guarded loops), sort/argmax dtype slots (real only), getitem returns owned copies (kills the ADR-0031/#259 aliasing class; reductions genuinely safe); float64 boundary stays closed | Python bindings / numpy | diff --git a/python/bertini/__init__.py b/python/bertini/__init__.py index 7de4c87d2..840e73809 100644 --- a/python/bertini/__init__.py +++ b/python/bertini/__init__.py @@ -122,9 +122,6 @@ from . import _slice_ops as _slice_ops _slice_ops.install(nag_algorithm.Slice) -from . import operators # `from bertini.operators import *` -> just the math ops - - # --- sympy interop (#295) ---------------------------------------------------------------------- # Make every function-tree node auto-convert to sympy (the `_sympy_` protocol), so sympy.sympify(node), # sympy.Matrix(array_of_nodes), and sympy.det(J) work directly. sympy is an optional dependency; the @@ -146,21 +143,39 @@ def __getattr__(name): raise AttributeError("module {!r} has no attribute {!r}".format(__name__, name)) -# --- numpy-compatible elementwise helpers for the mp types (#298, #301) ------------------------ -# Vectorized real/imag/abs/conj/round/sum/norm/is_real that stay mp-native (numpy's ufuncs/attrs on -# the custom dtypes are unreliable -- see docs/source/known_gotchas.rst). These are attributes of the +# --- numpy-friendly elementwise helpers for the mp types (#298, #301) -------------------------- +# Vectorized real/imag/abs/conj/round/sum/norm/is_real that stay mp-native, riding the native numpy +# ufunc loops on mp-dtype arrays and working element-wise on lists/mixed input (see the +# "Multiprecision numbers and NumPy" docs page, docs/source/numpy.rst). These are attributes of the # top-level module (bertini.abs, bertini.real, ...). The builtin-shadowing names (abs, round, sum) are # deliberately kept OUT of __all__, so `from bertini import *` never clobbers the Python builtins. from . import _numpy_helpers as _numpy_helpers + +# make np.real/np.imag/np.angle raise (instead of silently returning wrong values) +# on plain mp-complex arrays -- numpy has no user-dtype hook for component access, +# and a crash is better than incorrect values. See _numpy_guard for the story. +from . import _numpy_guard as _numpy_guard +_numpy_guard.install() + real = _numpy_helpers.real imag = _numpy_helpers.imag conj = _numpy_helpers.conj +arg = _numpy_helpers.arg norm = _numpy_helpers.norm is_real = _numpy_helpers.is_real abs = _numpy_helpers.abs # noqa: A001 (bertini.abs; not exported via *) round = _numpy_helpers.round # noqa: A001 sum = _numpy_helpers.sum # noqa: A001 +# --- the one-stop math vocabulary --------------------------------------------------------------- +# `from bertini.operators import *` gives sin/cos/.../abs/arg/real/imag/... that work on symbolic +# expressions AND numbers AND numpy containers alike, dispatching per argument. The top-level +# elementary functions are rebound to the polymorphic versions (a strict superset of the symbolic +# ones bound above): bertini.sin(x) works for a Variable, a real_mp, or an array. +from . import operators +from .operators import (sin, cos, tan, asin, acos, atan, exp, log, sqrt, # noqa: F811 + sinh, cosh, tanh, asinh, acosh, atanh) + # https://stackoverflow.com/questions/44834/what-does-all-mean-in-python @@ -170,7 +185,7 @@ def __getattr__(name): 'jacobian','random_matrix','random_vector','random_real','random_complex','coefficient','coefficients', 'complex_mp','real_mp','int_mp','rational_mp', 'nag_algorithm','default_precision','is_distinct_up_to', - 'real','imag','conj','norm','is_real', + 'real','imag','conj','arg','norm','is_real', 'tracking','endgame','logging','symbolics','parse','multiprec','random','parallel', 'operators', # everyday classes hoisted to the top level @@ -181,6 +196,7 @@ def __getattr__(name): # symbolic constants 'E','Pi','I', 'sin','cos','tan','asin','acos','atan','exp','log','sqrt', + 'sinh','cosh','tanh','asinh','acosh','atanh', 'canonicalize','monomial_order'] diff --git a/python/bertini/_numpy_guard.py b/python/bertini/_numpy_guard.py new file mode 100644 index 000000000..e5b47b24e --- /dev/null +++ b/python/bertini/_numpy_guard.py @@ -0,0 +1,111 @@ +# This file is part of Bertini 2. +# +# python/bertini/_numpy_guard.py is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# python/bertini/_numpy_guard.py is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this file. If not, see . +# +# Copyright(C) Bertini2 Development Team + +"""Make numpy's component accessors fail LOUDLY on multiprecision complex arrays. + +numpy's ``ndarray.real`` / ``.imag`` are hardwired to its three built-in complex types +(``PyArray_ISCOMPLEX`` in ``getset.c``): on any other dtype, ``.real`` returns the array +itself and ``.imag`` returns zeros -- **silently wrong** for ``complex_mp``, and there is +no user-dtype hook to fix or even detect it at the numpy level. ``np.real`` / ``np.imag`` +/ ``np.angle`` are thin wrappers over those attributes and inherit the lie. + +The attributes themselves are C-level and untouchable, but the module *functions* are +plain Python -- so importing bertini wraps them: called on a plain ndarray of the +multiprecision complex dtype they raise ``TypeError`` (wrong-by-construction is converted +to loud), and every other input passes straight through to the original numpy functions. +A crash is better than incorrect values. + +Correct spellings, always: + +* ``bertini.real(x)`` / ``bertini.imag(x)`` / ``bertini.multiprec.arg(x)`` +* solution points (:class:`bertini.records.Solution`) override ``.real`` / ``.imag`` + at the subclass level and are simply correct -- the guard lets them through. + +The remaining untouchable spelling is the raw ``.real`` / ``.imag`` attribute on a plain +ndarray you built yourself; see the "Multiprecision numbers and NumPy" docs page. +""" + +import functools as _functools + +import numpy as _np + +from bertini.multiprec import complex_mp as _complex_mp + +_CPLX_MP = _np.dtype(_complex_mp) + +_MSG = ( + "numpy.{name}() cannot work on an array of the multiprecision complex dtype: " + "numpy's component access is hardwired to its built-in complex types and would " + "silently return WRONG values for user dtypes (there is no hook for bertini to fix " + "it). Use bertini.real(x) / bertini.imag(x) / bertini.multiprec.arg(x) instead -- " + "or hold a bertini Solution, whose .real/.imag are correct." +) + + +def _is_plain_complex_mp_array(x): + # exact-type check: subclasses (e.g. bertini's Solution) override .real/.imag + # correctly and must pass through + return type(x) is _np.ndarray and x.dtype == _CPLX_MP + + +def _would_lie(val): + if _is_plain_complex_mp_array(val): + return True + # a list/tuple of complex_mp converts to a plain mp-dtype array INSIDE numpy's + # real()/imag(), landing on the same wrong attribute path + if isinstance(val, (list, tuple)): + try: + return _is_plain_complex_mp_array(_np.asanyarray(val)) + except Exception: + return False + return False + + +def _is_complex_mp_anything(val): + # np.angle never consults .real/.imag attributes: it branches on the DTYPE and + # would die in arctan2 with a cryptic error for every mp-complex input -- + # Solutions included, since the subclass property cannot help it. Catch them + # all and say what to use instead. + if isinstance(val, _complex_mp): + return True + try: + return _np.asanyarray(val).dtype == _CPLX_MP + except Exception: + return False + + +def _guarded(orig, name, applies): + @_functools.wraps(orig) + def wrapper(val, *args, **kwargs): + if applies(val): + raise TypeError(_MSG.format(name=name)) + return orig(val, *args, **kwargs) + + wrapper._bertini_guarded = True + wrapper._bertini_original = orig + return wrapper + + +def install(): + """Wrap ``np.real`` / ``np.imag`` / ``np.angle``. Idempotent.""" + for name, applies in (("real", _would_lie), + ("imag", _would_lie), + ("angle", _is_complex_mp_anything)): + orig = getattr(_np, name) + if getattr(orig, "_bertini_guarded", False): + continue + setattr(_np, name, _guarded(orig, name, applies)) diff --git a/python/bertini/_numpy_helpers.py b/python/bertini/_numpy_helpers.py index 6ea7b81d9..70d8aad81 100644 --- a/python/bertini/_numpy_helpers.py +++ b/python/bertini/_numpy_helpers.py @@ -15,30 +15,45 @@ # # Copyright(C) Bertini2 Development Team -"""Vectorized element-wise helpers for the multiprecision types (issues #298, #301). +"""Vectorized helpers for the multiprecision types (issues #298, #301). -NumPy's ``.real`` / ``np.abs`` / ``np.round`` / identity-seeded reductions are unreliable on the -custom ``real_mp`` / ``complex_mp`` dtypes -- that boundary cannot be patched in the bindings (see -``docs/source/known_gotchas.rst``). These helpers do the elementwise work themselves, over a scalar, -a list, or a numpy array, and return **mp-native** results (a numpy object array for array input), so -your values stay arbitrary-precision instead of collapsing to float64. +The mp dtypes have full native numpy support (ufuncs, reductions, sorting -- see the +"Multiprecision numbers and NumPy" docs page); these helpers are the everyday sugar on +top. They accept a scalar, a list, or a numpy array, and return **mp-native** results, +so values stay arbitrary-precision instead of collapsing to float64. On an mp-dtype +ndarray they ride the native numpy loops; on lists/tuples and mixed input they work +element-wise. bertini.real(pt) # real parts, as real_mp bertini.abs(pt) # magnitudes, as real_mp bertini.is_real(pt) # is every coordinate real (imag within tol)? + +They also paper over the one true numpy limitation for user dtypes: the ndarray +``.real`` / ``.imag`` attributes return silently wrong values on complex_mp arrays +(numpy cannot know a user dtype is complex-like), so component access must go through +these helpers or ``bertini.multiprec.real/imag/arg``. """ +import builtins as _builtins import numpy as _np from decimal import Decimal as _Decimal, ROUND_HALF_EVEN as _ROUND_HALF_EVEN +import bertini.multiprec as _mp from bertini.multiprec import real_mp as _real_mp, complex_mp as _complex_mp from bertini.multiprec import abs as _mp_abs, conj as _mp_conj +_MP_DTYPES = (_np.dtype(_real_mp), _np.dtype(_complex_mp)) + def _is_container(x): return isinstance(x, (list, tuple)) or (isinstance(x, _np.ndarray) and x.ndim > 0) +def _is_mp_array(x): + """A numpy array of one of the mp dtypes -- eligible for the native ufunc loops.""" + return isinstance(x, _np.ndarray) and x.ndim > 0 and x.dtype in _MP_DTYPES + + def _elementwise(fn, x): """Apply scalar ``fn`` over a scalar / list / numpy array, returning mp-native results (a numpy object array, preserving shape, for container input).""" @@ -74,7 +89,9 @@ def _imag_scalar(v): def _abs_scalar(v): if isinstance(v, (_complex_mp, _real_mp)): return _mp_abs(v) - return abs(v) + # explicitly the BUILTIN: this module's own `abs` shadows it at module scope, + # and the bare name recursed infinitely for plain python input + return _builtins.abs(v) def _conj_scalar(v): @@ -96,31 +113,69 @@ def _round_scalar(v, decimals): return _complex_mp(_round_real_mp(v.real, decimals), _round_real_mp(v.imag, decimals)) if isinstance(v, _real_mp): return _round_real_mp(v, decimals) - return round(v, decimals) + return _builtins.round(v, decimals) # explicitly the builtin (module `round` shadows it) # --- the public helpers ----------------------------------------------------------------------- +# each takes the native numpy path when handed an mp-dtype array (the C++ loops -- +# fast, and the result is a proper mp-dtype array that keeps working with sort, +# reductions, and friends), and falls back to element-wise work for lists and +# mixed input. def real(x): - """Real part(s), as ``real_mp`` -- over a scalar / list / array (replaces numpy ``.real``).""" + """Real part(s), as ``real_mp`` -- over a scalar / list / array (replaces numpy ``.real``, + which returns silently wrong values on complex_mp arrays).""" + if _is_mp_array(x) and x.ndim == 1: + if x.dtype == _np.dtype(_complex_mp): + return _mp.real(x) + return x.copy() return _elementwise(_real_scalar, x) def imag(x): - """Imaginary part(s), as ``real_mp`` -- over a scalar / list / array (replaces numpy ``.imag``).""" + """Imaginary part(s), as ``real_mp`` -- over a scalar / list / array (replaces numpy ``.imag``, + which returns silently wrong values on complex_mp arrays).""" + if _is_mp_array(x) and x.ndim == 1: + if x.dtype == _np.dtype(_complex_mp): + return _mp.imag(x) + return _np.zeros(x.shape, dtype=_real_mp) return _elementwise(_imag_scalar, x) def abs(x): - """Magnitude(s), as ``real_mp`` -- over a scalar / list / array (replaces ``np.abs``).""" + """Magnitude(s), as ``real_mp`` -- over a scalar / list / array (same as ``np.abs``).""" + if _is_mp_array(x): + return _np.abs(x) return _elementwise(_abs_scalar, x) def conj(x): - """Complex conjugate(s) -- over a scalar / list / array.""" + """Complex conjugate(s) -- over a scalar / list / array (same as ``np.conj``).""" + if _is_mp_array(x): + return _np.conj(x) return _elementwise(_conj_scalar, x) +def _arg_scalar(v): + if isinstance(v, _complex_mp): + return _mp.arg(v) + if isinstance(v, _real_mp): + return _mp.arg(_complex_mp(v)) + import cmath + return cmath.phase(complex(v)) + + +def arg(x): + """Argument(s) -- the angle from 0 -- as ``real_mp``, over a scalar / list / array + (the ``np.angle`` replacement; numpy's own cannot work on mp dtypes). Beware the + branch cut.""" + if _is_mp_array(x) and x.ndim == 1: + if x.dtype == _np.dtype(_real_mp): + x = x.astype(_np.dtype(_complex_mp)) # registered safe cast, exact + return _mp.arg(x) + return _elementwise(_arg_scalar, x) + + def round(x, decimals=0): """Round to ``decimals`` places, staying mp-native -- over a scalar / list / array (``np.round``).""" return _elementwise(lambda v: _round_scalar(v, decimals), x) @@ -133,12 +188,19 @@ def is_real(point, tol=1e-10): just_real = [pt for pt in solutions if bertini.is_real(pt)] """ + if _is_mp_array(point) and point.ndim == 1: + if point.dtype == _np.dtype(_real_mp): + return True + # mp-vs-float tolerance orderings are exact and native + return bool(_np.all(_np.abs(_mp.imag(point)) < tol)) flat = _np.asarray(point, dtype=object).reshape(-1) return all(float(_abs_scalar(_imag_scalar(c))) < tol for c in flat) def sum(x): - """Sum of a 1-D collection, staying mp-native (sidesteps numpy's identity-reduction gotcha).""" + """Sum of a collection, staying mp-native.""" + if _is_mp_array(x) and x.size > 0: + return _np.sum(x) flat = list(_np.asarray(x, dtype=object).reshape(-1)) if not flat: return 0 @@ -150,6 +212,9 @@ def sum(x): def norm(x): """Euclidean (2-)norm of a 1-D collection, as ``real_mp`` -- ``sqrt(sum |x_i|^2)``.""" + if _is_mp_array(x) and x.ndim == 1 and x.size > 0: + a = _np.abs(x) # magnitudes, as real_mp + return _mp.sqrt(_np.dot(a, a)) # dot slot + scalar sqrt, all mp flat = _np.asarray(x, dtype=object).reshape(-1) acc = None for v in flat: diff --git a/python/bertini/multiprec/__init__.py b/python/bertini/multiprec/__init__.py index 4504c670b..39136cd1a 100644 --- a/python/bertini/multiprec/__init__.py +++ b/python/bertini/multiprec/__init__.py @@ -33,12 +33,12 @@ """ Multiprecision types, and functions that operate on them. -Numeric types exposed are +Numeric types exposed are -* Complex (Boost.Multiprecision mpc) -* Float (Boost.Multiprecision mpfr) -* Int (Boost.Multiprecision mpz) -* Rational (Boost.Multiprecision.mpq) +* complex_mp (Boost.Multiprecision mpc) +* real_mp (Boost.Multiprecision mpfr) +* int_mp (Boost.Multiprecision mpz) +* rational_mp (Boost.Multiprecision.mpq) This namespace also includes the mathematical operators, like `cos`, etc. """ @@ -48,7 +48,11 @@ from bertini._pybertini.multiprec import * # (no Vector helper: eigenpy makes the mp number types work as numpy dtypes directly, so a plain -# numpy array -- e.g. np.zeros(n, dtype=bertini.complex_mp) -- is the vector.) +# numpy array -- e.g. np.zeros(n, dtype=bertini.complex_mp) -- is the vector. numpy ufuncs +# (np.abs, np.exp, np.sum, ...) work on such arrays; see the "Multiprecision numbers and NumPy" +# docs page. For the real/imaginary parts or argument of a COMPLEX ARRAY use this module's +# real()/imag()/arg() -- the ndarray .real/.imag attributes and np.angle return silently wrong +# values for user-defined dtypes, a numpy limitation.) __all__ = dir(_pybmp) diff --git a/python/bertini/operators.py b/python/bertini/operators.py index af44d7853..bbdb5e646 100644 --- a/python/bertini/operators.py +++ b/python/bertini/operators.py @@ -19,21 +19,112 @@ # as well as COPYING. Bertini2 is provided with permitted # additional terms in the b2/licenses/ directory. -"""The symbolic math vocabulary, gathered for ``from bertini.operators import *``. +"""The whole math vocabulary in one namespace -- symbols and numbers alike. -These are the elementary functions (``sin``, ``cos``, ...) and constants (``E``, ``Pi``, ``I``) used -to *build* symbolic systems on :class:`~bertini.Variable`\\ s. They also live at the top level -(``bertini.sin``, ``bertini.Pi``, ...); this module exists only so you can pull the math vocabulary -into your namespace without importing the rest of ``bertini``:: +One import gives you functions that work on *everything*: a symbolic +:class:`~bertini.Variable`/expression, a multiprecision number, a numpy array of them, +or a plain python number:: from bertini.operators import * - f = sin(x) + Pi*y - E -These are the *symbolic* operators; the numeric elementary functions (acting on multiprecision -numbers rather than expression nodes) live in :mod:`bertini.multiprec`. + f = sin(x) + Pi*y - E # symbolic (x, y Variables -> an expression) + v = sin(real_mp('0.5')) # numeric (full precision) + m = abs(solutions[0]) # numpy arrays (multiprecision dtypes included) + t = arg(complex_mp(1, 1)) # components (arg/real/imag/conj) + +Dispatch is by argument: a function-tree node builds a symbolic node; everything else +takes the numeric path (multiprecision scalars and mp-dtype numpy arrays go through the +native precision-preserving loops; python numbers and float arrays are plain numpy). +You never have to remember whether a name lives in ``bertini.multiprec`` or at the top +level -- it is here. + +The functions with no symbolic counterpart (``abs``, ``arg``, ``real``, ``imag``, +``conj``, ``round``, ``sum``, ``norm``, ``is_real``, and the hyperbolics) raise a clear +``TypeError`` when handed a symbolic expression. + +``abs``, ``round``, and ``sum`` shadow the python builtins **within your namespace** +when you star-import this module -- that is the point (they fall back to builtin +behavior on plain python input), but it is opt-in: ``from bertini import *`` never +shadows builtins. """ -# Re-exported from the top-level bertini package (defined there before this module is imported). -from . import sin, cos, tan, asin, acos, atan, exp, log, sqrt, E, Pi, I # noqa: F401 +import numpy as _np + +from bertini._pybertini.function_tree import AbstractNode as _AbstractNode + +from . import symbolics as _sym +from . import _numpy_helpers as _nh + +# the symbolic constants, ready to drop into expressions +from . import E, Pi, I # noqa: F401 + + +def _polymorphic(name, sym_fn, num_fn, doc): + def f(x): + if isinstance(x, _AbstractNode): + return sym_fn(x) + return num_fn(x) + f.__name__ = name + f.__qualname__ = name + f.__doc__ = doc + ("\n\nPolymorphic: builds a symbolic node for a function-tree " + "argument, computes numerically (precision-preserving, numpy " + "containers included) for everything else.") + return f + + +def _numeric_only(name, num_fn, doc): + def f(x, *args, **kwargs): + if isinstance(x, _AbstractNode): + raise TypeError( + f"{name}() is not defined for symbolic expressions -- it is a numeric " + "operation. Evaluate the expression first, or use the symbolic " + "functions (sin, cos, exp, ...) to build systems.") + return num_fn(x, *args, **kwargs) + f.__name__ = name + f.__qualname__ = name + f.__doc__ = doc + ("\n\nNumeric: multiprecision scalars, numpy arrays (mp dtypes " + "included), lists, and plain python numbers.") + return f + + +# --- the elementary functions with symbolic twins: full polymorphic dispatch ------------------- + +sin = _polymorphic('sin', _sym.sin, _np.sin, "Sine.") +cos = _polymorphic('cos', _sym.cos, _np.cos, "Cosine.") +tan = _polymorphic('tan', _sym.tan, _np.tan, "Tangent.") +asin = _polymorphic('asin', _sym.asin, _np.arcsin, "Arcsine.") +acos = _polymorphic('acos', _sym.acos, _np.arccos, "Arccosine.") +atan = _polymorphic('atan', _sym.atan, _np.arctan, "Arctangent.") +exp = _polymorphic('exp', _sym.exp, _np.exp, "Exponential, base e.") +log = _polymorphic('log', _sym.log, _np.log, "Natural logarithm.") +sqrt = _polymorphic('sqrt', _sym.sqrt, _np.sqrt, "Square root.") + +# --- numeric-only elementary functions (no symbolic node exists) ------------------------------- + +sinh = _numeric_only('sinh', _np.sinh, "Hyperbolic sine.") +cosh = _numeric_only('cosh', _np.cosh, "Hyperbolic cosine.") +tanh = _numeric_only('tanh', _np.tanh, "Hyperbolic tangent.") +asinh = _numeric_only('asinh', _np.arcsinh, "Hyperbolic arcsine.") +acosh = _numeric_only('acosh', _np.arccosh, "Hyperbolic arccosine.") +atanh = _numeric_only('atanh', _np.arctanh, "Hyperbolic arctangent.") + +# --- components, magnitudes, and friends (numeric-only) ---------------------------------------- + +abs = _numeric_only('abs', _nh.abs, "Magnitude(s), as real_mp for mp input.") # noqa: A001 +arg = _numeric_only('arg', _nh.arg, "Argument(s) (angle from 0), as real_mp. Beware the branch cut.") +real = _numeric_only('real', _nh.real, "Real part(s), as real_mp for mp input.") +imag = _numeric_only('imag', _nh.imag, "Imaginary part(s), as real_mp for mp input.") +conj = _numeric_only('conj', _nh.conj, "Complex conjugate(s).") +round = _numeric_only('round', _nh.round, "Round to N DECIMAL digits, staying mp-native.") # noqa: A001 +sum = _numeric_only('sum', _nh.sum, "Sum of a collection, staying mp-native.") # noqa: A001 +norm = _numeric_only('norm', _nh.norm, "Euclidean (2-)norm, as real_mp.") +is_real = _numeric_only('is_real', _nh.is_real, "Is every coordinate real (|imag| < tol)?") + -__all__ = ['sin', 'cos', 'tan', 'asin', 'acos', 'atan', 'exp', 'log', 'sqrt', 'E', 'Pi', 'I'] +__all__ = [ + 'sin', 'cos', 'tan', 'asin', 'acos', 'atan', + 'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh', + 'exp', 'log', 'sqrt', + 'abs', 'arg', 'real', 'imag', 'conj', 'round', 'sum', 'norm', 'is_real', + 'E', 'Pi', 'I', +] diff --git a/python/bertini/records.py b/python/bertini/records.py index 492c7340e..bd24e4037 100644 --- a/python/bertini/records.py +++ b/python/bertini/records.py @@ -144,6 +144,32 @@ def __array_wrap__(self, out_arr, context=None, return_scalar=False): result.annotations = {} return result + # numpy's ndarray .real/.imag are hardwired to its built-in complex types: on a + # multiprecision-complex array the base attributes return silently WRONG values + # (.real gives the complex values, .imag gives zeros), with no hook for user + # dtypes. A python subclass property CAN shadow the C-level attribute, so + # solution points -- the mp-complex arrays users actually hold -- are correct. + # Plain ndarrays are covered by the guarded np.real/np.imag/np.angle (see + # bertini._numpy_guard) and by bertini.real/imag. + + @property + def real(self): + """The real parts -- correct also for the multiprecision complex dtype.""" + from bertini.multiprec import complex_mp + if self.dtype == _np.dtype(complex_mp): + from bertini import _numpy_helpers as _nh + return _nh.real(_np.asarray(self)) + return _np.ndarray.real.__get__(self) + + @property + def imag(self): + """The imaginary parts -- correct also for the multiprecision complex dtype.""" + from bertini.multiprec import complex_mp + if self.dtype == _np.dtype(complex_mp): + from bertini import _numpy_helpers as _nh + return _nh.imag(_np.asarray(self)) + return _np.ndarray.imag.__get__(self) + class SolveResult: """What ``solve`` returns: the solutions plus a claim ticket on the recorded run. diff --git a/python/docs/source/index.rst b/python/docs/source/index.rst index 2b6078af2..3479634fd 100644 --- a/python/docs/source/index.rst +++ b/python/docs/source/index.rst @@ -23,7 +23,7 @@ The Python bindings for Bertini 2 .. toctree:: :maxdepth: 2 - known_gotchas + numpy detailed/api zbib diff --git a/python/docs/source/known_gotchas.rst b/python/docs/source/known_gotchas.rst deleted file mode 100644 index b0b7c2033..000000000 --- a/python/docs/source/known_gotchas.rst +++ /dev/null @@ -1,106 +0,0 @@ -⚠️ Known gotchas -***************************** - -A few sharp edges fall out of how Bertini 2's multiprecision numbers are exposed to -NumPy. They are collected here with the idiom that works and the idiom that bites. - -.. _gotcha-numpy-reductions: - -NumPy reductions over multiprecision arrays -============================================ - -Bertini 2 exposes :class:`~bertini.real_mp` (variable-precision real) and -:class:`~bertini.complex_mp` (variable-precision complex) as **custom NumPy -dtypes** (via eigenpy). Element-wise math, indexing, ``@`` / :func:`numpy.dot`, -:func:`numpy.linalg.norm`, :func:`numpy.cumsum` and friends all work on arrays of these -dtypes. - -The sharp edge is the **identity-seeded reductions** -- :func:`numpy.sum`, -:func:`numpy.prod`, :func:`numpy.mean`, and a bare ``ufunc.reduce`` with no ``initial=``: - -.. code-block:: python - - >>> import numpy as np - >>> from bertini.multiprec import complex_mp - >>> v = np.array([complex_mp(3), complex_mp(4)], dtype=complex_mp) - >>> np.sum(v) # ✗ DON'T -- may crash - Traceback (most recent call last): - ... - SystemError: returned NULL without setting an exception - -Why this happens ----------------- - -To reduce an array, NumPy needs the reduction's **identity element** -- ``0`` for a sum, -``1`` for a product -- materialized *in the array's dtype*. For these custom dtypes, -NumPy (2.x) tries to build that identity from a Python ``int`` and fails *inside its own -machinery*, before ever calling Bertini's dtype code, returning ``NULL`` without setting a -Python exception (hence the bare ``SystemError``). It is a limitation of the NumPy ↔ -eigenpy boundary for legacy user-defined dtypes, **not** something Bertini can patch in its -bindings, and whether it triggers depends on the exact NumPy / eigenpy / Eigen build -- so -it may "work" on one machine and crash on another. Treat it as always-unsupported. - -What to do instead ------------------- - -Give the reduction a value to start from, or use an operation that seeds itself from the -data (``dot`` / ``norm`` / a plain Python loop). All of these are stable across builds: - -.. doctest:: - - >>> import numpy as np - >>> from bertini.multiprec import real_mp, complex_mp - >>> v = np.array([complex_mp(3), complex_mp(4)], dtype=complex_mp) - >>> w = np.array([real_mp(1), real_mp(2), real_mp(3)], dtype=real_mp) - - >>> # ✓ sum: hand ufunc.reduce an explicit, correctly-typed identity - >>> complex(np.add.reduce(v, initial=complex_mp(0))) - (7+0j) - >>> float(np.add.reduce(w, initial=real_mp(0))) - 6.0 - - >>> # ✓ sum: or just use Python's built-in sum() - >>> float(sum(w)) - 6.0 - - >>> # ✓ Euclidean norm works directly (it routes through the dtype's dot slot) - >>> float(np.linalg.norm(v)) - 5.0 - - >>> # ✓ sum of squares: np.dot does not conjugate, so dot(v, v) == sum(v**2) - >>> complex(np.dot(v, v)) - (25+0j) - - >>> # ✓ mean: reduce with an identity, then divide by the count - >>> float(np.add.reduce(w, initial=real_mp(0)) / w.size) - 2.0 - -In short: anywhere you would reach for ``np.sum(a)`` / ``np.prod(a)`` / ``np.mean(a)`` on a -``real_mp`` or ``complex_mp`` array, reach for ``np.add.reduce(a, initial=...)`` (or -``np.multiply.reduce(a, initial=...)``), ``np.dot`` / ``np.linalg.norm``, or a plain Python -``sum`` instead. - -The bertini helpers do this for you ------------------------------------ - -The same boundary makes ``arr.real`` / ``np.abs(arr)`` / ``np.round(arr)`` unreliable on these -dtypes. So bertini ships elementwise helpers that operate over a scalar, a list, or a numpy array and -return **mp-native** results (never a float64 collapse): - -.. code-block:: python - - import bertini - - bertini.real(pt) # real parts, as real_mp (replaces arr.real) - bertini.imag(pt) # imaginary parts, as real_mp - bertini.abs(pt) # magnitudes, as real_mp (replaces np.abs) - bertini.conj(pt) # complex conjugates - bertini.round(pt, 8) # rounded, staying mp (replaces np.round) - bertini.sum(pt) # sum, staying mp (sidesteps the reduction gotcha above) - bertini.norm(pt) # Euclidean norm, as real_mp - bertini.is_real(pt) # is every coordinate real (|imag| < tol)? -> bool - -They live at the top level (``bertini.abs``, ...); the builtin-shadowing names (``abs``, ``round``, -``sum``) are deliberately kept out of ``from bertini import *``, so a star-import never clobbers the -Python builtins. For the tolerance point comparison behind de-duplication, see -:func:`bertini.is_distinct_up_to`. diff --git a/python/docs/source/numpy.rst b/python/docs/source/numpy.rst new file mode 100644 index 000000000..a587534b1 --- /dev/null +++ b/python/docs/source/numpy.rst @@ -0,0 +1,164 @@ +Multiprecision numbers and NumPy +******************************** + +Bertini 2 exposes :class:`~bertini.real_mp` (variable-precision real) and +:class:`~bertini.complex_mp` (variable-precision complex) as **native NumPy dtypes** +(via eigenpy). A plain numpy array *is* the multiprecision vector -- no wrapper types, +no conversion step -- and the numpy you already know works on it: + +* element-wise arithmetic, comparisons, ``@`` / :func:`numpy.dot` / + :func:`numpy.linalg.norm`, +* the ufunc family: :func:`numpy.abs`, :func:`numpy.conj`, ``exp`` / ``log`` / + the trigonometric and hyperbolic functions and their inverses, ``power``, ``sign``, + ``sqrt``, the rounding family (``floor`` / ``ceil`` / ``trunc`` / ``rint``, with + ``rint`` / ``round`` on complex rounding each component), ``minimum`` / ``maximum``, + ``isnan`` / ``isinf`` / ``isfinite``, and friends, +* sorting and order statistics on the real type: :func:`numpy.sort`, + :func:`numpy.argsort`, :func:`numpy.searchsorted`, :func:`numpy.argmax` / + :func:`numpy.argmin`, :func:`numpy.median`, +* reductions: :func:`numpy.sum`, :func:`numpy.prod`, :func:`numpy.mean`, + :func:`numpy.cumsum`, ``ufunc.reduce``. + +The contract behind all of it: every ufunc loop calls the same multiprecision function +that the scalar :mod:`bertini.multiprec` function of the same name binds, at the +operands' precision. ``np.exp(a)[i]`` is exactly ``mp.exp(a[i])`` -- numpy is a fast +way to say the same thing, never a detour through float64. + +.. _numpy-float64-boundary: + +Your digits are protected: the float64 boundary +================================================ + +A Python ``float`` is a 53-bit binary number: ``0.1`` as a float is *not* the decimal +``0.1`` you typed. So Bertini draws a deliberate line -- a float64 **value** never +silently enters a multiprecision computation. Values come from strings +(``real_mp('0.1')``) or integers (exact), and high-precision results only drop to +double when you explicitly ask (``float(x)``, ``complex(z)``, ``arr.astype(float)``, +``arr.astype(complex)``). + +**Comparisons are a different story** -- checking a residual against a tolerance never +pollutes anything, so it just works, floats and all: + +.. doctest:: + + >>> import numpy as np + >>> from bertini.multiprec import real_mp, complex_mp + >>> a = np.array([complex_mp(3), complex_mp(4)]) + >>> b = np.array([complex_mp(3), complex_mp(4)]) + + >>> # the closeness check: compare against a double tolerance, get bools back + >>> bool(np.all(np.abs(a - b) < 1e-10)) + True + + >>> # the comparison is exact -- 1e-22 is not "rounded away" against 1e-30 + >>> bool(np.all(np.array([real_mp('1e-22')]) < 1e-30)) + False + +This mirrors the C++ solvers, whose tolerances are doubles too. + +Where the boundary shows as an error, that is it doing its job: + +* ``a + 0.1`` raises (``ufunc ... not supported``) -- a float value would have entered + the computation. Say ``a + real_mp('0.1')``. +* ``a == 0.1`` raises -- exact equality against a float literal is the trap the + boundary exists for. Compare against ``real_mp('0.1')``, or use a tolerance. +* ``np.isclose`` / ``np.allclose`` raise ``DTypePromotionError`` -- they *compute* + ``atol + rtol*np.abs(b)`` with float64 internally. Use the one-liner above. +* ``np.round(a, decimals)`` with nonzero ``decimals`` fails (it scales by a float + power of ten internally); plain ``np.round(a)`` / :func:`numpy.rint` work, and + :func:`bertini.round` does decimal-digit rounding at full precision. + +.. _numpy-complex-components: + +Component access on complex arrays +=================================== + +The one place numpy itself cannot be taught about a user-defined dtype: the ndarray +attributes ``.real`` and ``.imag``. numpy hardwires them to its own built-in complex +types -- on any other dtype ``.real`` returns the array itself and ``.imag`` returns +zeros, with no hook for the bindings to fix or even detect it. + +Bertini defends every spelling it can reach: + +* **Solution points are simply correct.** The arrays returned by ``bertini.solve`` + override ``.real`` / ``.imag`` at the subclass level, so ``sol.real`` / ``sol.imag`` + give the true components at full precision. +* **The numpy functions fail loudly.** Importing bertini wraps :func:`numpy.real` / + :func:`numpy.imag` / :func:`numpy.angle`: on a plain mp-complex array they raise a + ``TypeError`` naming the right tool, instead of silently returning wrong values. + All other inputs pass straight through to numpy. + +.. warning:: + + The raw ``.real`` / ``.imag`` **attributes on a plain ndarray** are the one + spelling nothing can reach: on a ``complex_mp`` array you built yourself (not a + Solution), they return wrong values silently. Complex *scalars* are fine -- + ``w[0].imag`` is exact. + +The component accessors -- same results, full precision, array-capable: + +.. doctest:: + + >>> import numpy as np + >>> import bertini.multiprec as mp + >>> from bertini.multiprec import complex_mp + >>> w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + + >>> [float(x) for x in mp.real(w)] + [1.0, 3.0] + >>> [float(x) for x in mp.imag(w)] + [2.0, 4.0] + + >>> # the np.angle equivalent + >>> [round(float(x), 4) for x in mp.arg(w)] + [1.1071, 0.9273] + +Convenience helpers +==================== + +For everyday work with solution points, the top-level helpers accept a scalar, a list, +or a numpy array, and return multiprecision-native results: + +.. code-block:: python + + import bertini + + bertini.real(pt) # real parts, as real_mp + bertini.imag(pt) # imaginary parts, as real_mp + bertini.abs(pt) # magnitudes, as real_mp + bertini.conj(pt) # complex conjugates + bertini.round(pt, 8) # rounded to 8 DECIMAL digits, staying mp + bertini.sum(pt) # sum, staying mp + bertini.norm(pt) # Euclidean norm, as real_mp + bertini.is_real(pt) # is every coordinate real (|imag| < tol)? -> bool + +On multiprecision-dtype arrays these ride the native numpy loops; on lists and mixed +input they work element-wise. The builtin-shadowing names (``abs``, ``round``, +``sum``) are deliberately kept out of ``from bertini import *``, so a star-import +never clobbers the Python builtins. For the tolerance point comparison behind +de-duplication, see :func:`bertini.is_distinct_up_to`. + +And for the whole vocabulary in one go -- the elementary functions *and* these +helpers, working on symbolic expressions, multiprecision numbers, and numpy +containers alike, dispatched per argument:: + + from bertini.operators import * + + f = sin(x) + Pi*y # symbolic (x, y Variables) + v = sin(real_mp('0.5')) # numeric, full precision + m = abs(solutions[0]) # numpy containers, mp dtypes included + t = arg(complex_mp(1, 1)) # components: arg/real/imag/conj + +This star-import *does* shadow ``abs``/``round``/``sum`` in your namespace -- that is +its point (they fall back to builtin behavior on plain python input), and it is +opt-in. + +Boundaries by design +===================== + +* **Complex numbers are unordered.** ``complex_mp`` arrays do not support + :func:`numpy.sort`, :func:`numpy.argmax`, ``minimum`` / ``maximum``, or ``<``. + (numpy's ``complex128`` sorts lexicographically for historical reasons; Bertini + does not reproduce that.) Sort a derived real quantity: ``np.argsort(np.abs(w))``. +* **No casts to integer types.** ``arr.astype(int)`` fails; go through double + explicitly if you truly want it. diff --git a/python/docs/source/tutorials/settings_and_precision/precision_models/index.rst b/python/docs/source/tutorials/settings_and_precision/precision_models/index.rst index 9bb6bff2c..5feb89a33 100644 --- a/python/docs/source/tutorials/settings_and_precision/precision_models/index.rst +++ b/python/docs/source/tutorials/settings_and_precision/precision_models/index.rst @@ -41,8 +41,8 @@ Reading solutions: convert with ``complex()`` The model changes the **type** of the numbers you get back. A double solve returns NumPy ``complex128``; a multiprecision solve returns arrays of :class:`bertini.complex_mp` (NumPy -arrays with ``dtype`` ``Complex``). The portable habit -- works in every model -- is to convert each -coordinate with :func:`complex`: +arrays with ``dtype`` ``complex_mp``). The portable habit -- works in every model -- is to convert +each coordinate with :func:`complex`: .. testcode:: diff --git a/python/test/classes/numpy_ufuncs_test.py b/python/test/classes/numpy_ufuncs_test.py new file mode 100644 index 000000000..75f72bbe4 --- /dev/null +++ b/python/test/classes/numpy_ufuncs_test.py @@ -0,0 +1,574 @@ +# This file is part of Bertini 2. +# +# python/test/classes/numpy_ufuncs_test.py is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# python/test/classes/numpy_ufuncs_test.py is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with python/test/classes/numpy_ufuncs_test.py. If not, see . +# +# Copyright(C) Bertini2 Development Team +# +# See for a copy of the license, +# as well as COPYING. Bertini2 is provided with permitted +# additional terms in the b2/licenses/ directory. + +# individual authors of this file include: +# +# silviana amethyst +# summer 2026 + +""" +Tests for the numpy ufunc coverage of the multiprecision dtypes. + +Every registered loop calls the same boost::multiprecision free function that +the ``bertini.multiprec`` scalar function of the same name binds, so a ufunc +applied to an array must agree element-wise with the scalar function -- +exactly, not approximately, since both run identical code at identical +precision. + +Also pins: +- the identity-seeded reductions (np.sum/np.prod/np.mean) that older docs + declared unsupported (they work on current numpy; regression-guard them), +- the sort/argsort/searchsorted/argmax family (dtype compare/arg slots), +- numpy semantics corners: rint half-to-even, remainder sign-of-divisor, + nan propagation in minimum/maximum vs fmin/fmax, +- precision preservation through every loop shape, +- unwritten-slot safety (np.empty) per the ADR-0006 guard doctrine, +- the mp.imag copy-paste regression (returned the real part), and the array + overloads of multiprec real/imag/arg that replace the silently-wrong + ndarray .real/.imag attributes and np.angle. +""" + +import numpy as np +import pytest + +import bertini.multiprec as mp +from bertini.multiprec import complex_mp, real_mp + + +@pytest.fixture(params=[real_mp, complex_mp], ids=["real_mp", "complex_mp"]) +def dtype(request): + return request.param + + +def _sample(dtype): + """A small array of values safe for every domain-restricted function.""" + if dtype is real_mp: + return np.array([real_mp('0.25'), real_mp('0.5'), real_mp('0.75')]) + return np.array([complex_mp('0.25', '0.125'), + complex_mp('0.5', '-0.25'), + complex_mp('-0.75', '0.375')]) + + +# ufuncs defined for both dtypes, with the matching multiprec scalar function +UFUNCS_BOTH = [ + (np.exp, mp.exp), (np.log, mp.log), (np.sqrt, mp.sqrt), + (np.sin, mp.sin), (np.cos, mp.cos), (np.tan, mp.tan), + (np.arcsin, mp.asin), (np.arccos, mp.acos), (np.arctan, mp.atan), + (np.sinh, mp.sinh), (np.cosh, mp.cosh), (np.tanh, mp.tanh), + (np.arcsinh, mp.asinh), (np.arccosh, None), (np.arctanh, mp.atanh), +] + + +class TestElementwiseAgreesWithScalarFunctions: + """np.f(arr)[i] must equal mp.f(arr[i]) exactly (identical code path).""" + + @pytest.mark.parametrize( + "ufunc,scalar", UFUNCS_BOTH, + ids=[u.__name__ for u, _ in UFUNCS_BOTH]) + def test_transcendental(self, dtype, ufunc, scalar): + if ufunc is np.arccosh: + # acosh needs |x| >= 1 on the real line + v = (np.array([real_mp('1.5'), real_mp(2)]) if dtype is real_mp + else np.array([complex_mp('1.5', '0.5')])) + scalar = mp.acosh + else: + v = _sample(dtype) + out = ufunc(v) + assert out.dtype == np.dtype(dtype) + for got, x in zip(out, v): + assert got == scalar(x) + + def test_absolute(self, dtype): + v = _sample(dtype) + out = np.abs(v) + # output is ALWAYS real, also for complex input (the magnitude) + assert out.dtype == np.dtype(real_mp) + for got, x in zip(out, v): + assert got == mp.abs(x) + + def test_conjugate(self, dtype): + v = _sample(dtype) + out = np.conj(v) + assert out.dtype == np.dtype(dtype) + for got, x in zip(out, v): + assert got == (mp.conj(x) if dtype is complex_mp else x) + + def test_power(self, dtype): + v = _sample(dtype) + out = np.power(v, v) + for got, x in zip(out, v): + assert got == x ** x + + def test_reciprocal(self, dtype): + v = _sample(dtype) + out = np.reciprocal(v) + for got, x in zip(out, v): + assert got == dtype(1) / x + + def test_square_negative_positive(self, dtype): + v = _sample(dtype) + assert all(np.square(v)[i] == v[i] * v[i] for i in range(len(v))) + assert all(np.negative(v)[i] == -v[i] for i in range(len(v))) + assert all(np.positive(v)[i] == v[i] for i in range(len(v))) + + def test_real_only_transcendentals(self): + v = _sample(real_mp) + pairs = [(np.log10, mp.log), (np.exp2, None), (np.log2, None), + (np.expm1, None), (np.log1p, None), (np.cbrt, None)] + # spot-check values against the double versions loosely; the exact + # contract (same boost call as a scalar) has no bound scalar twin for + # these, so compare against float64 at double precision. + for ufunc, _ in pairs: + got = ufunc(v) + want = ufunc(np.array([float(x) for x in v])) + for g, w in zip(got, want): + assert abs(float(g) - w) < 1e-14, ufunc.__name__ + + def test_sign_real(self): + v = np.array([real_mp(-3), real_mp(0), real_mp(2)]) + assert [str(s) for s in np.sign(v)] == ['-1', '0', '1'] + + def test_sign_complex_is_unit_modulus(self): + # numpy-2 semantics: sign(z) = z/|z|, 0 at 0 + z = complex_mp(3, 4) + s = np.sign(np.array([z, complex_mp(0)])) + assert s[0] == complex_mp('0.6', '0.8') + assert s[1] == complex_mp(0) + + def test_arctan2_hypot_copysign(self): + # atan2/hypot use dedicated algorithms, so agreement with the composed + # formulas is to the last ulp, not bit-exact + tol = real_mp('1e-25') + a = np.array([real_mp(1), real_mp(-2)]) + b = np.array([real_mp(3), real_mp(4)]) + assert mp.abs(np.arctan2(a, b)[0] - mp.atan(a[0] / b[0])) < tol + assert mp.abs(np.hypot(a, b)[1] - mp.sqrt(a[1] * a[1] + b[1] * b[1])) < tol + got = np.copysign(a, np.array([real_mp(-1), real_mp(1)])) + assert [str(x) for x in got] == ['-1', '2'] + + +class TestNumpySemanticsCorners: + """Corners where numpy's semantics differ from naive C/boost calls.""" + + def test_rint_rounds_half_to_even(self): + # regression: boost's rint rounds half AWAY from zero; numpy (and the + # loop, via mpfr_rint in MPFR_RNDN) rounds half to even + v = np.array([real_mp('0.5'), real_mp('1.5'), real_mp('2.5'), + real_mp('-0.5'), real_mp('-2.5')]) + assert [str(x) for x in np.rint(v)] == ['0', '2', '2', '-0', '-2'] + + def test_rint_and_round_on_complex(self): + # numpy defines rint (and hence np.round) for complex, component-wise + w = np.array([complex_mp('1.5', '2.5'), complex_mp('-0.5', '3.4')]) + want = np.rint(np.array([1.5 + 2.5j, -0.5 + 3.4j])) + for got, ref in zip(np.rint(w), want): + assert complex(got) == ref + for got, ref in zip(np.round(w), want): + assert complex(got) == ref + + def test_floor_ceil_trunc(self): + v = np.array([real_mp('1.7'), real_mp('-1.7')]) + assert [str(x) for x in np.floor(v)] == ['1', '-2'] + assert [str(x) for x in np.ceil(v)] == ['2', '-1'] + assert [str(x) for x in np.trunc(v)] == ['1', '-1'] + + def test_remainder_takes_sign_of_divisor(self): + a = np.array([real_mp(7), real_mp(-7), real_mp(7), real_mp(-7)]) + b = np.array([real_mp(3), real_mp(3), real_mp(-3), real_mp(-3)]) + assert [str(x) for x in np.mod(a, b)] == ['1', '2', '-2', '-1'] + # fmod keeps C semantics (sign of dividend) + assert [str(x) for x in np.fmod(a, b)] == ['1', '-1', '1', '-1'] + assert [str(x) for x in np.floor_divide(a, b)] == ['2', '-3', '-3', '2'] + + def test_minimum_maximum_propagate_nan_fmin_fmax_ignore_it(self): + nan = real_mp('nan') + one = np.array([real_mp(1)]) + nans = np.array([nan]) + assert mp.abs(np.fmax(nans, one)[0] - real_mp(1)) == 0 + assert mp.abs(np.fmin(nans, one)[0] - real_mp(1)) == 0 + assert np.isnan(np.maximum(nans, one))[0] + assert np.isnan(np.minimum(nans, one))[0] + + def test_minimum_maximum_values(self): + a = np.array([real_mp(1), real_mp(5)]) + b = np.array([real_mp(3), real_mp(2)]) + assert [str(x) for x in np.minimum(a, b)] == ['1', '2'] + assert [str(x) for x in np.maximum(a, b)] == ['3', '5'] + + def test_predicates(self, dtype): + good = _sample(dtype) + assert not np.isnan(good).any() + assert np.isfinite(good).all() + assert not np.isinf(good).any() + assert np.isnan(good).dtype == np.dtype(bool) + if dtype is real_mp: + bad = np.array([real_mp('nan'), real_mp('inf'), real_mp(1)]) + else: + bad = np.array([complex_mp('nan', '0'), complex_mp('0', 'inf'), + complex_mp(1, 1)]) + assert list(np.isnan(bad)) == [True, False, False] + assert list(np.isinf(bad)) == [False, True, False] + assert list(np.isfinite(bad)) == [False, False, True] + + def test_signbit(self): + v = np.array([real_mp(-2), real_mp(0), real_mp(3)]) + assert list(np.signbit(v)) == [True, False, False] + + +class TestSortingAndArgExtrema: + """The dtype compare/argmax/argmin slots (real only -- complex is unordered).""" + + def test_sort_and_argsort(self): + v = np.array([real_mp(3), real_mp(1), real_mp(2)]) + assert [str(x) for x in np.sort(v)] == ['1', '2', '3'] + assert list(np.argsort(v)) == [1, 2, 0] + + def test_argmax_argmin_max_min(self): + v = np.array([real_mp(3), real_mp(1), real_mp(7), real_mp(2)]) + assert np.argmax(v) == 2 + assert np.argmin(v) == 1 + assert str(np.max(v)) == '7' + assert str(np.min(v)) == '1' + + def test_identityless_reduce_does_not_corrupt_input(self): + # numpy >= 2.5 initializes the accumulator of an identityless reduce + # (np.min/np.max) by BITWISE copy of element 0, so the accumulator and + # v[0] share one mpfr allocation. The loops must never free or write + # through an output slot's existing allocation (slot_write): before + # that rule, np.max silently rewrote v[0] and freed its storage, and + # the next reduce crashed the interpreter (use-after-free -> corrupted + # allocator). Values AND the input array must survive, repeatedly. + v = np.array([real_mp(3), real_mp(1), real_mp(7), real_mp(2)]) + for _ in range(3): + assert str(np.max(v)) == '7' + assert [str(x) for x in v] == ['3', '1', '7', '2'] + assert str(np.min(v)) == '1' + assert [str(x) for x in v] == ['3', '1', '7', '2'] + + def test_argmax_nan_wins(self): + # numpy float semantics: the first nan is the arg-extremum + v = np.array([real_mp(1), real_mp('nan'), real_mp(3)]) + assert np.argmax(v) == 1 + assert np.argmin(v) == 1 + + def test_searchsorted_and_median(self): + v = np.array([real_mp(1), real_mp(2), real_mp(4)]) + assert np.searchsorted(v, real_mp(3)) == 2 + assert str(np.median(v)) == '2' + + def test_complex_stays_unordered(self): + w = np.array([complex_mp(1, 2), complex_mp(0, 1)]) + with pytest.raises(TypeError): + np.sort(w) + + +class TestReductions: + """Identity-seeded reductions -- previously documented as crashing. + + They work on current numpy; these tests exist so any numpy/eigenpy + combination that breaks them again fails loudly here instead of in + user code. + """ + + def test_sum_prod_mean_real(self): + v = np.array([real_mp(1), real_mp(2), real_mp(3)]) + assert np.sum(v) == real_mp(6) + assert np.prod(v) == real_mp(6) + assert np.mean(v) == real_mp(2) + + def test_sum_prod_mean_complex(self): + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + assert np.sum(w) == complex_mp(4, 6) + assert np.prod(w) == complex_mp(-5, 10) + assert np.mean(w) == complex_mp(2, 3) + + def test_bare_reduce_and_cumsum(self, dtype): + v = np.array([dtype(1), dtype(2), dtype(3)]) + assert np.add.reduce(v) == dtype(6) + assert list(np.cumsum(v)) == [dtype(1), dtype(3), dtype(6)] + + def test_reduce_with_explicit_initial_still_works(self, dtype): + # the old explicit-initial workaround must keep working (portable to old numpy) + v = np.array([dtype(1), dtype(2)]) + assert np.add.reduce(v, initial=dtype(10)) == dtype(13) + + +class TestCloseness: + """The float64 boundary: tolerance ORDERINGS against a double are allowed + (a comparison yields a bool -- no float flows into a multiprecision value; + this mirrors the C++ solvers' double ToleranceT and the scalar + GreatLessVisitor). Everything that would let a float VALUE into + an mp computation stays closed: mixed equality, mixed arithmetic, + np.isclose's internal float tolerances.""" + + def test_tolerance_comparison_with_float(self, dtype): + v = _sample(dtype) + w = v + dtype('1e-20') + assert np.all(np.abs(v - w) <= 1e-10) + assert not np.all(np.abs(v - (w + dtype(1))) <= 1e-10) + # both operand orders, and float64 arrays as well as scalars + assert np.all(1e-10 >= np.abs(v - w)) + assert np.all(np.abs(v - w) < np.full(len(v), 1e-10)) + + def test_tolerance_comparison_is_exact_not_sloppy(self): + # the double is compared exactly (boost mixed compare), not by rounding + # the mp value down to double first + tiny = real_mp('1e-22') + assert np.all(np.array([tiny]) < 1e-10) + assert not np.any(np.array([tiny]) < 1e-30) + + def test_allclose_idiom_all_mp_still_works(self, dtype): + v = _sample(dtype) + w = v + dtype('1e-20') + assert np.all(np.abs(v - w) <= real_mp('1e-10')) + + def test_float_equality_stays_blocked(self, dtype): + # exact equality against a float literal is the 0.1-intent trap; it is + # not bound at the scalar level either + v = _sample(dtype) + with pytest.raises(TypeError): + v == 0.1 + + def test_float_arithmetic_stays_blocked(self, dtype): + v = _sample(dtype) + with pytest.raises(TypeError): + v + 0.1 + + def test_isclose_itself_still_raises(self, dtype): + # its internal float64 rtol/atol cannot promote; if this ever starts + # passing, numpy grew user-dtype promotion -- revisit the numpy docs page + v = _sample(dtype) + with pytest.raises(TypeError): + np.isclose(v, v) + + +class TestExplicitDownConversion: + """astype is the explicit, conscious truncation to double precision.""" + + def test_astype_float_and_complex(self): + v = np.array([real_mp('1.5'), real_mp(2)]) + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + assert list(v.astype(float)) == [1.5, 2.0] + assert list(v.astype(complex)) == [1.5 + 0j, 2.0 + 0j] + assert list(w.astype(complex)) == [1 + 2j, 3 + 4j] + + def test_astype_int_stays_forbidden(self, dtype): + v = np.array([dtype(1)]) + with pytest.raises(TypeError): + v.astype(np.int64) + + +class TestPrecisionPreservation: + """Outputs carry the operands' precision, not the ambient default -- + including through the mixed real/complex division inside sign and + reciprocal (the known boost precision-mis-tagging hazard).""" + + HIGH = 50 + + def _high_precision_sample(self, dtype): + mp.default_precision(self.HIGH) + v = (np.array([real_mp('1.5')]) if dtype is real_mp + else np.array([complex_mp('1.5', '2.5')])) + mp.default_precision(30) + assert v[0].precision == self.HIGH + return v + + @pytest.mark.parametrize("ufunc", [np.exp, np.sqrt, np.abs, np.sign, + np.reciprocal, np.conj, np.rint], + ids=lambda u: u.__name__) + def test_unary_output_precision(self, dtype, ufunc): + if dtype is complex_mp and ufunc is np.rint: + pytest.skip("rint is real-only") + v = self._high_precision_sample(dtype) + assert ufunc(v)[0].precision == self.HIGH + + def test_binary_output_precision(self, dtype): + v = self._high_precision_sample(dtype) + assert np.power(v, v)[0].precision == self.HIGH + + +class TestUnwrittenSlotSafety: + """Every new loop shape must survive never-written np.empty slots + (which hold the all-zero BMP sentinel) -- the ADR-0006 doctrine.""" + + def test_unary_loops_on_empty(self, dtype): + e = np.empty(3, dtype=dtype) + for ufunc in (np.exp, np.sin, np.conj, np.sign, np.abs, + np.isnan, np.isfinite): + ufunc(e) # must not crash + + def test_binary_loops_on_empty(self, dtype): + e = np.empty(3, dtype=dtype) + np.power(e, e) + if dtype is real_mp: + np.minimum(e, e) + np.arctan2(e, e) + np.mod(e, e) + + def test_sort_and_argmax_on_empty(self): + e = np.empty(4, dtype=real_mp) + np.sort(e) + np.argmax(e) + np.argmin(e) + + def test_reductions_on_zeros(self, dtype): + z = np.zeros(3, dtype=dtype) + assert np.sum(z) == dtype(0) + + +class TestScalarsAreOwnedCopies: + """Regression tests for the getitem aliasing fix. + + getitem used to return boost::ref into the numpy buffer (as stock eigenpy + does), so a scalar extracted from a temporary array -- most visibly the + result of np.sum/np.mean -- dangled once the array was freed: reading it + later gave zeros/garbage or SIGABRT inside mpfr (the ADR-0031 / #259 + hazard class). getitem now returns an owned copy. + """ + + def test_reduce_scalar_survives_its_array(self, dtype): + s = np.sum(np.array([dtype(1), dtype(2), dtype(3)])) + # the source (temporary) array is gone; s must still be intact + assert s == dtype(6) + assert str(s) is not None # printing used to MPFR-assert on the corpse + assert s / dtype(3) == dtype(2) + + def test_mean_is_correct_not_silently_zero(self, dtype): + # np.mean's internal divide ran on a dangling extraction and returned 0 + v = np.array([dtype(1), dtype(2), dtype(3)]) + assert np.mean(v) == dtype(2) + + def test_stored_indexed_elements_stay_distinct(self): + # the ADR-0031 shape, now safe at the binding level (still copy in + # Python code by convention) + pts = np.array([complex_mp(1, 1), complex_mp(2, 2), complex_mp(3, 3)]) + kept = [pts[i] for i in range(3)] + del pts + assert [str(k) for k in kept] == ['(1,1)', '(2,2)', '(3,3)'] + + def test_mutating_an_extracted_scalar_leaves_the_array_alone(self): + v = np.array([real_mp(1), real_mp(2)]) + x = v[0] + x += real_mp(10) + assert v[0] == real_mp(1) + + +class TestComponentAccessors: + """The multiprec real/imag/arg array overloads, and the scalar imag + regression.""" + + def test_scalar_imag_returns_imaginary_part(self): + # regression: mp.imag was bound to boost::multiprecision::real by a + # copy-paste error, so it returned the REAL part + z = complex_mp(1, 2) + assert mp.real(z) == real_mp(1) + assert mp.imag(z) == real_mp(2) + + def test_array_real_imag(self): + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + r, i = mp.real(w), mp.imag(w) + assert r.dtype == np.dtype(real_mp) and i.dtype == np.dtype(real_mp) + assert [str(x) for x in r] == ['1', '3'] + assert [str(x) for x in i] == ['2', '4'] + + def test_array_arg_replaces_np_angle(self): + w = np.array([complex_mp(1, 1), complex_mp(-1, 0)]) + a = mp.arg(w) + assert a.dtype == np.dtype(real_mp) + assert a[0] == mp.arg(w[0]) + assert a[1] == mp.arg(w[1]) + + def test_ndarray_real_imag_attributes_are_untrustworthy(self): + # documenting-by-test: numpy cannot know a legacy user dtype is + # complex-like, so ndarray .real returns the complex values themselves + # and .imag returns zeros. If numpy ever fixes this, the accessors + # above become optional and the numpy docs page should be updated. + w = np.array([complex_mp(1, 2)]) + assert w.real.dtype == np.dtype(complex_mp) # not real_mp! + assert w.imag[0] == complex_mp(0) # wrong value, by numpy + + +class TestGuardedNumpyComponentFunctions: + """np.real/np.imag/np.angle raise on plain mp-complex arrays instead of + silently returning wrong values (bertini._numpy_guard) -- a crash is better + than incorrect values. Solution points override .real/.imag at the subclass + level and pass through correct.""" + + def test_np_real_imag_raise_on_plain_complex_mp_array(self): + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + with pytest.raises(TypeError, match="bertini.real"): + np.real(w) + with pytest.raises(TypeError, match="bertini.real"): + np.imag(w) + with pytest.raises(TypeError, match="bertini.real"): + np.angle(w) + + def test_np_real_imag_raise_on_lists_of_complex_mp(self): + # a list converts to a plain mp array inside numpy, same wrong path + with pytest.raises(TypeError): + np.imag([complex_mp(1, 2)]) + + def test_guard_passes_everything_else_through(self): + # ordinary numpy is untouched + z = np.array([1 + 2j, 3 + 4j]) + assert list(np.real(z)) == [1.0, 3.0] + assert list(np.imag(z)) == [2.0, 4.0] + assert np.angle(np.array([1j]))[0] == pytest.approx(np.pi / 2) + # real_mp arrays are not complex: base semantics are already correct + v = np.array([real_mp(1), real_mp(2)]) + assert list(np.real(v)) == [real_mp(1), real_mp(2)] + assert list(np.imag(v)) == [real_mp(0), real_mp(0)] + # mp-complex SCALARS go through the (correct) scalar properties + assert np.real(complex_mp(1, 2)) == real_mp(1) + assert np.imag(complex_mp(1, 2)) == real_mp(2) + + def test_guard_is_idempotent(self): + import bertini._numpy_guard as guard + before = np.real + guard.install() + assert np.real is before + + def test_solution_real_imag_are_correct(self): + from bertini.records import Solution + s = Solution(np.array([complex_mp(1, 2), complex_mp(3, 4)])) + assert [str(x) for x in s.real] == ['1', '3'] + assert [str(x) for x in s.imag] == ['2', '4'] + assert s.real.dtype == np.dtype(real_mp) + # np.real/np.imag on a Solution route through the subclass property + assert [str(x) for x in np.real(s)] == ['1', '3'] + assert [str(x) for x in np.imag(s)] == ['2', '4'] + + def test_solution_real_imag_correct_for_double_solves_too(self): + from bertini.records import Solution + s = Solution(np.array([1 + 2j, 3 + 4j])) + assert list(s.real) == [1.0, 3.0] + assert list(s.imag) == [2.0, 4.0] + + def test_np_angle_raises_helpfully_for_all_mp_complex(self): + # np.angle branches on the DTYPE (never the .real/.imag attributes), so + # not even the Solution subclass can make it work -- the guard turns the + # cryptic arctan2 failure into a pointer at mp.arg, for every spelling + from bertini.records import Solution + for val in (np.array([complex_mp(1, 1)]), + Solution(np.array([complex_mp(1, 1)])), + complex_mp(1, 1)): + with pytest.raises(TypeError, match="arg"): + np.angle(val) diff --git a/python/test/classes/operators_test.py b/python/test/classes/operators_test.py new file mode 100644 index 000000000..5c54eff96 --- /dev/null +++ b/python/test/classes/operators_test.py @@ -0,0 +1,153 @@ +# This file is part of Bertini 2. +# +# python/test/classes/operators_test.py is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# python/test/classes/operators_test.py is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this file. If not, see . +# +# Copyright(C) Bertini2 Development Team + +"""bertini.operators: one star-import, functions that work on symbols AND numbers AND +numpy containers alike, dispatched per argument.""" + +import builtins + +import numpy as np +import pytest + +import bertini as pb +import bertini.multiprec as mp +import bertini.operators as ops +from bertini.multiprec import complex_mp, real_mp +from bertini._pybertini.function_tree import AbstractNode + + +class TestPolymorphicDispatch: + """sin & friends: symbolic on expressions, numeric on everything else.""" + + def test_symbolic_on_variables(self): + x = pb.Variable('x') + f = ops.sin(x) + ops.Pi * x - ops.E + assert isinstance(ops.sin(x), AbstractNode) + assert isinstance(f, AbstractNode) + + def test_numeric_on_mp_scalars(self): + v = real_mp('0.5') + assert ops.sin(v) == mp.sin(v) + assert ops.exp(v) == mp.exp(v) + z = complex_mp('0.5', '0.25') + assert ops.sqrt(z) == mp.sqrt(z) + + def test_numeric_on_mp_arrays(self): + v = np.array([real_mp('0.25'), real_mp('0.5')]) + out = ops.cos(v) + assert out.dtype == np.dtype(real_mp) + assert out[0] == mp.cos(v[0]) + + def test_numeric_on_lists_and_python_numbers(self): + assert ops.sin(0.0) == 0.0 + out = ops.tan([real_mp('0.25'), real_mp('0.5')]) + assert out[1] == mp.tan(real_mp('0.5')) + + def test_asin_maps_to_arcsin(self): + v = real_mp('0.5') + assert ops.asin(v) == mp.asin(v) + assert ops.acos(v) == mp.acos(v) + assert ops.atan(v) == mp.atan(v) + + def test_hyperbolics_numeric(self): + v = real_mp('0.5') + assert ops.sinh(v) == mp.sinh(v) + assert ops.atanh(v) == mp.atanh(v) + + def test_hyperbolics_reject_symbols(self): + x = pb.Variable('x') + with pytest.raises(TypeError, match="symbolic"): + ops.sinh(x) + + +class TestComponentsAndFriends: + """abs/arg/real/imag/conj/round/sum/norm/is_real, all in the same namespace.""" + + def test_components_on_arrays(self): + w = np.array([complex_mp(1, 2), complex_mp(3, 4)]) + assert [str(t) for t in ops.real(w)] == ['1', '3'] + assert [str(t) for t in ops.imag(w)] == ['2', '4'] + assert str(ops.abs(np.array([complex_mp(3, 4)]))[0]) == '5' + assert ops.arg(w)[0] == mp.arg(w[0]) + assert complex(ops.conj(w)[1]) == complex(3, -4) + + def test_arg_on_scalars_and_reals(self): + assert ops.arg(complex_mp(0, 1)) == mp.arg(complex_mp(0, 1)) + # arg of a negative real is pi + assert mp.abs(ops.arg(np.array([real_mp(-2)]))[0] - mp.arg(complex_mp(-2))) == 0 + # plain python numbers give floats + assert ops.arg(1j) == pytest.approx(np.pi / 2) + + def test_sum_norm_is_real(self): + v = np.array([real_mp(3), real_mp(4)]) + assert ops.sum(v) == real_mp(7) + assert ops.norm(v) == real_mp(5) + assert ops.is_real(np.array([complex_mp(1)])) is True + + def test_round_stays_decimal(self): + assert str(ops.round(real_mp('2.34567'), 2)) == '2.35' + + def test_numeric_only_reject_symbols(self): + x = pb.Variable('x') + for fn in (ops.abs, ops.arg, ops.real, ops.imag, ops.conj, + ops.round, ops.sum, ops.norm, ops.is_real): + with pytest.raises(TypeError, match="symbolic"): + fn(x) + + def test_builtin_fallback_on_plain_python(self): + # the shadowing names still behave sanely on plain python input + assert ops.abs(-3) == 3 + assert ops.sum([1, 2, 3]) == 6 + assert ops.round(2.345, 1) == builtins.round(2.345, 1) + + +class TestStarImportSurface: + def test_star_import_gives_the_whole_vocabulary(self): + ns = {} + exec("from bertini.operators import *", ns) + for name in ('sin', 'cos', 'tan', 'asin', 'acos', 'atan', + 'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh', + 'exp', 'log', 'sqrt', + 'abs', 'arg', 'real', 'imag', 'conj', + 'round', 'sum', 'norm', 'is_real', + 'E', 'Pi', 'I'): + assert name in ns, name + + def test_one_import_covers_symbols_and_numbers(self): + # the point of the module, as a single flow + ns = {} + exec("from bertini.operators import *", ns) + x = pb.Variable('x') + f = ns['sin'](x) # symbolic + assert isinstance(f, AbstractNode) + val = ns['sin'](real_mp('0.5')) # mp scalar + assert val == mp.sin(real_mp('0.5')) + w = np.array([complex_mp(1, 2)]) # numpy container + assert ns['imag'](w)[0] == real_mp(2) + + def test_top_level_functions_are_polymorphic_too(self): + x = pb.Variable('x') + assert isinstance(pb.sin(x), AbstractNode) + assert pb.sin(real_mp('0.5')) == mp.sin(real_mp('0.5')) + assert pb.arg(complex_mp(1, 1)) == mp.arg(complex_mp(1, 1)) + + def test_from_bertini_star_still_never_shadows_builtins(self): + ns = {} + exec("from bertini import *", ns) + assert 'abs' not in ns + assert 'round' not in ns + assert 'sum' not in ns diff --git a/python/test/tracking/endgame_test.py b/python/test/tracking/endgame_test.py index 2b2d9c917..8b35215aa 100644 --- a/python/test/tracking/endgame_test.py +++ b/python/test/tracking/endgame_test.py @@ -137,9 +137,8 @@ def test_using_total_degree_ss(): for soln in dehomogenized_solns: diff = exact_soln - soln - # NB: np.sum / np.prod / np.mean over multiprecision (mpfr/mpc) arrays can raise - # SystemError on some numpy + eigenpy builds -- numpy cannot construct the reduction - # identity element for these custom dtypes. See the "Known gotchas" page in the docs. - # np.dot(diff, diff) == sum(diff_i**2) (numpy's dot does not conjugate) and goes through - # the dtype's dot slot, which works everywhere; it preserves this assertion exactly. + # np.dot(diff, diff) == sum(diff_i**2) (numpy's dot does not conjugate); it goes + # through the dtype's dot slot. np.sum would work too on current numpy, but dot + # stays portable to older numpy builds where the identity-seeded reduce raised + # SystemError for user dtypes. assert mp.abs(np.sqrt(np.dot(diff, diff))) < 1e-10 diff --git a/python_bindings/CMakeLists.txt b/python_bindings/CMakeLists.txt index 8252361e5..a20f30372 100644 --- a/python_bindings/CMakeLists.txt +++ b/python_bindings/CMakeLists.txt @@ -293,6 +293,13 @@ endif() # from scratch, causing peak-RAM spikes during parallel builds and inflating log output. target_precompile_headers(_pybertini PRIVATE "include/python_common.hpp") +# clang rejects a PCH whose input file's mtime changed, even with identical content -- +# which happens when the wheel build re-configures into a shared bld/. See the note +# beside bertini2_exe's PCH in core/CMakeLists.txt. +if(CMAKE_CXX_COMPILER_ID MATCHES "Clang") + target_compile_options(_pybertini PRIVATE "$<$:SHELL:-Xclang -fno-pch-timestamp>") +endif() + cmake_print_variables(SKBUILD) if(${SKBUILD}) # see https://stackoverflow.com/questions/1242904/finding-python-site-packages-directory-with-cmake diff --git a/python_bindings/include/eigenpy_interaction.hpp b/python_bindings/include/eigenpy_interaction.hpp index bcd0455b1..00d8219fc 100644 --- a/python_bindings/include/eigenpy_interaction.hpp +++ b/python_bindings/include/eigenpy_interaction.hpp @@ -7,7 +7,9 @@ #include "python_common.hpp" +#include #include +#include #include #include @@ -76,7 +78,15 @@ namespace eigenpy } - // template specialization for real numbers + // template specialization for real numbers. + // + // NB: eigenpy's stock getitem (and an earlier version of this one) returns + // boost::ref(slot) — a Python object ALIASING the numpy buffer. That is + // the root of the ADR-0031 hazard family (#259: stored elements silently + // collapse or SIGABRT once the buffer is reused/freed), and it makes + // scalars extracted from temporary arrays — np.sum/np.mean results, most + // visibly — dangle outright. Returning an owned COPY kills the whole + // class: an indexed element is a durable value, as numpy users expect. template <> struct getitem { @@ -90,13 +100,14 @@ namespace eigenpy { mpfr_scalar = NumT(0); } - boost::python::object m(boost::ref(mpfr_scalar)); + boost::python::object m(mpfr_scalar); // owned copy — never boost::ref (see above) Py_INCREF(m.ptr()); return m.ptr(); } }; - // a template specialization for complex numbers + // a template specialization for complex numbers; see the real one for the + // copy-not-ref rationale. template <> struct getitem { @@ -110,7 +121,7 @@ namespace eigenpy { mpfr_scalar = NumT(0); } - boost::python::object m(boost::ref(mpfr_scalar)); + boost::python::object m(mpfr_scalar); // owned copy — never boost::ref (see above) Py_INCREF(m.ptr()); return m.ptr(); } @@ -156,19 +167,45 @@ namespace eigenpy // These replace eigenpy's EIGENPY_REGISTER_{BINARY,UNARY}_UFUNC loop bodies // (and its gufunc_matrix_multiply), which read input slots unguarded and // segfault inside libmpfr/libmpc on never-written np.zeros/np.empty slots. - // Writes into the output slot go through BMP operator=, which initializes - // a zeroed destination itself. + // + // Writing an output slot: NEVER through BMP operator= on the slot's + // existing value. numpy may hand a loop an output slot that bitwise- + // ALIASES another slot's mpfr allocation -- numpy 2.5 initializes the + // accumulator of an identityless reduce (np.min/np.max) by memcpy of + // element 0, so the accumulator and v[0] share one set of limbs. A + // plain assignment move-frees the slot's old limbs (freeing v[0]'s + // storage out from under it: use-after-free, double-free, corrupted + // allocator, SIGSEGV a few calls later) and writes through shared + // storage (silently mutating v[0]'s value). slot_write below is the + // only sanctioned store: compute the value FIRST (the slot may also be + // an input), memset the slot to BMP's uninitialized sentinel, then + // move the fresh value in -- no existing allocation is ever freed or + // written through. If the slot held a uniquely-owned value, its + // allocation leaks (numpy never destructs user-dtype elements anyway); + // same crash-into-bounded-leak trade as HardenSetitem / ADR-0006. + + // store `val` into a numpy-managed mp slot without freeing or writing + // through the slot's existing (possibly aliased) allocation. + template + inline void slot_write(T& slot, T val) + { + std::memset(static_cast(&slot), 0, sizeof(T)); + slot = std::move(val); // move into the sentinel: steals val's limbs, frees nothing + } struct op_add { template static T apply(T const& x, T const& y) { return T(x + y); } }; struct op_subtract { template static T apply(T const& x, T const& y) { return T(x - y); } }; struct op_multiply { template static T apply(T const& x, T const& y) { return T(x * y); } }; struct op_divide { template static T apply(T const& x, T const& y) { return T(x / y); } }; - struct op_equal { template static bool apply(T const& x, T const& y) { return x == y; } }; - struct op_not_equal { template static bool apply(T const& x, T const& y) { return x != y; } }; - struct op_greater { template static bool apply(T const& x, T const& y) { return x > y; } }; - struct op_less { template static bool apply(T const& x, T const& y) { return x < y; } }; - struct op_greater_equal { template static bool apply(T const& x, T const& y) { return x >= y; } }; - struct op_less_equal { template static bool apply(T const& x, T const& y) { return x <= y; } }; + // comparison functors are heterogeneous (two type parameters) so the same + // functor serves the mp-vs-mp loops AND the mp-vs-double tolerance loops + // (boost::multiprecision compares a number against a double exactly). + struct op_equal { template static bool apply(T const& x, U const& y) { return x == y; } }; + struct op_not_equal { template static bool apply(T const& x, U const& y) { return x != y; } }; + struct op_greater { template static bool apply(T const& x, U const& y) { return x > y; } }; + struct op_less { template static bool apply(T const& x, U const& y) { return x < y; } }; + struct op_greater_equal { template static bool apply(T const& x, U const& y) { return x >= y; } }; + struct op_less_equal { template static bool apply(T const& x, U const& y) { return x <= y; } }; struct op_negative { template static T apply(T const& x) { return T(-x); } }; struct op_square { template static T apply(T const& x) { return T(x * x); } }; @@ -182,6 +219,232 @@ namespace eigenpy } }; + // trait: is this scalar the complex mp type? several ops (conjugate, sign, + // the isnan/isinf/isfinite predicates) need a different body for complex. + template struct is_complex_mp : std::false_type {}; + template <> struct is_complex_mp : std::true_type {}; + + // re-tag `val` to carry `ref`'s precision. guards against boost's mixed + // real/complex arithmetic occasionally mis-tagging the result's precision + // (division is the known offender); .precision(n) preserves the value. + template + inline T at_precision_of(T val, T const& ref) + { + if (val.precision() != ref.precision()) + val.precision(ref.precision()); + return val; + } + + // ----- unary ops, same-type output ------------------------------------ + // bodies call the same boost::multiprecision free functions the multiprec + // module binds as scalar functions, so np.exp(arr)[i] == mp.exp(arr[i]). + + struct op_positive { template static T apply(T const& x) { return x; } }; + struct op_reciprocal + { + template static T apply(T const& x) + { + T one(1); + one.precision(x.precision()); + return at_precision_of(T(one / x), x); + } + }; + struct op_conjugate + { + template static T apply(T const& x) + { + if constexpr (is_complex_mp::value) + return T(conj(x)); + else + return x; + } + }; + // numpy-2 sign semantics: real -> -1/0/+1; complex -> z/|z| (0 at 0). + struct op_sign + { + template static T apply(T const& x) + { + if constexpr (is_complex_mp::value) + { + if (x == 0) + return at_precision_of(T(0), x); + bertini::real_mp const mag(abs(x)); + return at_precision_of(T(bertini::real_mp(x.real() / mag), + bertini::real_mp(x.imag() / mag)), x); + } + else + { + T res(x > 0 ? 1 : (x < 0 ? -1 : 0)); + return at_precision_of(std::move(res), x); + } + } + }; + + struct op_exp { template static T apply(T const& x) { return T(exp(x)); } }; + struct op_log { template static T apply(T const& x) { return T(log(x)); } }; + struct op_log10 { template static T apply(T const& x) { return T(log10(x)); } }; + struct op_exp2 { template static T apply(T const& x) { return T(exp2(x)); } }; + struct op_log2 { template static T apply(T const& x) { return T(log2(x)); } }; + struct op_expm1 { template static T apply(T const& x) { return T(expm1(x)); } }; + struct op_log1p { template static T apply(T const& x) { return T(log1p(x)); } }; + struct op_cbrt { template static T apply(T const& x) { return T(cbrt(x)); } }; + + struct op_sin { template static T apply(T const& x) { return T(sin(x)); } }; + struct op_cos { template static T apply(T const& x) { return T(cos(x)); } }; + struct op_tan { template static T apply(T const& x) { return T(tan(x)); } }; + struct op_arcsin { template static T apply(T const& x) { return T(asin(x)); } }; + struct op_arccos { template static T apply(T const& x) { return T(acos(x)); } }; + struct op_arctan { template static T apply(T const& x) { return T(atan(x)); } }; + struct op_sinh { template static T apply(T const& x) { return T(sinh(x)); } }; + struct op_cosh { template static T apply(T const& x) { return T(cosh(x)); } }; + struct op_tanh { template static T apply(T const& x) { return T(tanh(x)); } }; + struct op_arcsinh { template static T apply(T const& x) { return T(asinh(x)); } }; + struct op_arccosh { template static T apply(T const& x) { return T(acosh(x)); } }; + struct op_arctanh { template static T apply(T const& x) { return T(atanh(x)); } }; + + // real-only rounding family + struct op_floor { template static T apply(T const& x) { return T(floor(x)); } }; + struct op_ceil { template static T apply(T const& x) { return T(ceil(x)); } }; + struct op_trunc { template static T apply(T const& x) { return T(trunc(x)); } }; + // numpy's rint is round-half-to-EVEN; boost's rint rounds half away from + // zero, so call mpfr directly in MPFR_RNDN (nearest, ties to even). + // Unlike the rest of the rounding family, numpy defines rint for complex + // (component-wise) — np.round on a complex array goes through it. + struct op_rint + { + static bertini::real_mp rint_one(bertini::real_mp const& x) + { + bertini::real_mp out(0); + out.precision(x.precision()); + mpfr_rint(out.backend().data(), x.backend().data(), MPFR_RNDN); + return out; + } + + template static T apply(T const& x) + { + if constexpr (is_complex_mp::value) + return at_precision_of(T(rint_one(x.real()), rint_one(x.imag())), x); + else + return rint_one(x); + } + }; + + // ----- unary ops, cross-type output ----------------------------------- + + // absolute: real -> real, complex -> real (the magnitude). + struct op_absolute + { + template static bertini::real_mp apply(T const& x) + { + return bertini::real_mp(abs(x)); + } + }; + struct op_fabs { template static T apply(T const& x) { return T(fabs(x)); } }; + + // predicates -> bool. the isnan/isinf/isfinite family are function-like + // macros in C , so call the boost versions qualified. + struct op_isnan + { + template static bool apply(T const& x) + { + if constexpr (is_complex_mp::value) + return boost::multiprecision::isnan(x.real()) || boost::multiprecision::isnan(x.imag()); + else + return boost::multiprecision::isnan(x); + } + }; + struct op_isinf + { + template static bool apply(T const& x) + { + if constexpr (is_complex_mp::value) + return boost::multiprecision::isinf(x.real()) || boost::multiprecision::isinf(x.imag()); + else + return boost::multiprecision::isinf(x); + } + }; + struct op_isfinite + { + template static bool apply(T const& x) + { + if constexpr (is_complex_mp::value) + return boost::multiprecision::isfinite(x.real()) && boost::multiprecision::isfinite(x.imag()); + else + return boost::multiprecision::isfinite(x); + } + }; + struct op_signbit + { + template static bool apply(T const& x) + { + return boost::multiprecision::signbit(x); + } + }; + + // ----- binary ops ------------------------------------------------------ + + struct op_power { template static T apply(T const& x, T const& y) { return T(pow(x, y)); } }; + struct op_arctan2 { template static T apply(T const& x, T const& y) { return T(atan2(x, y)); } }; + struct op_hypot { template static T apply(T const& x, T const& y) { return T(hypot(x, y)); } }; + struct op_copysign { template static T apply(T const& x, T const& y) { return T(copysign(x, y)); } }; + // numpy fmod keeps C fmod's sign-of-dividend semantics + struct op_fmod { template static T apply(T const& x, T const& y) { return T(fmod(x, y)); } }; + // numpy remainder/mod takes the sign of the DIVISOR (python % semantics) + struct op_remainder + { + template static T apply(T const& x, T const& y) + { + T r(fmod(x, y)); + if (r != 0 && ((r < 0) != (y < 0))) + r += y; + return r; + } + }; + struct op_floor_divide + { + template static T apply(T const& x, T const& y) + { + return T(floor(x / y)); + } + }; + // minimum/maximum propagate nan (numpy semantics); fmin/fmax ignore it + struct op_minimum + { + template static T apply(T const& x, T const& y) + { + if (boost::multiprecision::isnan(x)) return x; + if (boost::multiprecision::isnan(y)) return y; + return y < x ? y : x; + } + }; + struct op_maximum + { + template static T apply(T const& x, T const& y) + { + if (boost::multiprecision::isnan(x)) return x; + if (boost::multiprecision::isnan(y)) return y; + return x < y ? y : x; + } + }; + struct op_fmin + { + template static T apply(T const& x, T const& y) + { + if (boost::multiprecision::isnan(x)) return y; + if (boost::multiprecision::isnan(y)) return x; + return y < x ? y : x; + } + }; + struct op_fmax + { + template static T apply(T const& x, T const& y) + { + if (boost::multiprecision::isnan(x)) return y; + if (boost::multiprecision::isnan(y)) return x; + return x < y ? y : x; + } + }; + template void guarded_binary_op( char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, @@ -195,7 +458,7 @@ namespace eigenpy T const& x = value_or_zero(*reinterpret_cast(i0), zero); T const& y = value_or_zero(*reinterpret_cast(i1), zero); T& res = *reinterpret_cast(o); - res = Op::apply(x, y); + slot_write(res, Op::apply(x, y)); // never plain operator= -- see slot_write i0 += is0; i1 += is1; o += os; @@ -222,6 +485,44 @@ namespace eigenpy } } + // mixed mp-vs-float64 ORDERING comparison (Reversed swaps operand order: + // false = (mp, double), true = (double, mp)). Comparisons against a + // double tolerance -- np.abs(a - b) < 1e-10 -- are safe: boost compares a + // number against a double exactly, and the result is a bool, so no float + // ever flows INTO a multiprecision value. This mirrors the scalar + // bindings (GreatLessVisitor) and the C++ solvers' double + // ToleranceT. Deliberately orderings-only: mixed EQUALITY with a float + // literal is the 0.1-intent trap the unsafe double->mp cast exists to + // block, and it is not bound at the scalar level either. + template + void guarded_mixed_compare_op( + char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, + EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *steps, void * /*data*/) + { + npy_intp is0 = steps[0], is1 = steps[1], os = steps[2], n = *dimensions; + char *i0 = args[0], *i1 = args[1], *o = args[2]; + const T zero(0); + for (npy_intp k = 0; k < n; ++k) + { + bool& res = *reinterpret_cast(o); + if constexpr (Reversed) + { + double const& x = *reinterpret_cast(i0); + T const& y = value_or_zero(*reinterpret_cast(i1), zero); + res = Op::apply(x, y); + } + else + { + T const& x = value_or_zero(*reinterpret_cast(i0), zero); + double const& y = *reinterpret_cast(i1); + res = Op::apply(x, y); + } + i0 += is0; + i1 += is1; + o += os; + } + } + template void guarded_unary_op( char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, @@ -234,7 +535,32 @@ namespace eigenpy { T const& x = value_or_zero(*reinterpret_cast(i), zero); T& res = *reinterpret_cast(o); - res = Op::apply(x); + slot_write(res, Op::apply(x)); // never plain operator= -- see slot_write + i += is; + o += os; + } + } + + // unary loop with an output type different from the input type + // (absolute: complex -> real; the isnan family: T -> bool). Writes into + // mp-typed output slots go through BMP operator=, which initializes a + // zeroed destination itself; bool slots are plain bytes. + template + void guarded_unary_op_out( + char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, + EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *steps, void * /*data*/) + { + npy_intp is = steps[0], os = steps[1], n = *dimensions; + char *i = args[0], *o = args[1]; + const T zero(0); + for (npy_intp k = 0; k < n; ++k) + { + T const& x = value_or_zero(*reinterpret_cast(i), zero); + OutT& res = *reinterpret_cast(o); + if constexpr (std::is_trivially_copyable_v) + res = Op::apply(x); + else + slot_write(res, Op::apply(x)); // never plain operator= -- see slot_write i += is; o += os; } @@ -267,7 +593,7 @@ namespace eigenpy b += is2_n; } T& res = *reinterpret_cast(op); - res = sum; + slot_write(res, std::move(sum)); // never plain operator= -- see slot_write ip2 += is2_p; op += os_p; } @@ -323,7 +649,54 @@ namespace eigenpy p0 += is0; p1 += is1; } - *reinterpret_cast(op) = acc; + slot_write(*reinterpret_cast(op), std::move(acc)); // never plain operator= -- see slot_write + } + + // guarded element comparison for the PyArray_ArrFuncs `compare` slot + // (np.sort / argsort / searchsorted / unique). eigenpy leaves this slot + // empty for user dtypes ("type does not have compare function"). Only + // installed for the real type — complex has no ordering. nan compares + // false both ways (weak-ordering violation, same as C doubles): sorting + // arrays containing nan gives an unspecified nan position, not a crash. + template + int guarded_compare(const void *a, const void *b, void * /*arr*/) + { + const T zero(0); + T const& x = value_or_zero(*static_cast(a), zero); + T const& y = value_or_zero(*static_cast(b), zero); + if (x < y) return -1; + if (y < x) return 1; + return 0; + } + + // guarded argmax/argmin for the PyArray_ArrFuncs slots (np.argmax / + // np.argmin / np.max / np.min dispatch through these for user dtypes on + // some numpy paths). Mirrors numpy's float semantics: a nan wins + // immediately (first nan is the arg-extremum). numpy hands these a + // contiguous buffer. + template + int guarded_argminmax(void *data, npy_intp n, npy_intp *extremum_ind, void * /*arr*/) + { + const T zero(0); + T const* p = static_cast(data); + *extremum_ind = 0; + if (n == 0) + return 0; + T best = value_or_zero(p[0], zero); + if (boost::multiprecision::isnan(best)) + return 0; + for (npy_intp k = 1; k < n; ++k) + { + T const& v = value_or_zero(p[k], zero); + if (boost::multiprecision::isnan(v) || (Max ? best < v : v < best)) + { + *extremum_ind = k; + if (boost::multiprecision::isnan(v)) + return 0; + best = v; + } + } + return 0; } } // namespace internal @@ -354,6 +727,32 @@ namespace eigenpy } }; + // mp -> complex128, for the explicit down-conversion arr.astype(complex) + // (registered unsafe, like mp -> double: you consciously truncate). + // boost mp numbers have no conversion operator to std::complex, so go + // through the components. + template <> + struct cast> + { + static std::complex run(bertini::real_mp const& from) + { + if (internal::mpfr_slot::uninitialized(from)) + return {0.0, 0.0}; + return {from.convert_to(), 0.0}; + } + }; + + template <> + struct cast> + { + static std::complex run(bertini::complex_mp const& from) + { + if (internal::mpfr_slot::uninitialized(from)) + return {0.0, 0.0}; + return {from.real().convert_to(), from.imag().convert_to()}; + } + }; + // Install the zero-initialization setitem guard for an MPFR-backed dtype. // Call immediately after eigenpy::registerNewType(), before any arrays @@ -378,6 +777,29 @@ namespace eigenpy funcs->dotfunc = reinterpret_cast(&internal::guarded_dotfunc); } + // Fill the element-comparison slot (empty in eigenpy's registration), enabling + // np.sort / np.argsort / np.searchsorted / np.unique. Real type only — + // complex has no ordering. Call immediately after eigenpy::registerNewType. + template + void HardenCompare() + { + PyArray_Descr *descr = Register::getPyArrayDescr(); + PyArray_ArrFuncs *funcs = PyDataType_GetArrFuncs(descr); + funcs->compare = reinterpret_cast(&internal::guarded_compare); + } + + // Fill the argmax/argmin slots (empty in eigenpy's registration), enabling + // np.argmax / np.argmin ("data type not ordered" otherwise). Real type only. + // Call immediately after eigenpy::registerNewType. + template + void HardenArgMinMax() + { + PyArray_Descr *descr = Register::getPyArrayDescr(); + PyArray_ArrFuncs *funcs = PyDataType_GetArrFuncs(descr); + funcs->argmax = reinterpret_cast(&internal::guarded_argminmax); + funcs->argmin = reinterpret_cast(&internal::guarded_argminmax); + } + // register a single guarded loop on the named numpy ufunc, mirroring the // error handling of eigenpy's EIGENPY_REGISTER_*_UFUNC macros. inline void registerGuardedLoop(PyObject *numpy, char const *ufunc_name, @@ -414,13 +836,20 @@ namespace eigenpy // i lifted this from EigenPy and adapted it: all loops are the guarded // versions from internal:: above (eigenpy's read input slots unguarded — - // see the header comment), and the ordering comparitors are a compile-time - // option because they are NOT defined for complex types (instantiating - // them for complex_mp would be a hard error). + // see the header comment), and the ordering-dependent set is a compile-time + // option because ordering is NOT defined for complex types (instantiating + // those functors for complex_mp would be a hard error). Coverage beyond + // eigenpy's arithmetic core (absolute, conjugate, the transcendental family, + // rounding, min/max, the isnan predicates) closes the documented + // "ufunc not supported" gotchas — every loop body calls the same + // boost::multiprecision free function the multiprec module binds as the + // scalar function of the same name. template void registerGuardedUfunct() { const int type_code = Register::getTypeCode(); + const int bool_code = Register::getTypeCode(); + const int real_code = Register::getTypeCode(); PyObject *numpy_str; #if PY_MAJOR_VERSION >= 3 @@ -434,6 +863,19 @@ namespace eigenpy import_ufunc(); + // registration helpers: (in...) -> out signatures. the types array is + // copied by PyUFunc_RegisterLoopForType, so stack storage is fine. + auto unary = [&](char const* name, PyUFuncGenericFunction loop, int out_code) + { + int types[2] = {type_code, out_code}; + registerGuardedLoop(numpy, name, type_code, loop, types, 2); + }; + auto binary = [&](char const* name, PyUFuncGenericFunction loop, int out_code) + { + int types[3] = {type_code, type_code, out_code}; + registerGuardedLoop(numpy, name, type_code, loop, types, 3); + }; + // Matrix multiply { int types[3] = {type_code, type_code, type_code}; @@ -442,49 +884,109 @@ namespace eigenpy types, 3); } - // Binary operators - { - int types[3] = {type_code, type_code, type_code}; - registerGuardedLoop(numpy, "add", type_code, - &internal::guarded_binary_op, types, 3); - registerGuardedLoop(numpy, "subtract", type_code, - &internal::guarded_binary_op, types, 3); - registerGuardedLoop(numpy, "multiply", type_code, - &internal::guarded_binary_op, types, 3); - registerGuardedLoop(numpy, "divide", type_code, - &internal::guarded_binary_op, types, 3); - } - - // Comparison operators + // Binary arithmetic + binary("add", &internal::guarded_binary_op, type_code); + binary("subtract", &internal::guarded_binary_op, type_code); + binary("multiply", &internal::guarded_binary_op, type_code); + binary("divide", &internal::guarded_binary_op, type_code); + binary("power", &internal::guarded_binary_op, type_code); + + // Equality comparisons (defined for real and complex alike) + binary("equal", &internal::guarded_compare_op, bool_code); + binary("not_equal", &internal::guarded_compare_op, bool_code); + + // Unary, same-type output + unary("negative", &internal::guarded_unary_op, type_code); + unary("positive", &internal::guarded_unary_op, type_code); + unary("square", &internal::guarded_unary_op, type_code); + unary("sqrt", &internal::guarded_unary_op, type_code); + unary("reciprocal", &internal::guarded_unary_op, type_code); + unary("conjugate", &internal::guarded_unary_op, type_code); + unary("sign", &internal::guarded_unary_op, type_code); + unary("exp", &internal::guarded_unary_op, type_code); + unary("log", &internal::guarded_unary_op, type_code); + unary("log10", &internal::guarded_unary_op, type_code); + unary("sin", &internal::guarded_unary_op, type_code); + unary("cos", &internal::guarded_unary_op, type_code); + unary("tan", &internal::guarded_unary_op, type_code); + unary("arcsin", &internal::guarded_unary_op, type_code); + unary("arccos", &internal::guarded_unary_op, type_code); + unary("arctan", &internal::guarded_unary_op, type_code); + unary("sinh", &internal::guarded_unary_op, type_code); + unary("cosh", &internal::guarded_unary_op, type_code); + unary("tanh", &internal::guarded_unary_op, type_code); + unary("arcsinh", &internal::guarded_unary_op, type_code); + unary("arccosh", &internal::guarded_unary_op, type_code); + unary("arctanh", &internal::guarded_unary_op, type_code); + + // absolute: real -> real, complex -> real (magnitude) + unary("absolute", &internal::guarded_unary_op_out, real_code); + + // predicates -> bool + unary("isnan", &internal::guarded_unary_op_out, bool_code); + unary("isinf", &internal::guarded_unary_op_out, bool_code); + unary("isfinite", &internal::guarded_unary_op_out, bool_code); + + // rint is the one rounding ufunc numpy defines for complex too + // (component-wise) — np.round dispatches through it + unary("rint", &internal::guarded_unary_op, type_code); + + if constexpr (WithOrderingComparitors) // the ordering-dependent set; NOT defined for complex types { - int types[3] = {type_code, type_code, Register::getTypeCode()}; - registerGuardedLoop(numpy, "equal", type_code, - &internal::guarded_compare_op, types, 3); - registerGuardedLoop(numpy, "not_equal", type_code, - &internal::guarded_compare_op, types, 3); - - if constexpr (WithOrderingComparitors) // NOT defined for complex types + binary("greater", &internal::guarded_compare_op, bool_code); + binary("less", &internal::guarded_compare_op, bool_code); + binary("greater_equal", &internal::guarded_compare_op, bool_code); + binary("less_equal", &internal::guarded_compare_op, bool_code); + + // mixed mp-vs-float64 orderings, both operand orders: the tolerance + // idiom `np.abs(a - b) < 1e-10`. Orderings ONLY -- see + // guarded_mixed_compare_op for why equality stays mp-vs-mp. + auto mixed_ordering = [&](char const* name, PyUFuncGenericFunction fwd, + PyUFuncGenericFunction rev) { - registerGuardedLoop(numpy, "greater", type_code, - &internal::guarded_compare_op, types, 3); - registerGuardedLoop(numpy, "less", type_code, - &internal::guarded_compare_op, types, 3); - registerGuardedLoop(numpy, "greater_equal", type_code, - &internal::guarded_compare_op, types, 3); - registerGuardedLoop(numpy, "less_equal", type_code, - &internal::guarded_compare_op, types, 3); - } - } - - // Unary operators - { - int types[2] = {type_code, type_code}; - registerGuardedLoop(numpy, "negative", type_code, - &internal::guarded_unary_op, types, 2); - registerGuardedLoop(numpy, "square", type_code, - &internal::guarded_unary_op, types, 2); - registerGuardedLoop(numpy, "sqrt", type_code, - &internal::guarded_unary_op, types, 2); + int types_td[3] = {type_code, NPY_DOUBLE, bool_code}; + int types_dt[3] = {NPY_DOUBLE, type_code, bool_code}; + registerGuardedLoop(numpy, name, type_code, fwd, types_td, 3); + registerGuardedLoop(numpy, name, type_code, rev, types_dt, 3); + }; + mixed_ordering("greater", + &internal::guarded_mixed_compare_op, + &internal::guarded_mixed_compare_op); + mixed_ordering("less", + &internal::guarded_mixed_compare_op, + &internal::guarded_mixed_compare_op); + mixed_ordering("greater_equal", + &internal::guarded_mixed_compare_op, + &internal::guarded_mixed_compare_op); + mixed_ordering("less_equal", + &internal::guarded_mixed_compare_op, + &internal::guarded_mixed_compare_op); + + // real-only unary: rounding family (sans rint, registered for both + // above), fabs, real-only transcendentals + unary("floor", &internal::guarded_unary_op, type_code); + unary("ceil", &internal::guarded_unary_op, type_code); + unary("trunc", &internal::guarded_unary_op, type_code); + unary("fabs", &internal::guarded_unary_op, type_code); + unary("exp2", &internal::guarded_unary_op, type_code); + unary("log2", &internal::guarded_unary_op, type_code); + unary("expm1", &internal::guarded_unary_op, type_code); + unary("log1p", &internal::guarded_unary_op, type_code); + unary("cbrt", &internal::guarded_unary_op, type_code); + + unary("signbit", &internal::guarded_unary_op_out, bool_code); + + // real-only binary + binary("arctan2", &internal::guarded_binary_op, type_code); + binary("hypot", &internal::guarded_binary_op, type_code); + binary("copysign", &internal::guarded_binary_op, type_code); + binary("fmod", &internal::guarded_binary_op, type_code); + binary("remainder", &internal::guarded_binary_op, type_code); + binary("floor_divide", &internal::guarded_binary_op, type_code); + binary("minimum", &internal::guarded_binary_op, type_code); + binary("maximum", &internal::guarded_binary_op, type_code); + binary("fmin", &internal::guarded_binary_op, type_code); + binary("fmax", &internal::guarded_binary_op, type_code); } Py_DECREF(numpy); diff --git a/python_bindings/src/mpfr_export.cpp b/python_bindings/src/mpfr_export.cpp index e93350c0c..aeb567adc 100644 --- a/python_bindings/src/mpfr_export.cpp +++ b/python_bindings/src/mpfr_export.cpp @@ -260,10 +260,30 @@ namespace bertini{ using boost::multiprecision::imag; real_mp (*reeeal)(const T&) = &boost::multiprecision::real; - real_mp (*imaaag)(const T&) = &boost::multiprecision::real; + // regression note: this was `&boost::multiprecision::real` (copy-paste), + // so mp.imag(z) returned the REAL part. + real_mp (*imaaag)(const T&) = &boost::multiprecision::imag; def("real",reeeal, (arg("val")), "get the real part"); //,return_value_policy() def("imag",imaaag, (arg("val")), "get the imaginary part"); //,return_value_policy() + // array overloads: numpy's .real/.imag ndarray attributes (and np.real / + // np.imag / np.angle) return silently WRONG values for legacy user + // dtypes -- numpy does not know complex_mp is complex-like, so .real + // returns the complex values themselves and .imag returns zeros. These + // element-wise overloads are the sanctioned array component accessors. + def("real", +[](Vec const& v) { + Vec r(v.size()); + for (Eigen::Index i = 0; i < v.size(); ++i) + r(i) = v(i).real(); + return r; + }, (arg("val")), "get the real parts of an array of complex numbers, as an array of real_mp"); + def("imag", +[](Vec const& v) { + Vec r(v.size()); + for (Eigen::Index i = 0; i < v.size(); ++i) + r(i) = v(i).imag(); + return r; + }, (arg("val")), "get the imaginary parts of an array of complex numbers, as an array of real_mp"); + // and then a few more free functions // def("abs2",&T::abs2); @@ -277,6 +297,14 @@ namespace bertini{ real_mp (*aaaarg)(const T&) = &boost::multiprecision::arg; def("arg",aaaarg, "the argument, or the angle from 0. beware the branch cut."); + // array overload: the np.angle replacement (np.angle goes through the + // broken .imag/.real ndarray attributes -- see the note at real/imag). + def("arg", +[](Vec const& v) { + Vec r(v.size()); + for (Eigen::Index i = 0; i < v.size(); ++i) + r(i) = boost::multiprecision::arg(v(i)); + return r; + }, (arg("val")), "the arguments (angles from 0) of an array of complex numbers, as an array of real_mp. beware the branch cut."); // def("square",&square); // def("cube",&cube); @@ -478,9 +506,14 @@ namespace bertini{ eigenpy::registerNewType(); eigenpy::HardenSetitem(); // zero slots before assignment — see eigenpy_interaction.hpp & ADR-0003 eigenpy::HardenDotfunc(); // np.dot/np.inner guard — see eigenpy_interaction.hpp - // guarded loops (real type — orderings included); eigenpy's registerCommonUfunc - // loops read input slots unguarded and crash on never-written np.zeros/np.empty slots. - eigenpy::registerGuardedUfunct(); + eigenpy::HardenCompare(); // element compare slot: np.sort/argsort/searchsorted + eigenpy::HardenArgMinMax(); // argmax/argmin slots + + // casts must be registered BEFORE the ufunc loops: registering the + // mixed mp-vs-double ordering loops makes numpy query the mp<->double + // casts, and a cast first registered after it has been queried is + // ignored (numpy RuntimeWarning "registered/modified ... after the + // cast had been used"). // you can convert from integer types with no fear eigenpy::registerCast(true); @@ -497,6 +530,13 @@ namespace bertini{ // both directions are scary. eigenpy::registerCast(false); eigenpy::registerCast(false); + // explicit down-conversion arr.astype(complex) — unsafe, you consciously truncate + eigenpy::registerCast>(false); + + // guarded loops (real type — orderings included, plus the mp-vs-double + // tolerance orderings); eigenpy's registerCommonUfunc loops read input + // slots unguarded and crash on never-written np.zeros/np.empty slots. + eigenpy::registerGuardedUfunct(); IMPLICITLY_CONVERTIBLE(int,T); @@ -573,8 +613,8 @@ namespace bertini{ eigenpy::registerNewType(); eigenpy::HardenSetitem(); // zero slots before assignment — see eigenpy_interaction.hpp & ADR-0003 eigenpy::HardenDotfunc(); // np.dot/np.inner guard — see eigenpy_interaction.hpp - eigenpy::registerUfunct_without_comparitors(); + // casts before ufunc loops — see the note in ExposeFloat. // you can safely convert from integer types to Complex's, there's no loss possible eigenpy::registerCast(true); @@ -590,10 +630,14 @@ namespace bertini{ // these conversions are unsafe. you can, but you probably shouldn't eigenpy::registerCast(false); eigenpy::registerCast(false); + // explicit down-conversion arr.astype(complex) — unsafe, you consciously truncate + eigenpy::registerCast>(false); // it's ok to convert from variable precision Float to Complex, that's ok! eigenpy::registerCast(true); + eigenpy::registerUfunct_without_comparitors(); + IMPLICITLY_CONVERTIBLE(int,T); IMPLICITLY_CONVERTIBLE(long,T); IMPLICITLY_CONVERTIBLE(int64_t,T);