Coverage for src/chebpy/singfun.py: 100%
132 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-01 13:43 +0000
1"""Implementation of :class:`Singfun` for functions with endpoint singularities.
3A :class:`Singfun` represents a function on a bounded interval ``[a, b]``
4that is analytic on the open interval ``(a, b)`` but may have algebraic or
5logarithmic branch-type singularities at one or both endpoints. It is a
6sibling of :class:`~chebpy.bndfun.Bndfun` and
7:class:`~chebpy.compactfun.CompactFun` under
8:class:`~chebpy.classicfun.Classicfun`: the only structural novelty is that
9the bijective map between the storage variable ``t in [-1, 1]`` and the
10logical variable ``x in [a, b]`` is a non-affine, endpoint-clustering
11exponential transform (see :mod:`chebpy.maps`) rather than the affine
12:class:`~chebpy.utilities.Interval` map.
14For a function ``f`` with branch-type endpoint behaviour, the composition
15``f(m(t))`` is analytic in a Bernstein ellipse around ``[-1, 1]`` and is
16therefore resolved by ordinary Chebyshev interpolation in ``t``. All the
17existing :class:`~chebpy.classicfun.Classicfun` plumbing (``__call__``,
18``roots``, the binary operators) is reused unchanged via the ``map``
19property override; only the calculus operations (``sum``, ``cumsum``,
20``diff``) need bespoke implementations because the affine-Jacobian
21shortcuts in :class:`~chebpy.classicfun.Classicfun` no longer apply.
23This is a v1 implementation (Phase 3 of the singfun plan); see
24``docs/plans/03-singfun-mapped-integration.md`` for the broader design and
25the remaining closure / fallback work.
26"""
28from __future__ import annotations
30from typing import Any
32from .bndfun import Bndfun
33from .classicfun import Classicfun, techdict
34from .exceptions import InvalidSingularitySide, NotSubinterval
35from .maps import DoubleSlitMap, MapParams, SingleSlitMap
36from .settings import _preferences as prefs
37from .utilities import Interval, IntervalMap
40def _build_map(a: float, b: float, sing: str, params: MapParams) -> IntervalMap:
41 """Construct the appropriate non-affine map for the requested singularity pattern.
43 Args:
44 a: Left endpoint of the logical interval.
45 b: Right endpoint of the logical interval.
46 sing: One of ``"left"``, ``"right"``, or ``"both"``.
47 params: A :class:`MapParams` instance carrying ``(L, alpha)``.
49 Returns:
50 IntervalMap: A :class:`SingleSlitMap` (for ``"left"`` / ``"right"``)
51 or :class:`DoubleSlitMap` (for ``"both"``).
53 Raises:
54 InvalidSingularitySide: If ``sing`` is not one of the recognised values.
55 """
56 if sing in ("left", "right"):
57 return SingleSlitMap(a, b, params, side=sing)
58 if sing == "both":
59 return DoubleSlitMap(a, b, params)
60 msg = f"sing must be 'left', 'right', or 'both'; got {sing!r}"
61 raise InvalidSingularitySide(msg)
64class Singfun(Classicfun):
65 """Functions with branch-type endpoint singularities on a bounded interval.
67 A :class:`Singfun` stores:
69 - ``self.onefun`` (inherited): a standard :class:`~chebpy.onefun.Onefun`
70 (typically a :class:`~chebpy.chebtech.Chebtech`) on ``[-1, 1]``
71 representing ``f(m(t))``, which is analytic by construction.
72 - ``self._interval`` (inherited): an :class:`~chebpy.utilities.Interval`
73 ``Interval(a, b)`` carrying the logical support endpoints.
74 - ``self._map``: a non-affine :class:`~chebpy.utilities.IntervalMap`
75 (a :class:`~chebpy.maps.SingleSlitMap` or
76 :class:`~chebpy.maps.DoubleSlitMap`) — the actual bijection between
77 the reference and logical variables. Returned by the
78 :attr:`map` override so that
79 :meth:`~chebpy.classicfun.Classicfun.__call__` and
80 :meth:`~chebpy.classicfun.Classicfun.roots` route through the
81 non-affine map without further changes.
83 Examples:
84 ``sqrt`` has a branch-point singularity at the left endpoint, which a
85 plain polynomial approximation resolves only slowly. Declaring the
86 side clusters the points there instead:
88 >>> import numpy as np
89 >>> f = Singfun.initfun_adaptive(np.sqrt, [0.0, 1.0], sing="left")
90 >>> f.map.side
91 'left'
93 The integral of ``sqrt`` over [0, 1] is 2/3:
95 >>> bool(abs(f.sum() - 2.0 / 3.0) < 1e-10)
96 True
98 The map is non-affine, unlike the plain
99 :class:`~chebpy.utilities.Interval` a :class:`~chebpy.bndfun.Bndfun`
100 would carry:
102 >>> from chebpy.maps import SingleSlitMap
103 >>> isinstance(f.map, SingleSlitMap)
104 True
105 """
107 # Mixed-subclass binary ops (Singfun + Bndfun, etc.) reconstruct on the
108 # operand with the highest priority. Singfun outranks Bndfun/CompactFun
109 # because the singularity must be preserved in the result.
110 _singularity_priority: int = 10
112 def __init__(self, onefun: Any, interval: Any, map_: IntervalMap) -> None:
113 """Create a new :class:`Singfun`.
115 Args:
116 onefun: The :class:`~chebpy.onefun.Onefun` on ``[-1, 1]``
117 representing ``f(m(t))``.
118 interval: The logical support :class:`~chebpy.utilities.Interval`
119 ``Interval(a, b)`` (always finite for v1).
120 map_: The non-affine :class:`~chebpy.utilities.IntervalMap` between
121 ``[-1, 1]`` and ``[a, b]``.
122 """
123 super().__init__(onefun, interval)
124 self._map = map_
126 def _rebuild(self, onefun: Any) -> Singfun:
127 """Construct a new :class:`Singfun` preserving the map.
129 Used by every operation in :class:`~chebpy.classicfun.Classicfun`
130 that produces a new instance from a replacement ``onefun``
131 (``copy``, ``simplify``, the unary operators, the binary operators
132 between two same-map :class:`Singfun` instances).
133 """
134 return type(self)(onefun, self._interval, self._map)
136 def _can_share_onefun_with(self, other: Any) -> bool:
137 """Two :class:`Singfun` instances share a t-grid only when their maps match.
139 The ``Onefun`` coefficients of a :class:`Singfun` represent ``f(m(t))``
140 sampled at Chebyshev nodes in ``t``-space; if the maps differ, those
141 nodes encode different logical points and onefun-level arithmetic is
142 no longer correct. In that case the parent class falls back to
143 rebuilding the result adaptively on the dominant operand's map.
144 """
145 if not super()._can_share_onefun_with(other):
146 return False
147 return self._maps_equal(self._map, other._map)
149 def _rebuild_from_callable(self, f: Any) -> Singfun:
150 """Adaptively rebuild a :class:`Singfun` evaluating callable ``f`` on this map."""
151 m = self._map
152 if isinstance(m, SingleSlitMap):
153 return type(self).initfun_adaptive(f, self._interval, sing=m.side, params=m.params)
154 if isinstance(m, DoubleSlitMap):
155 return type(self).initfun_adaptive(f, self._interval, sing="both", params=m.params)
156 # Defensive: a Singfun map is always a SingleSlitMap or DoubleSlitMap.
157 msg = "Singfun._rebuild_from_callable: unknown map type" # pragma: no cover - defensive
158 raise NotImplementedError(msg) # pragma: no cover - defensive
160 @staticmethod
161 def _maps_equal(m1: IntervalMap, m2: IntervalMap) -> bool:
162 """Structural equality check for the maps used by :class:`Singfun`."""
163 if isinstance(m1, SingleSlitMap) and isinstance(m2, SingleSlitMap):
164 return m1.side == m2.side and m1.params == m2.params and m1.support == m2.support
165 if isinstance(m1, DoubleSlitMap) and isinstance(m2, DoubleSlitMap):
166 return m1.params == m2.params and m1.support == m2.support
167 return False
169 # ------------
170 # properties
171 # ------------
172 @property
173 def map(self) -> IntervalMap:
174 """Return the non-affine clustering map used by this :class:`Singfun`."""
175 return self._map
177 # --------------------------
178 # alternative constructors
179 # --------------------------
180 @classmethod
181 def initempty(cls) -> Singfun:
182 """Initialise an empty :class:`Singfun` with a default left-clustered map."""
183 iv = Interval(-1.0, 1.0)
184 m = SingleSlitMap(-1.0, 1.0, side="left")
185 onefun = techdict[prefs.tech].initempty(interval=iv)
186 return cls(onefun, iv, m)
188 @classmethod
189 def initconst(
190 cls,
191 c: Any,
192 interval: Any,
193 *,
194 sing: str = "left",
195 params: MapParams | None = None,
196 ) -> Singfun:
197 """Initialise a constant :class:`Singfun`.
199 Args:
200 c: The constant value.
201 interval: The bounded logical interval.
202 sing: Which endpoint(s) to cluster (``"left"`` / ``"right"`` / ``"both"``).
203 Default ``"left"``.
204 params: Slit-strip map parameters; if ``None``, :class:`MapParams`
205 defaults are used.
206 """
207 a, b = float(interval[0]), float(interval[1])
208 iv = Interval(a, b)
209 m = _build_map(a, b, sing, params if params is not None else MapParams())
210 onefun = techdict[prefs.tech].initconst(c, interval=iv)
211 return cls(onefun, iv, m)
213 @classmethod
214 def initidentity(
215 cls,
216 interval: Any,
217 *,
218 sing: str = "left",
219 params: MapParams | None = None,
220 ) -> Singfun:
221 """Initialise the identity ``f(x) = x`` as a :class:`Singfun`.
223 Note that ``f(x) = x`` is itself analytic on ``[a, b]`` and does not
224 require a clustering map; this constructor exists primarily for
225 symmetry with the :class:`Bndfun` API and for testing.
226 """
227 a, b = float(interval[0]), float(interval[1])
228 iv = Interval(a, b)
229 m = _build_map(a, b, sing, params if params is not None else MapParams())
230 # Sample x = m(t) at the t-Chebyshev nodes so the Onefun encodes f(m(t)) = m(t).
231 onefun = techdict[prefs.tech].initfun(lambda t: m.formap(t), interval=iv)
232 return cls(onefun, iv, m)
234 @classmethod
235 def initfun_adaptive(
236 cls,
237 f: Any,
238 interval: Any,
239 *,
240 sing: str = "left",
241 params: MapParams | None = None,
242 ) -> Singfun:
243 """Adaptive constructor for a :class:`Singfun`.
245 Builds the underlying :class:`~chebpy.chebtech.Chebtech` (or
246 :class:`~chebpy.trigtech.Trigtech`) by adaptively sampling ``f``
247 composed with the chosen clustering map.
249 Args:
250 f: Callable accepting a NumPy array of logical points and returning
251 an array of function values.
252 interval: The bounded logical interval ``(a, b)``.
253 sing: Which endpoint(s) of the interval carry a branch-type
254 singularity. One of ``"left"``, ``"right"``, ``"both"``.
255 params: Slit-strip map parameters; if ``None``, :class:`MapParams`
256 defaults are used.
258 Returns:
259 Singfun: The newly constructed :class:`Singfun`.
260 """
261 a, b = float(interval[0]), float(interval[1])
262 iv = Interval(a, b)
263 m = _build_map(a, b, sing, params if params is not None else MapParams())
264 onefun = techdict[prefs.tech].initfun(lambda t: f(m.formap(t)), interval=iv)
265 return cls(onefun, iv, m)
267 @classmethod
268 def initfun_fixedlen(
269 cls,
270 f: Any,
271 interval: Any,
272 n: int,
273 *,
274 sing: str = "left",
275 params: MapParams | None = None,
276 ) -> Singfun:
277 """Fixed-length constructor for a :class:`Singfun` (``n`` Chebyshev coefficients)."""
278 a, b = float(interval[0]), float(interval[1])
279 iv = Interval(a, b)
280 m = _build_map(a, b, sing, params if params is not None else MapParams())
281 onefun = techdict[prefs.tech].initfun(lambda t: f(m.formap(t)), n, interval=iv)
282 return cls(onefun, iv, m)
284 # ----------
285 # calculus
286 # ----------
287 def sum(self) -> Any:
288 r"""Definite integral of the function over ``[a, b]``.
290 Computed in the reference variable via the change-of-variables
292 .. math::
294 \int_a^b f(x)\,dx \;=\; \int_{-1}^{1} (f \circ m)(t)\, m'(t)\, dt.
296 The integrand ``onefun(t) * m'(t)`` is built adaptively as a
297 standard :class:`~chebpy.chebtech.Chebtech` on ``[-1, 1]``; ``m'(t)``
298 vanishes at the clustered endpoint(s), exactly absorbing the
299 integrable singularity of ``f``.
300 """
301 if self.onefun.isempty:
302 return 0.0
303 iv = Interval(-1.0, 1.0)
304 m = self._map
305 onefun = self.onefun
306 integrand = techdict[prefs.tech].initfun(lambda t: onefun(t) * m.drvmap(t), interval=iv)
307 return integrand.sum()
309 def cumsum(self) -> Singfun:
310 r"""Indefinite integral ``F(x) = \int_a^x f(s)\,ds`` as a :class:`Singfun`.
312 The chain rule gives ``F(m(t)) = \int_{-1}^t (f \circ m)(s)\, m'(s)\, ds``,
313 which is the cumulative sum of ``onefun(t) * m'(t)`` and is even more
314 regular than ``f`` itself. The result re-uses the same map.
315 """
316 if self.onefun.isempty:
317 return self._rebuild(self.onefun.copy())
318 iv = Interval(-1.0, 1.0)
319 m = self._map
320 onefun = self.onefun
321 integrand = techdict[prefs.tech].initfun(lambda t: onefun(t) * m.drvmap(t), interval=iv)
322 return self._rebuild(integrand.cumsum())
324 def diff(self) -> Singfun:
325 r"""Differentiation is not yet implemented for :class:`Singfun` (Phase 3 v1).
327 ``f'(x) = (f \\circ m)'(t) / m'(t)`` introduces a stronger
328 endpoint singularity (``1/m'`` blows up at the clustered endpoint),
329 which the current map cannot resolve. See plan 03 for the planned
330 treatment.
331 """
332 msg = (
333 "Singfun.diff is not implemented yet; differentiation introduces a "
334 "stronger endpoint singularity that the current map cannot resolve."
335 )
336 raise NotImplementedError(msg)
338 # -----------
339 # utilities
340 # -----------
341 def restrict(self, subinterval: Any) -> Classicfun:
342 """Restrict to a subinterval.
344 Behaviour depends on the relationship between ``subinterval`` and the
345 clustered endpoint(s):
347 * Trivial restriction (``subinterval == self.interval``) returns
348 ``self`` unchanged.
349 * A subinterval that **shares** the clustered endpoint (e.g. for
350 ``sing="left"`` a sub-range ``[a, c]``, or for ``sing="both"``
351 ``[a, c]`` / ``[c, b]``) is rebuilt as a :class:`Singfun` with a
352 rescaled map that retains the singular endpoint. Two-sided maps
353 are restricted to a one-sided map of the appropriate side.
354 * A purely interior subinterval (one that excludes the clustered
355 endpoint(s)) is returned as a :class:`~chebpy.bndfun.Bndfun`,
356 since the function is analytic there and the affine map suffices.
358 This is the closure-fallback described in plan 03 phase 4: the
359 result remains a usable :class:`~chebpy.classicfun.Classicfun` but
360 may change subclass.
361 """
362 if subinterval not in self.interval:
363 raise NotSubinterval(self.interval, subinterval)
364 a, b = float(self._interval[0]), float(self._interval[1])
365 sa, sb = float(subinterval[0]), float(subinterval[1])
366 if sa == a and sb == b:
367 return self
369 m = self._map
370 new_iv = Interval(sa, sb)
371 # Which clustered endpoints (if any) the subinterval still touches.
372 touches_left = sa == a
373 touches_right = sb == b
375 if isinstance(m, SingleSlitMap):
376 return self._restrict_single_slit(m, new_iv, touches_left=touches_left, touches_right=touches_right)
377 if isinstance(m, DoubleSlitMap):
378 return self._restrict_double_slit(m, new_iv, touches_left=touches_left, touches_right=touches_right)
380 # Unknown map type — conservative fallback.
381 return Bndfun.initfun_adaptive(self, new_iv) # pragma: no cover - defensive
383 def _restrict_single_slit(
384 self, m: SingleSlitMap, new_iv: Interval, *, touches_left: bool, touches_right: bool
385 ) -> Classicfun:
386 """Restrict a single-slit :class:`Singfun`.
388 Retains the :class:`Singfun` representation when the subinterval still
389 touches the clustered endpoint; otherwise the function is analytic on
390 the subinterval and drops to a :class:`~chebpy.bndfun.Bndfun`.
391 """
392 if (m.side == "left" and touches_left) or (m.side == "right" and touches_right):
393 return type(self).initfun_adaptive(self, new_iv, sing=m.side, params=m.params)
394 return Bndfun.initfun_adaptive(self, new_iv)
396 def _restrict_double_slit(
397 self, m: DoubleSlitMap, new_iv: Interval, *, touches_left: bool, touches_right: bool
398 ) -> Classicfun:
399 """Restrict a double-slit :class:`Singfun`.
401 A subinterval that still touches one clustered endpoint becomes a
402 one-sided :class:`Singfun`; a purely interior subinterval drops to a
403 :class:`~chebpy.bndfun.Bndfun`.
404 """
405 if touches_left and touches_right:
406 # Subinterval == self.interval handled by the caller; this branch is
407 # therefore unreachable in normal usage.
408 return self # pragma: no cover - defensive, see above
409 if touches_left:
410 return type(self).initfun_adaptive(self, new_iv, sing="left", params=m.params)
411 if touches_right:
412 return type(self).initfun_adaptive(self, new_iv, sing="right", params=m.params)
413 return Bndfun.initfun_adaptive(self, new_iv)
415 def translate(self, c: float) -> Singfun:
416 """Translate the function: ``g(x) = f(x - c)``.
418 The map is rebuilt for the shifted support; the underlying
419 ``onefun`` (which lives on ``[-1, 1]``) is unchanged.
420 """
421 a, b = float(self._interval[0]) + float(c), float(self._interval[1]) + float(c)
422 new_iv = Interval(a, b)
423 new_map = self._rebuild_map_for(a, b)
424 return type(self)(self.onefun, new_iv, new_map)
426 # -----------------
427 # internal helpers
428 # -----------------
429 def _rebuild_map_for(self, a: float, b: float) -> IntervalMap:
430 """Return a copy of ``self._map`` rescaled to a new logical interval."""
431 m = self._map
432 if isinstance(m, SingleSlitMap):
433 return SingleSlitMap(a, b, m.params, side=m.side)
434 if isinstance(m, DoubleSlitMap):
435 return DoubleSlitMap(a, b, m.params)
436 # Unknown map type — fall back to leaving the map unchanged; callers
437 # of translate that need a non-trivial rebuild should override.
438 return m # pragma: no cover - defensive