Coverage for src/chebpy/maps.py: 100%
184 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-01 13:43 +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, cast
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.
54 Examples:
55 >>> p = MapParams()
56 >>> p.L, p.alpha
57 (8.0, 1.0)
59 Both parameters must be strictly positive:
61 >>> MapParams(L=-1.0)
62 Traceback (most recent call last):
63 ...
64 ValueError: require L > 0
65 >>> MapParams(alpha=0.0)
66 Traceback (most recent call last):
67 ...
68 ValueError: require alpha > 0
70 A larger ``L`` shrinks the gap between the map's image and the
71 clustered endpoint:
73 >>> from chebpy.maps import SingleSlitMap
74 >>> wide = SingleSlitMap(0.0, 1.0, MapParams(L=4.0), side="left")
75 >>> tight = SingleSlitMap(0.0, 1.0, MapParams(L=8.0), side="left")
76 >>> bool(tight.gap < wide.gap)
77 True
78 """
80 L: float = 8.0
81 alpha: float = 1.0
83 def __post_init__(self) -> None:
84 """Validate that ``L`` and ``alpha`` are strictly positive."""
85 if not (self.L > 0):
86 msg = "require L > 0"
87 raise ValueError(msg)
88 if not (self.alpha > 0):
89 msg = "require alpha > 0"
90 raise ValueError(msg)
93def _as_array(x: float | np.ndarray) -> tuple[np.ndarray, bool]:
94 """Coerce ``x`` to a numpy array, remembering whether the input was scalar.
96 Returns:
97 tuple[numpy.ndarray, bool]: A ``(array, was_scalar)`` pair. ``was_scalar``
98 is ``True`` if ``x`` is a 0-d / Python scalar input, in which case
99 callers should ``.item()`` the result before returning.
100 """
101 arr = np.asarray(x, dtype=float)
102 return arr, arr.ndim == 0
105class SingleSlitMap:
106 """Paper-faithful semi-infinite slit-strip map ``phi_S``.
108 Maps the reference interval ``[-1, 1]`` to (approximately) the
109 logical interval ``[a, b]`` via the composition of an affine scaling
110 and the inverse semi-infinite slit-strip map of Adcock & Richardson
111 (arXiv:1305.2643). The map's derivative vanishes
112 super-algebraically at the chosen clustered endpoint, so a function
113 with an algebraic or logarithmic singularity at that endpoint becomes
114 analytic in the reference variable ``t``.
116 Mathematical form (with ``side='left'``)::
118 s(y) = L * (y - 1) / 2 # [-1, 1] -> [-L, 0]
119 gamma = (alpha / pi) * log(exp(pi / alpha) - 1)
120 u(s) = (alpha / pi) * log(1 + exp(pi * (s + gamma) / alpha))
121 x = a + (b - a) * u(s(y))
123 The shift ``gamma`` is chosen so that ``u(0) = 1``; ``u(-L)`` is a
124 small positive number :attr:`gap_unit`, equal to the fraction of the
125 interval ``[a, b]`` that is *not* covered by ``formap([-1, 1])``.
127 For ``side='right'`` the analogous reflection ``x = b - (b - a) *
128 u(-s(y))`` is used so the cluster lies at ``b``.
130 Args:
131 a: Left endpoint of the logical interval.
132 b: Right endpoint of the logical interval (must satisfy ``b > a``).
133 params: A :class:`MapParams` instance carrying ``(L, alpha)``.
134 If ``None``, ``MapParams()`` defaults are used.
135 side: Either ``"left"`` or ``"right"``.
137 Raises:
138 ValueError: If ``b <= a`` or ``side`` is not one of the supported
139 values.
140 """
142 def __init__(
143 self,
144 a: float,
145 b: float,
146 params: MapParams | None = None,
147 *,
148 side: str = "left",
149 ) -> None:
150 """Initialise a semi-infinite slit-strip clustering map.
152 See the class docstring for parameter descriptions.
153 """
154 if not (b > a):
155 msg = "require b > a"
156 raise ValueError(msg)
157 if side not in ("left", "right"):
158 msg = "side must be 'left' or 'right'"
159 raise ValueError(msg)
160 self.a = float(a)
161 self.b = float(b)
162 self.params = params if params is not None else MapParams()
163 self.side = side
164 # Pre-compute the shift so that u(0) = 1 exactly.
165 a_p = self.params.alpha
166 # gamma = (alpha/pi) * log(exp(pi/alpha) - 1); log1p(-exp(-pi/alpha)) is the
167 # numerically stable form of log(exp(pi/alpha) - 1) - pi/alpha.
168 pi_over_alpha = np.pi / a_p
169 self._gamma = (a_p / np.pi) * (pi_over_alpha + np.log1p(-np.exp(-pi_over_alpha)))
171 # Conveniences ------------------------------------------------------
172 @property
173 def alpha(self) -> float:
174 """Strip half-width ``alpha`` (equal to ``self.params.alpha``)."""
175 return self.params.alpha
177 @property
178 def L(self) -> float:
179 """Truncation length ``L`` (equal to ``self.params.L``)."""
180 return self.params.L
182 @property
183 def support(self) -> tuple[float, float]:
184 """Return the *nominal* logical support ``(a, b)`` as a plain tuple.
186 This is the interval the user requested; the actual image of
187 :meth:`formap` falls short by :attr:`gap` at the clustered
188 endpoint (negligible at the default ``L = 8``).
189 """
190 return (self.a, self.b)
192 @property
193 def gap_unit(self) -> float:
194 """Fraction of ``[0, 1]`` not covered by the underlying truncated map.
196 Equal to ``u(-L) = (alpha / pi) * log(1 + exp(pi * (-L + gamma) / alpha))``;
197 this is the unit-interval shortfall before the affine scaling to
198 ``[a, b]``.
199 """
200 a_p = self.params.alpha
201 z = np.pi * (-self.params.L + self._gamma) / a_p
202 return float((a_p / np.pi) * np.logaddexp(0.0, z))
204 @property
205 def gap(self) -> float:
206 """Distance between :meth:`formap` ``(-1)`` (or ``(+1)``) and the clustered endpoint."""
207 return (self.b - self.a) * self.gap_unit
209 # Internal helpers --------------------------------------------------
210 def _u_of_s(self, s: np.ndarray) -> np.ndarray:
211 """Compute ``u(s) = (alpha/pi) * log(1 + exp(pi(s+gamma)/alpha))`` stably."""
212 a_p = self.params.alpha
213 z = np.pi * (s + self._gamma) / a_p
214 return cast(np.ndarray, (a_p / np.pi) * np.logaddexp(0.0, z))
216 def _du_ds(self, s: np.ndarray) -> np.ndarray:
217 """Derivative ``du/ds = sigmoid(pi(s+gamma)/alpha)``."""
218 a_p = self.params.alpha
219 z = np.pi * (s + self._gamma) / a_p
220 # 1 / (1 + exp(-z)), numerically stable.
221 return cast(np.ndarray, np.where(z >= 0.0, 1.0 / (1.0 + np.exp(-z)), np.exp(z) / (1.0 + np.exp(z))))
223 def _s_of_u(self, u: np.ndarray) -> np.ndarray:
224 """Inverse of :meth:`_u_of_s`. Maps ``u in (0, 1]`` to ``s in (-inf, 0]``.
226 Uses the identity ``log(exp(w) - 1) = w + log1p(-exp(-w))`` for
227 ``w = pi * u / alpha`` to retain precision near ``u = 1``.
228 """
229 a_p = self.params.alpha
230 w = np.pi * u / a_p
231 # log(exp(w) - 1) = w + log1p(-exp(-w)).
232 log_em1 = w + np.log1p(-np.exp(-w))
233 return cast(np.ndarray, (a_p / np.pi) * log_em1 - self._gamma)
235 # IntervalMap protocol ---------------------------------------------
236 def formap(self, y: float | np.ndarray) -> Any:
237 """Map ``y in [-1, 1]`` to ``x`` clustered near the chosen endpoint."""
238 t, scalar = _as_array(y)
239 L = self.params.L
240 if self.side == "left":
241 s = L * (t - 1.0) * 0.5 # [-1, 1] -> [-L, 0]
242 u = self._u_of_s(s)
243 x = self.a + (self.b - self.a) * u
244 else: # side == "right": reflect.
245 s = L * (-t - 1.0) * 0.5 # [-1, 1] -> [0, -L]
246 u = self._u_of_s(s)
247 x = self.b - (self.b - self.a) * u
248 if scalar:
249 return float(x)
250 return x
252 def invmap(self, x: float | np.ndarray) -> Any:
253 """Map ``x in [a, b]`` back to ``y in [-1, 1]`` (analytical inverse)."""
254 xa, scalar = _as_array(x)
255 L = self.params.L
256 gap_u = self.gap_unit
257 if self.side == "left":
258 u = (xa - self.a) / (self.b - self.a)
259 # Points within the (tiny) gap near the clustered endpoint are
260 # mapped to t = -1; points at/beyond b map to t = +1. This keeps
261 # evaluation inside [-1, 1] so the onefun is not extrapolated.
262 u_safe = np.clip(u, gap_u, 1.0)
263 s = self._s_of_u(u_safe)
264 t = 2.0 * s / L + 1.0
265 else: # side == "right"
266 u = (self.b - xa) / (self.b - self.a)
267 u_safe = np.clip(u, gap_u, 1.0)
268 s = self._s_of_u(u_safe)
269 t = -(2.0 * s / L + 1.0)
270 if scalar:
271 return float(t)
272 return t
274 def drvmap(self, y: float | np.ndarray) -> Any:
275 """Return ``dx/dy`` of :meth:`formap` evaluated at ``y``."""
276 t, scalar = _as_array(y)
277 L = self.params.L
278 scale = self.b - self.a
279 if self.side == "left":
280 s = L * (t - 1.0) * 0.5
281 dxdy = scale * self._du_ds(s) * (L * 0.5)
282 else:
283 s = L * (-t - 1.0) * 0.5
284 # Reflection introduces a sign flip in ds/dy, plus the outer x = b - (b-a)*u
285 # introduces another, so they cancel.
286 dxdy = scale * self._du_ds(s) * (L * 0.5)
287 if scalar:
288 return float(dxdy)
289 return dxdy
291 def __repr__(self) -> str:
292 """Return a developer-friendly representation."""
293 return f"SingleSlitMap(a={self.a!r}, b={self.b!r}, params={self.params!r}, side={self.side!r})"
296class DoubleSlitMap:
297 """Paper-faithful infinite two-slit-strip map ``psi_S``.
299 Maps the reference interval ``[-1, 1]`` to (approximately) the
300 logical interval ``[a, b]`` via the composition of an affine scaling
301 and the inverse infinite two-slit-strip map of Adcock & Richardson
302 (arXiv:1305.2643). The map's derivative vanishes
303 super-algebraically at *both* endpoints ``t = ±1``, so functions
304 with simultaneous endpoint singularities at ``a`` and ``b`` (e.g.
305 ``sqrt((x - a)(b - x))``) become analytic in ``t``.
307 Mathematical form::
309 s(y) = L * y # [-1, 1] -> [-L, L]
310 u(s) = (alpha/pi) * [logaddexp(0, pi(s+1/2)/alpha)
311 - logaddexp(0, pi(s-1/2)/alpha)]
312 x = a + (b - a) * u(s(y))
314 The construction satisfies ``u(0) = 1/2`` and ``u(±inf) = (1±1)/2``;
315 with finite ``L`` the image of ``[-1, 1]`` is short of both endpoints
316 by :attr:`gap` (negligible at the default ``L = 8``).
318 Args:
319 a: Left endpoint of the logical interval.
320 b: Right endpoint of the logical interval (must satisfy ``b > a``).
321 params: A :class:`MapParams` instance carrying ``(L, alpha)``.
322 If ``None``, ``MapParams()`` defaults are used.
324 Raises:
325 ValueError: If ``b <= a``.
326 """
328 def __init__(
329 self,
330 a: float,
331 b: float,
332 params: MapParams | None = None,
333 ) -> None:
334 """Initialise a symmetric two-slit-strip clustering map.
336 See the class docstring for parameter descriptions.
337 """
338 if not (b > a):
339 msg = "require b > a"
340 raise ValueError(msg)
341 self.a = float(a)
342 self.b = float(b)
343 self.params = params if params is not None else MapParams()
345 # Conveniences ------------------------------------------------------
346 @property
347 def alpha(self) -> float:
348 """Strip half-width ``alpha`` (equal to ``self.params.alpha``)."""
349 return self.params.alpha
351 @property
352 def L(self) -> float:
353 """Truncation length ``L`` (equal to ``self.params.L``)."""
354 return self.params.L
356 @property
357 def support(self) -> tuple[float, float]:
358 """Return the *nominal* logical support ``(a, b)`` as a plain tuple."""
359 return (self.a, self.b)
361 @property
362 def gap_unit(self) -> float:
363 """Unit-interval shortfall at each clustered endpoint (equal at both ends by symmetry)."""
364 a_p = self.params.alpha
365 L = self.params.L
366 # u(-L) = (alpha/pi)*[logaddexp(0, pi(-L+1/2)/alpha) - logaddexp(0, pi(-L-1/2)/alpha)]
367 z_plus = np.pi * (-L + 0.5) / a_p
368 z_minus = np.pi * (-L - 0.5) / a_p
369 return float((a_p / np.pi) * (np.logaddexp(0.0, z_plus) - np.logaddexp(0.0, z_minus)))
371 @property
372 def gap(self) -> float:
373 """Distance between :meth:`formap` ``(-1)`` and ``a`` (and symmetrically at ``b``)."""
374 return (self.b - self.a) * self.gap_unit
376 # Internal helpers --------------------------------------------------
377 def _u_of_s(self, s: np.ndarray) -> np.ndarray:
378 """Compute ``u(s)`` stably via ``logaddexp``."""
379 a_p = self.params.alpha
380 z_plus = np.pi * (s + 0.5) / a_p
381 z_minus = np.pi * (s - 0.5) / a_p
382 return (a_p / np.pi) * (np.logaddexp(0.0, z_plus) - np.logaddexp(0.0, z_minus))
384 def _du_ds(self, s: np.ndarray) -> np.ndarray:
385 """Derivative ``du/ds = sigmoid(z_plus) - sigmoid(z_minus)``."""
386 a_p = self.params.alpha
387 z_plus = np.pi * (s + 0.5) / a_p
388 z_minus = np.pi * (s - 0.5) / a_p
389 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)))
390 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)))
391 return sig_plus - sig_minus
393 # IntervalMap protocol ---------------------------------------------
394 def formap(self, y: float | np.ndarray) -> Any:
395 """Map ``y in [-1, 1]`` to ``x`` clustered near both endpoints."""
396 t, scalar = _as_array(y)
397 s = self.params.L * t
398 u = self._u_of_s(s)
399 x = self.a + (self.b - self.a) * u
400 if scalar:
401 return float(x)
402 return x
404 def invmap(self, x: float | np.ndarray) -> Any:
405 """Map ``x in [a, b]`` back to ``y in [-1, 1]`` by Newton iteration.
407 ``u(s)`` has no closed-form inverse in elementary functions, so we
408 solve ``u(s) = u_target`` by a few Newton steps starting from a
409 well-conditioned initial guess based on the asymptotic behaviour
410 ``u(s) ~ 1/2 + s/(2*alpha)`` near ``s = 0`` and the saturating
411 single-slit limit elsewhere.
412 """
413 xa, scalar = _as_array(x)
414 u_target = (xa - self.a) / (self.b - self.a)
415 u_target = np.clip(u_target, np.finfo(float).tiny, 1.0 - np.finfo(float).tiny)
416 # Initial guess: invert the dominant single-slit branch.
417 # For u < 1/2 use the left slit; for u > 1/2 use the right slit's reflection.
418 a_p = self.params.alpha
419 with np.errstate(divide="ignore", invalid="ignore"):
420 # Left-branch guess: u ≈ (alpha/pi)*log(1 + exp(pi*(s+1/2)/alpha)) for s << 1/2.
421 # Solve: s_guess_L = (alpha/pi)*log(exp(pi*u/alpha) - 1) - 1/2.
422 w = np.pi * u_target / a_p
423 s_guess_L = (a_p / np.pi) * (w + np.log1p(-np.exp(-np.clip(w, 1e-30, None)))) - 0.5
424 # Right-branch guess: by symmetry u(s) = 1 - u(-s), so s_guess_R = -s_for_(1-u).
425 w_r = np.pi * (1.0 - u_target) / a_p
426 s_guess_R = -((a_p / np.pi) * (w_r + np.log1p(-np.exp(-np.clip(w_r, 1e-30, None)))) - 0.5)
427 s = np.where(u_target < 0.5, s_guess_L, s_guess_R)
428 # A handful of Newton iterations are enough for double precision across [-L, L].
429 for _ in range(40):
430 f_val = self._u_of_s(s) - u_target
431 df = self._du_ds(s)
432 # Guard against vanishing derivative at extreme s; these points won't move further.
433 step = np.where(df > 0.0, f_val / np.where(df > 0.0, df, 1.0), 0.0)
434 s = s - step
435 if np.all(np.abs(step) < 1e-15 * (1.0 + np.abs(s))):
436 break
437 t = s / self.params.L
438 if scalar:
439 return float(t)
440 return t
442 def drvmap(self, y: float | np.ndarray) -> Any:
443 """Return ``dx/dy`` of :meth:`formap` evaluated at ``y``."""
444 t, scalar = _as_array(y)
445 L = self.params.L
446 s = L * t
447 dxdy = (self.b - self.a) * self._du_ds(s) * L
448 if scalar:
449 return float(dxdy)
450 return dxdy
452 def __repr__(self) -> str:
453 """Return a developer-friendly representation."""
454 return f"DoubleSlitMap(a={self.a!r}, b={self.b!r}, params={self.params!r})"