|
| 1 | +############################################################################### |
| 2 | +# galpy.backend._coerce: backend DATA-coercion helpers. |
| 3 | +############################################################################### |
| 4 | +"""Backend data-coercion helpers: the single home for bringing numpy/Python |
| 5 | +data onto the active jax/torch backend. |
| 6 | +
|
| 7 | +PURPOSE |
| 8 | +------- |
| 9 | +This module does two related jobs: |
| 10 | +
|
| 11 | + * it brings numpy/Python *coordinate* data onto the active backend's array |
| 12 | + type (so e.g. ``torch.sqrt`` -- which rejects ``numpy.float64`` -- and |
| 13 | + ``Tensor`` arithmetic see real backend arrays), and |
| 14 | + * it anchors *stored numpy constants* (rotation matrices, lookup tables, a |
| 15 | + zero reference coordinate) onto the dtype/device of an input array, so the |
| 16 | + constant joins the computation as a same-dtype/same-device backend array. |
| 17 | +
|
| 18 | +It is the single home for data-coercion. Namespace *resolution* (which backend |
| 19 | +a call dispatches to) and the dtype/device *primitives* it builds on live in |
| 20 | +``_namespaces.py`` (``is_backend_array``, ``device_of``, ``asarray_on_device``, |
| 21 | +``match_input_dtype``); this module only consumes those primitives. |
| 22 | +
|
| 23 | +THE CORE INVARIANT |
| 24 | +------------------ |
| 25 | +Every function here is a STRICT PASS-THROUGH when ``xp is numpy``: it returns |
| 26 | +its inputs OBJECT-IDENTICALLY (no asarray, no copy, no dtype touch). This is |
| 27 | +what keeps the numpy code path BYTE-IDENTICAL to galpy's historical behaviour. |
| 28 | +Any new coercion helper added to this module MUST preserve this invariant -- |
| 29 | +guard the work behind ``if xp is numpy: return <inputs unchanged>`` first. |
| 30 | +
|
| 31 | +WHEN TO USE EACH |
| 32 | +---------------- |
| 33 | + * ``coerce_coords(xp, *coords)`` -- at the PUBLIC INPUT BOUNDARY (the |
| 34 | + ``@potential_physical_input`` decorator) to bring coordinate arguments onto |
| 35 | + the backend: plain Python/int scalars become float64 (galpy's interior |
| 36 | + precision), float arrays keep their dtype (so the float32 exit-cast policy |
| 37 | + still applies), and ``None`` passes through. |
| 38 | + * ``promote_scalars(xp, *vals)`` -- INSIDE coordinate transforms to promote |
| 39 | + plain Python scalars sitting alongside array arguments, anchored on the |
| 40 | + dtype/device of the first array, so mixed scalar/array inputs work on a |
| 41 | + backend whose functions require arrays. |
| 42 | + * ``as_backend_constant(xp, value, ref)`` -- to anchor a single STORED numpy |
| 43 | + constant (a rotation matrix, an offset, a table) on a backend ``ref`` array |
| 44 | + derived from the coordinate inputs. |
| 45 | + * ``zeros_like_backend(xp, R)`` -- for a backend ZERO reference coordinate |
| 46 | + (e.g. the ``z = 0`` plane a spherical-in-disguise wrapper feeds its wrapped |
| 47 | + potential). |
| 48 | +
|
| 49 | +WHY float64-INTERIOR / DEVICE-ANCHORING |
| 50 | +--------------------------------------- |
| 51 | +galpy computes in float64 internally: a bare ``asarray`` of a Python float |
| 52 | +yields torch float32 and silently misses galpy's tolerances, so plain scalars |
| 53 | +are lifted to ``xp.float64`` while genuine float arrays keep their own dtype. |
| 54 | +Anchoring constants and promoted scalars on an input array's dtype/device keeps |
| 55 | +the whole computation on one device and at one precision, which is required for |
| 56 | +torch (cross-device / mixed-dtype ops raise) and correct for jax. |
| 57 | +""" |
| 58 | + |
| 59 | +import numpy |
| 60 | + |
| 61 | +from ._namespaces import ( |
| 62 | + _is_floating_dtype, |
| 63 | + asarray_on_device, |
| 64 | + device_of, |
| 65 | +) |
| 66 | + |
| 67 | + |
| 68 | +def coerce_coords(xp, *coords): |
| 69 | + """Bring coordinate inputs onto the active backend's array type. |
| 70 | +
|
| 71 | + The dominant non-numpy failure mode is "the namespace resolved to a backend |
| 72 | + (forced harness, or a user mixing a backend tensor with a numpy/python arg) |
| 73 | + but a coordinate is still numpy/python", which torch rejects strictly |
| 74 | + (``torch.sqrt(numpy.float64)`` raises; ``numpy.ndarray * Tensor`` raises). |
| 75 | + Coercing the coordinates to backend arrays once, at the public input |
| 76 | + boundary, fixes it for every potential at once. |
| 77 | +
|
| 78 | + Rules (applied only when the backend is NOT numpy): |
| 79 | + * ``None`` is passed through (axisymmetric ``phi=None`` etc.). |
| 80 | + * a coordinate that already carries a *floating* dtype (a numpy/backend |
| 81 | + float32/float64 array or scalar) is moved onto the backend with its |
| 82 | + dtype PRESERVED, so the float32/exit-cast policy (``match_input_dtype``) |
| 83 | + still applies. |
| 84 | + * a plain Python scalar (``1.0``/``1``) or an integer array is brought to |
| 85 | + the backend's float64 -- galpy's interior precision; a bare ``asarray`` |
| 86 | + of a Python float would give torch float32 and miss the tolerances. |
| 87 | +
|
| 88 | + The numpy backend is a strict pass-through (``coords`` returned object- |
| 89 | + identical) -> the numpy path stays byte-identical. |
| 90 | + """ |
| 91 | + if xp is numpy: |
| 92 | + return coords |
| 93 | + dev = device_of(*coords) |
| 94 | + out = [] |
| 95 | + for c in coords: |
| 96 | + if c is None: |
| 97 | + out.append(c) |
| 98 | + continue |
| 99 | + dt = getattr(c, "dtype", None) |
| 100 | + if dt is not None and _is_floating_dtype(dt): |
| 101 | + out.append(asarray_on_device(xp, c, dev)) # preserve float dtype |
| 102 | + else: |
| 103 | + out.append(asarray_on_device(xp, c, dev, dtype=xp.float64)) |
| 104 | + return tuple(out) |
| 105 | + |
| 106 | + |
| 107 | +def promote_scalars(xp, *vals): |
| 108 | + """Promote plain Python scalars among ``vals`` to the active non-numpy |
| 109 | + namespace, anchored on the dtype/device of the first array argument, so |
| 110 | + that e.g. torch functions -- which require Tensors -- accept the mixed |
| 111 | + scalar/array inputs that the numpy path has always supported. The numpy |
| 112 | + path passes everything through untouched (byte-identical).""" |
| 113 | + if xp is numpy: |
| 114 | + return vals |
| 115 | + ref = next((v for v in vals if hasattr(v, "ndim")), None) |
| 116 | + if ref is None: |
| 117 | + # Nothing to anchor on (e.g. all-scalar inputs under a forced backend |
| 118 | + # default): pass through, the namespace's functions handle scalars |
| 119 | + return vals |
| 120 | + dtype = getattr(ref, "dtype", None) |
| 121 | + device = getattr(ref, "device", None) |
| 122 | + |
| 123 | + def _promote(v): |
| 124 | + if hasattr(v, "ndim"): |
| 125 | + return v |
| 126 | + try: |
| 127 | + return xp.asarray(v, dtype=dtype, device=device) |
| 128 | + except (TypeError, ValueError): |
| 129 | + # namespace without a device kwarg, or one that rejects the device |
| 130 | + # value (array-api jax exposes .device as the string 'cpu', which |
| 131 | + # jnp.asarray(device=...) refuses with ValueError). jax tracks device |
| 132 | + # automatically, so a plain asarray is correct; torch's .device is a |
| 133 | + # real object and never hits this fallback. |
| 134 | + return xp.asarray(v, dtype=dtype) |
| 135 | + |
| 136 | + return tuple(_promote(v) for v in vals) |
| 137 | + |
| 138 | + |
| 139 | +def as_backend_constant(xp, value, ref): |
| 140 | + """Bring a stored numpy constant (rotation matrix / offset) into the active |
| 141 | + namespace, anchored on the dtype/device of ``ref`` (a backend array derived |
| 142 | + from the coordinate inputs). The numpy path passes the stored array through |
| 143 | + untouched (byte-identical).""" |
| 144 | + if xp is numpy: |
| 145 | + return value |
| 146 | + dtype = getattr(ref, "dtype", None) |
| 147 | + device = getattr(ref, "device", None) |
| 148 | + try: |
| 149 | + return xp.asarray(value, dtype=dtype, device=device) |
| 150 | + except TypeError: # pragma: no cover - namespace without device= kwarg |
| 151 | + return xp.asarray(value, dtype=dtype) |
| 152 | + |
| 153 | + |
| 154 | +def zeros_like_backend(xp, R): |
| 155 | + """The numpy path passes the plain scalar through untouched |
| 156 | + (byte-identical); on a non-numpy backend the z = 0 reference |
| 157 | + coordinate is anchored on the inputs so the wrapped potential sees a |
| 158 | + backend array (torch functions require Tensors) on the right |
| 159 | + device/dtype.""" |
| 160 | + return 0.0 if xp is numpy else xp.zeros_like(R) |
0 commit comments