Coverage for src/chebpy/compactfun.py: 100%
287 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-01 13:43 +0000
1"""Implementation of functions on (semi-)infinite intervals via numerical-support truncation.
3This module provides the :class:`CompactFun` class, which sits next to
4:class:`~chebpy.bndfun.Bndfun` under :class:`~chebpy.classicfun.Classicfun`.
5It represents functions whose user-facing logical interval has one or both
6endpoints at ``±inf`` but whose **numerical support** — the set of points
7where the function differs from its asymptotic limit by more than a
8configured tolerance — is finite. Internally, a :class:`CompactFun` stores
9a standard :class:`~chebpy.onefun.Onefun` (Chebtech) on the discovered
10finite storage interval; outside that interval the function is reported as
11the corresponding asymptotic constant (``tail_left`` or ``tail_right``,
12default ``0``).
14This approach is a deliberate departure from MATLAB Chebfun's ``@unbndfun``
15(which uses a rational change of variables to map ``(-inf, inf)`` onto
16``[-1, 1]``). See ``docs/plans/02-compactfun-integration.md`` for the
17zero-tail design and ``docs/plans/02b-compactfun-tail-constants.md`` for
18the non-zero asymptote extension.
19"""
21from __future__ import annotations
23from typing import Any, cast
25import numpy as np
27from .bndfun import Bndfun
28from .classicfun import Classicfun, techdict
29from .exceptions import CompactFunConstructionError, DivergentIntegralError
30from .plotting import plotfun
31from .settings import _preferences as prefs
32from .utilities import Interval
35def _ensure_endpoints(interval: Any) -> tuple[float, float]:
36 """Return ``(a, b)`` floats from any 2-element interval-like object.
38 Accepts :class:`Interval`, ``numpy.ndarray``, list, or tuple. Both
39 endpoints may be ``±inf``.
40 """
41 a, b = interval[0], interval[1]
42 return float(a), float(b)
45def _discover_one_side(
46 f: Any, anchor: float, sign: int, tol: float, max_width: float, max_probes: int
47) -> tuple[float, float, float]:
48 """Discover the numerical-support boundary on one infinite side.
50 Probes ``f`` at ``anchor + sign * 2**k`` for ``k = 0, 1, 2, ...`` up to
51 the configured budget. Detects the asymptotic limit ``L`` of ``f`` on
52 this side (which may be zero or non-zero) and returns the smallest
53 finite boundary beyond which ``|f - L| < tol * scale``.
55 Args:
56 f: Callable being approximated.
57 anchor: Finite anchor point (the bounded endpoint of a semi-infinite
58 interval, or ``0.0`` for the doubly-infinite case).
59 sign: ``+1`` for the rightward (toward ``+inf``) side, ``-1`` for the
60 leftward side.
61 tol: Relative tolerance threshold for both convergence detection and
62 boundary placement.
63 max_width: Maximum permitted boundary distance from ``anchor``.
64 max_probes: Maximum number of geometric probes.
66 Returns:
67 Tuple ``(boundary, tail, vscale)`` where ``boundary`` is the finite
68 boundary, ``tail`` is the detected asymptotic constant (``0.0`` if
69 the function decays to zero), and ``vscale`` is the largest
70 absolute probed value on this side.
72 Raises:
73 CompactFunConstructionError: If ``f`` does not converge to a
74 constant within the probing budget or ``max_width``.
75 """
76 radii, values = _probe_side(f, anchor, sign, max_width, max_probes)
77 if not radii:
78 return anchor + sign * 1.0, 0.0, 0.0
80 vscale = max(abs(v) for v in values)
81 tail = _detect_tail(radii, values, anchor, sign, tol, vscale)
83 # Find the largest radius at which f is still "active" (above threshold
84 # relative to the tail).
85 threshold = tol * max(abs(tail), vscale, 1.0)
86 active_r = 0.0
87 for ri, vi in zip(radii, values, strict=False):
88 if abs(vi - tail) > threshold:
89 active_r = ri
91 boundary_r = max(2.0 * active_r, 1.0)
92 # Defensive: the last-3 convergence window forces active_r <= r_{N-3}, so
93 # boundary_r = 2*active_r stays below the largest probed radius (<= max_width).
94 if boundary_r > max_width: # pragma: no cover - defensive, see above
95 raise CompactFunConstructionError( # noqa: TRY003
96 f"Discovered numerical support exceeds max_width = {max_width:g}; "
97 f"heavy-tailed inputs are not supported in this release."
98 )
99 return anchor + sign * boundary_r, tail, vscale
102def _probe_side(f: Any, anchor: float, sign: int, max_width: float, max_probes: int) -> tuple[list[float], list[float]]:
103 """Geometrically probe ``f`` on one side, returning ``(radii, values)``.
105 Samples ``f`` at ``anchor + sign * 2**k`` for ``k = 0, 1, ...`` while the
106 radius stays within ``max_width`` and the probe budget is not exhausted.
108 Raises:
109 CompactFunConstructionError: If ``f`` returns a non-finite value.
110 """
111 radii: list[float] = []
112 values: list[float] = [] # signed values
113 r = 1.0
114 for _ in range(max_probes):
115 if r > max_width:
116 break
117 x = anchor + sign * r
118 try:
119 v = float(f(x))
120 except (FloatingPointError, OverflowError, ZeroDivisionError) as err:
121 raise CompactFunConstructionError( # noqa: TRY003
122 f"Could not evaluate f at probe x = {x:g} during numerical-support discovery"
123 ) from err
124 if not np.isfinite(v):
125 raise CompactFunConstructionError( # noqa: TRY003
126 f"f returned non-finite value {v} at probe x = {x:g}; CompactFun "
127 f"requires the function to be finite at all sampled points."
128 )
129 radii.append(r)
130 values.append(v)
131 r *= 2.0
132 return radii, values
135def _detect_tail(radii: list[float], values: list[float], anchor: float, sign: int, tol: float, vscale: float) -> float:
136 """Detect the asymptotic constant of ``f`` from its probed ``values``.
138 Requires at least three probes and that the last three agree to within
139 ``tol * max(vscale, 1)``; a tail below that threshold is reported as ``0``.
141 Raises:
142 CompactFunConstructionError: If there are too few probes or the last
143 few do not settle to a constant (heavy-tailed / oscillating input).
144 """
145 last_n = 3
146 if len(values) < last_n:
147 raise CompactFunConstructionError( # noqa: TRY003
148 f"Too few probes ({len(values)}) to determine the asymptotic "
149 f"behaviour of f near {'+' if sign > 0 else '-'}inf; "
150 f"increase numsupp_max_probes or numsupp_max_width."
151 )
153 # Convergence test: the last few signed probes must agree to tol*scale.
154 tail_window = values[-last_n:]
155 conv_threshold = tol * max(vscale, 1.0)
156 spread = max(tail_window) - min(tail_window)
157 if spread > conv_threshold:
158 # Function does not settle to a constant — heavy tail or oscillation.
159 raise CompactFunConstructionError( # noqa: TRY003
160 f"Function does not converge to a constant within "
161 f"{radii[-1]:g} of anchor {anchor:g} on the "
162 f"{'+' if sign > 0 else '-'}inf side (last {last_n} probes "
163 f"spread by {spread:g} > {conv_threshold:g}); heavy-tailed or "
164 f"oscillating inputs are not supported in this release."
165 )
166 tail = float(np.mean(tail_window))
167 if abs(tail) < conv_threshold:
168 tail = 0.0
169 return tail
172def _discover_numsupp(
173 f: Any, a: float, b: float, tol: float, max_width: float, max_probes: int
174) -> tuple[float, float, float, float]:
175 """Discover the storage interval and tail constants for ``f``.
177 Args:
178 f: Callable being approximated.
179 a: Left endpoint of the logical interval (may be ``-inf``).
180 b: Right endpoint of the logical interval (may be ``+inf``).
181 tol: Relative tolerance for support detection.
182 max_width: Maximum permitted storage interval width.
183 max_probes: Maximum probes per unbounded side.
185 Returns:
186 Tuple ``(a', b', tail_left, tail_right)`` where ``a' < b'`` are
187 finite floats and the tails are the detected asymptotic constants
188 (``0.0`` on any side whose logical endpoint is finite).
190 Raises:
191 CompactFunConstructionError: If support cannot be discovered.
192 """
193 left_inf = not np.isfinite(a)
194 right_inf = not np.isfinite(b)
196 if not (left_inf or right_inf):
197 return a, b, 0.0, 0.0
199 # Anchor: the finite endpoint of a semi-infinite interval, else 0.
200 if left_inf and right_inf:
201 anchor = 0.0
202 elif left_inf:
203 anchor = b
204 else:
205 anchor = a
207 if left_inf:
208 a_storage, tail_left, _ = _discover_one_side(f, anchor, -1, tol, max_width, max_probes)
209 else:
210 a_storage, tail_left = a, 0.0
212 if right_inf:
213 b_storage, tail_right, _ = _discover_one_side(f, anchor, +1, tol, max_width, max_probes)
214 else:
215 b_storage, tail_right = b, 0.0
217 if b_storage - a_storage > max_width:
218 raise CompactFunConstructionError( # noqa: TRY003
219 f"Discovered numerical support [{a_storage:g}, {b_storage:g}] exceeds "
220 f"max_width = {max_width:g}; heavy-tailed inputs are not supported "
221 f"in this release."
222 )
223 if b_storage <= a_storage: # pragma: no cover - defensive, see below
224 # Defensive: each discovered boundary is >= 1 from the anchor, so
225 # b_storage > a_storage always holds; kept as a safety net.
226 a_storage, b_storage = anchor - 1.0, anchor + 1.0
227 return a_storage, b_storage, tail_left, tail_right
230class CompactFun(Classicfun):
231 """Functions on (semi-)infinite intervals with finite numerical support.
233 A :class:`CompactFun` represents a function whose user-facing logical
234 interval has one or both endpoints at ``±inf`` but whose numerical
235 support — the set where the function differs from its asymptotic limit
236 by more than a configured tolerance — is finite. Internally it
237 inherits from :class:`Classicfun` and stores a standard
238 :class:`Onefun` on the discovered finite storage interval; outside that
239 interval the function is reported as the corresponding asymptotic
240 constant ``tail_left`` or ``tail_right`` (default ``0.0``).
242 Two intervals are tracked:
244 - ``self._interval`` (inherited): the finite storage interval where the
245 underlying ``Onefun`` lives.
246 - ``self._logical_interval``: the user-facing interval, which may have
247 ``±inf`` endpoints; returned by :attr:`support`.
249 Two scalar tail constants are tracked:
251 - ``tail_left``: the value reported for ``x < a_storage`` when the
252 logical-left endpoint is ``-inf``.
253 - ``tail_right``: the value reported for ``x > b_storage`` when the
254 logical-right endpoint is ``+inf``.
256 For finite logical intervals the storage and logical intervals coincide
257 and the tails are ignored, so a :class:`CompactFun` behaves identically
258 to :class:`~chebpy.bndfun.Bndfun`.
260 Attributes:
261 onefun: Inherited; the standard :class:`Onefun` on ``[-1, 1]``.
262 support: The logical interval (possibly with ``±inf`` endpoints).
263 numerical_support: The finite storage interval.
264 tail_left: Asymptotic value at ``-inf`` (``0.0`` if logical-left is finite).
265 tail_right: Asymptotic value at ``+inf`` (``0.0`` if logical-right is finite).
267 Examples:
268 A Gaussian on the whole real line. The logical support is infinite,
269 but the discovered numerical support is finite:
271 >>> import numpy as np
272 >>> f = CompactFun.initfun_adaptive(lambda x: np.exp(-(x**2)), (-np.inf, np.inf))
273 >>> f.support.tolist()
274 [-inf, inf]
275 >>> a, b = f.numerical_support
276 >>> bool(np.isfinite(a) and np.isfinite(b))
277 True
279 Integration over the infinite interval recovers ``sqrt(pi)``:
281 >>> bool(abs(f.sum() - np.sqrt(np.pi)) < 1e-12)
282 True
284 Outside the numerical support the function reports its asymptotic
285 limit rather than extrapolating the polynomial:
287 >>> f.tail_left, f.tail_right
288 (0.0, 0.0)
289 >>> float(f(1e6))
290 0.0
291 """
293 def __init__(
294 self,
295 onefun: Any,
296 interval: Any,
297 logical_interval: Any = None,
298 tail_left: float = 0.0,
299 tail_right: float = 0.0,
300 ) -> None:
301 """Create a new :class:`CompactFun` instance.
303 Args:
304 onefun: The :class:`Onefun` representing the function on ``[-1, 1]``.
305 interval: The finite storage :class:`Interval` (always finite).
306 logical_interval: The user-facing interval (possibly with ``±inf``
307 endpoints). Defaults to ``interval`` if omitted.
308 tail_left: Asymptotic value at ``-inf``. Default ``0.0``.
309 tail_right: Asymptotic value at ``+inf``. Default ``0.0``.
310 """
311 super().__init__(onefun, interval)
312 if logical_interval is None:
313 self._logical_interval = np.asarray(interval, dtype=float)
314 else:
315 self._logical_interval = np.asarray((float(logical_interval[0]), float(logical_interval[1])), dtype=float)
316 self._tail_left = float(tail_left)
317 self._tail_right = float(tail_right)
319 def _rebuild(self, onefun: Any, *, tail_left: float | None = None, tail_right: float | None = None) -> CompactFun:
320 """Construct a new :class:`CompactFun` preserving logical interval and tails.
322 Args:
323 onefun: Replacement :class:`Onefun` for the new instance.
324 tail_left: Optional override for the new instance's left tail.
325 Defaults to ``self.tail_left``.
326 tail_right: Optional override for the new instance's right tail.
327 Defaults to ``self.tail_right``.
328 """
329 new_tl = self._tail_left if tail_left is None else float(tail_left)
330 new_tr = self._tail_right if tail_right is None else float(tail_right)
331 return self.__class__(
332 onefun,
333 self._interval,
334 logical_interval=self._logical_interval,
335 tail_left=new_tl,
336 tail_right=new_tr,
337 )
339 # --------------------------
340 # alternative constructors
341 # --------------------------
342 @classmethod
343 def initempty(cls) -> CompactFun:
344 """Initialise an empty CompactFun on ``(-inf, +inf)``."""
345 storage = Interval(-1.0, 1.0)
346 onefun = techdict[prefs.tech].initempty(interval=storage)
347 return cls(onefun, storage, logical_interval=(-np.inf, np.inf))
349 @classmethod
350 def initconst(cls, c: Any, interval: Any) -> CompactFun:
351 """Initialise a constant function.
353 On an unbounded interval the constant ``c`` becomes the asymptotic
354 value on each unbounded side: ``tail_left = tail_right = c``. This
355 makes ``initconst`` total — every constant is representable on every
356 interval — but note that integrating a non-zero constant over an
357 unbounded logical interval will (correctly) raise
358 :class:`~chebpy.exceptions.DivergentIntegralError`.
359 """
360 a, b = _ensure_endpoints(interval)
361 c_val = float(c)
362 if not np.isfinite(a) and not np.isfinite(b):
363 storage = Interval(-1.0, 1.0)
364 elif not np.isfinite(a):
365 storage = Interval(b - 1.0, b)
366 elif not np.isfinite(b):
367 storage = Interval(a, a + 1.0)
368 else:
369 storage = Interval(a, b)
370 onefun = techdict[prefs.tech].initconst(c_val, interval=storage)
371 tail_left = c_val if not np.isfinite(a) else 0.0
372 tail_right = c_val if not np.isfinite(b) else 0.0
373 return cls(onefun, storage, logical_interval=(a, b), tail_left=tail_left, tail_right=tail_right)
375 @classmethod
376 def initidentity(cls, interval: Any) -> CompactFun:
377 """Initialise the identity function ``f(x) = x``.
379 The identity function is unbounded and so cannot be represented as a
380 :class:`CompactFun` on an unbounded interval. This method is provided
381 only for completeness and refuses any infinite endpoint.
382 """
383 a, b = _ensure_endpoints(interval)
384 if not (np.isfinite(a) and np.isfinite(b)):
385 raise CompactFunConstructionError( # noqa: TRY003
386 "The identity function f(x) = x cannot be represented as a CompactFun on an unbounded interval."
387 )
388 storage = Interval(a, b)
389 onefun = techdict[prefs.tech].initvalues(np.asarray(storage), interval=storage)
390 return cls(onefun, storage, logical_interval=(a, b))
392 @classmethod
393 def initfun_adaptive(cls, f: Any, interval: Any) -> CompactFun:
394 """Initialise from a callable using adaptive sampling.
396 Discovers the numerical support and asymptotic tail constants of
397 ``f`` on the (possibly unbounded) logical interval, then builds a
398 standard adaptive :class:`Onefun` on that finite storage interval.
400 Raises:
401 CompactFunConstructionError: If the numerical support cannot be
402 discovered within the configured tolerance and width budget,
403 or if ``f`` does not converge to a constant at ``±inf``.
404 """
405 a, b = _ensure_endpoints(interval)
406 a_s, b_s, tl, tr = _discover_numsupp(
407 f,
408 a,
409 b,
410 prefs.numsupp_tol,
411 prefs.numsupp_max_width,
412 prefs.numsupp_max_probes,
413 )
414 storage = Interval(a_s, b_s)
415 onefun = techdict[prefs.tech].initfun(lambda y: f(storage(y)), interval=storage)
416 return cls(onefun, storage, logical_interval=(a, b), tail_left=tl, tail_right=tr)
418 @classmethod
419 def initfun_fixedlen(cls, f: Any, interval: Any, n: int) -> CompactFun:
420 """Initialise from a callable using a fixed number of points.
422 Discovers numerical support and tails as in :meth:`initfun_adaptive`,
423 then builds a fixed-length :class:`Onefun` on the storage interval.
424 """
425 a, b = _ensure_endpoints(interval)
426 a_s, b_s, tl, tr = _discover_numsupp(
427 f,
428 a,
429 b,
430 prefs.numsupp_tol,
431 prefs.numsupp_max_width,
432 prefs.numsupp_max_probes,
433 )
434 storage = Interval(a_s, b_s)
435 onefun = techdict[prefs.tech].initfun(lambda y: f(storage(y)), n, interval=storage)
436 return cls(onefun, storage, logical_interval=(a, b), tail_left=tl, tail_right=tr)
438 # -------------------
439 # evaluation
440 # -------------------
441 def __call__(self, x: Any, how: str = "clenshaw") -> Any:
442 """Evaluate the function at ``x``.
444 Outside the storage interval, returns the corresponding tail constant
445 when the matching logical endpoint is ``±inf`` (default ``0.0``), or
446 ``0.0`` when the logical endpoint is finite.
447 """
448 scalar_input = np.isscalar(x) or np.ndim(x) == 0
449 x_arr = np.atleast_1d(np.asarray(x))
450 is_complex = bool(getattr(self.onefun, "iscomplex", False))
451 result = np.zeros(x_arr.shape, dtype=complex if is_complex else float)
452 a_s, b_s = self._interval
453 a_log, b_log = float(self._logical_interval[0]), float(self._logical_interval[1])
454 # Outside-storage values: tail constants where the logical edge is ±inf.
455 left_mask = x_arr < a_s
456 right_mask = x_arr > b_s
457 if not np.isfinite(a_log) and self._tail_left != 0.0:
458 result[left_mask] = self._tail_left
459 if not np.isfinite(b_log) and self._tail_right != 0.0:
460 result[right_mask] = self._tail_right
461 # Inside-storage values: standard onefun evaluation.
462 mask = (x_arr >= a_s) & (x_arr <= b_s)
463 if mask.any():
464 y = self._interval.invmap(x_arr[mask])
465 result[mask] = self.onefun(y, how)
466 if scalar_input:
467 return result.item()
468 return result
470 # ------------
471 # properties
472 # ------------
473 @property
474 def support(self) -> Any:
475 """Return the logical (user-facing) interval, possibly with ``±inf`` endpoints."""
476 return self._logical_interval
478 @property
479 def numerical_support(self) -> Any:
480 """Return the finite storage interval ``[a, b]`` discovered at construction."""
481 return np.asarray(self._interval)
483 @property
484 def tail_left(self) -> float:
485 """Asymptotic value of the function as ``x → -inf``.
487 Always ``0.0`` when the logical-left endpoint is finite.
488 """
489 return self._tail_left
491 @property
492 def tail_right(self) -> float:
493 """Asymptotic value of the function as ``x → +inf``.
495 Always ``0.0`` when the logical-right endpoint is finite.
496 """
497 return self._tail_right
499 @property
500 def endvalues(self) -> Any:
501 """Return values at the logical endpoints; tails at any ``±inf`` endpoint."""
502 a_log, b_log = float(self._logical_interval[0]), float(self._logical_interval[1])
503 yl = self._tail_left if not np.isfinite(a_log) else self.__call__(a_log)
504 yr = self._tail_right if not np.isfinite(b_log) else self.__call__(b_log)
505 return np.array([yl, yr])
507 def __repr__(self) -> str:
508 """Return a string representation showing the logical interval, size, and tails."""
509 a_log, b_log = self._logical_interval
510 if self._tail_left != 0.0 or self._tail_right != 0.0:
511 return (
512 f"{self.__class__.__name__}([{a_log}, {b_log}], {self.size}, "
513 f"tails=({self._tail_left}, {self._tail_right}))"
514 )
515 return f"{self.__class__.__name__}([{a_log}, {b_log}], {self.size})"
517 # ----------
518 # calculus
519 # ----------
520 def sum(self) -> Any:
521 """Compute the definite integral over the logical interval.
523 Raises:
524 DivergentIntegralError: If the logical interval is unbounded on
525 a side where the corresponding tail is non-zero (the integral
526 of a non-decaying function over a half-line diverges).
527 """
528 a_log, b_log = float(self._logical_interval[0]), float(self._logical_interval[1])
529 if (not np.isfinite(a_log)) and self._tail_left != 0.0:
530 raise DivergentIntegralError( # noqa: TRY003
531 f"Integrand has non-zero left asymptote tail_left={self._tail_left}; "
532 f"integral over (-inf, ...) diverges."
533 )
534 if (not np.isfinite(b_log)) and self._tail_right != 0.0:
535 raise DivergentIntegralError( # noqa: TRY003
536 f"Integrand has non-zero right asymptote tail_right={self._tail_right}; "
537 f"integral over (..., +inf) diverges."
538 )
539 return super().sum()
541 def cumsum(self) -> CompactFun:
542 """Compute the indefinite integral.
544 For a :class:`CompactFun` with zero asymptote on the unbounded
545 left/right side, the antiderivative is well-defined; it is itself a
546 :class:`CompactFun` whose right-tail equals ``∫f`` and whose
547 left-tail is ``0`` (anchored so ``F(-inf) = 0``).
549 Raises:
550 DivergentIntegralError: If the logical interval is unbounded on
551 a side where the corresponding tail is non-zero, in which
552 case the antiderivative diverges.
553 """
554 a_log, b_log = float(self._logical_interval[0]), float(self._logical_interval[1])
555 if (not np.isfinite(a_log)) and self._tail_left != 0.0:
556 raise DivergentIntegralError( # noqa: TRY003
557 f"Antiderivative diverges at -inf because tail_left={self._tail_left} != 0."
558 )
559 if (not np.isfinite(b_log)) and self._tail_right != 0.0:
560 raise DivergentIntegralError( # noqa: TRY003
561 f"Antiderivative diverges at +inf because tail_right={self._tail_right} != 0."
562 )
563 # Standard cumsum on the storage interval anchors F(a_storage) = 0.
564 # When logical-left is -inf with tail_left=0, this approximates
565 # F(-inf) = 0 (since f is below tolerance below a_storage).
566 inner = super().cumsum()
567 # The right-tail of F is the total integral.
568 total = float(super().sum())
569 # The left-tail is 0 when logical-left is -inf (anchor at -inf).
570 new_tail_left = 0.0
571 new_tail_right = total
572 return self.__class__(
573 inner.onefun,
574 inner._interval,
575 logical_interval=self._logical_interval,
576 tail_left=new_tail_left,
577 tail_right=new_tail_right,
578 )
580 def diff(self) -> CompactFun:
581 """Compute the derivative.
583 The derivative of a function with constant asymptotic limits has
584 zero asymptotes, so the result has ``tail_left = tail_right = 0``.
585 """
586 result = cast(CompactFun, super().diff())
587 result._tail_left = 0.0
588 result._tail_right = 0.0
589 return result
591 # -------------
592 # rootfinding
593 # -------------
594 def roots(self) -> Any:
595 """Find the roots, filtering out spurious roots in numerical-noise regions.
597 The underlying polynomial approximation can produce many spurious
598 roots in regions where the function has decayed to numerical noise
599 (typically near the boundary of the storage interval). We keep a
600 candidate root ``r`` only if both:
602 - ``f(r - δ)`` and ``f(r + δ)`` have opposite signs (the function
603 actually crosses zero), **and**
604 - ``max(|f(r - δ)|, |f(r + δ)|)`` exceeds ``numsupp_tol * vscale``
605 (the values are above numerical noise).
607 Here ``delta = 1e-3 * storage_width``. This heuristic does not preserve
608 double roots; that is a documented limitation since double roots are
609 uncommon in the decay-to-zero functions that :class:`CompactFun` is
610 designed for.
611 """
612 raw = super().roots()
613 if raw.size == 0:
614 return raw
615 a_s, b_s = float(self._interval[0]), float(self._interval[1])
616 vals = np.abs(np.atleast_1d(self.onefun.values()))
617 vscale = float(vals.max()) if vals.size else 1.0
618 threshold = prefs.numsupp_tol * max(vscale, 1.0)
619 delta = 1e-3 * (b_s - a_s)
620 left = np.clip(raw - delta, a_s, b_s)
621 right = np.clip(raw + delta, a_s, b_s)
622 f_left = np.atleast_1d(self.__call__(left))
623 f_right = np.atleast_1d(self.__call__(right))
624 sign_flip = np.sign(f_left) != np.sign(f_right)
625 above_noise = np.maximum(np.abs(f_left), np.abs(f_right)) > threshold
626 keep = sign_flip & above_noise
627 return np.sort(np.unique(raw[keep]))
629 # -----------
630 # utilities
631 # -----------
632 def restrict(self, subinterval: Any) -> Any:
633 """Restrict to a finite subinterval, returning a :class:`Bndfun`."""
634 sub_a, sub_b = _ensure_endpoints(subinterval)
635 if not (np.isfinite(sub_a) and np.isfinite(sub_b)):
636 raise NotImplementedError(
637 "CompactFun.restrict() requires a finite subinterval; "
638 "restriction to unbounded subintervals is not supported."
639 )
640 return Bndfun.initfun_adaptive(self, Interval(sub_a, sub_b))
642 def translate(self, c: float) -> CompactFun:
643 """Translate by ``c`` along the real line, preserving both intervals and tails."""
644 new_storage = Interval(float(self._interval[0]) + c, float(self._interval[1]) + c)
645 a_log, b_log = float(self._logical_interval[0]), float(self._logical_interval[1])
646 new_logical = (a_log + c, b_log + c)
647 return self.__class__(
648 self.onefun,
649 new_storage,
650 logical_interval=new_logical,
651 tail_left=self._tail_left,
652 tail_right=self._tail_right,
653 )
655 # ------------
656 # arithmetic
657 # ------------
658 def __neg__(self) -> CompactFun:
659 """Return ``-f``; negates both tail constants."""
660 result = cast(CompactFun, super().__neg__())
661 result._tail_left = -self._tail_left
662 result._tail_right = -self._tail_right
663 return result
665 def __add__(self, other: Any) -> Any:
666 """Pointwise addition; combines tail constants additively."""
667 result = super().__add__(other)
668 if isinstance(result, CompactFun):
669 other_tl, other_tr = self._other_tails(other)
670 result._tail_left = self._tail_left + other_tl
671 result._tail_right = self._tail_right + other_tr
672 return result
674 def __radd__(self, other: Any) -> Any:
675 """Right-hand addition for scalar + CompactFun."""
676 result = super().__radd__(other)
677 if isinstance(result, CompactFun):
678 other_tl, other_tr = self._other_tails(other)
679 result._tail_left = self._tail_left + other_tl
680 result._tail_right = self._tail_right + other_tr
681 return result
683 def __sub__(self, other: Any) -> Any:
684 """Pointwise subtraction; combines tail constants additively."""
685 result = super().__sub__(other)
686 if isinstance(result, CompactFun):
687 other_tl, other_tr = self._other_tails(other)
688 result._tail_left = self._tail_left - other_tl
689 result._tail_right = self._tail_right - other_tr
690 return result
692 def __rsub__(self, other: Any) -> Any:
693 """Right-hand subtraction for scalar - CompactFun."""
694 result = super().__rsub__(other)
695 if isinstance(result, CompactFun):
696 other_tl, other_tr = self._other_tails(other)
697 result._tail_left = other_tl - self._tail_left
698 result._tail_right = other_tr - self._tail_right
699 return result
701 def __mul__(self, other: Any) -> Any:
702 """Pointwise multiplication; combines tail constants multiplicatively."""
703 result = super().__mul__(other)
704 if isinstance(result, CompactFun):
705 other_tl, other_tr = self._other_tails(other)
706 result._tail_left = self._tail_left * other_tl
707 result._tail_right = self._tail_right * other_tr
708 return result
710 def __rmul__(self, other: Any) -> Any:
711 """Right-hand multiplication for scalar * CompactFun."""
712 result = super().__rmul__(other)
713 if isinstance(result, CompactFun):
714 other_tl, other_tr = self._other_tails(other)
715 result._tail_left = self._tail_left * other_tl
716 result._tail_right = self._tail_right * other_tr
717 return result
719 def _other_tails(self, other: Any) -> tuple[float, float]:
720 """Extract ``(tail_left, tail_right)`` from a binary-op operand.
722 For a :class:`CompactFun` operand, returns its tail attributes; for
723 a scalar, returns ``(scalar, scalar)``.
724 """
725 if isinstance(other, CompactFun):
726 return other._tail_left, other._tail_right
727 if np.isscalar(other):
728 v = float(cast(Any, other))
729 return v, v
730 # Anything else (e.g. a different Classicfun subclass) is treated as
731 # zero-tailed; tail propagation may be inexact in that case.
732 return 0.0, 0.0
734 # ----------
735 # plotting
736 # ----------
737 @property
738 def plot_support(self) -> tuple[float, float]:
739 """Return a finite ``[a, b]`` plotting window.
741 Replaces any ``±inf`` logical endpoint with the corresponding
742 numerical-support endpoint padded by 10% of the storage width
743 (minimum padding of 1.0) so the decay-to-zero region is visible.
744 """
745 a_s, b_s = float(self._interval[0]), float(self._interval[1])
746 a_log, b_log = float(self._logical_interval[0]), float(self._logical_interval[1])
747 pad = max(0.1 * (b_s - a_s), 1.0)
748 a = a_log if np.isfinite(a_log) else a_s - pad
749 b = b_log if np.isfinite(b_log) else b_s + pad
750 return (a, b)
752 def plot(self, ax: Any = None, **kwds: Any) -> Any:
753 """Plot the function over a finite window derived from its numerical support.
755 For doubly- or singly-infinite logical intervals, the plotting window
756 defaults to the numerical-support interval padded by 10% on each
757 unbounded side. Pass an explicit ``support=(a, b)`` keyword to override.
758 """
759 support = kwds.pop("support", self.plot_support)
760 return plotfun(self, support, ax=ax, **kwds)