solprop.stepper.schedule#

Fixed step schedules: discover, calibrate, replay, and step landing.

Adaptive stepping does not remove the convergence study, because rtol is a local tolerance and the global error it buys is problem dependent. Fixed stepping is not less accurate, and it wins outright for ensembles (identical steps make the numerical error a common systematic that cancels rather than a per-member variance that adds), for reproducibility, and on a device. So both are first class:

discover – one adaptive run records h(z) as a FixedSchedule;
calibrate – calibrate() replays it at h, h/r and h/r**2 and measures the observed order, with a gate that refuses an order far above the scheme’s;
replay – replay() rescales the schedule to the step size that reaches a target error, and production runs, ensemble members and scan points execute that.

Separately, landing() is the per-step rule that clips a trial step so it ends exactly on the next save plane or hard kink, which is what replaces dense output (see docs/dev/decisions.md). See docs/dev/stepper.md sections 3 and 8.

Module Attributes

`ORDER_MARGIN`	How far the observed order may exceed the scheme's before it is refused.
`CONCLUSIVE_FACTOR`	How far above the round-off floor the coarsest difference must sit to be believed.

Functions

`calibrate`(schedule, evaluate, order, requested)	Replay a schedule at three levels and measure the order it converges at.
`discover`(steps[, z_start])	Record the step sizes an adaptive run chose, as a schedule that can be replayed.
`landing`(z, dz, targets)	Clip a trial step to the caps and the next position it must land on.
`observed_order`(coarse, fine[, refinement])	Return the convergence order implied by two successive differences.
`replay`(schedule, calibration, target_error)	Rescale a discovered schedule to the step size that reaches a target error.

Classes

`Calibration`(observed_order, declared_order, ...)	What a three-level replay measured, and what it licenses.
`FixedSchedule`(steps[, z_start])	A recorded sequence of step lengths covering one propagation.
`Landing`(dz, target, z_target[, stretched])	The step a landing rule allows, and what bound it.
`LandingTarget`(*values)	What bound the length of a step.
`StepTargets`(save_planes, hard_kinks, ...[, ...])	Everything that can shorten a trial step.

solprop.stepper.schedule.CONCLUSIVE_FACTOR = 1000.0#

How far above the round-off floor the coarsest difference must sit to be believed.

Below the floor the differences measure accumulated rounding rather than truncation and the ratio means nothing; the same protocol has returned an observed order of -0.83 there and would have been read as healthy.

class solprop.stepper.schedule.Calibration(observed_order, declared_order, refinement, differences, achieved_error, requested, conclusive, message)[source]#

Bases: object

What a three-level replay measured, and what it licenses.

Attributes:

observed_orderfloat: Order implied by the two differences, dimensionless; NaN when none could be read.
declared_orderint: Order the scheme claims, dimensionless.
refinementint: Ratio between successive step sizes of the ladder.
differencestuple: (|Q0 - Q1|, |Q1 - Q2|) for the observable Q at the three levels, in the units of the observable.
achieved_errorfloat: Richardson estimate of the error of the base schedule, same units. This is the global error the run actually delivered.
requestedfloat: Local tolerance the discovering run asked for, dimensionless.
conclusivebool: Whether the coarsest difference sits far enough above the round-off floor for the order to mean anything.
messagestr: One line saying what was measured, for logs and for the output metadata.

achieved_error: float#

conclusive: bool#

declared_order: int#

differences: tuple[float, float]#

message: str#

observed_order: float#

property ratio: float[source]#

Return the achieved global error as a multiple of the requested tolerance.

Returns:

float: achieved_error / requested. It is not 1 and it is not constant: the map from a local tolerance to the global error is problem dependent and misses by one to three orders in either direction, which is the whole reason this is measured rather than assumed.

refinement: int#

requested: float#

scale_for(target_error)[source]#

Return the step-size scale that would reach target_error.

Parameters:

target_errorfloat: Global error wanted, in the units of the observable, positive.

Returns:

float: Factor to pass to FixedSchedule.rescaled(), dimensionless.

Raises:

StepperError: If target_error is not positive.
ConvergenceError: If the calibration was inconclusive or read no order, so there is no measured exponent to extrapolate along. Assuming one is exactly the mistake calibration exists to prevent.

Examples

>>> from solprop.stepper import Calibration
>>> calibration = Calibration(
...     observed_order=5.0,
...     declared_order=5,
...     refinement=2,
...     differences=(1e-6, 3.1e-8),
...     achieved_error=1e-6,
...     requested=1e-8,
...     conclusive=True,
...     message="",
... )
>>> calibration.ratio
100.0
>>> f"{calibration.scale_for(1e-8):.4f}"
'0.3981'

class solprop.stepper.schedule.FixedSchedule(steps, z_start=0.0)[source]#

Bases: object

A recorded sequence of step lengths covering one propagation.

Equality is by identity (eq=False) because steps is an array.

Parameters:

stepsnumpy.ndarray: Step lengths of shape (n_steps,), m, strictly positive and finite.
z_startfloat, optional: Position the first step starts from, m.

Raises:

StepperError: If the schedule is empty or any step is not finite and positive.

Examples

>>> import numpy as np
>>> from solprop.stepper import FixedSchedule
>>> schedule = FixedSchedule(np.array([0.4, 0.4, 0.2]))
>>> schedule.n_steps, schedule.z_end
(3, 1.0)
>>> schedule.refined(2).n_steps
6

property edges: Float[ndarray, 'n_edges'][source]#

Return the positions the steps start and end at.

Returns:

numpy.ndarray: Positions of shape (n_steps + 1,), m, accumulated in step order so that they match what a driver marching the schedule would reach.

property max_step: float[source]#

Return the longest step.

Returns:

float: Length, m.

property min_step: float[source]#

Return the shortest step.

Returns:

float: Length, m.

property n_steps: int[source]#

Return the number of steps.

Returns:

int: Number of steps.

refined(factor)[source]#

Return the schedule with every step split into factor equal pieces.

This is the refinement the three-level ladder uses. It preserves the shape of the schedule, including wherever it landed on a kink or a save plane, so the levels are nested and the ratio of their differences measures truncation and nothing else.

Parameters:

factorint: Number of pieces per step, at least 1. A power of two divides a double exactly, so each piece is exactly h / factor; the accumulated endpoint can still move by an ulp or two, because the cumulative sum now has more terms in it.

Returns:

FixedSchedule: The refined schedule, covering the same span.

Raises:

StepperError: If factor is below 1.

Examples

>>> import numpy as np
>>> from solprop.stepper import FixedSchedule
>>> FixedSchedule(np.array([1.0, 0.5])).refined(2).steps.tolist()
[0.5, 0.5, 0.25, 0.25]

rescaled(factor)[source]#

Return a schedule of the same shape and span with steps scaled.

The step-size profile transfers while the endpoint stays put: the walk takes factor times the local step of the original schedule at each position, so a schedule discovered at one tolerance can be replayed at another without losing where it chose to go slowly. When the final piece comes out shorter than half its predecessor, the last two are evened up instead of leaving a sliver.

Parameters:

factorfloat: Scale applied to the local step, dimensionless and positive. Below 1 refines, above 1 coarsens.

Returns:

FixedSchedule: The rescaled schedule, spanning the same interval to within the rounding of the accumulated sum.

Raises:

StepperError: If factor is not finite and positive.

Examples

>>> import numpy as np
>>> from solprop.stepper import FixedSchedule
>>> schedule = FixedSchedule(np.array([0.4, 0.4, 0.2]))
>>> coarse = schedule.rescaled(2.0)
>>> coarse.steps.tolist(), coarse.z_end
([0.5, 0.5], 1.0)

steps: Float[ndarray, 'n_steps']#

property total_length: float[source]#

Return the propagation length the schedule covers.

Returns:

float: Length, m.

classmethod uniform(z_start, z_end, dz, *, save_planes=(), hard_kinks=())[source]#

Build a plain fixed-step schedule that lands on every plane and kink.

A fixed-step run with no adaptivity at all. The run is cut into segments at every save plane and hard kink, and each segment is covered by the fewest equal steps of at most dz, so no step exceeds dz and every target is an edge of the schedule (to rounding in the accumulated edges, which the driver snaps away as it does for a landing). There is no dense output, so a fixed schedule must contain the planes; this is the constructor that guarantees it.

Parameters:

z_startfloat: Where the run starts, m.
z_endfloat: Where it ends, m, after z_start.
dzfloat: Largest step, m, positive. Steps are shortened, never lengthened, to fit the segments.
save_planescollections.abc.Sequence or numpy.ndarray, optional: Positions the driver will save at, m, within [z_start, z_end]. Positions at either end are dropped.
hard_kinkscollections.abc.Sequence or numpy.ndarray, optional: Positions where the linear operator jumps or kinks, m, within [z_start, z_end]. Genuine kinks only, as for StepTargets.

Returns:

FixedSchedule: The schedule.

Raises:

StepperError: If z_end is not after z_start, dz is not finite and positive, or a position is non-finite or outside the run.

Examples

>>> from solprop.stepper import FixedSchedule
>>> FixedSchedule.uniform(0.0, 1.0, 0.3, save_planes=[0.5]).steps
array([0.25, 0.25, 0.25, 0.25])

property z_end: float[source]#

Return the position the last step ends at.

Returns:

float: Final position, m.

z_start: float = 0.0#

class solprop.stepper.schedule.Landing(dz, target, z_target, stretched=False)[source]#

Bases: object

The step a landing rule allows, and what bound it.

Attributes:

dzfloat: Step to take, m, positive and no longer than the trial step.
targetLandingTarget: What bound it.
z_targetfloat or None: The position the step ends on for a landing, m; None for a cap, which shortens the step without aiming it anywhere. A driver snaps its position to this value after the step rather than accumulating z + dz, so that a save plane is hit exactly however many steps have gone before.
stretchedbool: True if the step was lengthened beyond the caps to reach z_target.

dz: float#

stretched: bool = False#

target: LandingTarget#

z_target: float | None#

class solprop.stepper.schedule.LandingTarget(*values)[source]#

Bases: StrEnum

What bound the length of a step.

NONE means the controller’s own proposal stood. MAX_STEP and INTERVAL_WIDTH are caps: they shorten the step but the step ends nowhere in particular. HARD_KINK, SAVE_PLANE and END are landings: the step ends exactly on that position. When several coincide the later member of this list wins, so a save plane that sits on a kink is reported as the save plane.

END = 'end'#

HARD_KINK = 'hard_kink'#

INTERVAL_WIDTH = 'interval_width'#

MAX_STEP = 'max_step'#

NONE = 'none'#

SAVE_PLANE = 'save_plane'#

solprop.stepper.schedule.ORDER_MARGIN = 3#

How far the observed order may exceed the scheme’s before it is refused.

Three, because the two failure modes have to be told apart. A legitimate excursion – a strongly nonlinear stretch where the asymptotic constant is still moving – reads about half an order high. A schedule that is off the convergence curve altogether reads nine orders high: refining removed something that was not truncation error at all. A two-level extrapolant cannot distinguish them, which is why the gate needs three levels (docs/dev/stepper.md section 10, trap 4).

class solprop.stepper.schedule.StepTargets(save_planes, hard_kinks, breakpoints, dz_max, z_end, stretch=1.1)[source]#

Bases: object

Everything that can shorten a trial step.

Equality is by identity (eq=False) because the fields are arrays.

Parameters:

save_planesnumpy.ndarray: Positions of shape (n_z,) the driver will save at, m, strictly increasing. Steps land on them exactly, which is what replaces dense output.
hard_kinksnumpy.ndarray: Positions of shape (n_kinks,) where the linear operator jumps or kinks – pressure-stage boundaries, taper corners – m, strictly increasing. Only genuine kinks belong here, never the edges the adaptive Chebyshev split introduces to resolve smooth variation: a typical fibre has thirty intervals and two kinks, so stopping at every edge would force fifteen times more boundaries than the physics needs.
breakpointsnumpy.ndarray: Chebyshev interval edges of shape (n_breakpoints,), m, strictly increasing, at least two.
dz_maxfloat: Largest step allowed anywhere, m, positive.
z_endfloat: Where the propagation stops, m. The last step lands on it exactly rather than stepping past it.
stretchfloat, optional: How far a step may be lengthened, as a factor on the capped trial step, to land on a target just beyond it instead of taking the step and then a sliver. The default, 1.1, is the controller’s growth cap: a stretch the controller already permits every step, raising the local error by at most 1.1**5. Pass 1 to switch stretching off. A stretch never exceeds dz_max.

Raises:

StepperError: If an array is not strictly increasing, breakpoints has fewer than two entries, dz_max is not positive, or stretch is below 1.

breakpoints: Float[ndarray, 'n_breakpoints']#

dz_max: float#

hard_kinks: Float[ndarray, 'n_kinks']#

save_planes: Float[ndarray, 'n_z']#

stretch: float = 1.1#

z_end: float#

solprop.stepper.schedule.calibrate(schedule, evaluate, order, requested, refinement=2, floor=0.0)[source]#

Replay a schedule at three levels and measure the order it converges at.

The second step, calibration. evaluate is asked for a scalar observable at h, h/r and h/r**2; the two differences give the observed order and, by Richardson extrapolation, the error the base schedule actually carried.

The observable is a scalar by design: this layer never touches an array of the field shape (docs/dev/architecture.md section 3). The driver supplies whatever functional it wants calibrated – an energy, a peak power, a band’s content.

Parameters:

scheduleFixedSchedule: The base schedule, usually one that discover() recorded.
evaluatecollections.abc.Callable: Runs a schedule and returns a scalar observable.
orderint: Order the scheme claims, at least 1.
requestedfloat: Local tolerance the discovering run asked for, dimensionless and positive.
refinementint, optional: Ratio between successive step sizes, at least 2.
floorfloat, optional: Round-off floor of the observable, same units, non-negative. When the coarsest difference is within CONCLUSIVE_FACTOR of it, the result is marked inconclusive and no order is pronounced. Left at 0 every ladder is believed, which is only safe when the caller knows the differences are above the floor.

Returns:

Calibration: The measurement.

Raises:

StepperError: If order, requested or refinement is out of range.
ConvergenceError: If the observed order exceeds order + ORDER_MARGIN. Refining removed something that is not truncation error, so the base schedule is off the convergence curve rather than merely large, and extrapolating from it would be meaningless.

Examples

>>> import numpy as np
>>> from solprop.stepper import FixedSchedule, calibrate
>>> base = FixedSchedule(np.full(8, 0.125))
>>> # A synthetic order-5 observable: Q(h) = 1 + h**5.
>>> quantity = lambda s: 1.0 + s.max_step**5
>>> calibration = calibrate(base, quantity, order=5, requested=1e-8)
>>> round(calibration.observed_order, 6)
5.0

solprop.stepper.schedule.discover(steps, z_start=0.0)[source]#

Record the step sizes an adaptive run chose, as a schedule that can be replayed.

This is the first step, discovery. The driver collects the length of every accepted step; this turns that record into a schedule, which calibrate() then measures and replay() rescales.

Parameters:

stepscollections.abc.Sequence: Length of every accepted step, m, in the order they were taken.
z_startfloat, optional: Position the run started from, m.

Returns:

FixedSchedule: The recorded schedule.

Raises:

StepperError: If the record is empty or holds a step that is not finite and positive.

Examples

>>> from solprop.stepper import discover
>>> discover([0.1, 0.2, 0.15]).n_steps
3

solprop.stepper.schedule.landing(z, dz, targets)[source]#

Clip a trial step to the caps and the next position it must land on.

Two kinds of bound apply (docs/dev/stepper.md section 3):

caps – dz_max, and the width of the Chebyshev interval containing z. The containing interval, never the narrowest one anywhere: on a real fibre the finest local scale is four orders below the median, and capping globally by it would turn a metre that is smooth almost everywhere into over a million steps (section 10, trap 2).
landings – the next hard kink, the next save plane and the end of the run, each of which the step ends exactly on.

The nearest landing wins when it lies within the capped step. When it lies just beyond, within targets.stretch times the capped step and no further than dz_max, the step is stretched to reach it rather than taking the step and then a sliver; Hairer’s DOPRI5 stretches its final step by one percent and MATLAB’s ode45 by ten. The floor dz_min is never consulted here: a landing may be as short as the geometry demands (solprop.stepper.controller.default_dz_min()).

Parameters:

zfloat: Current position, m, inside the range the breakpoints cover and before z_end.
dzfloat: Step the controller proposed, m, positive.
targetsStepTargets: The caps and landing positions.

Returns:

Landing: The allowed step and what bound it.

Raises:

StepperError: If dz is not positive, or z is outside the breakpoints or at or past z_end.

Examples

>>> import numpy as np
>>> from solprop.stepper import StepTargets, landing
>>> targets = StepTargets(
...     save_planes=np.array([0.25, 0.5, 0.75, 1.0]),
...     hard_kinks=np.array([0.4]),
...     breakpoints=np.array([0.0, 1.0]),
...     dz_max=0.3,
...     z_end=1.0,
... )
>>> stop = landing(0.0, 1.0, targets)
>>> stop.dz, stop.target.value, stop.z_target
(0.25, 'save_plane', 0.25)
>>> landing(0.25, 1.0, targets).z_target      # the kink, before the next plane
0.4
>>> landing(0.9, 0.05, targets).target        # nothing binds
<LandingTarget.NONE: 'none'>
>>> stop = landing(0.0, 0.24, targets)        # 0.25 lies within 1.1 x 0.24
>>> stop.z_target, stop.stretched
(0.25, True)

solprop.stepper.schedule.observed_order(coarse, fine, refinement=2)[source]#

Return the convergence order implied by two successive differences.

Parameters:

coarsefloat: Difference between the two coarsest levels, in the units of the observable.
finefloat: Difference between the two finest levels, same units.
refinementint, optional: Ratio between successive step sizes.

Returns:

float: log(coarse / fine) / log(refinement), or NaN when either difference is not positive, in which case no order can be read.

Examples

>>> from solprop.stepper import observed_order
>>> observed_order(1.6e-3, 1e-4)
4.0

solprop.stepper.schedule.replay(schedule, calibration, target_error)[source]#

Rescale a discovered schedule to the step size that reaches a target error.

The third step, replay. Production runs, ensemble members and scan points execute the schedule this returns, so their numerical error is a common systematic rather than a per-run variance.

Parameters:

scheduleFixedSchedule: The schedule that was calibrated.
calibrationCalibration: The measurement from calibrate().
target_errorfloat: Global error wanted, in the units of the calibrated observable, positive.

Returns:

FixedSchedule: The schedule to run.

Raises:

StepperError: If target_error is not positive.
ConvergenceError: If the calibration was inconclusive, so no exponent was measured.

Examples

>>> import numpy as np
>>> from solprop.stepper import FixedSchedule, calibrate, replay
>>> base = FixedSchedule(np.full(8, 0.125))
>>> quantity = lambda s: 1.0 + s.max_step**5
>>> calibration = calibrate(base, quantity, order=5, requested=1e-8)
>>> replay(base, calibration, calibration.achieved_error / 32).n_steps
16