solprop.stepper.controller#

Step-size control: the error-norm scaling and the PI controller.

The backend reduces the embedded pair’s error estimate to a scalar; everything here is scalar logic on that number (docs/dev/stepper.md section 3). Two pieces make up the control loop:

a norm scale, the denominator the backend divides the error by, so that a scaled norm of 1 means the step sat exactly at the requested tolerance;
the controller, which turns a history of scaled norms into a step-size factor.

The controller is the proportional-integral rule of Soederlind and Wang [1] rather than the elementary rule of [2], because the integral term damps the response to fluctuations in the error estimate. It is a frozen dataclass and PIController.propose() is a pure function of the norm history, so replaying a run’s norms reproduces its step sizes bitwise. That is what makes a fixed schedule a faithful replay rather than an approximation.

References#

[1]

G. Soederlind and L. Wang, “Adaptive time-stepping and computational stability”, J. Comput. Appl. Math. 185, 225 (2006). doi:10.1016/j.cam.2005.03.008

[2]

E. Hairer, S. P. Noersett and G. Wanner, Solving Ordinary Differential Equations I, 2nd ed., Springer (1993), section II.4.

Module Attributes

`START_GROWTH`	Largest factor the start-up rule may take the first step's successor by.
`DZ_MIN_RELATIVE`	Default step floor as a fraction of the propagation length.

Functions

`default_dz_min`(length)	Return the default step floor for a run of the given length.
`limit_step`(dz, factor, dz_min, dz_max)	Apply a step-size factor and clamp the result to the allowed range.
`weak_norm_scale`(norm_y, norm_y_new, rtol, atol)	Return the weak error-norm scale, one global number for the whole state.
`weighted_norm_scale`(norm_y, norm_y_new, ...)	Return a per-bin error-norm scale that tightens chosen parts of the spectrum.

Classes

`ControlRecord`(requested_rtol, ...)	What a completed run achieved against what it asked for.
`PIController`(rtol, atol[, safety, ...])	Proportional-integral step-size controller with lag and a reject-free mode.

class solprop.stepper.controller.ControlRecord(requested_rtol, requested_atol, n_accepted, n_rejected, max_norm, mean_norm, worst_step)[source]#

Bases: object

What a completed run achieved against what it asked for.

The scaled error norm is achieved-over-requested for the local error: 1 means the step sat exactly at the requested tolerance. The map from that local tolerance to the global error at the output is a different matter entirely – it is problem dependent and misses by one to three orders in either direction, which is why every run records what it achieved instead of assuming (docs/dev/stepper.md section 3). The global ratio comes from solprop.stepper.schedule.Calibration.

Attributes:

requested_rtolfloat: Relative tolerance the run asked for, dimensionless.
requested_atolfloat: Absolute tolerance the run asked for, in state units.
n_acceptedint: Number of accepted steps.
n_rejectedint: Number of rejected attempts.
max_normfloat: Largest scaled error norm over the accepted steps, dimensionless; 0 for a run with no accepted steps.
mean_normfloat: Mean scaled error norm over the accepted steps, dimensionless; 0 for a run with no accepted steps.
worst_stepint: Index, among the accepted steps, of the one with the largest norm; -1 for a run with no accepted steps. The step history says where it was.

classmethod from_norms(controller, accepted_norms, n_rejected=0)[source]#

Summarise a run from the scaled norms of its accepted steps.

Parameters:

controllerPIController: The controller the run used, which carries the requested tolerances.
accepted_normscollections.abc.Sequence: Scaled error norm of every accepted step, dimensionless.
n_rejectedint, optional: Number of rejected attempts.

Returns:

ControlRecord: The summary.

Examples

>>> from solprop.stepper import ControlRecord, PIController
>>> control = PIController(rtol=1e-8, atol=0.0)
>>> record = ControlRecord.from_norms(control, [0.2, 0.8, 0.5], n_rejected=1)
>>> record.tolerance_ratio, f"{record.rejection_rate:.2f}"
(0.8, '0.25')

max_norm: float#

mean_norm: float#

property met_tolerance: bool[source]#

Return whether every accepted step sat within the requested tolerance.

A fixed schedule is replayed without control, so this is the check that a fixed run was not too coarse: the embedded estimate is computed regardless, and a largest norm above 1 means a step exceeded the tolerance the norm is scaled by. An adaptive run in the rejecting mode meets it by construction; a reject-free run may not, by design.

Returns:

bool: True if the largest scaled norm is at most 1.

n_accepted: int#

n_rejected: int#

property rejection_rate: float[source]#

Return the fraction of attempts that were rejected.

Returns:

float: Between 0 and 1; 0 for a run with no attempts.

requested_atol: float#

requested_rtol: float#

property tolerance_ratio: float[source]#

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

Returns:

float: The largest scaled norm over the accepted steps. Below 1 the run stayed inside its local tolerance everywhere.

worst_step: int#

solprop.stepper.controller.DZ_MIN_RELATIVE = 1e-09#

Default step floor as a fraction of the propagation length.

The floor exists to fail fast, not to protect accuracy: a run whose controller wants a step below a billionth of the length would need over a billion steps and never finish, so raising at once turns hours of silence into an immediate ConvergenceError. The finest legitimate scale on a real gradient fibre sits about four orders below the median step (docs/dev/stepper.md section 10, trap 2), five orders above this floor, so it never binds a healthy run.

class solprop.stepper.controller.PIController(rtol, atol, safety=0.8, max_factor=1.1, min_factor=0.2, beta1=0.12, beta2=-0.04, lag=1, reject_free=False, weights=None, estimator_order=4, region='active')[source]#

Bases: object

Proportional-integral step-size controller with lag and a reject-free mode.

The factor applied after a step whose scaled norms are e_n and e_{n-1} is

factor = clip(safety * (1/e_n)**beta1 * (1/e_{n-1})**beta2, min_factor, max_factor),

the Lund (PI.4.2) rule of [1]. A pure proportional controller reacts to every fluctuation of the estimate; the integral term damps that, which is what keeps the step size smooth enough for the step history to be worth replaying.

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

Parameters:

rtolfloat: Relative tolerance, dimensionless and positive.
atolfloat: Absolute tolerance, in state units, non-negative.
safetyfloat: Factor applied to every proposal, dimensionless, in (0, 1]. The default 0.8 is the value measured with the reject-free mode below.
max_factorfloat: Largest growth per step, dimensionless and above 1. The default 1.1 is a tight cap by the usual standards; it is what makes the reject-free mode safe, and SEP 14 measured a looser one to cost accuracy in weak spectral bands, which the norm does not see, for little saving once the start-up rule (start_factor()) has removed the slow growth from dz_init.
min_factorfloat: Smallest shrink per step, dimensionless, in (0, 1).
beta1float: Proportional gain, dimensionless and positive.
beta2float: Integral gain, dimensionless. The defaults 3/25 and -1/25 are Soederlind’s (0.6, -0.2) already divided by k = estimator_order + 1 = 5. The exponent follows the estimator’s order, not the scheme’s: an order-5 pair with an order-4 embedded companion has k = 5. Use for_estimator_order() for any other estimator.
lagint: Number of steps by which the update is deferred, so the factor computed from step n is applied to step n + 1 + lag. This lets the error norm come back from a device asynchronously. Lag 1 was measured to cost nothing (error x1.00-1.05, cost x0.99-1.00); lag 2 drifted to x2.45 on one configuration.
reject_freebool: Never reject a step whose error norm is finite, letting the lagged update absorb the overshoot. Rejections happen 0.0-0.4 % of the time in practice, so this abolishes something that was barely happening, and it removes the per-step device synchronisation. A non-finite norm is still rejected: no lagged update can absorb a poisoned state.
weightsnumpy.ndarray or None: Per-bin weights of shape (n_omega,) selecting weighted_norm_scale(); None selects weak_norm_scale().
estimator_orderint: Order q of the error estimate, which sets the exponent of the shrink applied to a rejected step. 4 for an embedded order-5(4) pair.
regionstr: The spectral bins the norm scores. "active" (the default) scores every bin the propagation evolves, omega_min / 2 < omega < 1.1 omega_max: the requested band and the margin around it, where generated light propagates and the spectral collar absorbs it. "band" scores only the requested wavelength limits, the norm before this choice: light in the margin is then integrated without control, and can be wrong by as much as it holds while the band’s error stays within tolerance (docs/dev/stepper.md).

Raises:

StepperError: If a tolerance, gain or factor is outside its valid range, or the region is not "active" or "band".

References

[1]

G. Soederlind and L. Wang, “Adaptive time-stepping and computational stability”, J. Comput. Appl. Math. 185, 225 (2006). doi:10.1016/j.cam.2005.03.008

Examples

>>> from solprop.stepper import PIController
>>> control = PIController(rtol=1e-8, atol=1e-30)
>>> control.accepts(0.5), control.accepts(2.0)
(True, False)
>>> f"{control.propose([0.5, 0.5]):.4f}"
'0.8694'

accepts(norm)[source]#

Return whether a step with this scaled error norm is accepted.

Parameters:

normfloat: Scaled error norm of the trial step, dimensionless.

Returns:

bool: True to keep the step. reject_free keeps every step with a finite norm; a NaN or infinite norm is rejected whatever the mode, because the state is already poisoned and no later update can undo that.

Examples

>>> import math
>>> from solprop.stepper import PIController
>>> free = PIController(rtol=1e-8, atol=0.0, reject_free=True)
>>> free.accepts(1e6), free.accepts(math.nan)
(True, False)

atol: float#

beta1: float = 0.12#

beta2: float = -0.04#

equivalent_norm(norm, dz_taken, dz_proposed)[source]#

Rescale the norm of a clipped or stretched step to the step proposed.

Landing shortens a step onto a save plane or a kink, or stretches it a little to reach one; the controller must not learn from that step as if it had chosen it. Its state is the proposal: the driver bases the next step on the proposal, not on the executed step, and feeds the history the norm the proposal would have scored under the asymptotic error model, norm * (dz_proposed / dz_taken) ** (q + 1) with q the estimator’s order. Without this a step onto a plane restarts the controller from the sliver and it regrows at max_factor per step; a thousandfold recovery at 1.1 takes over seventy steps, and over two hundred planes that is more steps than the run needs.

The same rescaling serves a rejection: a clipped step that failed means the longer proposal would have failed harder, and rejection_factor applied to the equivalent norm shrinks the proposal accordingly.

Parameters:

normfloat: Scaled error norm of the step actually taken, dimensionless. A non-finite norm is returned as it is: it must still be rejected.
dz_takenfloat: Length of the step taken, m, positive.
dz_proposedfloat: Length the controller proposed before landing, m, positive.

Returns:

float: The norm to enter into the history, dimensionless.

Raises:

StepperError: If either length is not finite and positive.

Examples

>>> from solprop.stepper import PIController
>>> control = PIController(rtol=1e-8, atol=0.0)
>>> control.equivalent_norm(0.3, dz_taken=1e-3, dz_proposed=1e-3)
0.3
>>> f"{control.equivalent_norm(1e-6, dz_taken=1e-4, dz_proposed=1e-3):.2f}"
'0.10'

estimator_order: int = 4#

classmethod for_estimator_order(estimator_order, rtol, atol)[source]#

Build a controller whose PI gains match an error estimate of order q.

The Lund gains are beta1 = 0.6 / k and beta2 = -0.2 / k with k = q + 1. Using the scheme’s order in place of the estimator’s is a standing mistake: it makes the controller oscillate rather than converge.

Parameters:

estimator_orderint: Order q for which the error estimate behaves as C h**(q+1), at least 1.
rtolfloat: Relative tolerance, dimensionless and positive.
atolfloat: Absolute tolerance, in state units, non-negative.

Returns:

PIController: A controller with every other setting at its default.

Examples

>>> from solprop.stepper import PIController
>>> control = PIController.for_estimator_order(4, rtol=1e-8, atol=0.0)
>>> control.beta1 == 3 / 25 and control.beta2 == -1 / 25
True

lag: int = 1#

max_factor: float = 1.1#

min_factor: float = 0.2#

propose(norm_history)[source]#

Return the step-size factor for the next step.

A pure function of the history, so the same norms always give the same factors, bit for bit. That is what lets a discovered schedule be replayed.

Parameters:

norm_historycollections.abc.Sequence: Scaled error norms of the accepted steps so far, oldest first, dimensionless. A value of 1 means the step sat exactly at tolerance.

Returns:

float: Factor to multiply the step size by, dimensionless, in [min_factor, max_factor].

Examples

>>> from solprop.stepper import PIController
>>> control = PIController(rtol=1e-8, atol=0.0, lag=1)
>>> control.propose([])          # nothing to go on yet
1.0
>>> control.propose([0.1])       # the lag has not filled yet either
1.0
>>> f"{control.propose([0.1, 0.2]):.4f}"
'1.0546'
>>> control.propose([1e-6, 0.2])   # growth is capped
1.1

region: Literal['active', 'band'] = 'active'#

reject_free: bool = False#

rejection_factor(norm)[source]#

Return the factor to shrink by after a rejected step.

Applied at once, without the lag: a rejected step has to be retried now. The exponent is the deadbeat 1/(q+1) of the estimator’s order, and the result can never exceed 1, so a rejection cannot grow the step.

Parameters:

normfloat: Scaled error norm of the rejected step, dimensionless.

Returns:

float: Factor in [min_factor, 1], dimensionless.

Examples

>>> from solprop.stepper import PIController
>>> control = PIController(rtol=1e-8, atol=0.0)
>>> f"{control.rejection_factor(4.0):.4f}"
'0.6063'

rtol: float#

safety: float = 0.8#

scale(norm_y, norm_y_new)[source]#

Return the scale the backend divides the error estimate by.

Parameters:

norm_y, norm_y_newfloat: Euclidean norms of the state before and after the trial step, in state units, non-negative.

Returns:

float or numpy.ndarray: A single scale for the weak norm, or one per spectral bin when weights were given.

start_factor(norm)[source]#

Return the factor the start-up rule takes the first accepted step by.

The first step of a run is a guess, dz_init. Its own scaled norm says how far it is from the tolerance, so the step after it is the elementary deadbeat one, safety * norm**(-1/(q+1)) with q the estimator’s order, at most START_GROWTH times the first and never below min_factor. A run then starts near its natural step instead of growing to it at the cap of 1.1, which from the usual dz_init takes tens of steps (SEP 14 section 4.2).

Parameters:

normfloat: Scaled error norm of the first accepted step at its proposed length, dimensionless.

Returns:

float: Factor in [min_factor, START_GROWTH], dimensionless.

Examples

>>> from solprop.stepper import PIController
>>> control = PIController(rtol=1e-6, atol=0.0)
>>> f"{control.start_factor(1e-5):.2f}"
'8.00'
>>> control.start_factor(0.0), control.start_factor(1e-30)
(100.0, 100.0)

weights: Float[ndarray, 'n_omega'] | None = None#

solprop.stepper.controller.START_GROWTH = 100.0#: Largest factor the start-up rule may take the first step’s successor by.

solprop.stepper.controller.default_dz_min(length)[source]#

Return the default step floor for a run of the given length.

DZ_MIN_RELATIVE times the length. The floor applies to the controller’s proposal, never to a step that landing clipped: a short step onto a save plane or a kink is a landing, not a failure to meet the tolerance, and solprop.stepper.schedule.landing() does not consult the floor.

Parameters:

lengthfloat: Propagation length, m, finite and positive.

Returns:

float: The floor, m.

Raises:

StepperError: If length is not finite and positive.

Examples

>>> from solprop.stepper import default_dz_min
>>> default_dz_min(2.0)
2e-09

solprop.stepper.controller.limit_step(dz, factor, dz_min, dz_max)[source]#

Apply a step-size factor and clamp the result to the allowed range.

Falling below dz_min is a failure, not something to clamp away. Clamping to the minimum step and accepting it would let a run that cannot meet its tolerance silently abandon it with no warning, no counter and no record in the output, so this raises instead.

Parameters:

dzfloat: The step just taken, m, positive.
factorfloat: Step-size factor from the controller, dimensionless and positive.
dz_minfloat: Smallest step the run may take, m, non-negative.
dz_maxfloat: Largest step the run may take, m, at least dz_min.

Returns:

float: The next step, m, in [dz_min, dz_max].

Raises:

StepperError: If dz_max is smaller than dz_min.
ConvergenceError: If the requested step falls below dz_min.

Examples

>>> from solprop.stepper import limit_step
>>> limit_step(1e-3, 1.1, dz_min=1e-9, dz_max=1e-2)
0.0011
>>> limit_step(1e-3, 1000.0, dz_min=1e-9, dz_max=1e-2)
0.01

solprop.stepper.controller.weak_norm_scale(norm_y, norm_y_new, rtol, atol)[source]#

Return the weak error-norm scale, one global number for the whole state.

The scaled norm the controller consumes is ‖e‖ / scale with

scale = rtol * max(‖y‖, ‖y_new‖) + atol,

reduced over the in-band mask. A single global scale is deliberate. Scaling per element instead is catastrophic here: the in-band dynamic range of a hollow-fibre spectrum reaches 1e17 and roughly a third of in-band bins carry no light, so a 5e-15 error in an empty bin scores several thousand and the step collapses to zero before a single step is accepted (docs/dev/stepper.md section 10, trap 1).

Note that atol is additive, not a floor on the state norm: a floor would never engage at realistic field magnitudes.

Parameters:

norm_y, norm_y_newfloat: Euclidean norms of the state before and after the trial step, in state units (sqrt(W) for modal amplitudes), non-negative.
rtolfloat: Relative tolerance, dimensionless and positive.
atolfloat: Absolute tolerance, in state units, non-negative.

Returns:

float: The scale to divide the error norm by, in state units.

Examples

>>> from solprop.stepper import weak_norm_scale
>>> f"{weak_norm_scale(2.0, 3.0, rtol=1e-8, atol=0.0):.3e}"
'3.000e-08'
>>> weak_norm_scale(0.0, 0.0, rtol=1e-8, atol=1e-30)
1e-30

solprop.stepper.controller.weighted_norm_scale(norm_y, norm_y_new, rtol, atol, weights)[source]#

Return a per-bin error-norm scale that tightens chosen parts of the spectrum.

scale_i = rtol * w_i * max(‖y‖, ‖y_new‖) + atol. The magnitude stays global, so this cannot reproduce the collapse that per-bin magnitudes cause (see weak_norm_scale()); the weights only redistribute the tolerance. A weight below 1 tightens that bin, and all-ones reproduces the weak scale exactly. This is the option for observables five to eight orders below the pump, such as a dispersive wave, which a whole-field norm leaves with no error control at all.

Note

Provisional. Which weighting actually serves weak observables best – fixed per-band weights, per-mode weights or a rolling per-bin scale – is the first experiment of docs/dev/stepper.md section 9 and is not yet settled.

Parameters:

norm_y, norm_y_newfloat: Euclidean norms of the state before and after the trial step, in state units, non-negative.
rtolfloat: Relative tolerance, dimensionless and positive.
atolfloat: Absolute tolerance, in state units, non-negative.
weightsnumpy.ndarray: Positive, dimensionless weights of shape (n_omega,), one per spectral bin.

Returns:

numpy.ndarray: Per-bin scales of shape (n_omega,), in state units.

Examples

>>> import numpy as np
>>> from solprop.stepper import weak_norm_scale, weighted_norm_scale
>>> weights = np.array([1.0, 0.1])
>>> weighted_norm_scale(2.0, 3.0, 1e-8, 0.0, weights)
array([3.e-08, 3.e-09])
>>> ones = np.ones(2)
>>> bool(
...     np.array_equal(
...         weighted_norm_scale(2.0, 3.0, 1e-8, 0.0, ones),
...         np.full(2, weak_norm_scale(2.0, 3.0, 1e-8, 0.0)),
...     )
... )
True