API reference¶

Generated from the docstrings in the source, so this page cannot drift from the code.

What each subpackage answers¶

One line each. These are the first sentence of each subpackage's own docstring, so they are also what help() shows you:

Subpackage	It answers
`ppsds.ambiguity`	What am I willing to admit about the system? — the belief, written as a set
`ppsds.dynamics`	The system itself — known dynamics and mechanical models, plus identifying a nominal model from rollouts
`ppsds.probabilistic_predictor`	A prediction is a distribution, not a point. The kinds sit side by side: `drpp/` derives one from an ambiguity set, `pmf/` takes one you supply
`ppsds.scoring_rule`	How well did a prediction do? — one number
`ppsds.predictability`	How predictable is the system at all? — the ceiling, as a reference line
`ppsds.environment`	Run it as an environment — for a system you declared an ambiguity set for
`ppsds.numerics`	Small numerical tools, and the weighted-atom pmf
`ppsds.control_policy`	A control policy, for the experiments that need one
`ppsds.experiment`	Running experiments at scale: the batch driver and its validated cache
`ppsds.visualization`	Drawing the results: score curves, covariance ellipses, trajectory overlays

The module path then states the scope layer by layer: ambiguity/additive_drift/conic_moment reads as assumption (additive drift) / family (conic moment), which is why the class names need not repeat it. Swapping the family means adding a module at that level, not editing this one.

Two layers inside one package¶

It matters which one you are importing:

	modules	dependencies
core	`ambiguity`, `dynamics`, `probabilistic_predictor`, `environment`, `numerics`, `scoring_rule`, `control_policy`	numpy, scipy, cvxpy, gymnasium
periphery	`experiment`, `visualization`	additionally matplotlib, pandas, seaborn, tqdm (`pip install "ppsds[experiments]"`)

A core module never imports the periphery, so pip install ppsds does not pull matplotlib in.

Predictors¶

The DRPP over conic-moment ambiguity sets:

conic_moment ¶

DRPP predictors over an additive-drift + conic-moment ambiguity set.

The path probabilistic_predictor / drpp / additive_drift / conic_moment states the scope layer by layer: DRPP family / additive-drift assumption / conic-moment ambiguity set. Another ambiguity-set family gets a module on the same level -- :func:eigvec_fixed_minmax solves equation (19) of the paper and holds only for the conic-moment ambiguity set; a different family needs a predictor of its own.

Three predictors

All three share the same private skeleton :class:_ConicMomentPredictor; their only substantive difference is "where the predictive covariance comes from":

================================= ================================================== :class:NominalPredictor Nominal: Sigma_pred = Sigma_bar :class:NoiseDrpp Noise-robust only: Sigma_pred = gamma[2] * Sigma_bar :class:EigDrpp The paper's Eig-DRPP: solves :func:eigvec_fixed_minmax each step ================================= ==================================================

.. note:: Noise-DRPP now has a class of its own. It used to have none and was stood in for by "handing the nominal predictor an inflated gamma[2] * Sigma_bar" -- two roles mixed together, and the name gave nothing away. This class internally computes that product first and then hands it to the skeleton (the same expression as the old code), so the numbers are bit-for-bit unchanged.

Core optimization

:func:eigvec_fixed_minmax -- solves equation (19) of the paper: it turns the functional maximin into a Euclidean cone program (CVXPY + SCS) and returns (optimal eigenvalues, optimal value, worst direction j).

Relation to the predictor contract

All three classes satisfy :class:~ppsds.probabilistic_predictor.base.ProbabilisticPredictor (the three methods predict / update / reset, and the three attributes x_pred / Sigma_pred / lamb_pred) but do not inherit from it -- that is a typing.Protocol, with zero runtime overhead and recognition by structure. This repository has tripped over "providing an interface" before: the earlier ABC wrote the attributes as @property (data descriptors) while the implementing classes used plain attribute assignment, so to actually inherit it would raise AttributeError at construction time, and ABCMeta.register() hid that incompatibility as well. If you give an interface you have to be able to keep it, so here we deliberately keep a structural contract and add no inheritance.

Usage

::

from ppsds.probabilistic_predictor.drpp.additive_drift.conic_moment import (
    NominalPredictor, NoiseDrpp, EigDrpp,
)

pred = EigDrpp(f_bar, mu_bar, Sigma_bar, gamma)
pred.update(x, u)      # solve one cone program
pred.predict()         # this step's predictive distribution
pred.lamb_pred         # predictive covariance eigenvalues -- the adversarial
                       # dynamics selects the worst direction from these

NominalPredictor ¶

NominalPredictor(f_bar, mu_bar, Sigma_bar, multivariate_rv=multivariate_normal)

Bases: _ConicMomentPredictor

The nominal predictor: no robustification at all.

The predictive mean uses the nominal dynamics f_bar, and the predictive covariance is simply the nominal covariance Sigma_bar.

Parameters:

Name	Type	Description	Default
`f_bar`	`callable`	The nominal dynamics `f(x, u) -> x_next`.	required
`mu_bar`	`(x_dim,) ndarray`	Nominal noise mean.	required
`Sigma_bar`	`(x_dim, x_dim) ndarray`	Nominal noise covariance.	required
`multivariate_rv`	`scipy distribution class`	The predictive distribution family (default `multivariate_normal`).	`multivariate_normal`

NoiseDrpp ¶

NoiseDrpp(f_bar, mu_bar, Sigma_bar, gamma, multivariate_rv=multivariate_normal)

Bases: _ConicMomentPredictor

Noise-DRPP: robust to noise only -- inflates the nominal covariance by gamma[2].

The prediction law is word-for-word the same as :class:NominalPredictor; the only difference is "which covariance it acts on".

.. note:: This class internally computes gamma[2] * Sigma_bar first and then hands it to the skeleton, the same expression as the formulation that preceded it ("handing the nominal predictor an inflated covariance"), so the numbers are bit-for-bit unchanged. That order is not arbitrary: eig_real(γ₂·Σ̄) and γ₂ · eig_real(Σ̄) may differ in their last bits, and lamb_pred is fed to the adversarial dynamics, where a difference in the last bits gets amplified.

Parameters:

Name	Type	Description	Default
`f_bar`	`callable`	The nominal dynamics `f(x, u) -> x_next`.	required
`mu_bar`	`(x_dim,) ndarray`	Nominal noise mean.	required
`Sigma_bar`	`(x_dim, x_dim) ndarray`	The nominal covariance -- the inflation happens inside this class, so callers pass the nominal value.	required
`gamma`	`list of length 3`	The ambiguity-set radius; only `gamma[2]` (the covariance scaling factor) is used.	required
`multivariate_rv`	`scipy distribution class`	The predictive distribution family (default `multivariate_normal`).	`multivariate_normal`

EigDrpp ¶

EigDrpp(f_bar, mu_bar, Sigma_bar, gamma, multivariate_rv=multivariate_normal)

Bases: _ConicMomentPredictor

Eig-DRPP: the distributionally robust predictor with eigenvectors fixed.

It solves :func:eigvec_fixed_minmax once per step; only the eigenvalues of the predictive covariance are optimized, while the eigenvectors are fixed to those of Sigma_bar -- that is what "restricted" refers to. The predictive mean is the same as the nominal predictor's.

This is the suboptimal solution of DRPP; the paper bounds its gap to the global optimum.

Parameters:

Name	Type	Description	Default
`f_bar`	`callable`	The nominal dynamics `f(x, u) -> x_next`.	required
`mu_bar`	`(x_dim,) ndarray`	Nominal noise mean.	required
`Sigma_bar`	`(x_dim, x_dim) ndarray`	Nominal noise covariance (its eigenvectors are used directly as the fixed basis).	required
`gamma`	`list of length 3`	The ambiguity-set radius, in the same format as :func:`eigvec_fixed_minmax`.	required
`multivariate_rv`	`scipy distribution class`	The predictive distribution family (default `multivariate_normal`).	`multivariate_normal`

Attributes:

Name	Type	Description
`j`	`int or None`	The index of the worst eigen-direction selected at this step (diagnostic).
`val_wst`	`float or None`	The worst-case expected-log-score lower bound at this step (diagnostic).

reset ¶

reset()

Restore a reusable state (restore the nominal quantities and redo the eigendecomposition).

gamma is never rewritten, so it needs no restoring. See the contract notes in :mod:ppsds.probabilistic_predictor.base.

trajectory_prediction ¶

trajectory_prediction(scoring_rule, x_traj, u_traj, abg_params=None)

Run the full predict-update cycle over one trajectory.

.. deprecated:: Deprecated; please migrate to the Gymnasium layer in :mod:ppsds.environment.

Parameters:

Name	Type	Description	Default
`scoring_rule`	`ScoringRule`		required
`x_traj`	`(T+1, dim_x) ndarray`		required
`u_traj`	`(T+1, dim_u) ndarray`		required
`abg_params`	`optional`	Reserved for a future adaptive DRPP (currently unused).	`None`

Returns:

Type	Description
`dict`	`{'x_pred', 'j_pred', 'Sigma_pred', 'score'}`.

eigvec_fixed_minmax ¶

eigvec_fixed_minmax(lamb, gamma, x, u)

Solve the "eigenvectors fixed" DRPP cone program (equation 19 of the paper).

This is the central optimization of the DRPP framework. Given the eigenvalues lamb of the nominal covariance and the uncertainty parameters gamma, it finds the worst-case optimal eigenvalues of the predictive covariance.

The functional maximin over probability measures can be converted equivalently into a Euclidean cone program, solved with CVXPY + SCS.

Parameters:

Name	Type	Description	Default
`lamb`	`(dim_x,) ndarray`	The eigenvalues of the nominal covariance `Sigma_bar`.	required
`gamma`	`list of length 3`	The uncertainty parameters: `gamma[0]`: a callable `f(x, u)` returning the squared upper bound of the drift perturbation; `gamma[1]`: float, the bound on the noise-mean uncertainty; `gamma[2]`: float, the covariance scaling factor.	required
`x`	`(dim_x,) ndarray`	The current state.	required
`u`	`(dim_u,) ndarray`	The current control input.	required

Returns:

Name	Type	Description
`optimal_lamb`	`(dim_x,) ndarray`	The worst-case optimal eigenvalues of the predictive covariance.
`optimal_value`	`float`	The optimal objective value (the lower bound of the expected log score).
`optimal_j`	`int`	The index of the selected "worst" eigen-direction.

Notes

⚠️ The SCS parameters (max_iters=5000, eps=1e-9) were tuned for cone-program numerical stability, and changing them affects the discrete choice of j: a perturbation of 1e-14 between nearly tied candidates is enough to flip j (after the flip the score differs by only ~1e-9, but the result is then no longer bit-for-bit identical). Do not change them in passing.

The optimization sweeps over j (each possible worst direction) and takes the minimum.

Predictors whose law is a pmf you supply:

atom ¶

Predictors whose predictive law is a pmf you supply, over a set of atoms.

The second category of probabilistic predictor beside drpp/. Where a DRPP predictor derives a Gaussian from a nominal model and an ambiguity set, these take an :class:~ppsds.numerics.atom_pmf.AtomPmf as given and predict with it.

Why this is useful even without an optimizer

Because it is the piece that lets a prediction you chose by hand be evaluated as rigorously as one a model produced: same environment, same scoring rule, same statistics. That makes the interactive bet-placing picture and an experiment the same object, rather than two demos that happen to look alike.

It is deliberately not a search over pmfs. Nothing here picks the atoms or the weights; it evaluates the ones it is handed. Two ways to hand them over:

:class:StaticAtomPmfPredictor -- the atoms are in the state frame and never move. This is the shape of "here is what I think the next state is", one step at a time.
:class:DriftShiftedAtomPmfPredictor -- the atoms are in the noise frame, and each step shifts them by the nominal drift. This is the shape of a predictor that keeps up with a trajectory, and it is the same reduction the predictability paper uses for an i.i.d. process.

StaticAtomPmfPredictor ¶

StaticAtomPmfPredictor(atom_pmf)

Bases: _AtomPmfPredictor

Predicts the same pmf every step: the atoms you place are where you say they are.

This is the shape of a one-step prediction made by hand -- "the next state is near here, with these weights" -- repeated. It is also the shape a bet-placing game has, which is what makes the two comparable.

DriftShiftedAtomPmfPredictor ¶

DriftShiftedAtomPmfPredictor(atom_pmf, f_bar)

Bases: _AtomPmfPredictor

Shifts a noise-frame pmf by the nominal drift every step.

Predicts p_w(x_k - f(x_{k-1}, u_{k-1})): the atoms describe the noise, and moving them by the drift is what turns them into a prediction about the state. A predictor of this shape is what the predictability paper's convergence result assumes.

Parameters:

Name	Type	Description	Default
`atom_pmf`	`AtomPmf`	The pmf over the noise, in the drift's frame.	required
`f_bar`	`callable`	The nominal drift `f_bar(x, u) -> (x_dim,)`.	required

The contract both satisfy. Three methods: predict() returns the predictive distribution for the current step, update(x, u) advances the internal state after observing the state and control, and reset() restores it for reuse — clearing the outputs alone is not enough, because a predictor rewrites the nominal parameters in place. Three attributes must exist on every implementation — x_pred, Sigma_pred, lamb_pred — because the environment layer's adversarial mode reads them to pick the worst-case perturbation direction.

base ¶

Interface contract for probabilistic predictors (typing.Protocol).

Every DRPP predictor follows a stateful predict() -> update() cycle:

predict() returns the predictive distribution of :math:x_{k+1} given the current internal state;
update(x_observe, u_observe) advances the internal state once the true state and control have been observed.

The statefulness exists so that expensive computations (such as the cone program of Eig-DRPP) can be cached from step to step.

Why Protocol and not ABC

This module used to be an ABC plus ABCMeta.register() virtual subclasses, but that design had a trap:

The contract declared lamb_pred / x_pred / Sigma_pred as @property. @property is a data descriptor, whereas the concrete predictors assign them as plain instance attributes (self.x_pred = ...). The two are incompatible -- to actually inherit it would raise AttributeError at construction time.
Registering virtually through register() is precisely what hid that incompatibility: isinstance() returned True, while not one clause of the contract was actually enforced.

With Protocol: attributes are declared by annotation (a plain instance attribute satisfies them), methods are declared with ..., there is zero runtime overhead, and register() is no longer needed. Use mypy when you want static checking.

Protocol cannot provide default implementations, so implementers must supply their own reset().

Ambiguity sets¶

conic_moment ¶

Ambiguity sets based on conic moments.

This module covers only one branch of the model assumptions::

x_{k+1} = f(x_k, u_k) + delta_f + w_k

that is

additive drift: the drift error delta_f enters the state added to the nominal drift f;
additive noise: the noise w likewise enters the state added to it, and does not change with the state (homogeneous);
conic moment uncertainty: the norm of delta_f, the noise mean mu and the noise covariance Sigma each lie in a set parameterized by gamma.

.. note:: The ambiguity set itself does not assume Gaussianity. The conic moment constraints use only the first two moments of the noise. Gaussianity appears only in two downstream places: (a) the generative kernel used for sampling (:meth:~AdditiveDriftConicMomentSet.transition), and (b) the logarithmic scoring rule. Do not read "moment constraints" as "distributional assumptions".

.. note:: This module does not cover: state-dependent noise, multiplicative noise, w correlated with x, or drift errors that enter the state nonlinearly. Those require another base class -- and another predictor (the existing Eig-DRPP predictor requires the three things f / mu / Sigma).

Class hierarchy

:class:AdditiveDriftConicMomentSet is the base class, carrying the parts common to both families (the sampling-trajectory skeleton, the theoretical bounds, the selection of the worst-case direction, transition/validate); the two subclasses differ only in parameterization -- that is, in what sds_params looks like:

================================= ============================= ============================== Subclass Dynamics sds_params keys ================================= ============================= ============================== :class:LinearAmbiguitySet x' = A x + B u + w A / B / mu / Sigma :class:NonlinearAmbiguitySet x' = f(x,u) + delta_f + w delta_f / mu / Sigma ================================= ============================= ==============================

The difference between the two is only one of representation: the linear family writes the drift error as a matrix perturbation ((A_bar+dA)x + (B_bar+dB)u, which expands to f_bar(x,u) + delta_f) and folds it into the noise mean in the worst case; the nonlinear family gives the delta_f vector directly. Given the same gamma and the same prediction covariance, the two update to the same state.

On the class names: this module's path ambiguity/additive_drift/conic_moment already says it all -- "assumption (additive drift) + family (conic moment)" -- so the subclass names need not repeat it: they say only their own business (linear / nonlinear). The base class name AdditiveDriftConicMomentSet keeps the full name, because it is the reference point when adding a family (see below).

Adding another family

The landing spot is a sibling of this module: ambiguity/additive_drift/<family>.py. The additive_drift level is the assumption; conic_moment is one family under it, and a different family is a new module beside this one rather than a branch inside it.

Do not subclass :class:AdditiveDriftConicMomentSet -- the representation of the radius and the worst-case rule both differ. Only part of the base class is family-agnostic:

=========================================== ========================================== Family-agnostic conic moment specific =========================================== ========================================== x_dim / u_dim / f_bar / gamma (three-element radius) mu_bar / Sigma_bar param_keys / validate worst_case_params (worst direction + covariance shrink) transition (additive noise kernel) value_upper_bound_true_transition / sample_trajectory (by _traj_fields) value_lower_bound_eigvec_fixed mean() / sample() abstract interface =========================================== ==========================================

When a second family actually lands, then extract the left column into a shared base class; right now there is only one family, so extracting would be guessing.

AdditiveDriftConicMomentSet ¶

AdditiveDriftConicMomentSet(x_dim, u_dim, mu_bar, Sigma_bar, *, f_err_callable, gamma1, gamma2, f_err_max)

Base class for conic moment ambiguity sets: assumes [additive drift + additive noise].

"Additive" is not just an adjective in the name -- it is enforced by the base class: :meth:transition implements x' = mean(x, u; params) + w -- the drift is supplied by the subclass's :meth:mean, and the noise enters the state added to it. Subclasses cannot get around this.

Attributes:

Name	Type	Description
`x_dim, u_dim`	`int`	State / control dimensions. `dim_x` / `dim_u` are compatible aliases for them.
`mu_bar`	`ndarray of shape (x_dim,)`	Nominal noise mean.
`Sigma_bar`	`ndarray of shape (x_dim, x_dim)`	Nominal noise covariance.
`gamma1`	`float`	The squared radius of the noise mean: `\|mu - mu_bar\| <= sqrt(gamma1)`.
`gamma2`	`float`	Covariance scaling factor: `Sigma = gamma2 * Sigma_bar - ...`.
`f_err_max`	`float`	The un-squared hard upper bound on the additive drift norm `\|\|delta_f\|\|`.

dim_x `property` ¶

dim_x

Alias for x_dim.

dim_u `property` ¶

dim_u

Alias for u_dim.

gamma `property` ¶

gamma

The canonical format consumed by the predictor: [f_err, gamma1, gamma2].

The meaning of the three elements (all of them concern "squared", the place most easily filled in wrong):

gamma[0] -- a callable f(x, u) returning an upper bound on the squared norm of the drift ||delta_f||^2. Note this differs from f_err_max: the latter is the un-squared hard upper bound.
gamma[1] -- the squared radius of the noise mean.
gamma[2] -- the covariance scaling factor.

This list is consumed directly by the DRPP predictor (drpp_eigvec_fixed_minmax) and by the adversarial dynamics, and it is the only shape that must be aligned with them.

mean ¶

mean(x, u, params)

Nominal drift + per-parameter perturbation; does not include the noise mean mu.

sample ¶

sample()

Sample one SDS realization from the ambiguity set, returning a dict shaped like sds_params.

validate ¶

validate(params)

Validate that every key of sds_params is accepted by this ambiguity set.

A missing key is allowed -- it means that term takes its nominal value (on the first step of adversarial mode there is no prediction yet, so {} is returned). Only an unknown key raises: that is precisely where the old "guess the format via 'A' in params" used to swallow the error.

transition ¶

transition(x, u, params, rng)

Take one step: x' ~ N(mean(x,u;params) + mu, Sigma).

This line is the executable form of the "additive" assumption -- the drift and the noise enter the state added together, rather than through a state-dependent transform. Choosing a Gaussian kernel for the noise is a sampling-side choice (and also the kernel of the logarithmic scoring rule), not a distributional assumption of the ambiguity set.

rng is supplied by the caller (gym.Env.np_random) so that results are reproducible; when mu / Sigma are missing the nominal values are used.

sample_trajectory ¶

sample_trajectory(T, period=1)

Sample one piecewise-constant trajectory of SDS parameters.

Parameters:

Name	Type	Description	Default
`T`	`int`	Number of time steps.	required
`period`	`int`	How many steps each SDS realization persists before resampling: `1` means step-by-step i.i.d., `100` means slowly varying.	`1`

Returns:

Type	Description
`dict`	Keys are declared by :attr:`_traj_fields`; the leading dimension of each value is `T+1`.

value_upper_bound_true_transition ¶

value_upper_bound_true_transition()

Upper bound: the tight value of the expected log score when the true transition is known.

It is the closed form of the Gaussian log score -(d·log2π + d + log det(Σ))/2 evaluated when the prediction covariance equals the noise covariance (the largest one in the set, gamma[2]*Sigma_bar). It is a tight upper bound: no predictor can score higher than "true transition known". It is a pure property of the set and involves no relaxation.

value_lower_bound_eigvec_fixed ¶

value_lower_bound_eigvec_fixed()

Lower bound: the expected log score that Eig-DRPP guarantees at x=0, u=0.

Opposite to the upper bound: the upper bound assumes the true transition is known, whereas this lower bound assumes the worst case has already occurred and that the predictor may only use the eigenvectors fixed relaxation (fixed to the eigenvectors of Sigma_bar) -- that is, the design space of Eig-DRPP. It takes the optimal value of the minmax under that relaxation, i.e. "how much one can at least obtain in the worst case".

.. note:: This method calls :func:~ppsds.probabilistic_predictor.drpp.additive_drift.conic_moment.eigvec_fixed_minmax from probabilistic_predictor.drpp.additive_drift.conic_moment, so ambiguity depends on probabilistic_predictor. This is not a typo: the lower bound is defined by the Eig-DRPP minmax, and without the predictor it has no meaning.

worst_case_params ¶

worst_case_params(x, u, lamb_pred)

For a given predictor, construct the SDS parameters that maximize its score loss.

First the worst eigendirection j is selected, then the covariance is shrunk, and finally the subclass's :meth:_assemble_worst assembles this family's sds_params.

Parameters:

Name	Type	Description	Default
`x`	`ndarray of shape (x_dim,)`	Current state.	required
`u`	`ndarray of shape (u_dim,)`	Current control input.	required
`lamb_pred`	`ndarray of shape (x_dim,)`	Eigenvalues of the predictor's current prediction covariance.	required

Returns:

Type	Description
`dict`	This family's `sds_params`, plus the key `'j'` (index of the worst direction).

LinearAmbiguitySet ¶

LinearAmbiguitySet(A_bar, B_bar, mu_bar, Sigma_bar, alphas, f_err_max)

Bases: AdditiveDriftConicMomentSet

Linear time-invariant SDS: x_{k+1} = A x_k + B u_k + w_k.

The uncertainty is on (A, B, mu, Sigma), parameterized by alphas and f_err_max. The drift error still enters additively -- here it is written as a matrix perturbation delta_f = dA x + dB u, which expands to f_bar(x, u) + delta_f.

Parameters:

Name	Type	Description	Default
`A_bar`	`ndarray of shape (x_dim, x_dim)`	Nominal dynamics matrix.	required
`B_bar`	`ndarray of shape (x_dim, u_dim)`	Nominal input matrix.	required
`mu_bar`	`ndarray of shape (x_dim,)`	Nominal noise mean.	required
`Sigma_bar`	`ndarray of shape (x_dim, x_dim)`	Nominal noise covariance.	required
`alphas`	`dict`	The uncertainty radius of each component, with keys `'A'` / `'B'` / `'mu'` / `'Sigma'`: `alphas['A']`: the maximum perturbation of each superdiagonal element of `A` (un-squared magnitude; the sampler draws uniformly on `[-alphas['A'], alphas['A']]`). `alphas['B']`: the maximum perturbation of `B[0,-1]`, same convention. `alphas['mu']`: an upper bound on `\|\|mu_a\|\|^2` (the squared norm of the noise-mean deviation), i.e. `sum_i (mu_i - mu_bar_i)^2 <= alphas['mu']`. Hence the per-component box half-width is `sqrt(alphas['mu'] / x_dim)` -- that is what :meth:`sample` uses, so that `max\|\|mu_a\|\|^2` exactly exhausts the declared bound. `alphas['Sigma']`: the scaling factor of the nominal covariance.	required
`f_err_max`	`float`	The hard upper bound on the perturbation norm `\|\|delta_f\|\|`, un-squared. :meth:`f_err` compares it against the un-squared `norm_max` and returns the square of the smaller, so its return value is a squared quantity whose upper bound is `f_err_max**2`.	required

f_bar ¶

f_bar(x, u)

Nominal dynamics: A_bar @ x + B_bar @ u.

f_err ¶

f_err(x, u)

An upper bound on the squared additive drift norm, clipped by f_err_max.

Taking a_A||x|| + a_B||u|| (a sum) rather than the combined norm a * sqrt(||x||^2 + ||u||^2) is deliberate: the drift produced by the sampler is delta_f = dA x + dB u, and by the triangle inequality its tight upper bound is exactly this sum (||dA x|| <= a_A||x||, ||dB u|| <= a_B||u||, see :meth:sample). The combined-norm expression is always <= this one, so using it would place some samples outside the declared ambiguity set -- i.e. the set could not contain its own generator. "It contains it" is the criterion.

mean ¶

mean(x, u, params)

When A / B are given, replace the nominal matrices wholesale -- this is the existing semantics of the linear family's parameterization.

sample ¶

sample()

Sample one SDS realization from the ambiguity set.

A: add a U[-alphas['A'], alphas['A']] perturbation to the superdiagonal;
B: add a U[-alphas['B'], alphas['B']] perturbation to B[0,-1];
mu: per-component box half-width sqrt(alphas['mu'] / x_dim), so that max||mu_a||^2 exactly exhausts alphas['mu'] (the paper takes mu = 0.5 * alpha_3 with alpha_3 ~ U[-1,1]^d, whose half-width is exactly this value);
Sigma: alphas['Sigma'] * Sigma_bar - mu_a mu_a^T.

NonlinearAmbiguitySet ¶

NonlinearAmbiguitySet(f_bar, mu_bar, Sigma_bar, gamma, legacy_sampler=False)

Bases: AdditiveDriftConicMomentSet

Nonlinear SDS: x_{k+1} = f(x_k, u_k) + delta_f + w_k.

This is the most direct form of "additive drift" -- delta_f is added straight to the nominal drift. The NLink in the name comes from this repository's original use case (an n-link manipulator), but the class itself requires only that f_bar be a callable and is not tied to any particular system.

Parameters:

Name	Type	Description	Default
`f_bar`	`dict`	Must contain `'function'` (a callable `f(x, u)` returning the nominal next state), `'x_dim'`, and `'u_dim'`.	required
`mu_bar`	`ndarray of shape (x_dim,)`	Nominal noise mean.	required
`Sigma_bar`	`ndarray of shape (x_dim, x_dim)`	Nominal noise covariance.	required
`gamma`	`list of length 3`	`gamma[0]`: a dict containing `'function'` (a callable returning an upper bound on the squared norm of the drift at `(x,u)`), `'max'` (the hard upper bound on the un-squared norm `\|\|delta_f\|\|`), `'x_dim'`, and `'u_dim'`; `gamma[1]`: the squared radius of the noise mean (`\|mu - mu_bar\| <= sqrt(gamma[1])`); `gamma[2]`: the covariance scaling factor.	required
`legacy_sampler`	`bool`	Reproduce the old generator -- `delta_f ~ N(0, (f_err_max/3)^2 I)`. When `x_dim` is large that sampler violates its own declared bound (`/3` is a per-component 3-sigma rule of thumb that ignores the dimension: `E\|\|delta_f\|\|^2 = x_dim * (f_err_max/3)^2`, which already exceeds `f_err_max^2` once `x_dim > 9`). Kept only to reproduce old results.	`False`

mean ¶

mean(x, u, params)

The nominal drift f_bar(x, u) plus the additive delta_f (zero by default).

sample ¶

sample()

Sample one SDS realization from the ambiguity set.

mu uses rejection sampling: the construction Sigma = gamma[2]*Sigma_bar - mu_a mu_a^T may be non-positive-definite when mu approaches the box boundary, and the ambiguity set contains only valid covariances, so invalid (mu, Sigma) pairs are discarded and redrawn.

Scoring rules¶

logarithm_score ¶

Logarithmic scoring rule: score a predictive distribution against the true state.

score(predictive_dist, x_true) returns logpdf(x_true) — higher is better. DRPP's theoretical bounds (value_upper_bound_true_transition / value_lower_bound_eigvec_fixed) and the comparisons of predictors in the experiments all use it.

For multivariate_rv only logpdf is required (duck typing); it need not be a scipy distribution object.

LogarithmScore ¶

The standard logarithmic scoring rule.

score ¶

score(multivariate_rv, x)

Returns logpdf(x).

epsilon_logarithm_score ¶

The epsilon-logarithm score: the logarithm of the predictive probability near the outcome.

Given a neighbourhood radius eps, score a predictive law p on a realised outcome x by

.. math:: \mathcal L_\epsilon(p, x) = \log \int_{|s - x|_\infty \le \epsilon} p(s)\, ds ,

which is the ordinary logarithm score when eps = 0. Taking a neighbourhood is what lets the score be evaluated on a pmf, whose support need not contain the truth's -- the ordinary logarithm score is -inf for such a predictor almost surely, the epsilon version is not. \|\cdot\|_\infty is the norm the predictability paper uses, chosen for the convenience of integrating over a box.

Which path is taken

:meth:EpsilonLogarithmScore.score dispatches on the law it is handed, cheapest first:

eps == 0 -- logpdf, the ordinary logarithm score;
a law exposing box_integral(x, eps) -- the pmf path, e.g. :class:~ppsds.numerics.random_variables.AtomPmf;
a Gaussian with diagonal covariance -- the closed form prod_i [Phi((x_i + eps - mu_i) / sigma_i) - Phi((x_i - eps - mu_i) / sigma_i)], which is what makes a 10**5-trajectory sweep of a Gaussian predictor cheap;
anything else with a cdf -- inclusion-exclusion over the 2 ** d_x box corners. This is the general fallback and costs exponentially in d_x, so it suits small d_x.

Input is a single outcome (d_x,) or a batch (n, d_x), returning a scalar or (n,). A zero box probability gives -inf, silently: that is the score for a predictor that puts no mass anywhere near what happened.

.. warning:: The score is discontinuous at eps = 0, and this is a property of the definition rather than a numerical artefact. At exactly eps = 0 the score is log p(x); for eps > 0 the box integral is about (2 eps) ** d_x * p(x), so the score is about log p(x) + d_x log(2 eps) and rushes to -inf as eps -> 0+. Said the other way round, L_eps - d_x log(2 eps) -> log p(x): the ordinary logarithm score is recovered only after subtracting the box's log volume, and the I_{eps > 0} indicator in the predictability formulas is what encodes this. Consequences worth knowing: a pmf predictor scores -inf at eps = 0 for any finite eps too, only less so; and a slider sweeping eps down to zero shows a cliff, not a landing.

EpsilonLogarithmScore ¶

EpsilonLogarithmScore(eps)

The epsilon-logarithm scoring rule.

Parameters:

Name	Type	Description	Default
`eps`	`float`	Neighbourhood radius, `>= 0`. Zero degenerates to the ordinary logarithm score.	required

score ¶

score(multivariate_rv, x)

Score multivariate_rv on the realised outcome(s) x; higher is better.

Parameters:

Name	Type	Description	Default
`multivariate_rv`	`object`	The predictive law. Any of: one exposing `box_integral`; a Gaussian exposing `mean` / `cov`; or anything exposing `cdf`. With `eps == 0` it must expose `logpdf`.	required
`x`	`array_like`	A single outcome `(d_x,)` or a batch `(n, d_x)`.	required

Returns:

Type	Description
`float or ndarray`	The score, or the vector of scores for a batch.

bin_to_list ¶

bin_to_list(x, d)

Bit decomposition of x into d most-significant-first bits.

Predictability: the ceiling¶

bounds ¶

Two reference values for the expected epsilon-logarithm score.

Given a true one-step law q (written here in the state frame; for an additive-drift system it is the noise law translated by the drift, so it is the same shape every step), the epsilon-logarithm score of a predictive law p on a realised outcome x is

.. math:: \mathcal L_\epsilon(p, x) = \log \int_{|s - x|_\infty \le \epsilon} p(s)\, ds ,

and this module provides the two horizontal lines the predictability literature plots against it:

:func:predictability_upper_bound -- :math:d_x \log(2\epsilon) - h_d(q). An upper bound on the expected score of any predictor, so it is a property of the system alone: no predictor enters. This is the "probabilistic predictability" of the system.
:func:approximation_value -- :math:d_x \log(2\epsilon) - \mathrm H(q \| p). The O(epsilon) approximation of the expected score of the given predictor p.

.. note:: The upper bound is a bound, not a target: it is attained only by some systems, so a gap between a predictor's score and this line is not by itself evidence that the predictor is poor. Where the bound is attained, the gap closes from below.

Both formulas carry the same indicator: the d_x log(2 epsilon) volume term is present only for epsilon > 0. At epsilon = 0 the score degenerates to the ordinary logarithm score, the volume term vanishes, and the two lines become -h_d(q) and -H(q || p) respectively.

Where the predictability upper bound sits relative to DRPP

These two answer "how well could anything do, for a system known exactly" -- an upper bound on prediction quality. DRPP answers "what can be guaranteed for a system known only through an ambiguity set" -- a lower bound. They are complementary, and they are a different pair of bounds from value_upper_bound_true_transition / value_lower_bound_eigvec_fixed, which bracket the relaxation gap inside one ambiguity set.

predictability_upper_bound ¶

predictability_upper_bound(law_q, eps, d_x=None)

:math:d_x \log(2\epsilon) - h_d(q) -- the best expected score achievable, in the bound.

Parameters:

Name	Type	Description	Default
`law_q`	`object`	The true continuous one-step law; needs `entropy()`.	required
`eps`	`float`	Neighbourhood radius, `>= 0`. Zero degenerates to the ordinary logarithm score.	required
`d_x`	`int`	State dimension. Inferred from `law_q` when omitted, which is correct whenever the process noise spans the state -- true for every additive-drift system here. Pass it explicitly if the state is higher-dimensional than the noise.	`None`

Returns:

Type	Description
`float`	The bound, in nats.

approximation_value ¶

approximation_value(law_q, law_p, eps, d_x=None, *, bounds=None, n_quad=None)

:math:d_x \log(2\epsilon) - \mathrm H(q \| p) -- the O(epsilon) score approximation.

The true expected score differs from this by O(epsilon), so it is the line a continuous predictor's score converges to. It is not available for an atom/pmf predictor: :func:~ppsds.predictability.entropy.cross_entropy raises, because such a predictor has log p = -inf almost everywhere.

Parameters:

Name	Type	Description	Default
`law_q`	`object`	True and predictive laws. See :func:`~ppsds.predictability.entropy.cross_entropy`.	required
`law_p`	`object`	True and predictive laws. See :func:`~ppsds.predictability.entropy.cross_entropy`.	required
`eps`	`float`	Neighbourhood radius, `>= 0`.	required
`d_x`	`int`	State dimension; inferred from `law_q` when omitted.	`None`
`bounds`	`optional`	Passed through to :func:`~ppsds.predictability.entropy.cross_entropy`.	`None`
`n_quad`	`optional`	Passed through to :func:`~ppsds.predictability.entropy.cross_entropy`.	`None`

Returns:

Type	Description
`float`	The approximation, in nats.

entropy ¶

Differential entropy of the true law, and its cross-entropy against a predictive law.

Two per-step quantities, both expectations under the true one-step law q:

:func:differential_entropy -- :math:h_d(q);
:func:cross_entropy -- :math:\mathrm H(q \,\|\, p) = -\mathbb E_{x \sim q} \log p(x).

Why this module exists

The predictability upper bound is d_x log(2 eps) - h_d(q), and the O(eps) approximation of a given predictor's expected score is d_x log(2 eps) - H(q || p). Both need an entropy of a law, not of a sample, so both live here rather than being re-derived at each call site.

A "law" is duck-typed: :func:cross_entropy needs .pdf, :func:differential_entropy needs .entropy(). Both :class:~ppsds.numerics.random_variables.multivariate_uniform and the scipy.stats multivariate distributions qualify.

.. warning:: q here is the true continuous one-step law. The bound formula above does not apply when q is itself an atom/pmf law: differential_entropy would then return its Shannon entropy, which is a different quantity (non-negative, where a differential entropy may be negative). Predict with a pmf; take the entropy of the continuous truth.

differential_entropy ¶

differential_entropy(law_q)

Return :math:h_d(q), the differential entropy of the true one-step law.

Delegates to the law's own entropy() -- for :class:~ppsds.numerics.random_variables.multivariate_uniform that is the closed form log(prod(high - low)); for a scipy.stats distribution it is scipy's.

Parameters:

Name	Type	Description	Default
`law_q`	`object`	The true continuous one-step law; must expose `entropy()`.	required

Returns:

Type	Description
`float`	The entropy. May be negative -- a differential entropy is a log of a volume, and the volume can be below 1.

cross_entropy ¶

cross_entropy(law_q, law_p, *, bounds=None, n_quad=None)

Return :math:\mathrm H(q \,\|\, p) = -\mathbb E_{x \sim q} \log p(x).

The expectation is integrated over q's support by tensor-product trapezoid rule -- adequate for d_x <= 3, which is what the predictability work needs. For a box law (.low / .high, as :class:~ppsds.numerics.random_variables.multivariate_uniform has) the support is taken from the law; otherwise pass bounds=(low, high).

Parameters:

Name	Type	Description	Default
`law_q`	`object`	The true law; must expose `pdf`, and `low` / `high` unless `bounds` is given.	required
`law_p`	`object`	The predictive law; must expose `pdf`. A pmf/atom predictor does not -- see the `TypeError` below, which is deliberate.	required
`bounds`	`tuple of array_like`	`(low, high)`, each of length `d_x`. Defaults to the box `q` declares.	`None`
`n_quad`	`int`	Grid points per axis. Defaults to 2001 in 1-D, 201 in 2-D, 41 above -- roughly 4e4 quadrature points either way.	`None`

Returns:

Type	Description
`float`	The cross-entropy, in nats.

Raises:

Type Description

TypeError

If law_p has no pdf. This is the pmf case, and it is worth being loud about: an atom law assigns zero density almost everywhere, so log p(x) = -inf and the cross-entropy is +inf. That is exactly why the epsilon-logarithm score exists -- score a pmf with :class:~ppsds.scoring_rule.epsilon_logarithm_score.EpsilonLogarithmScore, and use :func:~ppsds.predictability.bounds.predictability_upper_bound as the reference line rather than the O(eps) approximation.

ValueError

If q puts mass where p is zero (the cross-entropy is then genuinely +inf).

trajectory ¶

Aggregating per-step scores into the curves the predictability figures plot.

A rollout produces a score per step per trajectory, so a (n_traj, T) array. Two shapes get drawn from it:

:func:cumulative_average -- the running mean along the step axis, giving one curve per trajectory. This is the curve that converges to the expected score, at rate O_p(T ** -0.5) when the process noise is i.i.d. -- which is why a single trajectory is enough to estimate the expectation, and why these figures do not need repeated samplings.
:func:across_trajectories -- per-step mean and quantiles over trajectories, giving the expected-score curve and its confidence band.

Both are numpy-only and live in the core. :func:cumulative_average is the same operation as ppsds.visualization.data_transform.cum_average, which now delegates here so that the core and the plotting periphery cannot drift apart.

cumulative_average ¶

cumulative_average(scores)

Running mean along the step axis of a (n_traj, T) score block.

Entry [n, k] is the mean of scores[n, 0:k+1]. The first step is taken as-is, which matters for a one-step prediction: there the score is the average.

across_trajectories ¶

across_trajectories(scores, quantiles=(0.025, 0.975))

Per-step mean and quantiles over trajectories.

Parameters:

Name	Type	Description	Default
`scores`	`array_like`	`(n_traj, T)` -- typically already passed through :func:`cumulative_average`.	required
`quantiles`	`sequence of float`	Levels to compute per step. The default gives the paper's 95% band.	`(0.025, 0.975)`

Returns:

Type	Description
`dict`	`{'mean': (T,), 'quantile_levels': tuple, 'low': (T,), 'high': (T,)}`. With a single trajectory the quantiles equal the mean, so the band collapses rather than raising -- the figures stay drawable while a rollout is still short.

The system itself¶

known_sds ¶

A stochastic dynamical system that is known exactly: a drift plus a known noise law.

The counterpart to :mod:~ppsds.ambiguity.additive_drift.conic_moment, which describes a system known only through an ambiguity set. Here there is no uncertainty about the model at all -- x' = f(x, u) + w with f given and w drawn from a law you name -- so the noise can be non-Gaussian, which the ambiguity sets' Gaussian transition cannot express.

This is the object behind the predictability questions ("how well could anything predict this system?"), and the truthful half of a simulation: an ambiguity set is what a predictor believes, while this is what actually generated the data.

It belongs in dynamics/ rather than ambiguity/ because it is the thing being predicted, not a belief about it: ambiguity/ is reserved for sets of possible models.

AdditiveDriftSds ¶

AdditiveDriftSds(f_bar, noise_law, x0=None)

x' = f(x, u) + w with f known and w from a named law.

Parameters:

Name	Type	Description	Default
`f_bar`	`callable`	The drift `f_bar(x, u) -> (x_dim,)`.	required
`noise_law`	`object`	The process noise. Must expose `rvs`, taking `size` and a keyword for the random source -- either `rng` (as :class:`~ppsds.numerics.random_variables.multivariate_uniform` and :class:`~ppsds.numerics.atom_pmf.AtomPmf` do) or `random_state` (as `scipy.stats` distributions do). Which one is detected once, at construction.	required
`x0`	`array_like`	The default initial state. Defaults to the origin of the noise's dimension.	`None`

Raises:

Type	Description
`ValueError`	If `x_dim` cannot be determined from `x0` or the noise law.

transition ¶

transition(x, u, rng=None)

One step: the drift, plus a fresh noise draw. (x_dim,).

Parameters:

Name	Type	Description	Default
`x`	`array_like`	State and control.	required
`u`	`array_like`	State and control.	required
`rng`	`Generator`	Source of randomness for the noise draw; `None` uses the global `np.random` state. Pass one for a reproducible rollout.	`None`

sample_trajectory ¶

sample_trajectory(T, rng=None, x0=None, u=None)

Roll out T steps, returning both the states and the noise that drove them.

The noise is returned alongside because it is often the quantity of interest: when the dynamics are known to the predictor, the one-step prediction reduces to a prediction of the noise, so scoring happens on w rather than on x.

Parameters:

Name	Type	Description	Default
`T`	`int`	Number of steps. The returned states cover `T + 1` time points, the noise `T`.	required
`rng`	`Generator`	Source of randomness; `None` uses the global state.	`None`
`x0`	`array_like`	Initial state; defaults to the one given at construction.	`None`
`u`	`array_like`	The control. Either constant, shaped `(u_dim,)`, or per step, `(T, u_dim)`. Defaults to the empty vector, which is what an uncontrolled system wants; a system whose drift does `B @ u` will then fail loudly on the shape rather than quietly simulate a control nobody asked for.	`None`

Returns:

Type	Description
`dict`	`{'x': (T + 1, x_dim), 'w': (T, x_dim)}`.

identification ¶

Identify a linear nominal model from rollouts, and diagnose the residual.

This is the first step of "connecting DRPP to your own environment": you usually cannot obtain the system's analytical model, but you can run a few trajectories and fit it approximately to

.. math::

x_{k+1} \approx A x_k + B u_k + c + w_k

This module does only least squares + diagnostics, not rigorous system identification — the goal is not to find the true values, but to give you a reasonable nominal model plus a set of residual statistics on which to base a tolerance declaration.

On "what noise actually is"

A model is merely a human's understanding of some physical phenomenon, and all models are wrong. For this repository's additive SDS:

:math:f is the deterministic law you believe you understand;
:math:w is everything else you have not yet understood.

Both can be wrong — that is just modeling. So the residual this module computes need not "look like noise" to deserve being called noise; whether it is, and how much uncertainty to assign it, is yours to declare.

.. warning:: There is a practical hurdle (not a theoretical one): DRPP evaluates predictions with a logarithmic score, and the logarithmic score requires the predictive distribution to have a density. If you use the residual covariance as is for :math:\Sigma, it is often numerically singular —

some state components may be identically linear (e.g. the exact Euler step x' = x + dt*xdot);
observations may lie on a low-dimensional manifold (e.g. the robotic-arm observation (cos θ, sin θ, θdot) satisfies cos²+sin²=1, a 2D cylinder in 3D space).

The support of a degenerate distribution is a null set, so any value even slightly off it scores -inf (measured: at a relative distance of only 1.4e-6 from the support, all 50 steps still gave -inf).

So it must be regularized to positive definite. That is not a stopgap — even dimensions that look deterministic can be assigned uncertainty, and that is precisely part of modeling. The magnitude is yours to choose, and it directly determines the magnitude of the score (see :func:regularize).

LinearFit `dataclass` ¶

LinearFit(A: ndarray, B: ndarray, c: ndarray)

The least-squares solution of x' ≈ A x + B u + c.

f_bar ¶

f_bar(x, u) -> np.ndarray

The nominal drift — usable directly as a DRPP predictor's f_bar.

predict ¶

predict(X: ndarray, U: ndarray) -> np.ndarray

Batch prediction: X_next ≈ X @ A.T + U @ B.T + c.

residuals ¶

residuals(X, U, X_next) -> np.ndarray

The residual w = x_next - f_bar(x, u), with the same shape as X.

encode_action ¶

encode_action(space, a) -> np.ndarray

Turn a policy's output into the action features used for regression.

Discrete(n) → one-hot. Note that the integer itself must not be used directly: the action indices are nominal (there is no ordering among 0/1/2), and treating them as ordinal numbers imposes a linear structure that does not exist, inflating the residual.
Box → flatten as is.

What is returned is the u inside f_bar(x, u) — so the predictor's update(x, u) must use the same encoding.

collect_rollouts ¶

collect_rollouts(env, policy, n_episodes, T, seed0=0, reset_seed=None)

Run several trajectories with a fixed policy; returns (X, U, X_next).

Parameters:

Name	Type	Description	Default
`env`	`Env`		required
`policy`	`callable`	`policy(obs) -> action`. The policy you use throughout. Its return value is passed as is to `env.step()` (a discrete action space needs an int, so do not convert it to an array first).	required
`n_episodes`	`int`	Number of trajectories.	required
`T`	`int`	Maximum number of steps per trajectory.	required
`seed0`	`int`	`reset(seed=seed0)` for the first trajectory; subsequent ones use +1 each (for reproducibility).	`0`
`reset_seed`	`int`	If given, every trajectory uses the same `seed` (for observing the variation from the same initial condition).	`None`

Returns:

Type	Description
`X, U, X_next : ndarray`	`U` is the action features encoded by :func:`encode_action` (one-hot for a discrete space), not the raw actions.

fit_linear ¶

fit_linear(X, U, X_next) -> LinearFit

Least-squares fit of x' ≈ A x + B u + c.

Uses lstsq rather than the normal equations: [X, U, 1] is often ill-conditioned (state components differ greatly in scale, and some components are identically linear), and the normal equations square that conditioning.

diagnose ¶

diagnose(residuals: ndarray, increments: ndarray) -> dict

Diagnose the residual — to inform modeling decisions, not as a pass/fail criterion.

Returns:

Type Description

dict

residual_ratio Mean residual norm / mean state-increment norm. How large counts as large depends on your application. lag1_autocorr Average of the per-dimension lag-1 autocorrelation. High means the residual still contains structure the model can explain (i.e. "the model is not there yet"); low means it is more like i.i.d. This is only a diagnostic — see the module docstring: whether the residual looks like noise does not affect whether it can be declared as noise. zero_directions Number of numerically zero eigenvalues. > 0 means the covariance is singular and must be regularized. eigenvalues / condition_number

regularize ¶

regularize(Sigma, floor: float) -> np.ndarray

Lift a singular covariance to positive definite: Sigma + floor * I.

This step is modeling, not a stopgap. See the module docstring: even dimensions that look deterministic can be assigned uncertainty. floor is exactly "how much uncertainty you are willing to admit in those directions".

⚠️ The magnitude of floor directly determines the magnitude of the score (because it determines the width of the predictive distribution), so it should not be an arbitrary number. Two reasonable sources:

Physics: perturbations you did not observe, actuator/sensor noise, unmodeled high-frequency dynamics;
Consistency: pick an initial value, run once, and see whether the resulting score magnitude matches your intuition about "how accurate the prediction is", then adjust.

A negative floor or an overly large floor makes the prediction meaningless, but there is no single correct answer — that is precisely what "modeling" means.

format_report ¶

format_report(diag: dict, fit: LinearFit | None = None) -> str

Lay out the result of :func:diagnose as a human-readable report.

mechanical ¶

Robotic-arm (n-link) dynamics models.

This module holds only physical models: the continuous-time equations and their discretization, with no uncertainty involved. The "nominal drift + additive perturbation" assumption and the ambiguity sets live in :mod:ppsds.ambiguity.additive_drift.conic_moment.

:class:TwoLink — two links, analytical form: the inertia matrix, Coriolis force, and gravity all have closed-form expressions, with a full set of physical parameters (mass, link length, center of mass, moment of inertia, damping).
:class:NLink — n links, numerically assembled. The stiffness matrix is omitted and the Coriolis force is approximated by empirical damping, so it is not an exact robotic-arm model; its purpose is to supply a nonlinear, strongly coupled nominal drift for DRPP to use. It also provides :meth:NLink.equilibrium_for_pose and :meth:NLink.linearize_at_discrete for linearizing around an equilibrium point.

Both provide f_continuous(x, u) and discretize(x, u, dt) (fourth-order Runge–Kutta).

NLink ¶

NLink(n=3, params=None)

A numerically assembled model of an n-link robotic arm (see the module docstring for the caveats).

Parameters:

Name	Type	Description	Default
`n`	`int`	Number of links. State dimension `2n` (the first n are joint angles, the last n are angular velocities), control dimension `n` (joint torques).	`3`
`params`	`dict`	Physical parameters shared by all links: `m` / `l` / `lc` / `I` / `d` / `g`.	`None`

equilibrium_for_pose ¶

equilibrium_for_pose(q_star=None, qd_star=None)

Returns the equilibrium state and the required torque at pose q_star (velocities default to 0).

linearize_at_discrete ¶

linearize_at_discrete(x_star, u_star, dt=0.02, eps=1e-06)

Numerically differentiate the discretized dynamics to get A / B at (x_star, u_star).

rk4_step ¶

rk4_step(f_cont, x, u, dt)

A single fourth-order Runge–Kutta step. f_cont(x, u) is the continuous-time derivative.

The environment layer¶

Serves the "simulate your own ambiguity" path, plus one exception: KnownSdsEnv steps a system that is known exactly. If you are applying a predictor to an environment you already have, none of this is needed.

base ¶

DRPP's Gymnasium environment (generic, not subclassed per system).

This subpackage serves scene 2 only -- "custom simulated ambiguity": you declare an ambiguity set by hand, then generate dynamics from it and simulate. Scene 1 (applying DRPP to your own gymnasium environment) does not need it, only :mod:ppsds.probabilistic_predictor + :mod:ppsds.scoring_rule.

Env is assembled from two things:

an ambiguity set (:class:~ppsds.ambiguity.additive_drift.conic_moment.AdditiveDriftConicMomentSet and its subclasses) -- it supplies the nominal model, the state/control dimensions, and "how to advance one step given sds_params";
a :class:~ppsds.environment.dynamics_mode.DynamicsMode -- it decides which set of sds_params each step uses.

The key set of sds_params is declared by the ambiguity set itself (param_keys) and checked at the step() entry point, which fails loudly on a mismatch -- rather than the old behaviour of guessing the format via 'A' in sds_params and silently taking a different branch when the guess was wrong.

Env ¶

Env(ambiguity_set, mode: DynamicsMode, x0: ndarray | None = None, max_steps: int = 50, bounds: tuple | None = None)

Bases: Env

Gymnasium environment for an SDS, assembled from an ambiguity set + a DynamicsMode.

Why it is called Env and not DRPPEnv: the qualifier "DRPP" carries no information (the whole package is DRPP) but does suggest adversarialism is intrinsic to it. In fact adversarialism is only one of four modes -- FixedDynamics / RandomDynamics / PreSampledDynamics are not adversarial at all. And ppsds.environment.Env already states where the class belongs.

Spaces

observation_space: Box(shape=(x_dim,), dtype=float64)
action_space: Box(shape=(u_dim,), dtype=float64)

step() returns

observation: next state :math:x_{k+1}
reward: always 0.0 (DRPP prediction quality is assessed externally)
terminated: True when state goes out of bounds
truncated: True after max_episode_steps
info: {'sds_params': {...}, 'mode': str}

Parameters:

Name	Type	Description	Default
`ambiguity_set`	`AdditiveDriftConicMomentSet`	Supplies the nominal model (`f_bar` / `mu_bar` / `Sigma_bar`), the dimensions, :meth:`transition` and :attr:`param_keys`.	required
`mode`	`DynamicsMode`	Strategy selecting `sds_params` each step.	required
`x0`	`ndarray of shape (x_dim,)`	Initial state. If `None`, defaults to zeros.	`None`
`max_steps`	`int`	Episode length (truncation horizon).	`50`
`bounds`	`tuple of (low, high) ndarrays`	State bounds for early termination. `None` means no bounds check.	`None`

ambiguity_set `property` ¶

ambiguity_set

The underlying ambiguity set object.

dynamics_mode `property` ¶

dynamics_mode: DynamicsMode

The active dynamics mode strategy.

reset ¶

reset(seed=None, options=None)

Start a new episode.

Parameters:

Name	Type	Description	Default
`seed`	`int`		`None`
`options`	`dict`	May contain `'x0'` to override the initial state.	`None`

Returns:

Type	Description
`(observation, info)`

step ¶

step(action)

Advance the simulation by one step.

Parameters:

Name	Type	Description	Default
`action`	`ndarray of shape (u_dim,)`	Control input.	required

Returns:

Type	Description
`(observation, reward, terminated, truncated, info)`

inject_eigvec_fixed_info ¶

inject_eigvec_fixed_info(lamb_pred, u, x=None) -> None

Feed the predictor's predicted eigenvalues into the current mode (for the "eigenvectors fixed" adversarial mode).

The eigvec_fixed in the name is deliberate: the only consumer of this information is :class:~ppsds.environment.dynamics_mode.EigvecFixedAdversarialDynamics -- it delivers the Eig-DRPP bound (value_lower_bound_eigvec_fixed). Calling it inject_adversarial_info would present this very narrow use as a general one.

This is the only public channel from the agent side to a :class:DynamicsMode. It is necessary because of the paper's semantics: the adversarial SDS must "pick the worst SDS from the predicted covariance after the predictor emits its pdf" -- i.e. a best response, not a fixed strategy independent of the predictor. So it is not an architectural defect.

For modes that do not consume predictor information (FixedDynamics / RandomDynamics / PreSampledDynamics) the base set_info is a no-op, so this can be called unconditionally.

Parameters:

Name	Type	Description	Default
`lamb_pred`	`ndarray of shape (x_dim,) or None`	Eigenvalues of the predictive covariance. `None` skips the injection -- on the first `step()` there is no prediction yet, and the mode falls back to nominal parameters internally.	required
`u`	`ndarray of shape (u_dim,)`	Control input for this step.	required
`x`	`ndarray of shape (x_dim,)`	State for this step; defaults to `self.x`. Passing it explicitly makes the call site self-explanatory and avoids depending on the implementation detail of when the state advances.	`None`

dynamics_mode ¶

Dynamics-mode strategies: which set of SDS parameters each step uses.

Everything a mode needs from an ambiguity set is its generation ability:

sample_trajectory (used by :class:RandomDynamics),
worst_case_params (used by :class:EigvecFixedAdversarialDynamics).

:class:SamplingAmbiguitySet pins that contract down as a Protocol: a mode does not know about "ambiguity sets" as a concept, only about objects that can hand out a parameter trajectory and a worst-case parameter. :class:PreSampledDynamics does not even need that -- its parameters were sampled outside.

Layout

Non-adversarial modes are flat in :mod:.non_adversarial; adversarial modes are filed by assumption / family, mirroring :mod:ppsds.ambiguity::

dynamics_mode/
├── base.py                     SamplingAmbiguitySet + DynamicsMode
├── non_adversarial.py          Fixed / Random / PreSampled
└── adversarial/
    └── additive_drift/
        └── conic_moment.py     EigvecFixedAdversarialDynamics

An adversarial mode's notion of "worst case" is defined entirely by the ambiguity-set family it targets, so a new family gets a new module here rather than a new flag on an existing mode. That is also why the one adversarial mode is named narrowly: see :class:EigvecFixedAdversarialDynamics.

Usage

::

mode = RandomDynamics(abg, period=100)
mode.reset(rng, T=50)
for step in range(T):
    params = mode.get_params(step)   # keys declared by abg.param_keys

DynamicsMode ¶

Bases: ABC

Strategy interface: pick one set of SDS parameters per step.

A concrete mode implements reset() (called when an episode starts) and get_params(step_idx) (called before every env.step()).

reset `abstractmethod` ¶

reset(rng: Generator, T: int) -> None

(Re)initialize, preparing for a new episode.

Parameters:

Name	Type	Description	Default
`rng`	`Generator`	Seeded generator obtained from `env.np_random`.	required
`T`	`int`	Number of steps in the episode.	required

get_params `abstractmethod` ¶

get_params(step_idx: int) -> dict

Return the SDS parameters to use at this step.

Parameters:

Name	Type	Description	Default
`step_idx`	`int`	Current step index (0-based).	required

Returns:

Type	Description
`dict`	The key set is declared by the ambiguity set's `param_keys`; a missing key means "take the nominal value for that item" (so `{}` means "all nominal").

set_info ¶

set_info(lamb_pred, x, u) -> None

Optional hook: inject predictor information before step().

By default it does nothing -- only :class:EigvecFixedAdversarialDynamics consumes it. It is defined on the base class so that the caller (Env.inject_eigvec_fixed_info) does not have to sniff with hasattr(mode, 'set_info').

SamplingAmbiguitySet ¶

Bases: Protocol

The minimal contract :class:DynamicsMode needs: it can generate SDS parameters.

:class:~ppsds.ambiguity.additive_drift.conic_moment.AdditiveDriftConicMomentSet satisfies it naturally, but the two have no inheritance relation -- a mode uses this structural contract only.

sample_trajectory ¶

sample_trajectory(T: int, period: int = 1) -> dict

Sample one piecewise-constant parameter trajectory; each value's first axis is T+1.

worst_case_params ¶

worst_case_params(x, u, lamb_pred) -> dict

Worst-case parameters given the predicted covariance eigenvalues (includes key 'j').

FixedDynamics ¶

Bases: DynamicsMode

Always use the nominal dynamics -- no perturbation at all.

The degenerate case: when the ambiguity set shrinks to a single point, simulating is equivalent to running that deterministic system directly.

PreSampledDynamics ¶

PreSampledDynamics(sds_params: dict)

Bases: DynamicsMode

Replay a parameter trajectory that was sampled outside.

The only difference from :class:RandomDynamics is who samples: here the caller samples it and passes it in. This suits the reproducible setting where one parameter trajectory is shared by every episode.

Parameters:

Name	Type	Description	Default
`sds_params`	`dict of ndarrays`	The output of `abg.sample_trajectory(T, period)`. Every value has first axis `T+1`.	required

RandomDynamics ¶

RandomDynamics(abg: SamplingAmbiguitySet, period: int = 1)

Bases: DynamicsMode

Re-sample SDS parameters from the ambiguity set every period steps.

reset() samples the whole parameter trajectory in one go via abg.sample_trajectory(T, period). Therefore:

same seed -> every episode sees the same sds_params ("sample once, run many");
different seeds -> episodes are independent.

Parameters:

Name	Type	Description	Default
`abg`	`SamplingAmbiguitySet`		required
`period`	`int`	Number of steps between two re-samples (1 = i.i.d., larger = closer to constant).	`1`

EigvecFixedAdversarialDynamics ¶

EigvecFixedAdversarialDynamics(abg: SamplingAmbiguitySet)

Bases: DynamicsMode

Pick SDS parameters adversarially so the predictor's score is worst.

Why the name is so narrow: it delivers the ambiguity set's :meth:~ppsds.ambiguity.additive_drift.conic_moment.AdditiveDriftConicMomentSet.value_lower_bound_eigvec_fixed -- the bound Eig-DRPP guarantees under the "eigenvectors fixed" relaxation. Calling it AdversarialDynamics would overclaim: this mode's "worst" is defined entirely by abg.worst_case_params(), so swapping in another ambiguity-set family changes the worst-case rule, and that family needs a mode class of its own.

It needs information from the predictor -- specifically lamb_pred (the eigenvalues of its current predictive covariance). set_info(lamb_pred, x, u) must be called before each env.step().

The actual worst-case computation is delegated to abg.worst_case_params(x, u, lamb_pred), which returns a dict shaped like sample()'s output.

Parameters:

Name	Type	Description	Default
`abg`	`SamplingAmbiguitySet`	Must implement `worst_case_params()`.	required

set_info ¶

set_info(lamb_pred: ndarray, x: ndarray, u: ndarray) -> None

Inject predictor information before the next step().

In adversarial mode this must be called once per step.

Parameters:

Name	Type	Default
`lamb_pred`	`ndarray of shape (x_dim,)`	required
`x`	`ndarray of shape (x_dim,)`	required
`u`	`ndarray of shape (u_dim,)`	required

known_sds ¶

A Gymnasium environment for a system that is known exactly.

Sibling of :mod:~ppsds.environment.base, and the third thing this subpackage serves: where Env composes an ambiguity set with a :class:~ppsds.environment.dynamics_mode.DynamicsMode, :class:KnownSdsEnv steps an :class:~ppsds.dynamics.known_sds.AdditiveDriftSds directly. No ambiguity set is involved, so a non-Gaussian noise law -- which the ambiguity sets' Gaussian transition cannot express -- is exactly what this environment is for.

It exists so that "set up an environment, predict, score" is one story across the package: the same reset(seed=...) / step(action) cycle the rest of the library uses, on a system whose truth is known, which is what the predictability questions are asked about.

KnownSdsEnv ¶

KnownSdsEnv(sds, x0=None, max_steps=50)

Bases: Env

Steps a known :class:~ppsds.dynamics.known_sds.AdditiveDriftSds.

Parameters:

Name	Type	Description	Default
`sds`	`AdditiveDriftSds`	The system. Its `noise_law` and `f_bar` stay reachable through this environment, so a predictor and the predictability reference lines can be built from the same object.	required
`x0`	`array_like`	Initial state; defaults to the system's own.	`None`
`max_steps`	`int`	Horizon after which `step` reports `truncated`. As in :class:`~ppsds.environment.base.Env`, nothing is ever `terminated`: these systems run forever, they are just being simulated for a while.	`50`

Attributes:

Name	Type	Description
`x`	`ndarray or None`	The current state, public as it is on :class:`~ppsds.environment.base.Env`. `None` before the first `reset`.

noise_law `property` ¶

noise_law

The process-noise law, for building scoring rules and predictability references.

f_bar `property` ¶

f_bar

The known drift.

x0 `property` ¶

x0

The initial state this environment resets to.

reset ¶

reset(seed=None, options=None)

Start a fresh episode; options['x0'] overrides the initial state.

step ¶

step(action)

Advance one step; returns (state, reward, terminated, truncated, info).

The reward is always 0.0: this is a prediction environment, and what a prediction earns is measured by a scoring rule, not by a reward signal.

Numerical helpers¶

atom_pmf ¶

A discrete predictive law: probability atoms with weights.

:class:AtomPmf is the "predictive law may be a pmf" object. Its support need not contain the support of the truth, and that is the point: under a neighbourhood scoring rule a small set of well-placed atoms is a perfectly usable prediction, whereas under the ordinary logarithm score it would be hopeless almost surely.

Which is what :meth:AtomPmf.logpdf is here to make concrete. It exists, and it returns -inf, because an atom law has no density: scoring it with :class:~ppsds.scoring_rule.logarithm_score.LogarithmScore gives -inf for every outcome, correctly. Score it with :class:~ppsds.scoring_rule.epsilon_logarithm_score.EpsilonLogarithmScore instead -- the box integral :meth:AtomPmf.box_integral is what that score evaluates.

AtomPmf ¶

AtomPmf(atoms, weights)

A weighted set of atoms: sum_i weights[i] * delta(atoms[i]).

Parameters:

Name	Type	Description	Default
`atoms`	`array_like`	`(n_atoms, d_x)`. `(d_x,)` for a single atom is promoted to one row.	required
`weights`	`array_like`	`(n_atoms,)`, non-negative. Normalised on construction, so raw counts or stakes are accepted as they are; `AtomPmf.weights` is what the law actually uses.	required

Raises:

Type	Description
`ValueError`	If the shapes disagree, a weight is negative, or the weights sum to zero.

dimensions `property` ¶

dimensions

The state dimension d_x, matching multivariate_uniform's attribute name.

mean `property` ¶

mean

The weighted mean of the atoms.

cov `property` ¶

cov

The weighted covariance of the atoms -- the moment summary of the law.

entropy ¶

entropy()

The Shannon entropy of the weights, -sum w log w.

Not a differential entropy, and not interchangeable with one: this is non-negative, while a differential entropy is a log of a volume and may be negative. The predictability upper bound takes the entropy of the continuous truth, so do not pass an :class:AtomPmf as that truth.

logpdf ¶

logpdf(x)

-inf everywhere -- an atom law has no density.

Provided rather than omitted so that scoring a pmf with the ordinary logarithm score runs and reports -inf, instead of failing with an attribute error that hides the reason. Use the epsilon-logarithm score.

rvs ¶

rvs(size=1, rng=None)

Draw size atoms according to the weights. (size, d_x).

Parameters:

Name	Type	Description	Default
`size`	`int`	Number of draws.	`1`
`rng`	`Generator`	Source of randomness; `None` uses the global `np.random` state (the same convention as :class:`~ppsds.numerics.random_variables.multivariate_uniform`).	`None`

box_integral ¶

box_integral(x, eps)

Total weight within the L-infinity box of half-width eps around each point.

Returns sum_i weights[i] over atoms with ||atom_i - x||_inf <= eps. This is the quantity the epsilon-logarithm score takes the logarithm of, and the reason a pmf is scoreable at all.

Parameters:

Name	Type	Description	Default
`x`	`array_like`	`(d_x,)` for one outcome, or `(n, d_x)` for a batch -- the shape a rollout produces.	required
`eps`	`float`	Box half-width, `>= 0`.	required

Returns:

Type	Description
`float or ndarray`	A scalar, or `(n,)`.

Notes

Exact, and vectorised over the outcome batch: every atom is tested against every outcome, in chunks that bound the intermediate boolean array. The cost is O(n_atoms * n_outcomes * d_x), which sounds worse than a spatial index but is not in practice -- measured at 0.9 s for the full 10**5 trajectories x 30 steps (3e6 outcomes) against 50 atoms, and it stays exact, which cell-summing over a lattice would not be: with cells two eps wide and a box one eps either side, a neighbouring cell only overlaps the box partially, so summing whole cells over-counts. If a pmf ever carries thousands of atoms, a cKDTree ball query with p = inf is the next step; at tens to hundreds it would be slower, not faster.

linalg ¶

Small numerical utilities that keep behavior consistent across numpy versions.

Why eig_real is needed

What dtype np.linalg.eig returns for a real matrix is not consistent across numpy builds: for the same real symmetric covariance matrix, numpy 2.3.5 returns float64 while 2.5.4 returns complex128 (imaginary parts exactly 0, real parts bit-for-bit identical to the former).

The downstream cone program requires real-valued input and raises The 'minimize' objective must be real valued on the latter — and the error point is very far from the real cause.

A covariance matrix is necessarily symmetric, so its eigenvalues are necessarily real, which makes taking the real part lossless; for return values that are already real arrays, .real returns them unchanged, bit for bit.

eig_real ¶

eig_real(matrix)

A real-valued version of np.linalg.eig: returns (eigenvalues, eigenvectors), both guaranteed to be real arrays.

Should only be used when the input is a symmetric matrix (e.g. a covariance) — in that case the eigenvalues are necessarily real, so taking the real part is lossless.

random_variables ¶

Small random-variable utilities.

:class:multivariate_uniform — a multivariate uniform distribution with independent components. Used by the core: when sampling an ambiguity set, the noise mean is drawn from a box (see :mod:ppsds.ambiguity.additive_drift.conic_moment).

.. note:: rvs() defaults to the global np.random and has no seed of its own, so reproducibility is the caller's responsibility (e.g. np.random.seed(...) before sampling) -- or pass an explicit rng, which is the form new code should use.

pdf() / cdf() accept a single point (d,) and return a scalar, or many points (n, d) and return (n,).

A multivariate_gmm (Gaussian mixture) class lived here until 2026-10-16. It was merged in with src/lab/ on 2026-10-12, had zero callers for its whole life, and shipped in the wheel as unexercised code. Removed rather than carried: a class nothing calls and nothing tests is a claim the package cannot back up.

multivariate_uniform ¶

multivariate_uniform(low=-np.ones(2), high=np.ones(2))

A multivariate uniform distribution with independent components; the box is given by low / high.

Initializes a multivariate uniform random variable.

Parameters: low (numpy.ndarray): A 1D array representing the lower bounds of the uniform distribution for each dimension. high (numpy.ndarray): A 1D array representing the upper bounds of the uniform distribution for each dimension.

rvs ¶

rvs(size=1, rng=None)

Generates a sample from the multivariate uniform distribution.

Parameters: size (int): The number of samples to generate. rng (numpy.random.Generator, optional): Source of randomness. None (the default) uses the global np.random state, which is what the existing callers rely on; pass a Generator for a reproducible draw.

Returns: numpy.ndarray: A 2D array of shape (size, dimensions) containing the samples.

pdf ¶

pdf(x)

Density at a point (d,), or at many points (n, d) -> (n,).

cdf ¶

cdf(x)

CDF at a point (d,) (scalar), or at many points (n, d) -> (n,).

Control policy¶

lqr ¶

LQR control policy.

The regulation policy used by the linear experiments: drive the state to a target x_d by infinite-horizon discrete-time LQR.

Experiments and figures¶

Both need pip install "ppsds[experiments]".

cache ¶

Validation and atomic writing for the experiment pickle cache.

Background

The experiment runner used to skip recomputation based on os.path.isfile alone, while the cache key was a directory path that did not encode -t (the trajectory count). Measured: the seed=3/period=100 directory held only 5 trajectories, and re-running with -t 100 was silently skipped -- so a figure claiming 100 trajectories was drawn from 5. And data/** is in .gitignore, so overwriting loses it permanently.

This module replaces "skip if the file exists" with "skip only if four gates pass":

Content validation -- the pickle reads back, the key set matches exactly, and every array's first axis is n_traj (per-trajectory) or T+1 (per-step parameters shared across trajectories). That is enough to catch truncated/half-written files and files with the wrong trajectory count, without a hardcoded shape table per file.
Fingerprint validation -- _cache_manifest.json in the output directory records fingerprint = cache format version + source hash + parameter snapshot. Any change invalidates.
Dependency validation -- predictor results depend on traj.pkl (which carries the shared sds_params). Saving records each dependency's sha256; loading compares them, so a changed dependency recomputes this file.
Atomic writing -- write .tmp then os.replace, which rules out half-written files at the root.

The directory layout is unchanged (still data/.../seed=N/period=P): the fingerprint lives only in the manifest, so the existing data/ tree keeps its names.

Usage

::

ctx = CacheContext(dir_out, snapshot={...}, force=args.force)
data = ctx.load('nominal.pkl', KEYS, n_traj, T, deps=['traj.pkl'])
if data is None:
    data = ...                      # actually recompute
    ctx.save(data, 'nominal.pkl', deps=['traj.pkl'])
print(ctx.summary())

CacheContext ¶

CacheContext(dir_out, snapshot, force=False)

Cache context for one experiment run.

Parameters:

Name	Type	Description	Default
`dir_out`	`str`	Output directory for this run (keeps the existing `data/...` naming; directory names are not changed).	required
`snapshot`	`dict`	Every parameter that determines the output content. The caller is responsible for listing them all -- especially in the NLink case, where constants are bound into lambda default arguments (`dt` / `gamma_0` etc. are not attributes of `abg`), so they must be listed explicitly or changing them will not invalidate.	required
`force`	`bool`	If true, recompute everything (the manifest is still updated on write).	`False`

load ¶

load(name, keys, n_traj, T, deps=())

Return the data if it is fresh, else None (recording the reason).

save ¶

save(data, name, deps=())

Write atomically, recording the fingerprint and each dependency's sha256 in the manifest.

source_fingerprint ¶

source_fingerprint()

Hash of the __main__ and ppsds.* sources loaded in this process.

It walks sys.modules rather than a hardcoded file list, so it is sensitive only to code that is actually used: editing a file nobody imported (dead code, say) does not needlessly invalidate the cache -- and that costs thousands of cone programs. It also works on a dirty working tree: a git hash does not change before you commit, whereas this does, and that is precisely when the cache must be invalidated.

⚠️ When to call: because it reads the "already loaded" set, it must be called after all imports are finished. CacheContext is constructed inside __main__, which satisfies that naturally. Calling it midway through import would make the same code produce different fingerprints.

build_fingerprint ¶

build_fingerprint(snapshot)

Combined fingerprint of format version + source hash + parameter snapshot.

verify ¶

verify(path, keys, n_traj, T)

Read and validate path; return the data if it passes, else None.

Every array's first axis must be n_traj (per-trajectory) or T+1 (per-step parameters shared across trajectories); and if the former has two or more dimensions, axis 1 must be T+1. Those two conditions are enough to catch "wrong trajectory count" and "half-written/truncated", without hardcoding a shape table per file.

predictability ¶

The predictability figure: per-step score curves against the two reference lines.

One panel shows, for a single neighbourhood radius:

a few individual trajectories' running averages (thin) -- the point being that they converge to the expectation within a couple of dozen steps, so one trajectory is enough;
the expected score over many trajectories (solid) with a confidence band;
the predictability upper bound (dashed), and optionally the O(eps) approximation for the given predictor (dashed).

Both reference values are passed in, never recomputed here. They come from :mod:ppsds.predictability.bounds, and computing them inside a plotting function would both hide that dependency and invite the figure to disagree with the number a test asserts.

The approximation line is optional on purpose: an atom/pmf predictor has no cross-entropy against a continuous truth (see :func:~ppsds.predictability.entropy.cross_entropy), so for those panels there is no approximation to draw, only the upper bound and the score itself.

Which trajectories get drawn

By default, three chosen at random per panel, not "the first three". This is not a detail: drawing the same three trajectories in every panel makes it look as though the individual curves were a property of the radius, when they are a property of the sample. The choice is seeded from eps so a figure is reproducible without being identical across panels.

choose_trajectories ¶

choose_trajectories(n_traj, eps, count=DEFAULT_N_CURVES)

count distinct trajectory indices, differing between radii but reproducible.

Seeded from eps so that every panel draws its own sample and a re-run draws the same figure. Fewer trajectories than requested degrades to all of them, in order -- a short rollout should thin the figure, not break it.

plot_predictability_panel ¶

plot_predictability_panel(ax, scores, *, eps, upper_bound, approximation=None, trajectory_indexes=None, draw_band=True, legend=True)

Draw one panel from a (n_traj, T) block of per-step scores.

Parameters:

Name	Type	Description	Default
`ax`	`Axes`	The axes to draw on.	required
`scores`	`array_like`	`(n_traj, T)` -- the per-step epsilon-logarithm scores, not yet averaged.	required
`eps`	`float`	The neighbourhood radius. It also seeds :func:`choose_trajectories` when `trajectory_indexes` is not given.	required
`upper_bound`	`float`	Predictability upper bound, drawn as the dashed red line.	required
`approximation`	`float`	The `O(eps)` approximation for this predictor, drawn as the dashed blue line. `None` omits it, which is the pmf case.	`None`
`trajectory_indexes`	`sequence of int`	Which individual trajectories to draw thin. Defaults to :func:`choose_trajectories` -- three at random for this radius, so panels do not share a sample. Indices are clamped into range, warning if changed, rather than raising: a short rollout should degrade the figure, not the run.	`None`
`draw_band`	`bool`	Whether to shade the confidence band between the two quantiles.	`True`
`legend`	`bool`	Whether to draw the legend, whose four entries are the published figure's.	`True`

Returns:

Type	Description
`dict`	`{'mean', 'low', 'high'}`, each `(T,)` -- the aggregates actually drawn, so a caller can print the numbers it just plotted instead of recomputing them.

Next¶

This page is a reference, not a route. If you arrived here from a search:

What the guarantee costs, before you rely on it: Research
Get something running first: Getting Started
How the pieces fit together: the map — the routes through the examples, and which subpackage each one uses
Every example, with its real output: Examples

API reference¶

What each subpackage answers¶

Two layers inside one package¶

Predictors¶

conic_moment ¶

NominalPredictor ¶

NoiseDrpp ¶

EigDrpp ¶

reset ¶

trajectory_prediction ¶

eigvec_fixed_minmax ¶

atom ¶

StaticAtomPmfPredictor ¶

DriftShiftedAtomPmfPredictor ¶

base ¶

Ambiguity sets¶

conic_moment ¶

AdditiveDriftConicMomentSet ¶

dim_x property ¶

dim_u property ¶

gamma property ¶

mean ¶

sample ¶

validate ¶

transition ¶

sample_trajectory ¶

value_upper_bound_true_transition ¶

value_lower_bound_eigvec_fixed ¶

worst_case_params ¶

LinearAmbiguitySet ¶

f_bar ¶

f_err ¶

mean ¶

sample ¶

NonlinearAmbiguitySet ¶

mean ¶

sample ¶

Scoring rules¶

logarithm_score ¶

LogarithmScore ¶

score ¶

epsilon_logarithm_score ¶

EpsilonLogarithmScore ¶

score ¶

bin_to_list ¶

Predictability: the ceiling¶

bounds ¶

predictability_upper_bound ¶

approximation_value ¶

entropy ¶

differential_entropy ¶

cross_entropy ¶

trajectory ¶

cumulative_average ¶

across_trajectories ¶

The system itself¶

known_sds ¶

AdditiveDriftSds ¶

transition ¶

sample_trajectory ¶

identification ¶

LinearFit dataclass ¶

f_bar ¶

predict ¶

residuals ¶

encode_action ¶

collect_rollouts ¶

fit_linear ¶

diagnose ¶

regularize ¶

format_report ¶

mechanical ¶

NLink ¶

equilibrium_for_pose ¶

linearize_at_discrete ¶

rk4_step ¶

The environment layer¶

base ¶

Env ¶

ambiguity_set property ¶

dim_x `property` ¶

dim_u `property` ¶

gamma `property` ¶

LinearFit `dataclass` ¶

ambiguity_set `property` ¶

dynamics_mode `property` ¶

reset `abstractmethod` ¶

get_params `abstractmethod` ¶

noise_law `property` ¶

f_bar `property` ¶

x0 `property` ¶

dimensions `property` ¶

mean `property` ¶

cov `property` ¶