solprop.translate.chebyshev#

Adaptive piecewise Chebyshev representation of functions of z.

The solver rests on the cumulative linear operator Phi(omega, z) = int D dz: because D is diagonal in omega it commutes with itself, so exp(Phi(z2) - Phi(z1)) is the exact linear flow over any step. This module builds the representation that carries it — a piecewise Chebyshev polynomial in z with array-valued coefficients, built by Chebfun-style adaptive splitting, integrated exactly in the Chebyshev basis and evaluated by Clenshaw recursion.

The design, the evidence behind each choice and the traps that are silent are in docs/dev/stepper.md section 4; the manual page is docs/guide/numerics/linear-operator.md.

Module Attributes

`DEFAULT_TOL`	Default tolerance: the bound on the uniform error relative to the value scale.
`DEFAULT_MAX_DEGREE`	Default highest Chebyshev degree tried on one interval before it is bisected.
`DEFAULT_MAX_INTERVALS`	Default interval budget; exceeding it sets the fallback flag rather than raising.
`MIN_WIDTH_FRACTION`	Default narrowest interval the split may produce, as a fraction of the total span.
`INITIAL_DEGREE`	Degree the adaptive search starts from, here and in the parameter expansion.
`ENVELOPE_FLOOR`	Comparison-scale floor of the acceptance test, as a fraction of the global envelope.
`SAMPLING_NOISE_FACTOR`	Margin of the local sampling-noise floor of the acceptance test (SEP 14 phase 5).

Functions

`build`(f, z0, z1, *[, tol, ...])	Build an adaptive piecewise Chebyshev representation of `f` on `[z0, z1]`.
`build_from_paths`(paths, f_of_params, z0, z1, *)	Resolve the scalar paths first, then the array-valued function on their union.
`clenshaw`(x, coefficients)	Evaluate one interval's Chebyshev series at one point, as NumPy's `chebval`.
`interleave_axis`(coarse, fine, axis)	Merge a nested grid's samples with the ones a doubling added, along one axis.
`lobatto_coefficients`(values)	Transform samples at the Chebyshev-Lobatto points into Chebyshev coefficients.
`lobatto_points`(degree, lo, hi)	Return the Chebyshev-Lobatto points of `degree` mapped onto `[lo, hi]`.
`refinement_points`(degree, lo, hi)	Return the points a degree doubling adds to `lobatto_points()`.
`resolve_axis`(values, axis, *, tol)	Transform one axis to Chebyshev coefficients and judge its tail.
`trim_axis`(coefficients, axis)	Drop the trailing coefficients of one axis that are transform noise.

Classes

PiecewiseChebyshev(coefficients, ...[, ...])

A piecewise Chebyshev polynomial in z with array-valued coefficients.

solprop.translate.chebyshev.DEFAULT_MAX_DEGREE: Final = 64#: Default highest Chebyshev degree tried on one interval before it is bisected.

solprop.translate.chebyshev.DEFAULT_MAX_INTERVALS: Final = 512#: Default interval budget; exceeding it sets the fallback flag rather than raising.

solprop.translate.chebyshev.DEFAULT_TOL: Final = 1e-12#: Default tolerance: the bound on the uniform error relative to the value scale.

solprop.translate.chebyshev.ENVELOPE_FLOOR: Final = 1e-10#

Comparison-scale floor of the acceptance test, as a fraction of the global envelope.

An element whose own values are negligible beside the rest must not drive the refinement, so its comparison scale is floored here. Exposed because the round-off floor solprop.translate.linop measures has to be measured in this same norm to be comparable with the tail it is meant to bound.

solprop.translate.chebyshev.INITIAL_DEGREE: Final = 8#: Degree the adaptive search starts from, here and in the parameter expansion.

solprop.translate.chebyshev.MIN_WIDTH_FRACTION: Final = 1e-06#: Default narrowest interval the split may produce, as a fraction of the total span.

class solprop.translate.chebyshev.PiecewiseChebyshev(coefficients, breakpoints, degrees, hard_kinks, budget_exceeded=False)[source]#

Bases: object

A piecewise Chebyshev polynomial in z with array-valued coefficients.

On interval i the value is sum_j coefficients[i, j] T_j(x) with x = (2 z - lo - hi) / (hi - lo) mapping [lo, hi] onto [-1, 1]. Intervals that need a lower degree carry zero coefficients in the high-order slots, which Clenshaw handles exactly and at no accuracy cost, so one vectorised pass evaluates every interval.

coefficients is listed first because jaxtyping resolves the symbolic axis n_intervals+1 of breakpoints only once n_intervals is bound by an earlier argument. Shape agreement between the fields is enforced by those annotations under the test session’s runtime checker; construction validates what shapes cannot express.

Attributes:

coefficientsnumpy.ndarray: Chebyshev coefficients, shape (n_intervals, n_deg, *value_shape), float64 or complex128. Units are those of the represented function.
breakpointsnumpy.ndarray: Interval edges in z, shape (n_intervals + 1,), m, strictly increasing.
degreesnumpy.ndarray: Degree actually used on each interval, shape (n_intervals,). Diagnostic: the zero padding makes it unnecessary for evaluation.
hard_kinksnumpy.ndarray: Positions in z, m, where the represented function itself jumps or kinks — stage boundaries, taper corners. Every entry is an interior breakpoint. These are tracked apart from the far more numerous edges the adaptive split introduces to resolve smooth variation, because steps must land on them and must not be capped by the others.
budget_exceededbool, optional: True when the split ran out of its interval budget and accepted unresolved intervals. The caller falls back to per-substep Gauss-Legendre on it.

Raises:

TranslationError: If the interval count is zero, the breakpoints are not finite and strictly increasing, the coefficient dtype is neither float64 nor complex128, or a hard kink is not an interior breakpoint.

Examples

>>> import numpy as np
>>> from solprop.translate.chebyshev import build
>>> phi = build(lambda z: 2.0 * z, 0.0, 1.0).integral()
>>> float(phi(0.5))
0.25

at_positions(z, *, helpers=Helpers(name='numba'))[source]#

Evaluate at several positions, one Clenshaw pass for all of them.

Each row is what calling the table returns at that position, bit for bit: the same interval choice (an interior breakpoint belongs to the interval starting there), the same map onto [-1, 1] and the same recursion, so a caller may batch the positions it knows without changing a value. Rows are kept as the last positions’ values, as a call of the table keeps them, so the next step’s anchor finds its weights already formed.

The pass is the caller’s backend’s (helpers, SEP 29 section 6.6): a stage’s rows are evaluated by the backend that runs the stage. The kept rows are shared whichever implementation filled them, because every implementation performs the operations of a one-point call in its order under strict arithmetic, so a kept row is the value any of them would form.

Parameters:

znumpy.ndarray: Positions, m, float64, shape (n_points,), each within the represented span (a few ulp outside is tolerated).
helperssolprop._helpers.Helpers, optional: The evaluating backend’s implementation; the shared Numba one by default, for callers outside a backend.

Returns:

numpy.ndarray: The values, shape (n_points, *value_shape), read-only.

Raises:

TranslationError: If a position lies outside the represented span.

Examples

>>> import numpy as np
>>> from solprop.translate.chebyshev import build
>>> f = build(lambda z: np.cos(z), 0.0, 1.0)
>>> rows = f.at_positions(np.array([0.3, 0.7]))
>>> bool(rows[0] == f(0.3)) and bool(rows[1] == f(0.7))
True

breakpoints: Float[ndarray, 'n_intervals+1']#

budget_exceeded: bool = False#

coefficients: Inexact[ndarray, 'n_intervals n_deg *value']#

degrees: Int[ndarray, 'n_intervals']#

derivative()[source]#

Return the derivative with respect to z.

Returns:

PiecewiseChebyshev: The derivative, with the same breakpoints and hard kinks and one fewer coefficient per interval (a single zero coefficient at the floor).

Examples

>>> from solprop.translate.chebyshev import build
>>> f = build(lambda z: z**2, 0.0, 1.0).derivative()
>>> bool(abs(f(0.5) - 1.0) < 1e-14)
True

hard_kinks: Float[ndarray, 'n_kinks']#

integral()[source]#

Return the antiderivative, continuous across intervals and zero at the start.

Integration is a coefficient recurrence in the Chebyshev basis, not a quadrature, so it is exact for the polynomial it is given: the only approximation left in the cumulative operator is the interpolation of the integrand. The constant of integration lives entirely in the T_0 coefficient and is fixed from the left, interval by interval, so the result is continuous by construction.

Returns:

PiecewiseChebyshev: The antiderivative with respect to z, vanishing at the first breakpoint, with the same breakpoints and hard kinks and one more coefficient per interval.

Examples

>>> from solprop.translate.chebyshev import build
>>> phi = build(lambda z: 3.0 * z**2, 0.0, 2.0).integral()
>>> bool(abs(phi(2.0) - 8.0) < 1e-14)
True

interval_containing(z)[source]#

Return the index of the interval that contains z.

A position sitting exactly on an interior breakpoint belongs to the interval starting there, and the right end of the span resolves to the last interval.

Parameters:

zfloat: Position, m. Positions up to 1e-10 of the span outside are tolerated, because step schedules land on the ends and rounding puts them a few ulp out.

Returns:

int: Interval index in [0, n_intervals - 1].

Raises:

TranslationError: If z lies further outside the span than that tolerance.

Examples

>>> from solprop.translate.chebyshev import build
>>> f = build(lambda z: z, 0.0, 1.0, mandatory_breakpoints=(0.5,))
>>> f.interval_containing(0.5)
1

interval_width(z)[source]#

Return the width of the interval containing z.

This is the local step cap of docs/dev/stepper.md section 3. The narrowest interval anywhere is not a cap: the split resolves a square-root entrance with intervals four orders of magnitude below the median, and capping globally on those would cost four orders of magnitude in step count.

Parameters:

zfloat: Position, m.

Returns:

float: Width of the containing interval, m.

Examples

>>> from solprop.translate.chebyshev import build
>>> f = build(lambda z: z, 0.0, 1.0, mandatory_breakpoints=(0.25,))
>>> f.interval_width(0.1)
0.25

next_hard_kink(z)[source]#

Return the first hard kink beyond z, or None if there is none.

Strictly beyond, with the same edge tolerance as interval_containing(), so a step starting on a kink clears it instead of being capped to zero length.

Parameters:

zfloat: Position, m.

Returns:

float or None: Position of the next hard kink, m, or None when none remains.

Examples

>>> from solprop.translate.chebyshev import build
>>> f = build(lambda z: z, 0.0, 1.0, mandatory_breakpoints=(0.4,))
>>> f.next_hard_kink(0.0), f.next_hard_kink(0.4)
(0.4, None)

property value_shape: tuple[int, ...][source]#

Return the shape of one value.

Returns:

tuple of int: Shape of one value; () for a scalar-valued function.

property z_span: tuple[float, float][source]#

Return the ends of the represented span.

Returns:

tuple of float: First and last breakpoint, m.

solprop.translate.chebyshev.SAMPLING_NOISE_FACTOR: Final = 4.0#

Margin of the local sampling-noise floor of the acceptance test (SEP 14 phase 5).

An interval is also accepted when twice its coefficient tail is within this factor of the relative change its interpolant makes across the rounding of its own node positions, eps |z| max|f'| / |f|: below that the sampled function carries no information about its own variation. Near a choked gas end, where the pressure goes as the square root of the distance to the exit, that floor rises above the requested tolerance within half a millimetre of a 4.2 m exit, and without it the split tiled that region at its width floor.

solprop.translate.chebyshev.build(f, z0, z1, *, tol=1e-12, mandatory_breakpoints=(), n_max=64, min_width=None, max_intervals=512)[source]#

Build an adaptive piecewise Chebyshev representation of f on [z0, z1].

The split is Chebfun-style: sample at Chebyshev-Lobatto points, accept when the coefficient tail is negligible, otherwise double the degree up to n_max, otherwise bisect. Declared breakpoints are mandatory edges and become the hard kinks, so a jump that is known in advance is never crossed; a jump that is not declared is located by the bisection down to min_width.

Parameters:

fcallable: Sampled vectorised: f(z) takes positions of shape (n_z,) in m and returns an array of shape (n_z, *value_shape), float64 or complex128. Sampling is vectorised because the operator is evaluated over a whole omega grid at once.
z0, z1float: Ends of the domain, m, with z1 > z0.
tolfloat, optional: Target bound on the uniform error relative to the value scale, > 0. An interval is accepted when twice the relative coefficient tail is at most tol; the factor two is the Chebyshev aliasing bound, so tol bounds the error itself. Use 1e-12 when the frame quantities are analytic, 1e-8 when differenced.
mandatory_breakpointsarray_like, optional: Positions in z, m, inside [z0, z1] that must be interval edges: stage boundaries, taper corners. The interior ones become the PiecewiseChebyshev hard kinks.
n_maxint, optional: Highest degree tried on one interval before bisecting, >= 3.
min_widthfloat or None, optional: Narrowest interval the split may produce, m. None means MIN_WIDTH_FRACTION * (z1 - z0). It is absolute and computed once from the whole span: derived per call it lets an already tiny segment recurse several further decades and exhaust the budget without getting closer to resolving anything.
max_intervalsint, optional: Interval budget, >= 1. Exceeding it stops the split, accepts what is left and sets the PiecewiseChebyshev fallback flag rather than raising. Mandatory breakpoints are honoured whatever the budget, so a domain declaring more of them than the budget allows still gets one interval per segment.

Returns:

PiecewiseChebyshev: The representation of f. Call PiecewiseChebyshev.integral() for its cumulative integral.

Raises:

TranslationError: If the domain, the tolerance, the degree or the budget is invalid, if a mandatory breakpoint lies outside the domain, or if f returns the wrong leading length, an empty value shape, or a non-finite sample.

Notes

What the build reproduces exactly, to round-off, whatever tol is:

any polynomial in z of degree at most n_max - 3 on each declared segment. The interpolation at the Chebyshev-Lobatto points is exact for it and the tail test then sees three vanishing coefficients. Degrees n_max - 2 to n_max are still represented exactly but are not recognised as resolved, so the split bisects them until it reaches min_width;
a constant or piecewise-constant function, as one interval of degree 0 per declared segment, whose integral is one interval of degree 1;
the round trip build(f).integral().derivative(), which returns the original coefficients up to the round-off of the two recurrences.

The other round trip, the one the interaction picture needs – exp(-dPhi) as the exact inverse of exp(+dPhi) – is the caller’s: negate one evaluated exponent rather than evaluating the table twice.

Examples

>>> import numpy as np
>>> from solprop.translate.chebyshev import build
>>> f = build(lambda z: np.exp(-z), 0.0, 1.0)
>>> int(f.degrees.size)
1
>>> bool(abs(f(0.25) - np.exp(-0.25)) < 1e-14)
True

solprop.translate.chebyshev.build_from_paths(paths, f_of_params, z0, z1, *, tol=1e-12, mandatory_breakpoints=(), n_max=64, min_width=None, max_intervals=512)[source]#

Resolve the scalar paths first, then the array-valued function on their union.

The scalar parameter paths — density, core radius, any other profile of z — are cheap to evaluate and carry the structure that makes the operator hard: they are resolved first, within each declared segment separately so that a declared kink is never crossed by a path fit, and the union of all the edges they produce becomes the mandatory split of the expensive array-valued function.

Parameters:

pathssequence of callable: Scalar paths, each mapping positions of shape (n_z,) in m to values of shape (n_z,).
f_of_paramscallable: Called as f_of_params(z, *values) with the positions and one array of shape (n_z,) per path, returning an array of shape (n_z, *value_shape). The positions come first so that an operator which also depends on z directly – a co-moving frame quantity, say – needs no path of its own, and so that a build with no paths at all is exactly build().
z0, z1float: Ends of the domain, m, with z1 > z0.
tolfloat, optional: Target bound on the uniform error relative to the value scale; see build().
mandatory_breakpointsarray_like, optional: Positions in z, m, that must be interval edges; the interior ones become the hard kinks.
n_maxint, optional: Highest degree tried on one interval before bisecting, >= 3.
min_widthfloat or None, optional: Narrowest interval the split may produce, m; see build().
max_intervalsint, optional: Interval budget for each stage of the build, >= 1.

Returns:

PiecewiseChebyshev: The representation of f_of_params evaluated along the paths.

Raises:

TranslationError: For the same reasons as build().

Examples

>>> import numpy as np
>>> from solprop.translate.chebyshev import build_from_paths
>>> density = lambda z: 1.0 + z
>>> f = build_from_paths([density], lambda z, rho: rho**2, 0.0, 1.0)
>>> bool(abs(f(0.5) - 2.25) < 1e-14)
True

solprop.translate.chebyshev.clenshaw(x, coefficients)[source]#

Evaluate one interval’s Chebyshev series at one point, as NumPy’s chebval.

The same values bit for bit, compiled: a degree-64 series of 18 values costs under a microsecond against NumPy’s 36 µs. Until numba is loaded NumPy evaluates instead: a run on the native backend needs numba for nothing else, and importing it costs more than a set-up’s thousand NumPy evaluations.

Parameters:

xfloat: Point in the interval’s own coordinate, in [-1, 1].
coefficientsnumpy.ndarray: The series’ coefficients, lowest degree first, shape (n_terms, *value), float64 or complex128.

Returns:

numpy.ndarray: The value, shape value, of the coefficients’ dtype.

Examples

>>> import numpy as np
>>> from solprop.translate.chebyshev import clenshaw
>>> float(clenshaw(0.5, np.array([1.0, 2.0, 3.0])))
0.5

solprop.translate.chebyshev.interleave_axis(coarse, fine, axis)[source]#

Merge a nested grid’s samples with the ones a doubling added, along one axis.

Parameters:

coarsenumpy.ndarray: Samples on the degree d grid, with d + 1 entries along axis.
finenumpy.ndarray: Samples at the d points refinement_points() added, same shape otherwise.
axisint: The axis the doubling ran along.

Returns:

numpy.ndarray: Samples on the degree 2 d grid, with 2 d + 1 entries along axis.

Examples

>>> import numpy as np
>>> from solprop.translate.chebyshev import interleave_axis
>>> coarse = np.array([[0.0, 2.0]])
>>> interleave_axis(coarse, np.array([[1.0]]), 1).tolist()
[[0.0, 1.0, 2.0]]

solprop.translate.chebyshev.lobatto_coefficients(values)[source]#

Transform samples at the Chebyshev-Lobatto points into Chebyshev coefficients.

The discrete Chebyshev transform is the type-I DCT of the samples up to the end weights; taking it as the FFT of the mirror extension lets complex-valued integrands go through the same path as real ones. The coefficients are those numpy.polynomial.chebyshev.chebval() and PiecewiseChebyshev evaluate, so one transform serves the adaptive builder here and the fixed-section expansions of solprop.translate.coupled.

Parameters:

valuesnumpy.ndarray: Samples at lobatto_points() in increasing position, shape (n_deg, *value_shape).

Returns:

numpy.ndarray: Coefficients of sum_j c_j T_j, shape (n_deg, *value_shape), real when the samples are real.

Examples

>>> import numpy as np
>>> from solprop.translate.chebyshev import lobatto_coefficients, lobatto_points
>>> nodes = lobatto_points(4, -1.0, 1.0)
>>> lobatto_coefficients(nodes**2).round(12) + 0.0
array([0.5, 0. , 0.5, 0. , 0. ])

solprop.translate.chebyshev.lobatto_points(degree, lo, hi)[source]#

Return the Chebyshev-Lobatto points of degree mapped onto [lo, hi].

Unlike the z sampler these are not inset at the ends: a parameter range has no discontinuity sitting on its own boundary, and its ends are values the paths actually reach, where the function must be evaluated as it stands.

Parameters:

degreeint: Polynomial degree, >= 1.
lo, hifloat: Ends of the range, in the parameter’s own units, with hi > lo.

Returns:

numpy.ndarray: The degree + 1 points, increasing from lo to hi.

Examples

>>> from solprop.translate.chebyshev import lobatto_points
>>> [round(value, 12) for value in lobatto_points(2, 0.0, 1.0).tolist()]
[0.0, 0.5, 1.0]

solprop.translate.chebyshev.refinement_points(degree, lo, hi)[source]#

Return the points a degree doubling adds to lobatto_points().

Parameters:

degreeint: Current degree; the result completes the grid of degree 2 * degree.
lo, hifloat: Ends of the range, in the parameter’s own units.

Returns:

numpy.ndarray: The degree new points, increasing, all interior.

Examples

>>> from solprop.translate.chebyshev import refinement_points
>>> [round(value, 6) for value in refinement_points(1, 0.0, 1.0).tolist()]
[0.5]

solprop.translate.chebyshev.resolve_axis(values, axis, *, tol)[source]#

Transform one axis to Chebyshev coefficients and judge its tail.

The samples along the other axes are left as they are, so the acceptance test measures the tail per element against a genuine function value, as it does for the state in the one-dimensional build, rather than against a coefficient that may itself be negligible.

Parameters:

valuesnumpy.ndarray: Samples at the Chebyshev-Lobatto points of axis, of any shape.
axisint: The axis to transform.
tolfloat: Target bound on the uniform error relative to the value scale, > 0. The criterion is build()’s: twice the relative coefficient tail, the factor two being the Chebyshev aliasing bound.

Returns:

numpy.ndarray: The same array with axis holding Chebyshev coefficients.
bool: Whether that axis is resolved to tol.

Examples

>>> import numpy as np
>>> from solprop.translate.chebyshev import lobatto_points, resolve_axis
>>> x = lobatto_points(8, -1.0, 1.0)
>>> coefficients, resolved = resolve_axis(2.0 * x**2 - 1.0, 0, tol=1e-12)
>>> resolved, [abs(round(value, 12)) for value in coefficients[:3].tolist()]
(True, [0.0, 0.0, 1.0])

solprop.translate.chebyshev.trim_axis(coefficients, axis)[source]#

Drop the trailing coefficients of one axis that are transform noise.

Parameters:

coefficientsnumpy.ndarray: Chebyshev coefficients along axis.
axisint: The axis to trim.

Returns:

numpy.ndarray: The same array with axis shortened to the coefficients that carry signal; at least one is always kept.

Examples

>>> import numpy as np
>>> from solprop.translate.chebyshev import trim_axis
>>> trim_axis(np.array([[1.0, 0.0, 0.0]]), 1).shape
(1, 1)