Coverage for src / chebpy / singfun.py: 90%
126 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-05 07:22 +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 .classicfun import Classicfun, techdict
33from .maps import DoubleSlitMap, MapParams, SingleSlitMap
34from .settings import _preferences as prefs
35from .utilities import Interval, IntervalMap
38def _build_map(a: float, b: float, sing: str, params: MapParams) -> IntervalMap:
39 """Construct the appropriate non-affine map for the requested singularity pattern.
41 Args:
42 a: Left endpoint of the logical interval.
43 b: Right endpoint of the logical interval.
44 sing: One of ``"left"``, ``"right"``, or ``"both"``.
45 params: A :class:`MapParams` instance carrying ``(L, alpha)``.
47 Returns:
48 IntervalMap: A :class:`SingleSlitMap` (for ``"left"`` / ``"right"``)
49 or :class:`DoubleSlitMap` (for ``"both"``).
51 Raises:
52 ValueError: If ``sing`` is not one of the recognised values.
53 """
54 if sing in ("left", "right"):
55 return SingleSlitMap(a, b, params, side=sing)
56 if sing == "both":
57 return DoubleSlitMap(a, b, params)
58 msg = f"sing must be 'left', 'right', or 'both'; got {sing!r}"
59 raise ValueError(msg)
62class Singfun(Classicfun):
63 """Functions with branch-type endpoint singularities on a bounded interval.
65 A :class:`Singfun` stores:
67 - ``self.onefun`` (inherited): a standard :class:`~chebpy.onefun.Onefun`
68 (typically a :class:`~chebpy.chebtech.Chebtech`) on ``[-1, 1]``
69 representing ``f(m(t))``, which is analytic by construction.
70 - ``self._interval`` (inherited): an :class:`~chebpy.utilities.Interval`
71 ``Interval(a, b)`` carrying the logical support endpoints.
72 - ``self._map``: a non-affine :class:`~chebpy.utilities.IntervalMap`
73 (a :class:`~chebpy.maps.SingleSlitMap` or
74 :class:`~chebpy.maps.DoubleSlitMap`) — the actual bijection between
75 the reference and logical variables. Returned by the
76 :attr:`map` override so that
77 :meth:`~chebpy.classicfun.Classicfun.__call__` and
78 :meth:`~chebpy.classicfun.Classicfun.roots` route through the
79 non-affine map without further changes.
80 """
82 # Mixed-subclass binary ops (Singfun + Bndfun, etc.) reconstruct on the
83 # operand with the highest priority. Singfun outranks Bndfun/CompactFun
84 # because the singularity must be preserved in the result.
85 _singularity_priority: int = 10
87 def __init__(self, onefun: Any, interval: Any, map_: IntervalMap) -> None:
88 """Create a new :class:`Singfun`.
90 Args:
91 onefun: The :class:`~chebpy.onefun.Onefun` on ``[-1, 1]``
92 representing ``f(m(t))``.
93 interval: The logical support :class:`~chebpy.utilities.Interval`
94 ``Interval(a, b)`` (always finite for v1).
95 map_: The non-affine :class:`~chebpy.utilities.IntervalMap` between
96 ``[-1, 1]`` and ``[a, b]``.
97 """
98 super().__init__(onefun, interval)
99 self._map = map_
101 def _rebuild(self, onefun: Any) -> Singfun:
102 """Construct a new :class:`Singfun` preserving the map.
104 Used by every operation in :class:`~chebpy.classicfun.Classicfun`
105 that produces a new instance from a replacement ``onefun``
106 (``copy``, ``simplify``, the unary operators, the binary operators
107 between two same-map :class:`Singfun` instances).
108 """
109 return type(self)(onefun, self._interval, self._map)
111 def _can_share_onefun_with(self, other: Any) -> bool:
112 """Two :class:`Singfun` instances share a t-grid only when their maps match.
114 The ``Onefun`` coefficients of a :class:`Singfun` represent ``f(m(t))``
115 sampled at Chebyshev nodes in ``t``-space; if the maps differ, those
116 nodes encode different logical points and onefun-level arithmetic is
117 no longer correct. In that case the parent class falls back to
118 rebuilding the result adaptively on the dominant operand's map.
119 """
120 if not super()._can_share_onefun_with(other):
121 return False
122 return self._maps_equal(self._map, other._map)
124 def _rebuild_from_callable(self, f: Any) -> Singfun:
125 """Adaptively rebuild a :class:`Singfun` evaluating callable ``f`` on this map."""
126 m = self._map
127 if isinstance(m, SingleSlitMap):
128 return type(self).initfun_adaptive(f, self._interval, sing=m.side, params=m.params)
129 if isinstance(m, DoubleSlitMap):
130 return type(self).initfun_adaptive(f, self._interval, sing="both", params=m.params)
131 msg = "Singfun._rebuild_from_callable: unknown map type"
132 raise NotImplementedError(msg) # pragma: no cover
134 @staticmethod
135 def _maps_equal(m1: IntervalMap, m2: IntervalMap) -> bool:
136 """Structural equality check for the maps used by :class:`Singfun`."""
137 if isinstance(m1, SingleSlitMap) and isinstance(m2, SingleSlitMap):
138 return m1.side == m2.side and m1.params == m2.params and m1.support == m2.support
139 if isinstance(m1, DoubleSlitMap) and isinstance(m2, DoubleSlitMap):
140 return m1.params == m2.params and m1.support == m2.support
141 return False
143 # ------------
144 # properties
145 # ------------
146 @property
147 def map(self) -> IntervalMap:
148 """Return the non-affine clustering map used by this :class:`Singfun`."""
149 return self._map
151 # --------------------------
152 # alternative constructors
153 # --------------------------
154 @classmethod
155 def initempty(cls) -> Singfun:
156 """Initialise an empty :class:`Singfun` with a default left-clustered map."""
157 iv = Interval(-1.0, 1.0)
158 m = SingleSlitMap(-1.0, 1.0, side="left")
159 onefun = techdict[prefs.tech].initempty(interval=iv)
160 return cls(onefun, iv, m)
162 @classmethod
163 def initconst(
164 cls,
165 c: Any,
166 interval: Any,
167 *,
168 sing: str = "left",
169 params: MapParams | None = None,
170 ) -> Singfun:
171 """Initialise a constant :class:`Singfun`.
173 Args:
174 c: The constant value.
175 interval: The bounded logical interval.
176 sing: Which endpoint(s) to cluster (``"left"`` / ``"right"`` / ``"both"``).
177 Default ``"left"``.
178 params: Slit-strip map parameters; if ``None``, :class:`MapParams`
179 defaults are used.
180 """
181 a, b = float(interval[0]), float(interval[1])
182 iv = Interval(a, b)
183 m = _build_map(a, b, sing, params if params is not None else MapParams())
184 onefun = techdict[prefs.tech].initconst(c, interval=iv)
185 return cls(onefun, iv, m)
187 @classmethod
188 def initidentity(
189 cls,
190 interval: Any,
191 *,
192 sing: str = "left",
193 params: MapParams | None = None,
194 ) -> Singfun:
195 """Initialise the identity ``f(x) = x`` as a :class:`Singfun`.
197 Note that ``f(x) = x`` is itself analytic on ``[a, b]`` and does not
198 require a clustering map; this constructor exists primarily for
199 symmetry with the :class:`Bndfun` API and for testing.
200 """
201 a, b = float(interval[0]), float(interval[1])
202 iv = Interval(a, b)
203 m = _build_map(a, b, sing, params if params is not None else MapParams())
204 # Sample x = m(t) at the t-Chebyshev nodes so the Onefun encodes f(m(t)) = m(t).
205 onefun = techdict[prefs.tech].initfun(lambda t: m.formap(t), interval=iv)
206 return cls(onefun, iv, m)
208 @classmethod
209 def initfun_adaptive(
210 cls,
211 f: Any,
212 interval: Any,
213 *,
214 sing: str = "left",
215 params: MapParams | None = None,
216 ) -> Singfun:
217 """Adaptive constructor for a :class:`Singfun`.
219 Builds the underlying :class:`~chebpy.chebtech.Chebtech` (or
220 :class:`~chebpy.trigtech.Trigtech`) by adaptively sampling ``f``
221 composed with the chosen clustering map.
223 Args:
224 f: Callable accepting a NumPy array of logical points and returning
225 an array of function values.
226 interval: The bounded logical interval ``(a, b)``.
227 sing: Which endpoint(s) of the interval carry a branch-type
228 singularity. One of ``"left"``, ``"right"``, ``"both"``.
229 params: Slit-strip map parameters; if ``None``, :class:`MapParams`
230 defaults are used.
232 Returns:
233 Singfun: The newly constructed :class:`Singfun`.
234 """
235 a, b = float(interval[0]), float(interval[1])
236 iv = Interval(a, b)
237 m = _build_map(a, b, sing, params if params is not None else MapParams())
238 onefun = techdict[prefs.tech].initfun(lambda t: f(m.formap(t)), interval=iv)
239 return cls(onefun, iv, m)
241 @classmethod
242 def initfun_fixedlen(
243 cls,
244 f: Any,
245 interval: Any,
246 n: int,
247 *,
248 sing: str = "left",
249 params: MapParams | None = None,
250 ) -> Singfun:
251 """Fixed-length constructor for a :class:`Singfun` (``n`` Chebyshev coefficients)."""
252 a, b = float(interval[0]), float(interval[1])
253 iv = Interval(a, b)
254 m = _build_map(a, b, sing, params if params is not None else MapParams())
255 onefun = techdict[prefs.tech].initfun(lambda t: f(m.formap(t)), n, interval=iv)
256 return cls(onefun, iv, m)
258 # ----------
259 # calculus
260 # ----------
261 def sum(self) -> Any:
262 r"""Definite integral of the function over ``[a, b]``.
264 Computed in the reference variable via the change-of-variables
266 .. math::
268 \int_a^b f(x)\,dx \;=\; \int_{-1}^{1} (f \circ m)(t)\, m'(t)\, dt.
270 The integrand ``onefun(t) * m'(t)`` is built adaptively as a
271 standard :class:`~chebpy.chebtech.Chebtech` on ``[-1, 1]``; ``m'(t)``
272 vanishes at the clustered endpoint(s), exactly absorbing the
273 integrable singularity of ``f``.
274 """
275 if self.onefun.isempty:
276 return 0.0
277 iv = Interval(-1.0, 1.0)
278 m = self._map
279 onefun = self.onefun
280 integrand = techdict[prefs.tech].initfun(lambda t: onefun(t) * m.drvmap(t), interval=iv)
281 return integrand.sum()
283 def cumsum(self) -> Singfun:
284 r"""Indefinite integral ``F(x) = \int_a^x f(s)\,ds`` as a :class:`Singfun`.
286 The chain rule gives ``F(m(t)) = \int_{-1}^t (f \circ m)(s)\, m'(s)\, ds``,
287 which is the cumulative sum of ``onefun(t) * m'(t)`` and is even more
288 regular than ``f`` itself. The result re-uses the same map.
289 """
290 if self.onefun.isempty:
291 return self._rebuild(self.onefun.copy())
292 iv = Interval(-1.0, 1.0)
293 m = self._map
294 onefun = self.onefun
295 integrand = techdict[prefs.tech].initfun(lambda t: onefun(t) * m.drvmap(t), interval=iv)
296 return self._rebuild(integrand.cumsum())
298 def diff(self) -> Singfun: # pragma: no cover - not yet implemented
299 r"""Differentiation is not yet implemented for :class:`Singfun` (Phase 3 v1).
301 ``f'(x) = (f \\circ m)'(t) / m'(t)`` introduces a stronger
302 endpoint singularity (``1/m'`` blows up at the clustered endpoint),
303 which the current map cannot resolve. See plan 03 for the planned
304 treatment.
305 """
306 msg = (
307 "Singfun.diff is not implemented yet; differentiation introduces a "
308 "stronger endpoint singularity that the current map cannot resolve."
309 )
310 raise NotImplementedError(msg)
312 # -----------
313 # utilities
314 # -----------
315 def restrict(self, subinterval: Any) -> Classicfun:
316 """Restrict to a subinterval.
318 Behaviour depends on the relationship between ``subinterval`` and the
319 clustered endpoint(s):
321 * Trivial restriction (``subinterval == self.interval``) returns
322 ``self`` unchanged.
323 * A subinterval that **shares** the clustered endpoint (e.g. for
324 ``sing="left"`` a sub-range ``[a, c]``, or for ``sing="both"``
325 ``[a, c]`` / ``[c, b]``) is rebuilt as a :class:`Singfun` with a
326 rescaled map that retains the singular endpoint. Two-sided maps
327 are restricted to a one-sided map of the appropriate side.
328 * A purely interior subinterval (one that excludes the clustered
329 endpoint(s)) is returned as a :class:`~chebpy.bndfun.Bndfun`,
330 since the function is analytic there and the affine map suffices.
332 This is the closure-fallback described in plan 03 phase 4: the
333 result remains a usable :class:`~chebpy.classicfun.Classicfun` but
334 may change subclass.
335 """
336 from .bndfun import Bndfun
338 if subinterval not in self.interval:
339 from .exceptions import NotSubinterval
341 raise NotSubinterval(self.interval, subinterval)
342 a, b = float(self._interval[0]), float(self._interval[1])
343 sa, sb = float(subinterval[0]), float(subinterval[1])
344 if sa == a and sb == b:
345 return self
347 m = self._map
348 new_iv = Interval(sa, sb)
350 # Decide which clustered endpoints (if any) the subinterval still touches.
351 touches_left = sa == a
352 touches_right = sb == b
354 if isinstance(m, SingleSlitMap):
355 if (m.side == "left" and touches_left) or (m.side == "right" and touches_right):
356 return type(self).initfun_adaptive(self, new_iv, sing=m.side, params=m.params)
357 # Interior or opposite-end restriction: drop to Bndfun.
358 return Bndfun.initfun_adaptive(self, new_iv)
360 if isinstance(m, DoubleSlitMap):
361 if touches_left and touches_right:
362 # Subinterval == self.interval handled above; this branch is
363 # therefore unreachable in normal usage.
364 return self # pragma: no cover
365 if touches_left:
366 return type(self).initfun_adaptive(self, new_iv, sing="left", params=m.params)
367 if touches_right:
368 return type(self).initfun_adaptive(self, new_iv, sing="right", params=m.params)
369 return Bndfun.initfun_adaptive(self, new_iv)
371 # Unknown map type — conservative fallback.
372 return Bndfun.initfun_adaptive(self, new_iv) # pragma: no cover
374 def translate(self, c: float) -> Singfun:
375 """Translate the function: ``g(x) = f(x - c)``.
377 The map is rebuilt for the shifted support; the underlying
378 ``onefun`` (which lives on ``[-1, 1]``) is unchanged.
379 """
380 a, b = float(self._interval[0]) + float(c), float(self._interval[1]) + float(c)
381 new_iv = Interval(a, b)
382 new_map = self._rebuild_map_for(a, b)
383 return type(self)(self.onefun, new_iv, new_map)
385 # -----------------
386 # internal helpers
387 # -----------------
388 def _rebuild_map_for(self, a: float, b: float) -> IntervalMap:
389 """Return a copy of ``self._map`` rescaled to a new logical interval."""
390 m = self._map
391 if isinstance(m, SingleSlitMap):
392 return SingleSlitMap(a, b, m.params, side=m.side)
393 if isinstance(m, DoubleSlitMap):
394 return DoubleSlitMap(a, b, m.params)
395 # Unknown map type — fall back to leaving the map unchanged; callers
396 # of translate that need a non-trivial rebuild should override.
397 return m # pragma: no cover