Coverage for src / chebpy / maps.py: 99%
180 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-05 07:22 +0000
1"""Non-affine bijective maps between [-1, 1] and a logical interval [a, b].
3These maps are concrete implementations of the
4:class:`~chebpy.utilities.IntervalMap` protocol that cluster grid points
5exponentially towards one or both endpoints, following the slit-strip
6constructions of Adcock & Richardson, *New exponential variable transform
7methods for functions with endpoint singularities*, SIAM J. Numer. Anal.
852(4), 1887–1912, 2014 (doi:10.1137/130920460; arXiv:1305.2643).
10Two families are provided:
12* :class:`SingleSlitMap` — semi-infinite slit-strip map ``phi_S`` of
13 Adcock & Richardson. Exponential clustering at one chosen endpoint
14 (``side='left'`` clusters near ``a``; ``side='right'`` near ``b``).
16* :class:`DoubleSlitMap` — infinite two-slit-strip map ``psi_S`` of
17 Adcock & Richardson. Exponential clustering at both endpoints.
19Each map is parameterised by a :class:`MapParams` ``(L, alpha)`` pair,
20where ``alpha > 0`` is the strip half-width and ``L > 0`` is the
21truncation length applied to the underlying conformal map. The forward
22map composes an affine scaling with the truncated paper map; with
23finite ``L`` the image of ``[-1, 1]`` under :meth:`formap` falls a
24distance ``gap`` short of the clustered endpoint(s). With the default
25``L = 8.0`` the gap is below ``1e-10`` and is invisible at working
26precision; smaller ``L`` (closer to the paper's empirical optimum
27``L ~ 1``) shrinks the resolved interval visibly but improves the
28convergence rate of the mapped Chebyshev expansion. See
29:attr:`SingleSlitMap.gap` / :attr:`DoubleSlitMap.gap` for the exact
30shortfall.
32Both classes are pure-Python, NumPy-vectorised, and stateless apart from
33their constructor parameters; they carry no ``Onefun`` payload.
34"""
36from dataclasses import dataclass
37from typing import Any
39import numpy as np
42@dataclass(frozen=True, slots=True)
43class MapParams:
44 """Parameters for the slit-strip clustering maps.
46 Args:
47 L: Positive truncation length. With finite ``L`` the map's image
48 falls short of the clustered endpoint(s) by a tiny ``gap``;
49 larger ``L`` shrinks the gap (default ``L = 8.0`` gives
50 ``gap < 1e-10`` at ``alpha = 1.0``).
51 alpha: Positive strip half-width. Controls the clustering
52 strength; smaller ``alpha`` clusters more aggressively.
53 """
55 L: float = 8.0
56 alpha: float = 1.0
58 def __post_init__(self) -> None:
59 """Validate that ``L`` and ``alpha`` are strictly positive."""
60 if not (self.L > 0):
61 msg = "require L > 0"
62 raise ValueError(msg)
63 if not (self.alpha > 0):
64 msg = "require alpha > 0"
65 raise ValueError(msg)
68def _as_array(x: float | np.ndarray) -> tuple[np.ndarray, bool]:
69 """Coerce ``x`` to a numpy array, remembering whether the input was scalar.
71 Returns:
72 tuple[numpy.ndarray, bool]: A ``(array, was_scalar)`` pair. ``was_scalar``
73 is ``True`` if ``x`` is a 0-d / Python scalar input, in which case
74 callers should ``.item()`` the result before returning.
75 """
76 arr = np.asarray(x, dtype=float)
77 return arr, arr.ndim == 0
80class SingleSlitMap:
81 """Paper-faithful semi-infinite slit-strip map ``phi_S``.
83 Maps the reference interval ``[-1, 1]`` to (approximately) the
84 logical interval ``[a, b]`` via the composition of an affine scaling
85 and the inverse semi-infinite slit-strip map of Adcock & Richardson
86 (arXiv:1305.2643). The map's derivative vanishes
87 super-algebraically at the chosen clustered endpoint, so a function
88 with an algebraic or logarithmic singularity at that endpoint becomes
89 analytic in the reference variable ``t``.
91 Mathematical form (with ``side='left'``)::
93 s(y) = L * (y - 1) / 2 # [-1, 1] -> [-L, 0]
94 gamma = (alpha / pi) * log(exp(pi / alpha) - 1)
95 u(s) = (alpha / pi) * log(1 + exp(pi * (s + gamma) / alpha))
96 x = a + (b - a) * u(s(y))
98 The shift ``gamma`` is chosen so that ``u(0) = 1``; ``u(-L)`` is a
99 small positive number :attr:`gap_unit`, equal to the fraction of the
100 interval ``[a, b]`` that is *not* covered by ``formap([-1, 1])``.
102 For ``side='right'`` the analogous reflection ``x = b - (b - a) *
103 u(-s(y))`` is used so the cluster lies at ``b``.
105 Args:
106 a: Left endpoint of the logical interval.
107 b: Right endpoint of the logical interval (must satisfy ``b > a``).
108 params: A :class:`MapParams` instance carrying ``(L, alpha)``.
109 If ``None``, ``MapParams()`` defaults are used.
110 side: Either ``"left"`` or ``"right"``.
112 Raises:
113 ValueError: If ``b <= a`` or ``side`` is not one of the supported
114 values.
115 """
117 def __init__(
118 self,
119 a: float,
120 b: float,
121 params: MapParams | None = None,
122 *,
123 side: str = "left",
124 ) -> None:
125 """Initialise a semi-infinite slit-strip clustering map.
127 See the class docstring for parameter descriptions.
128 """
129 if not (b > a):
130 msg = "require b > a"
131 raise ValueError(msg)
132 if side not in ("left", "right"):
133 msg = "side must be 'left' or 'right'"
134 raise ValueError(msg)
135 self.a = float(a)
136 self.b = float(b)
137 self.params = params if params is not None else MapParams()
138 self.side = side
139 # Pre-compute the shift so that u(0) = 1 exactly.
140 a_p = self.params.alpha
141 # gamma = (alpha/pi) * log(exp(pi/alpha) - 1); log1p(-exp(-pi/alpha)) is the
142 # numerically stable form of log(exp(pi/alpha) - 1) - pi/alpha.
143 pi_over_alpha = np.pi / a_p
144 self._gamma = (a_p / np.pi) * (pi_over_alpha + np.log1p(-np.exp(-pi_over_alpha)))
146 # Conveniences ------------------------------------------------------
147 @property
148 def alpha(self) -> float:
149 """Strip half-width ``alpha`` (equal to ``self.params.alpha``)."""
150 return self.params.alpha
152 @property
153 def L(self) -> float:
154 """Truncation length ``L`` (equal to ``self.params.L``)."""
155 return self.params.L
157 @property
158 def support(self) -> tuple[float, float]:
159 """Return the *nominal* logical support ``(a, b)`` as a plain tuple.
161 This is the interval the user requested; the actual image of
162 :meth:`formap` falls short by :attr:`gap` at the clustered
163 endpoint (negligible at the default ``L = 8``).
164 """
165 return (self.a, self.b)
167 @property
168 def gap_unit(self) -> float:
169 """Fraction of ``[0, 1]`` not covered by the underlying truncated map.
171 Equal to ``u(-L) = (alpha / pi) * log(1 + exp(pi * (-L + gamma) / alpha))``;
172 this is the unit-interval shortfall before the affine scaling to
173 ``[a, b]``.
174 """
175 a_p = self.params.alpha
176 z = np.pi * (-self.params.L + self._gamma) / a_p
177 return float((a_p / np.pi) * np.logaddexp(0.0, z))
179 @property
180 def gap(self) -> float:
181 """Distance between :meth:`formap` ``(-1)`` (or ``(+1)``) and the clustered endpoint."""
182 return (self.b - self.a) * self.gap_unit
184 # Internal helpers --------------------------------------------------
185 def _u_of_s(self, s: np.ndarray) -> np.ndarray:
186 """Compute ``u(s) = (alpha/pi) * log(1 + exp(pi(s+gamma)/alpha))`` stably."""
187 a_p = self.params.alpha
188 z = np.pi * (s + self._gamma) / a_p
189 return (a_p / np.pi) * np.logaddexp(0.0, z)
191 def _du_ds(self, s: np.ndarray) -> np.ndarray:
192 """Derivative ``du/ds = sigmoid(pi(s+gamma)/alpha)``."""
193 a_p = self.params.alpha
194 z = np.pi * (s + self._gamma) / a_p
195 # 1 / (1 + exp(-z)), numerically stable.
196 return np.where(z >= 0.0, 1.0 / (1.0 + np.exp(-z)), np.exp(z) / (1.0 + np.exp(z)))
198 def _s_of_u(self, u: np.ndarray) -> np.ndarray:
199 """Inverse of :meth:`_u_of_s`. Maps ``u in (0, 1]`` to ``s in (-inf, 0]``.
201 Uses the identity ``log(exp(w) - 1) = w + log1p(-exp(-w))`` for
202 ``w = pi * u / alpha`` to retain precision near ``u = 1``.
203 """
204 a_p = self.params.alpha
205 w = np.pi * u / a_p
206 # log(exp(w) - 1) = w + log1p(-exp(-w)).
207 log_em1 = w + np.log1p(-np.exp(-w))
208 return (a_p / np.pi) * log_em1 - self._gamma
210 # IntervalMap protocol ---------------------------------------------
211 def formap(self, y: float | np.ndarray) -> Any:
212 """Map ``y in [-1, 1]`` to ``x`` clustered near the chosen endpoint."""
213 t, scalar = _as_array(y)
214 L = self.params.L
215 if self.side == "left":
216 s = L * (t - 1.0) * 0.5 # [-1, 1] -> [-L, 0]
217 u = self._u_of_s(s)
218 x = self.a + (self.b - self.a) * u
219 else: # side == "right": reflect.
220 s = L * (-t - 1.0) * 0.5 # [-1, 1] -> [0, -L]
221 u = self._u_of_s(s)
222 x = self.b - (self.b - self.a) * u
223 if scalar:
224 return float(x)
225 return x
227 def invmap(self, x: float | np.ndarray) -> Any:
228 """Map ``x in [a, b]`` back to ``y in [-1, 1]`` (analytical inverse)."""
229 xa, scalar = _as_array(x)
230 L = self.params.L
231 gap_u = self.gap_unit
232 if self.side == "left":
233 u = (xa - self.a) / (self.b - self.a)
234 # Points within the (tiny) gap near the clustered endpoint are
235 # mapped to t = -1; points at/beyond b map to t = +1. This keeps
236 # evaluation inside [-1, 1] so the onefun is not extrapolated.
237 u_safe = np.clip(u, gap_u, 1.0)
238 s = self._s_of_u(u_safe)
239 t = 2.0 * s / L + 1.0
240 else: # side == "right"
241 u = (self.b - xa) / (self.b - self.a)
242 u_safe = np.clip(u, gap_u, 1.0)
243 s = self._s_of_u(u_safe)
244 t = -(2.0 * s / L + 1.0)
245 if scalar:
246 return float(t)
247 return t
249 def drvmap(self, y: float | np.ndarray) -> Any:
250 """Return ``dx/dy`` of :meth:`formap` evaluated at ``y``."""
251 t, scalar = _as_array(y)
252 L = self.params.L
253 scale = self.b - self.a
254 if self.side == "left":
255 s = L * (t - 1.0) * 0.5
256 dxdy = scale * self._du_ds(s) * (L * 0.5)
257 else:
258 s = L * (-t - 1.0) * 0.5
259 # Reflection introduces a sign flip in ds/dy, plus the outer x = b - (b-a)*u
260 # introduces another, so they cancel.
261 dxdy = scale * self._du_ds(s) * (L * 0.5)
262 if scalar:
263 return float(dxdy)
264 return dxdy
266 def __repr__(self) -> str: # pragma: no cover - trivial
267 """Return a developer-friendly representation."""
268 return f"SingleSlitMap(a={self.a!r}, b={self.b!r}, params={self.params!r}, side={self.side!r})"
271class DoubleSlitMap:
272 """Paper-faithful infinite two-slit-strip map ``psi_S``.
274 Maps the reference interval ``[-1, 1]`` to (approximately) the
275 logical interval ``[a, b]`` via the composition of an affine scaling
276 and the inverse infinite two-slit-strip map of Adcock & Richardson
277 (arXiv:1305.2643). The map's derivative vanishes
278 super-algebraically at *both* endpoints ``t = ±1``, so functions
279 with simultaneous endpoint singularities at ``a`` and ``b`` (e.g.
280 ``sqrt((x - a)(b - x))``) become analytic in ``t``.
282 Mathematical form::
284 s(y) = L * y # [-1, 1] -> [-L, L]
285 u(s) = (alpha/pi) * [logaddexp(0, pi(s+1/2)/alpha)
286 - logaddexp(0, pi(s-1/2)/alpha)]
287 x = a + (b - a) * u(s(y))
289 The construction satisfies ``u(0) = 1/2`` and ``u(±inf) = (1±1)/2``;
290 with finite ``L`` the image of ``[-1, 1]`` is short of both endpoints
291 by :attr:`gap` (negligible at the default ``L = 8``).
293 Args:
294 a: Left endpoint of the logical interval.
295 b: Right endpoint of the logical interval (must satisfy ``b > a``).
296 params: A :class:`MapParams` instance carrying ``(L, alpha)``.
297 If ``None``, ``MapParams()`` defaults are used.
299 Raises:
300 ValueError: If ``b <= a``.
301 """
303 def __init__(
304 self,
305 a: float,
306 b: float,
307 params: MapParams | None = None,
308 ) -> None:
309 """Initialise a symmetric two-slit-strip clustering map.
311 See the class docstring for parameter descriptions.
312 """
313 if not (b > a):
314 msg = "require b > a"
315 raise ValueError(msg)
316 self.a = float(a)
317 self.b = float(b)
318 self.params = params if params is not None else MapParams()
320 # Conveniences ------------------------------------------------------
321 @property
322 def alpha(self) -> float:
323 """Strip half-width ``alpha`` (equal to ``self.params.alpha``)."""
324 return self.params.alpha
326 @property
327 def L(self) -> float:
328 """Truncation length ``L`` (equal to ``self.params.L``)."""
329 return self.params.L
331 @property
332 def support(self) -> tuple[float, float]:
333 """Return the *nominal* logical support ``(a, b)`` as a plain tuple."""
334 return (self.a, self.b)
336 @property
337 def gap_unit(self) -> float:
338 """Unit-interval shortfall at each clustered endpoint (equal at both ends by symmetry)."""
339 a_p = self.params.alpha
340 L = self.params.L
341 # u(-L) = (alpha/pi)*[logaddexp(0, pi(-L+1/2)/alpha) - logaddexp(0, pi(-L-1/2)/alpha)]
342 z_plus = np.pi * (-L + 0.5) / a_p
343 z_minus = np.pi * (-L - 0.5) / a_p
344 return float((a_p / np.pi) * (np.logaddexp(0.0, z_plus) - np.logaddexp(0.0, z_minus)))
346 @property
347 def gap(self) -> float:
348 """Distance between :meth:`formap` ``(-1)`` and ``a`` (and symmetrically at ``b``)."""
349 return (self.b - self.a) * self.gap_unit
351 # Internal helpers --------------------------------------------------
352 def _u_of_s(self, s: np.ndarray) -> np.ndarray:
353 """Compute ``u(s)`` stably via ``logaddexp``."""
354 a_p = self.params.alpha
355 z_plus = np.pi * (s + 0.5) / a_p
356 z_minus = np.pi * (s - 0.5) / a_p
357 return (a_p / np.pi) * (np.logaddexp(0.0, z_plus) - np.logaddexp(0.0, z_minus))
359 def _du_ds(self, s: np.ndarray) -> np.ndarray:
360 """Derivative ``du/ds = sigmoid(z_plus) - sigmoid(z_minus)``."""
361 a_p = self.params.alpha
362 z_plus = np.pi * (s + 0.5) / a_p
363 z_minus = np.pi * (s - 0.5) / a_p
364 sig_plus = np.where(z_plus >= 0.0, 1.0 / (1.0 + np.exp(-z_plus)), np.exp(z_plus) / (1.0 + np.exp(z_plus)))
365 sig_minus = np.where(z_minus >= 0.0, 1.0 / (1.0 + np.exp(-z_minus)), np.exp(z_minus) / (1.0 + np.exp(z_minus)))
366 return sig_plus - sig_minus
368 # IntervalMap protocol ---------------------------------------------
369 def formap(self, y: float | np.ndarray) -> Any:
370 """Map ``y in [-1, 1]`` to ``x`` clustered near both endpoints."""
371 t, scalar = _as_array(y)
372 s = self.params.L * t
373 u = self._u_of_s(s)
374 x = self.a + (self.b - self.a) * u
375 if scalar:
376 return float(x)
377 return x
379 def invmap(self, x: float | np.ndarray) -> Any:
380 """Map ``x in [a, b]`` back to ``y in [-1, 1]`` by Newton iteration.
382 ``u(s)`` has no closed-form inverse in elementary functions, so we
383 solve ``u(s) = u_target`` by a few Newton steps starting from a
384 well-conditioned initial guess based on the asymptotic behaviour
385 ``u(s) ~ 1/2 + s/(2*alpha)`` near ``s = 0`` and the saturating
386 single-slit limit elsewhere.
387 """
388 xa, scalar = _as_array(x)
389 u_target = (xa - self.a) / (self.b - self.a)
390 u_target = np.clip(u_target, np.finfo(float).tiny, 1.0 - np.finfo(float).tiny)
391 # Initial guess: invert the dominant single-slit branch.
392 # For u < 1/2 use the left slit; for u > 1/2 use the right slit's reflection.
393 a_p = self.params.alpha
394 with np.errstate(divide="ignore", invalid="ignore"):
395 # Left-branch guess: u ≈ (alpha/pi)*log(1 + exp(pi*(s+1/2)/alpha)) for s << 1/2.
396 # Solve: s_guess_L = (alpha/pi)*log(exp(pi*u/alpha) - 1) - 1/2.
397 w = np.pi * u_target / a_p
398 s_guess_L = (a_p / np.pi) * (w + np.log1p(-np.exp(-np.clip(w, 1e-30, None)))) - 0.5
399 # Right-branch guess: by symmetry u(s) = 1 - u(-s), so s_guess_R = -s_for_(1-u).
400 w_r = np.pi * (1.0 - u_target) / a_p
401 s_guess_R = -((a_p / np.pi) * (w_r + np.log1p(-np.exp(-np.clip(w_r, 1e-30, None)))) - 0.5)
402 s = np.where(u_target < 0.5, s_guess_L, s_guess_R)
403 # A handful of Newton iterations are enough for double precision across [-L, L].
404 for _ in range(40):
405 f_val = self._u_of_s(s) - u_target
406 df = self._du_ds(s)
407 # Guard against vanishing derivative at extreme s; these points won't move further.
408 step = np.where(df > 0.0, f_val / np.where(df > 0.0, df, 1.0), 0.0)
409 s = s - step
410 if np.all(np.abs(step) < 1e-15 * (1.0 + np.abs(s))):
411 break
412 t = s / self.params.L
413 if scalar:
414 return float(t)
415 return t
417 def drvmap(self, y: float | np.ndarray) -> Any:
418 """Return ``dx/dy`` of :meth:`formap` evaluated at ``y``."""
419 t, scalar = _as_array(y)
420 L = self.params.L
421 s = L * t
422 dxdy = (self.b - self.a) * self._du_ds(s) * L
423 if scalar:
424 return float(dxdy)
425 return dxdy
427 def __repr__(self) -> str: # pragma: no cover - trivial
428 """Return a developer-friendly representation."""
429 return f"DoubleSlitMap(a={self.a!r}, b={self.b!r}, params={self.params!r})"