API Reference

prophys.symbolic.deg2rad(x)[source]#

Elementwise conversion of an angle expression from degrees to radians.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.rad2deg(x)[source]#

Elementwise conversion of an angle expression from radians to degrees.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.relu(x)[source]#

Elementwise rectified linear unit max(x, 0) of an expression.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.softplus(x)[source]#

Elementwise softplus log(1 + exp(x)) of an expression, a smooth positive ramp.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.sigmoid(x)[source]#

Elementwise logistic sigmoid 1 / (1 + exp(-x)) of an expression.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.tanh(x)[source]#

Elementwise hyperbolic tangent of an expression.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.clip(x, lo=None, hi=None)[source]#

Elementwise clamp of an expression to the interval [lo, hi].

Parameters:

x (Any) – Expression or value.
lo (float | None) – Lower bound, or None for no lower bound.
hi (float | None) – Upper bound, or None for no upper bound.

Return type:

prophys.symbolic.where(cond, a, b)[source]#

Elementwise selection: a where cond is true, otherwise b.

Both branches are evaluated, and gradients flow only through the selected one. Build cond with comparisons such as x > 0.

Parameters:

cond (Any)
a (Any)
b (Any)

Return type:

prophys.symbolic.minimum(a, b)[source]#

Elementwise minimum of two broadcastable expressions.

Parameters:

a (Any)
b (Any)

Return type:

prophys.symbolic.maximum(a, b)[source]#

Elementwise maximum of two broadcastable expressions.

Parameters:

a (Any)
b (Any)

Return type:

Reductions and array construction#

prophys.symbolic.sum_(x, axis=None)[source]#

Sum of an expression over axis (all elements when axis=None).

Parameters:

x (Any)
axis (int | None)

Return type:

prophys.symbolic.mean(x, axis=None)[source]#

Arithmetic mean of an expression over axis (all elements when axis=None).

Parameters:

x (Any)
axis (int | None)

Return type:

prophys.symbolic.prod(x, axis=None)[source]#

Product of an expression over axis (all elements when axis=None).

Parameters:

x (Any)
axis (int | None)

Return type:

prophys.symbolic.logsumexp(x, axis=None)[source]#

Numerically stable log(sum(exp(x))) over axis (all elements when axis=None).

Parameters:

x (Any)
axis (int | None)

Return type:

prophys.symbolic.softmin(x, axis=-1, tau=1.0)[source]#

Smooth, differentiable approximation of the minimum along axis.

Computes -tau * logsumexp(-x / tau), which lies slightly below the hard minimum and converges to it as tau -> 0.

Parameters:

x (Any) – Expression or array.
axis (int) – Axis to reduce.
tau (float) – Smoothing temperature, in the units of x.

Return type:

prophys.symbolic.softmax(x, axis=-1, tau=1.0)[source]#

Smooth, differentiable approximation of the maximum along axis.

Computes tau * logsumexp(x / tau), which lies slightly above the hard maximum and converges to it as tau -> 0. (This is the smooth maximum, not the normalized exponential “softmax” of machine learning.)

Parameters:

x (Any) – Expression or array.
axis (int) – Axis to reduce.
tau (float) – Smoothing temperature, in the units of x.

Return type:

prophys.symbolic.dot(a, b)[source]#

Dot product of two expressions, with jnp.dot semantics (inner product of vectors).

Parameters:

a (Any)
b (Any)

Return type:

prophys.symbolic.norm(x, axis=-1)[source]#

Euclidean (L2) norm of an expression along axis (the last axis by default).

Parameters:

x (Any)
axis (int)

Return type:

prophys.symbolic.stack(exprs, axis=0)[source]#

Stack expressions of the same shape along a new axis.

Parameters:

exprs (list[Any]) – Expressions or values of identical shape.
axis (int) – Position of the new axis in the result.

Return type:

Example:

xy = prp.stack([x, y], axis=-1)  # shape (..., 2)

prophys.symbolic.concatenate(exprs, axis=0)[source]#

Join expressions along an existing axis.

Parameters:

exprs (list[Any]) – Expressions or values whose shapes agree except along axis.
axis (int) – Axis along which to join.

Return type:

Special functions#

Evaluated by optimized native kernels with exact gradients.

prophys.symbolic.erf(x)#

Elementwise error function erf(x) = 2/sqrt(pi) * integral_0^x exp(-t^2) dt.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.erfc(x)#

Elementwise complementary error function 1 - erf(x), accurate far into the upper tail.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.erfcx(x)#

Elementwise scaled complementary error function exp(x^2) erfc(x), free of underflow for large x.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.ndtr(x)#

Elementwise standard normal cumulative distribution function Phi(x).

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.log_ndtr(x)#

Elementwise logarithm of the standard normal CDF, log Phi(x), accurate in the lower tail.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.exp1(x)#

Elementwise exponential integral E1(x) for x > 0 (e.g. the Theis well function).

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.i0e(x)#

Elementwise exponentially scaled modified Bessel function of the first kind, exp(-|x|) I0(x).

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.i1e(x)#

Elementwise exponentially scaled modified Bessel function of the first kind, exp(-|x|) I1(x).

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.k0e(x)#

Elementwise scaled modified Bessel function of the second kind, exp(x) K0(x) for x > 0.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.k1e(x)#

Elementwise scaled modified Bessel function of the second kind, exp(x) K1(x) for x > 0.

Parameters:: x (Any)
Return type:: Expr

prophys.symbolic.owens_t(h, a)[source]#

Elementwise Owen’s T function T(h, a).

T(h, a) is the probability mass of a standard bivariate normal over the wedge between x = h and y = a x. It is the building block of bivariate normal probabilities over wedges and polygons.

Parameters:

h (Any)
a (Any)

Return type:

Frames#

class prophys.frames.Frame(name, ndim=2, units='m', crs=None)[source]#

A named 2D or 3D coordinate system that spatial structures live in.

Structures declared in the same frame can be combined. Combining structures from different frames raises FrameError unless a FrameTransform between them is registered. Frames compare by value, so two Frame("site") instances are the same frame.

Parameters:

name (str) – Unique identifier used for equality and compatibility checks.
ndim (int) – Number of spatial dimensions, 2 or 3.
units (str) – Coordinate unit, which must be a length (e.g. "m" or "km"). It is parsed with prophys.units.parse() at construction and used for unit conversion between frames.
crs (str | None) – Optional EPSG or CRS tag. Purely informational, prophys does not read it.

Raises:

ValueError – If ndim is not 2 or 3.
prophys.units.UnitError – If units is not a length unit.

Example:

site = prp.Frame("site", ndim=2, units="m")

property unit: Unit#: The parsed coordinate unit of this frame.

check(other)[source]#

Raise FrameError unless other is the same frame as this one.

Parameters:: other (Frame) – Frame to compare against.
Return type:: None

class prophys.frames.FrameTransform(from_, to, matrix=None, offset=None)[source]#

A differentiable affine map between two frames: y = R @ x + t.

matrix has shape (to.ndim, from_.ndim), offset has shape (to.ndim,). Both are plain arrays (not Expr), because frame alignment is a fixed, known calibration, not an optimizable quantity. The map is differentiable with respect to the points.

Parameters:

from – Source frame.
to (Frame) – Target frame.
matrix (Any) – Linear part R. Defaults to the identity scaled by the unit conversion factor between the two frames (e.g. 1000 from km to m).
offset (Any) – Translation t in units of the target frame. Defaults to zero.
from_ (Frame)

property unit_scale: float#: The pure unit conversion factor between the two frames, ignoring any rotation or translation this transform also applies.

apply(points)[source]#

Map an (..., from_.ndim) array of points into the target frame.

Returns:: Array of shape (..., to.ndim).
Parameters:: points (Any)
Return type:: Any

inverse()[source]#

Return the inverse transform, mapping the target frame back to the source frame.

Requires a square, invertible matrix.

Return type:: FrameTransform

class prophys.frames.FrameRegistry[source]#

Registry of transforms between frames, looked up when points are aligned.

Structures consult their registry when combined with a structure from another frame, so related frames can be mixed without passing transforms through every call. default_registry is the instance for the whole process, used by structures created without a registry of their own.

Example:

default_registry.register(FrameTransform(local, utm, offset=[500000.0, 5.6e6]))
xy_utm = default_registry.align(xy_local, local, utm)

register(transform)[source]#

Register transform and its inverse, so points can be aligned in both directions.

A later registration for the same pair of frame names replaces the earlier one. Computing the inverse needs a square, invertible matrix.

Parameters:: transform (FrameTransform)
Return type:: None

get(from_, to)[source]#

Return the registered transform from from_ to to, or None if there is none.

Parameters:

from_ (Frame)
to (Frame)

Return type:

FrameTransform | None

align(points, from_, to)[source]#

Return points expressed in frame to, applying a registered transform if the frames differ.

Raises:

FrameError – If the frames differ and no transform is registered.

Parameters:

points (Any)
from_ (Frame)
to (Frame)

Return type:

exception prophys.frames.FrameError[source]#: Raised when two structures with incompatible frames are combined.

prophys.frames.default_registry: FrameRegistry#: Process-wide FrameRegistry used by every structure that is not given an explicit registry.

Physics primitives#

prophys.physics.superpose(receptors, sources, strengths, *, kind='steady', direction=0.0, travel=0.0, spread_across=(1.0, 1.0), spread_vertical=(1.0, 1.0), spread_along=(1.0, 1.0), decay=0.0, lid=0.0, anisotropy=1.0, radial_scale=1.0, vertical='none', images=0)[source]#

Superpose the Green’s function contributions of all sources at every receptor.

Computes sum_s strength_s * f(receptor, source_s), the building block of plume, puff, front and radial diffusion models.

Displacements are resolved along direction (radians): u along it, v across it. Spreads follow power laws sigma = a * t**b of a spread argument t. Kinds:

"steady": t = u, G(v; sigma_y) V(z) exp(-decay u) for u > 0 and zero otherwise, a steady advected plume.
"moving": t = travel, G(u - travel; sigma_x) G(v; sigma_y) V(z) exp(-decay travel), a released puff after it travelled.
"front": as "moving" with the along profile erfc((u - travel)/(sqrt(2) sigma_x))/2, a continuous release.
"radial_exp1": E1(rho**2 / radial_scale) exp(-decay rho) with rho**2 = u**2 + anisotropy v**2, radial transient diffusion (the Theis well function for anisotropy = 1).

G is the normal density. The vertical factor V is 1 (vertical="none"), G(z - H; sigma_z) ("free"), or the image series of a reflecting ground and a lid at height lid with images pairs ("reflect", no lid images when lid <= 0). Receptors are [R, 3], sources [S, 3] (the third coordinate is height). Every parameter may be scalar or carry leading sample dimensions, and the result is [..., R].

Parameters:

receptors (Any) – Receptor coordinates (x, y, height), shape [R, 3].
sources (Any) – Source coordinates (x, y, height), shape [S, 3].
strengths (Any) – Source strengths, shape [S].
kind (str) – "steady", "moving", "front" or "radial_exp1".
direction (Any) – Transport direction in radians.
travel (Any) – Travel distance of a puff or front ("moving", "front").
spread_across (tuple[Any, Any]) – (a, b) of the crosswind spread sigma_y = a t**b. Unused by "radial_exp1".
spread_vertical (tuple[Any, Any]) – (a, b) of the vertical spread sigma_z.
spread_along (tuple[Any, Any]) – (a, b) of the along wind spread sigma_x.
decay (Any) – First order decay rate per unit of u ("steady"), travel or radial distance.
lid (Any) – Mixing lid height for vertical="reflect", 0 or less for none.
anisotropy (Any) – Cross direction weight in "radial_exp1".
radial_scale (Any) – Denominator of the E1 argument in "radial_exp1".
vertical (str) – "none", "free" or "reflect".
images (int) – Number of image pairs for vertical="reflect".

Returns:

The superposed field at each receptor, shape [..., R].

Raises:

ProphysError – PHY102 for an unknown kind or vertical mode.

Return type:

prophys.physics.gaussian_polygon_probability(mean, cov_cholesky, vertices)[source]#

Compute the exact probability P(x in polygon) for a bivariate normal position x ~ N(mean, L L^T).

mean is [..., 2], cov_cholesky the lower Cholesky factor as [..., 3] = [l00, l10, l11] (or [..., 2, 2]), vertices a simple polygon [V, 2] in either orientation. Exact (Owen’s T edge decomposition, evaluated without cancellation), with exact gradients with respect to all three arguments.

Parameters:

mean (Any) – Mean position, shape [..., 2].
cov_cholesky (Any) – Lower Cholesky factor of the covariance, [..., 3] or [..., 2, 2].
vertices (Any) – Polygon vertices, shape [V, 2].

Returns:

The probability, shape [...].

Return type:

prophys.physics.overlap_area(region, polygon=None, *, disk=None)[source]#

Compute the exact overlap area of a polygonal region with a polygon or a disk.

region and polygon are simple polygons [..., V, 2] in either orientation, disk is [..., 3] = (center_x, center_y, radius). Computed as a boundary integral (Green’s theorem) with closed form crossings. Gradients with respect to every vertex, the centre and the radius are exact. At most one leading batch dimension per argument.

Parameters:

region (Any) – Region polygon, shape [..., V, 2].
polygon (Any) – Second polygon, shape [..., W, 2]. Pass exactly one of polygon and disk.
disk (Any) – Disk (center_x, center_y, radius), shape [..., 3].

Returns:

The area of region ∩ polygon or region ∩ disk.

Raises:

ProphysError – PHY101 unless exactly one of polygon and disk is given.

Return type:

prophys.physics.ray_integral(values, segments, *, origin=(0.0, 0.0), cell=(1.0, 1.0), reduce='sum', threshold=0.0)[source]#

Compute exact reductions of a raster along line segments (path integrals, maxima, lengths above a threshold).

values is [..., ny, nx] (rows along y), segments [M, 2, 2] (start and end points). Segments are clipped to the raster. reduce:

"sum": integral of v dl over cells of constant value, cell (r, c) spanning origin + (c, r) * cell to origin + (c+1, r+1) * cell (attenuation optical depth, travel cost, column density).
"bilinear": integral of f dl for the bilinear interpolant of values at the nodes origin + (c, r) * cell.
"max": the largest cell value crossed (line of sight, obstruction), 0 for a segment that misses the raster.
"length_above": the length inside cells whose value exceeds threshold.

Gradients flow to values and to the segment endpoints.

Parameters:

values (Any) – Raster values, shape [..., ny, nx] with rows along y.
segments (Any) – Segment end points, shape [M, 2, 2].
origin (tuple[float, float]) – World coordinate (x, y) of the raster origin.
cell (tuple[float, float]) – Cell size (dx, dy).
reduce (str) – "sum", "bilinear", "max" or "length_above".
threshold (float) – Value threshold for "length_above".

Returns:

One result per segment, shape [..., M].

Raises:

ProphysError – PHY102 for an unknown reduce mode.

Return type:

prophys.physics.propagate_linear(generator, duration, initial, source=None)[source]#

Compute the exact state after duration of the linear system dx/dt = generator @ x + source.

Uses the exact matrix exponential (scaling and squaring with Padé approximants, the step count follows from the norm of generator * duration) with constant sources through Van Loan’s augmentation, and exact gradients through the Fréchet derivative. Covers decay chains, compartment and box models, continuous time Markov chains (generator is the rate matrix transposed) and lumped thermal networks. generator is [..., d, d], duration [...], initial and source [..., d], with at most one leading batch dimension.

Parameters:

generator (Any) – System matrix, shape [d, d] or [B, d, d].
duration (Any) – Elapsed time, scalar or [B].
initial (Any) – Initial state, shape [d] or [B, d].
source (Any) – Optional constant source term, shape [d] or [B, d].

Returns:

The state at time duration, shape [d] ([B, d] when batched).

Return type:

prophys.physics.state_space_loglik(observations, *, transition, observation, process_cov, observation_cov, initial_mean, initial_cov)[source]#

Compute the exact log likelihood of an observation series under a linear Gaussian state space model.

Uses the Kalman filter prediction error decomposition.

x[t+1] = transition @ x[t] + N(0, process_cov), y[t] = observation @ x[t] + N(0, observation_cov), x[0] ~ N(initial_mean, initial_cov). observations is [T, p]. NaN entries are missing and handled exactly, and a step with nothing observed only predicts. Use it to calibrate process models against noisy sensor series. Gradients flow to every argument.

Parameters:

observations (Any) – Observation series, shape [T, p], NaN marks missing values.
transition (Any) – State transition matrix, shape [d, d].
observation (Any) – Observation matrix, shape [p, d].
process_cov (Any) – Process noise covariance, shape [d, d].
observation_cov (Any) – Observation noise covariance, shape [p, p].
initial_mean (Any) – Mean of the initial state, shape [d].
initial_cov (Any) – Covariance of the initial state, shape [d, d].

Returns:

The log likelihood, a scalar ([B] when an argument carries one leading batch dimension).

Return type:

prophys.physics.piecewise_polynomial(x, breaks, coeffs, *, input='identity', local=False, output='identity')[source]#

Evaluate a published piecewise polynomial fit, such as a log log scaling law, at x.

In bin k (breaks[k] <= x < breaks[k+1], the end bins extrapolate) the value is output(sum_j coeffs[k, j] * v**j) with v = input(x) - (input(breaks[k]) if local else 0). input is "identity", "ln" or "log10", output "identity", "exp" or "pow10". This covers log log scaling law fits (blast parameters versus scaled distance, empirical attenuation and comfort indices) without hard coding any of them. Gradients flow to x, the breaks (local form) and every coefficient.

Parameters:

x (Any) – Evaluation points, any shape.
breaks (Any) – Bin edges, increasing, shape [K+1].
coeffs (Any) – Polynomial coefficients per bin, lowest order first, shape [K, J].
input (str) – Transform of x before the polynomial: "identity", "ln" or "log10".
local (bool) – Whether the polynomial variable is measured from each bin’s left edge.
output (str) – Transform of the polynomial value: "identity", "exp" or "pow10".

Returns:

Values with the shape of x.

Raises:

ProphysError – PHY102 for an unknown input or output transform.

Return type:

prophys.physics.network_potentials(edges, conductance, injection, *, ground=0)[source]#

Solve a grounded linear network L(conductance) phi = injection for the node potentials.

edges is [m, 2] node pairs, conductance [..., m], injection [..., n] (net inflow per node), and phi[ground] = 0. Solved exactly by Cholesky factorization of the reduced Laplacian. A disconnected network or a conductance that is not positive gives NaN potentials. The same system describes DC power flow (susceptances, angles), resistor and heat conduction networks and linearized pipe networks. Edge flows are conductance * (phi[i] - phi[j]). Gradients flow to conductances and injections.

Parameters:

edges (Any) – Node index pairs, shape [m, 2].
conductance (Any) – Edge conductances, shape [m] or [B, m].
injection (Any) – Net inflow per node, shape [n] or [B, n].
ground (int) – Index of the reference node whose potential is 0.

Returns:

Node potentials, shape [n] ([B, n] when batched).

Return type:

prophys.physics.compound_distribution(severity_pmf, *, family='poisson', rate=None, shape=None, probability=None, trials=None, n_out)[source]#

Compute the exact pmf of a compound sum of a random number of severities on the lattice 0..n_out-1.

The sum is S = X_1 + ... + X_N with a random count N.

severity_pmf is the pmf of each X on 0..M (one lattice step per index). family selects the count distribution: "poisson" (rate), "negative_binomial" (shape r and probability read as the scale beta, mean r beta) or "binomial" (trials n, probability q). Panjer’s recursion runs on rescaled values, so it stays stable for large mean counts. Value at risk and tail expectations follow from cumulative sums of the result without Monte Carlo error. Gradients flow to the severity pmf and the count parameters, except the binomial trial count, which is discrete.

Parameters:

severity_pmf (Any) – Pmf of a single severity on 0..M, shape [..., M+1].
family (str) – Count distribution, "poisson", "negative_binomial" or "binomial".
rate (Any) – Poisson mean count.
shape (Any) – Negative binomial shape r.
probability (Any) – Negative binomial scale beta or binomial probability q.
trials (Any) – Binomial number of trials n.
n_out (int) – Number of lattice points of the result.

Returns:

The pmf of S, shape [n_out] ([B, n_out] with a batch dimension).

Raises:

ProphysError – PHY102 for an unknown family or missing count parameters.

Return type:

prophys.physics.spectral_random_field(points, modes, phases, *, length_scale, std)[source]#

Evaluate a stationary random field with a prescribed covariance at points (spectral randomization method).

The field is std * sqrt(2/N) * sum_k cos(modes[k] . x / length_scale + phases[k]) at points [P, d]. Draw modes [N, d] from the normalized spectral density of the desired covariance at unit length scale (standard normal for the squared exponential covariance, a scaled multivariate t for Matérn) and phases [N] uniformly on [0, 2 pi). The field then has that covariance with variance std**2. Because the draws are inputs, samples are exact closed form functions, differentiable in the locations, draws, length_scale and std (which may carry a sample dimension).

Parameters:

points (Any) – Evaluation locations, shape [P, d].
modes (Any) – Spectral wave vectors, shape [N, d].
phases (Any) – Random phases in [0, 2 pi), shape [N].
length_scale (Any) – Correlation length, in the units of points, scalar or [B].
std (Any) – Standard deviation of the field, scalar or [B].

Returns:

Field values, shape [P] ([B, P] when length_scale or std carry a sample dimension).

Return type:

prophys.physics.directional_sweep(positions, inflow, *, direction, strength, width, growth, combine='linear')[source]#

Compute node outputs reduced by the shadowing of upstream nodes along a direction, such as turbine wakes.

Nodes are processed in order of their projection on direction (radians), ties by index. Node i outputs inflow[i] minus the combined influence of every upstream node j: c_ij = delta_ij * y_j with delta_ij = strength * (width / s)**2 * exp(-v**2 / (2 s**2)), s = width + growth * u for downstream distance u > 0 and cross distance v (zero otherwise). combine is "linear" (sum) or "rss" (root sum of squares). This is the structure of wake shadowing in turbine arrays or sheltering in canopies, with the kernel coefficients left to the model. One pass, no iteration. Gradients flow to positions, inflow and all coefficients for a fixed node order.

Parameters:

positions (Any) – Node positions, shape [n, 2].
inflow (Any) – Undisturbed input at each node, shape [n] or [B, n].
direction (Any) – Flow direction in radians.
strength (Any) – Kernel amplitude.
width (Any) – Kernel width at the source node, in position units.
growth (Any) – Linear growth of the width with downstream distance.
combine (str) – "linear" or "rss" combination of upstream influences.

Returns:

Node outputs, shape [n] ([B, n] with a batch dimension).

Raises:

ProphysError – PHY102 for an unknown combine mode.

Return type:

prophys.physics.fatigue_damage(series, *, sn_exponent, sn_coefficient, ultimate_strength=0.0)[source]#

Compute the Palmgren Miner fatigue damage of a load or stress history by rainflow counting.

Cycles are counted by the rainflow method (ASTM E1049, residual ranges count as half cycles). Each cycle contributes count * S**m / K with S the cycle range, Goodman corrected to S / (1 - mean / U) when ultimate_strength U > 0. Counting is a single bounded stack pass. Gradients flow to the series samples at the cycle reversals and to m, K and U.

Parameters:

series (Any) – Load or stress history, shape [T] or [B, T].
sn_exponent (Any) – Exponent m of the stress life (S N) curve.
sn_coefficient (Any) – Coefficient K of the stress life (S N) curve (cycles to failure N = K / S**m).
ultimate_strength (Any) – Ultimate strength U for the Goodman mean stress correction. A value of 0 or less (the default is 0) disables it.

Returns:

Accumulated damage (failure at 1), a scalar or shape [B].

Return type:

Transformations#

The transformation classes are documented on Transformations; the generic wrapper is documented here.

class prophys.transformations.FunctionTransform(fn, name='fn')[source]#

Wrap a plain Expr -> Expr callable as a transform that is not bijective.

Parameters:

fn (Callable[[Any], Expr]) – Function of one expression returning an expression.
name (str) – Label used in representations.

Example:

cube = prp.transformations.FunctionTransform(lambda v: v**3, name="cube")

Distribution base classes and copula helpers#

The distribution families are documented on Distributions. Their common interface (log_prob, sample, cdf, quantile, mean, variance) is described once on Distribution.

class prophys.distributions.Distribution[source]#

Base class of all probability distributions.

A distribution is not itself an expression. It holds its parameters as expressions (so they can be Param, Input, RandomVariable or any other Expr) and builds symbolic nodes on demand: log_prob, sample, cdf, quantile and, where available in closed form, mean and variance. Wrap a distribution in a RandomVariable to use it inside a model.

The cumulative distribution cdf and its inverse quantile are mutual duals. A subclass provides whichever it has in closed form, and the other is derived from it by monotone root finding. A subclass that provides neither has both raise NotImplementedError, as do mean and variance where a family has no closed form.

Example:

d = prp.Gaussian(10.0, 2.0)
d.log_prob(12.0).evaluate()
d.quantile(0.99).evaluate()
draws = d.sample(jax.random.PRNGKey(0), (1000,)).evaluate()

log_prob(x)[source]#

Return the log density (or log mass for discrete families) at x, as an expression.

Parameters:: x (Any) – Values at which to evaluate, broadcast against the parameters.
Return type:: Expr

sample(key, shape=())[source]#

Return random draws of shape shape (broadcast with the parameter shape), as an expression.

Families with a quantile sample by pushing uniforms through it, so draws are reparameterized and differentiable in the parameters. Gamma and Beta use implicit reparameterization instead, and discrete families are not differentiable in their parameters.

Parameters:

key (Any) – A jax.random PRNG key.
shape (tuple[int, ...]) – Sample shape. It ends with the parameter batch shape, or is a
axes. (prefix of it when the parameters carry extra trailing)

Return type:

mean()[source]#

Return the closed form mean, as an expression.

Return type:: Expr

variance()[source]#

Return the closed form variance, as an expression.

Return type:: Expr

cdf(x)[source]#

Return the cumulative distribution function F(x) = P(X <= x), as an expression.

Families without a closed form derive it by numerically inverting quantile, with the result clamped to [1e-6, 1 - 1e-6].

Parameters:: x (Any)
Return type:: Expr

quantile(q)[source]#

Return the quantile function (inverse CDF) F^{-1}(q) for q in (0, 1), as an expression.

Families without a closed form derive it by numerically inverting cdf from the bracket [-1, 1], expanded as needed. A q the expanded bracket cannot reach gives NaN.

Parameters:: q (Any)
Return type:: Expr

params()[source]#

Collect all Param leaves reachable from the distribution’s parameters, keyed by name.

Return type:: dict[str, Expr]

class prophys.distributions.TransformedDistribution(base, bijector)[source]#

Distribution of f(X) for X ~ base and an invertible transform f.

log_prob uses the change of variables formula log p_Y(y) = log p_X(f^-1(y)) - log|det J_f(f^-1(y))|, and sample pushes base samples through the forward transform. The transform must provide inverse and log_det_jacobian, such as Linear. No cdf or quantile is provided.

Parameters:

base (Distribution) – Distribution of the untransformed variable.
bijector – Invertible Transform applied to it.

Example:

celsius = prp.TransformedDistribution(prp.Gaussian(0.0, 1.0), prp.Linear(a=3.0, b=15.0))

prophys.distributions.corr_cholesky_from_raw(raw, d)[source]#

Return the lower Cholesky factor of a correlation matrix built from an unconstrained vector.

Each raw entry is mapped through tanh to a partial correlation in (-1, 1), ordered row by row over the strictly lower triangle, so any real vector gives a valid correlation matrix (for d = 2, rho = tanh(raw[0])). This is the same construction GaussianCopula uses for its density, so a fitted raw denotes the same correlation matrix for density and sampling.

Parameters:

raw (Any) – Unconstrained vector of length d * (d - 1) / 2.
d (int) – Dimension of the correlation matrix.

Returns:

Lower triangular (d, d) array L with L @ L.T the correlation matrix.

Return type:

prophys.distributions.correlated_uniforms(z, raw, d)[source]#

Map independent standard normals to correlated uniforms of a Gaussian copula.

The normals are correlated through corr_cholesky_from_raw() and mapped to (0, 1) by the normal CDF (clipped away from 0 and 1). This is the sampling core of GaussianCopula, reusable to couple otherwise independent draws (such as process innovations).

Parameters:

z (Any) – Independent standard normals, shape (..., d).
raw (Any) – Unconstrained correlation vector of length d * (d - 1) / 2.
d (int) – Number of coupled variables.

Returns:

Uniforms of shape (..., d) with Gaussian copula dependence.

Return type:

Processes#

class prophys.Process(n_steps, dt=1.0)[source]#

Base class of distributions over discretely indexed paths (time series).

Samples carry the time index as their last axis. Subclasses implement path_from_uniforms(), a deterministic, differentiable map from (*batch, n_steps) i.i.d. uniforms to the path, so gradients flow to every process parameter and externally supplied (such as correlated) uniforms can drive the same map.

Parameters:

n_steps (int) – Number of time steps in a path, at least 1.
dt (float) – Time step, in the model’s time unit.

Raises:

ValueError – If n_steps is less than 1.

path_from_uniforms(u)[source]#

Map uniforms per step (*batch, n_steps) in (0, 1) to paths of the same shape.

The map is deterministic and differentiable in the process parameters.

Parameters:: u (Any)
Return type:: Expr

sample(key, shape=())[source]#

Draw paths of shape (*shape, n_steps) by pushing fresh uniforms through path_from_uniforms().

Parameters:

key (Any) – A jax.random PRNG key.
shape (tuple[int, ...]) – Number and layout of paths to draw.

Return type:

class prophys.Autoregressive(init, mean, phi, innovation, n_steps, dt=1.0, geometry=None)[source]#

First order autoregressive (AR(1)) process with a pluggable innovation distribution.

x_t = mean + phi (x_{t-1} - mean) + eps_t, x_0 = init, eps_t ~ innovation (i.i.d.). Paths hold x_1 to x_n_steps.

phi in (-1, 1) gives mean reversion with serial correlation, the minimal structure under which the past matters: consecutive steps are dependent, so path dependent quantities built downstream (running totals, drawdowns, excursions over a hard threshold) have nontrivial distributions, unlike under i.i.d. draws.

innovation is any Distribution exposing quantile (for sampling) and log_prob (for path density). Its parameters may be Param instances and are calibrated from observed paths like any other (fit_distribution works unchanged, since log_prob of a (batch, n_steps) array of paths has shape (batch,)). All process parameters may be Expr instances.

geometry sets the space the state lives in. The default LinearGeometry is the real line. A Circular geometry makes the mean reversion difference and the placement of each step operate modulo the circle’s period, so the process is correct for angular state (see circular_autoregressive). The innovation always lives in the Euclidean difference space, so any real valued innovation distribution composes with either geometry.

Parameters:

init (Any) – State before the first step, x_0.
mean (Any) – Long run level the process reverts to.
phi (Any) – Lag one coefficient.
innovation (Distribution) – Distribution of the innovations of each step.
n_steps (int) – Number of steps in a path.
dt (float) – Time step.
geometry (Geometry | None) – State space, LinearGeometry by default.

Example:

demand = prp.Autoregressive(init=100.0, mean=100.0, phi=prp.Param("phi", init=0.7, bounds=(-1.0, 1.0)),
                            innovation=prp.Gaussian(0.0, 5.0), n_steps=24)
path = prp.RandomVariable("demand", demand)  # (n_samples, 24) in a model

innovations(x)[source]#

Return the innovations eps_t implied by an observed path x, shape (..., n_steps).

Each residual is residual(x_t, predicted_t) in the state geometry, where predicted_t is the location predicted one step ahead from x_{t-1} (init for the first step).

Parameters:: x (Any)
Return type:: Expr

path_from_innovations(eps)[source]#

Return the path generated by given innovations eps of shape (*batch, n_steps).

Parameters:: eps (Any)
Return type:: Expr

path_from_uniforms(u)[source]#

Map uniforms per step to paths through the innovation quantile and the recursion.

Parameters:: u (Any)
Return type:: Expr

log_prob(x)[source]#

Return the log density of observed paths (..., n_steps), shape (...).

It is the innovation log density summed over the implied innovations (conditional on init).

Parameters:: x (Any)
Return type:: Expr

mean()[source]#

Return the expected path E[x_t] for t = 1..n_steps, shape (n_steps,).

It is the AR recursion applied to the innovation distribution’s mean, placed into the state geometry. For a Circular geometry this is the wrapped mean of the linear recursion, exact only in the concentrated regime where the state stays away from the wrap seam.

Return type:: Expr

class prophys.RandomWalk(init, innovation, n_steps, dt=1.0)[source]#

Random walk x_t = x_{t-1} + eps_t that accumulates its innovations.

The phi = 1 boundary of Autoregressive, where the level mean drops out. Drift and step scale live entirely in the innovation distribution, such as Gaussian(drift * dt, sigma * sqrt(dt)).

Parameters:

init (Any) – State before the first step.
innovation (Distribution) – Distribution of the increments of each step.
n_steps (int) – Number of steps in a path.
dt (float) – Time step.

prophys.ornstein_uhlenbeck(init, mean, theta, sigma, n_steps, dt=1.0)[source]#

Build an Ornstein Uhlenbeck process dx = theta (mean - x) dt + sigma dW on a regular time grid.

Uses the exact discretization: an Autoregressive with phi = exp(-theta dt) and zero mean Gaussian innovations of standard deviation sigma sqrt((1 - exp(-2 theta dt)) / (2 theta)), so the continuous time theta and sigma, which must be positive, stay meaningful under any dt.

A factory, not a class: the derived parameters are ordinary symbolic expressions of theta and sigma, so both remain leaves that can be calibrated or optimized like any other Param.

Parameters:

init (Any) – State before the first step.
mean (Any) – Long run mean.
theta (Any) – Mean reversion rate per unit time, positive.
sigma (Any) – Diffusion coefficient, positive.
n_steps (int) – Number of steps in a path.
dt (float) – Time step.

Returns:

An Autoregressive process.

Return type:

Autoregressive

prophys.circular_autoregressive(init, mean, phi, innovation, n_steps, dt=1.0, period=360.0)[source]#

Build an Autoregressive process on the circle, for angular state such as wind direction.

The circle has circumference period (default 360, that is degrees). Mean reversion follows the shortest signed arc and each state is folded into [0, period), so the process is correct for angular quantities where the linear recursion would revert the long way around and report spurious ramps across the wrap seam.

A factory, not a class: it is Autoregressive with a Circular geometry, so every parameter stays an Expr that can be calibrated or optimized, and it composes into a JointProcess like any other.

Parameters:

init (Any) – Angle before the first step.
mean (Any) – Mean direction the process reverts to.
phi (Any) – Lag one coefficient.
innovation (Distribution) – Distribution of the angular innovations of each step.
n_steps (int) – Number of steps in a path.
dt (float) – Time step.
period (float) – Circumference of the circle (360 for degrees, 2*pi for radians).

Returns:

An Autoregressive process with a Circular geometry.

Return type:

Autoregressive

class prophys.JointProcess(processes, corr_raw)[source]#

Several autoregressive processes whose innovations are coupled through a Gaussian copula.

sample(key, shape) returns (*shape, n_steps, k), the component axis last, matching the component slicing convention of JointRandomVariable, so joint["name"] yields a (*shape, n_steps) path.

corr_raw is an unconstrained k(k-1)/2 vector (may be a Param) with the same raw correlation parameterization as GaussianCopula. It always gives a valid correlation matrix under gradient updates, and one fitted value means the same matrix for density and sampling.

Parameters:

processes (Sequence[Autoregressive]) – At least two Autoregressive processes with equal n_steps.
corr_raw (Any) – Unconstrained correlation vector of length k * (k - 1) / 2.

Raises:

ValueError – If fewer than two processes are given or their n_steps differ.

Example:

rho_raw = prp.Param("rho_raw", shape=(1,), init=jnp.zeros(1))
joint = prp.JointProcess([wind_speed_ar, wind_dir_ar], corr_raw=rho_raw)
met = prp.JointRandomVariable("met", joint, components=["speed", "direction"])

sample(key, shape=())[source]#

Draw joint paths of shape (*shape, n_steps, k), component axis last.

Parameters:

key (Any)
shape (tuple[int, ...])

Return type:

log_prob(x)[source]#

Return the log density of joint paths (..., n_steps, k), shape (...).

Parameters:: x (Any)
Return type:: Expr

class prophys.Recurrence(name, step, init, drivers, params=(), output='path')[source]#

User defined, differentiable state that evolves step by step along the time axis.

A jax.lax.scan over the time axis as a symbolic node.

step(carry, drivers_t, params) -> (carry, output_t) is a pure function of jax.numpy values: drivers_t is a tuple with one entry per driver expression (each the (..., ) slice of that driver at one time step), and params a tuple of the evaluated params expressions (constant along the path). Drivers carry time as their last axis and are broadcast against each other, and init must broadcast against their common batch shape.

output="path" evaluates to the stacked outputs of all steps with time as the last axis, output="final" to the final carry. The scan is inside the same JAX trace as everything else, so jit, grad and vmap flow through the recursion, including into init, every driver, and every parameter.

Parameters:

name (str) – Label of the node.
step (Callable[..., Any]) – Pure function step(carry, drivers_t, params) -> (carry, output_t).
init (Any) – Initial carry.
drivers (Sequence[Any]) – Path valued expressions with time as the last axis (at least one).
params (Sequence[Any]) – Expressions constant along the path, passed to step.
output (str) – "path" for the outputs of all steps or "final" for the final carry.

Raises:

ValueError – If drivers is empty or output is not "path" or "final".

Example:

def reservoir(level, drivers_t, params):
    inflow, demand = drivers_t
    (capacity,) = params
    level = jnp.clip(level + inflow - demand, 0.0, capacity)
    return level, level

level = prp.Recurrence("level", reservoir, init=50.0, drivers=[inflow, demand], params=[capacity])

class prophys.Geometry[source]#

Base class of value spaces for process state: how to take differences and wrap values.

residual(x, ref)[source]#

Return the signed difference of x from ref (x - ref on the line).

Parameters:

x (Any)
ref (Any)

Return type:

wrap(x)[source]#

Return the canonical representative of x in the value space.

Parameters:: x (Any)
Return type:: Any

class prophys.LinearGeometry[source]#: The real line: difference is subtraction, placement is the identity.

class prophys.Circular(period=6.283185307179586)[source]#

The circle of circumference period, for angular state such as directions.

The difference is the shortest signed arc in [-period/2, period/2), and placement folds values into [0, period).

Parameters:: period (float) – Circumference, such as 360 for degrees or 2*pi for radians (the default).
Raises:: ValueError – If period is not positive.

prophys.over_time(x)[source]#

Append a trailing time axis so a scalar valued expression broadcasts against paths.

Maps shape (...,) to (..., 1).

Example:

excess = flow_path - prp.over_time(capacity)  # (n, T) - (n, 1)

Parameters:: x (Any)
Return type:: Expr

prophys.time_sum(x)[source]#

Sum a path valued expression along the time axis, (..., T) -> (...,).

Parameters:: x (Any)
Return type:: Expr

prophys.time_mean(x)[source]#

Average a path valued expression along the time axis, (..., T) -> (...,).

Parameters:: x (Any)
Return type:: Expr

prophys.time_diff(x, geometry=None)[source]#

Take the first difference along the time axis, (..., T) -> (..., T-1).

The building block for magnitudes of step changes (combine with abs_, relu, threshold kernels, …).

With a geometry the difference is taken in that space. A Circular geometry gives the shortest signed step for angular paths, so with period 360 a change from 359 deg to 1 deg reads as 2, not as minus 358.

Parameters:

x (Any)
geometry (Geometry | None)

Return type:

prophys.discounted_sum(x, rate, dt=1.0)[source]#

Sum a path along the time axis with exponential discounting, (..., T) -> (...,).

Computes sum_t x_t exp(-rate * dt * t) for steps t = 0..T-1, such as the net present value of a cash flow path.

Parameters:

x (Any) – Path valued expression with time as the last axis.
rate (Any) – Continuous discount rate per unit time, may be an Expr or Param.
dt (float) – Time step.

Return type:

Domain layer#

Observations and as_observations() are documented on Observations.

class prophys.domain.RandomVariable(name, distribution)[source]#

A named source of uncertainty inside an expression, backed by a Distribution.

A RandomVariable can appear anywhere in a model graph. The compiled model marginalizes it by Monte Carlo, drawing it once per sample and feeding the draw through every expression that uses it.

Evaluating a RandomVariable requires either:

a bound concrete value in env[name] (e.g. an observation), or
a drawn sample in env[“__rv_samples__”][name], which is how the engine’s Monte Carlo marginalization feeds samples through the graph while keeping everything inside a single JAX trace.

A bound value in env[name] takes precedence over a sample. Evaluating with neither raises KeyError.

The distribution may be a scalar family or a Process (a distribution over discretely indexed paths). The engine treats both identically, and a path valued sample simply carries the index or time axis as its last axis (the Monte Carlo sample axis stays axis 0).

Parameters:

name (str) – Unique name of the variable. Leaves with the same name share one draw.
distribution – A Distribution or Process the variable is drawn from.

Example:

wind = prp.RandomVariable("wind_speed", prp.Weibull(scale=8.0, concentration=2.0))
power = rated * prp.clip((wind - 3.0) / 9.0, 0.0, 1.0) ** 3

log_prob(x)[source]#

Return the log density of the variable’s distribution at x, as an expression.

Parameters:: x (Any)
Return type:: Expr

class prophys.domain.JointRandomVariable(name, distribution, components)[source]#

A vector valued random variable whose named components are jointly distributed.

The joint distribution is e.g. a GaussianCopula, a MultivariateGaussian or a coupled JointProcess.

The joint distribution’s sample must return the component axis as the last axis. Use component(name) (or rv[“name”]) to obtain scalar or path valued leaves for individual components. All components share one underlying draw, so their dependence survives Monte Carlo marginalization. Two separately declared RandomVariable instances never can, since the engine draws each named leaf independently.

Parameters:

name (str) – Unique name of the joint variable.
distribution – Joint distribution whose samples have the component axis last.
components (Sequence[str]) – Names of the components, in the order of that axis.

Raises:

ValueError – If the component names are not unique.

Example:

met = prp.JointRandomVariable("met", copula, components=["speed", "temperature"])
speed, temp = met["speed"], met["temperature"]

component(name)[source]#

Return the named component as a random variable that shares the joint draw.

Equivalent to rv[name].

Raises:: KeyError – If name is not one of the components.
Parameters:: name (str)
Return type:: JointComponent

class prophys.domain.JointComponent(parent, index, component_name)[source]#

One named component of a JointRandomVariable.

Deliberately carries the parent’s name and distribution: the engine’s upstream detection deduplicates RandomVariable leaves by name, so all components of one joint resolve to a single draw of the joint distribution. That shared draw is what carries the dependence structure. Evaluation then slices this component’s index off the last axis.

Parameters:

parent (JointRandomVariable)
index (int)
component_name (str)

class prophys.domain.UncertainAttribute(name, distribution, unit='', domain='')[source]#

A named, documented uncertain quantity, the primary unit of meaning in prophys.

An attribute is the model output a study cares about (a noise level, a flood depth, an annual energy yield): a distribution whose parameters are expressions of geometry, parameters and upstream random variables. It carries metadata (unit, domain description) used by diagnostics, plotting, unit checking and export, in addition to the plain RandomVariable graph behavior. Attributes are bundled into a ProbabilityModel.

Parameters:

name (str) – Unique name of the attribute.
distribution – Distribution of the attribute’s value.
unit (str) – Measurement unit, parsed with prophys.units.parse() (e.g. "dB", "m/s"). Empty for dimensionless.
domain (str) – Free text description of what the attribute measures.

Raises:

prophys.units.UnitError – If unit cannot be parsed.

Example:

level = prp.UncertainAttribute(
    "noise_level",
    prp.Gaussian(source_db - 20.0 * prp.log10(distance), prp.Param("sigma", init=2.0, transform="softplus")),
    unit="dB",
    domain="A-weighted sound pressure level at each house",
)

property measurement_unit#: The parsed Unit of this attribute.

convert(values, to)[source]#

Convert values of this attribute from its unit into unit to.

Parameters:

values – Values in the attribute’s unit.
to (str) – Target unit.

Raises:

prophys.units.UnitError – When to measures a different quantity, which is the check that makes an attribute declared in "kW" impossible to read as "kWh".

upstream_random_variables()[source]#

Return all RandomVariable leaves this attribute’s distribution depends on, keyed by name.

For example the wind speed and direction an attribute’s mean is derived from.

Return type:: dict[str, RandomVariable]

class prophys.domain.AttributeInteractions(attributes, coupling)[source]#

Joint dependence between several UncertainAttribute objects, for attributes that are not independent.

Parameters:

attributes (Sequence[UncertainAttribute]) – The coupled attributes, with unique names.
coupling (Any) – A GaussianCopula or Correlation over the attributes’ marginal distributions, in the same order.

log_prob(values)[source]#

Return the joint log density of one value per attribute, as an expression.

Parameters:: values (Sequence[Any])
Return type:: Expr

class prophys.domain.ProbabilityModel(*attributes, interactions=None)[source]#

A probabilistic model: a collection of UncertainAttribute objects, optionally coupled.

This is the standardized handoff object: build it once from symbolic pieces, then call compile() to obtain a CompiledModel with a stable log_prob/sample/expectation interface for calibration, design optimization, uncertainty analysis and export.

Parameters:

*attributes (UncertainAttribute) – One or more UncertainAttribute objects.
interactions (AttributeInteractions | None) – Optional AttributeInteractions coupling some of them.

Example:

model = prp.ProbabilityModel(noise_level, annual_energy)
compiled = model.compile(mode="opt", tau=0.5, n_samples=512)

attribute(name)[source]#

Return the attribute called name.

Raises:: KeyError – If the model has no such attribute.
Parameters:: name (str)
Return type:: UncertainAttribute

params()[source]#

Collect all Param leaves of all attributes, keyed by name.

Return type:: dict[str, Param]

random_variables()[source]#

Return all upstream RandomVariable leaves referenced anywhere in the model, keyed by name.

For example the wind speed and direction feeding a noise surrogate.

compile(mode='eval', tau=1.0, n_samples=256, seed=0)[source]#

Compile this model into a CompiledModel.

mode="eval" uses hard min/max/contains everywhere in the geometry graph (exact evaluation). mode="opt" uses the smooth (soft) approximations at temperature tau, suitable for gradient based calibration and design optimization. n_samples and seed set the default Monte Carlo marginalization sample size and PRNG seed for any upstream RandomVariable leaves. Compiling first verifies the license and checks the model against the free tier (see License). That limit applies to structural model size (points, polygon vertices, segments, …), never to n_samples.

Raises:

ValueError – If mode is not "eval" or "opt", or if mode="opt" with a tau that is not positive.

Parameters:

mode (str)
tau (float)
n_samples (int)
seed (int)

prophys.domain.find_random_variables(expr)[source]#

Collect all distinct RandomVariable leaves reachable from expr, keyed by name.

Parameters:: expr (Expr) – Any expression.
Returns:: Mapping from variable name to the RandomVariable leaf.
Return type:: dict[str, RandomVariable]

Fitting and inference#

The fitting, uncertainty, sensitivity and validation functions are documented on their narrative pages (see Documented elsewhere).

prophys.engine.compile_model(model, mode, tau, n_samples, seed)[source]#

Compile a probability model. This is the function behind ProbabilityModel.compile().

Prefer model.compile(mode=..., tau=..., n_samples=..., seed=...), which supplies defaults ("eval", 1.0, 256, 0). Before building the CompiledModel the model is checked by analyze_model(), which reads PROPHYS_LICENSE again and checks the whole model against the free tier.

Parameters:

model (ProbabilityModel) – The model to compile.
mode (str) – "eval" for exact (hard) geometry, "opt" for smooth approximations suited to gradient based fitting.
tau (float) – Smoothing temperature of "opt" mode. Must be positive there.
n_samples (int) – Default Monte Carlo sample count for marginalizing upstream random variables.
seed (int) – Default PRNG seed.

Returns:

A CompiledModel.

Raises:

ValueError – If mode is not "eval" or "opt", or if mode="opt" and tau is not positive.
ProphysError – LIC103 when the model needs a license that is not active.

Return type:

CompiledModel

prophys.engine.negative_log_likelihood(compiled, observations, *, weights=None, prior_neg_log_prob=None)[source]#

Build the negative log likelihood that finetune() minimizes.

This is the single definition of “the model’s fit to these observations”: finetune() optimizes it, parameter_uncertainty() differentiates it twice at the optimum, and potential_function() samples it. Sharing one builder makes the reported standard errors those of the fit that was actually performed, including its weights, prior term and observation model.

Parameters:

compiled (CompiledModel) – The model.
observations (Mapping[str, Any]) – Attribute name to a one dimensional array of exact observations or an Observations record.
weights (Mapping[str, float] | None) – Optional weight of each attribute’s log likelihood by name (default 1).
prior_neg_log_prob (Callable[[Mapping[str, Any]], Any] | None) – Optional f(constrained_params) -> scalar added to the loss.

Returns:

A function of a full dict of unconstrained parameter values (as CompiledModel.log_prob() takes them) returning the scalar loss -sum_a weight_a * sum_i log p(y_ai) + prior.

Raises:

KeyError – If an observed attribute name is not in the model.

Return type:

Callable[[Mapping[str, Any]], Any]

prophys.engine.confidence_interval(estimate, std_error, level=0.95)[source]#

Confidence interval estimate +/- z * std_error from normal theory.

Used by every report that pairs a value with a standard error.

Parameters:

estimate (float) – The point estimate.
std_error (float) – Its standard error.
level (float) – Two sided confidence level, strictly between 0 and 1.

Returns:

(lower, upper). Both bounds are nan when std_error is negative or not finite.

Raises:

ValueError – If level is not strictly between 0 and 1.

Return type:

tuple[float, float]

Example

>>> confidence_interval(10.0, 1.0, 0.95)
(8.040036015459968, 11.959963984540032)

Operational surrogates#

How to amortize a model into a guarded operational surrogate is explained on Probabilistic Training.

prophys.amortize(model, attribute, inputs, *, params=None, n_points=2000, n_draws=256, n_samples=None, stats=('mean', 'quantile:0.05', 'quantile:0.5', 'quantile:0.95', 'cvar:0.95'), head='histogram', bins=64, tails=None, loss=None, hidden=(128, 128, 128), activation='silu', members=3, n_steps=4000, learning_rate=0.002, batch_size=256, n_validation=64, reference_draws=None, thresholds=None, margin=2.0, audit_rate=0.01, disagreement=None, check_structure=True, normalize=True, seed=0)[source]#

Train an operational surrogate of one attribute’s predictive distribution over an envelope.

The full model is run by Monte Carlo at n_points Latin hypercube points of the envelope, n_draws draws each, with vectorised calls of sample() over batches of points (at most 32768 draws per call). The draws of each point become a distribution target (target_kind="samples") for a network with a HistogramHead (or another probabilistic head), trained with pretrain() as an anchored ensemble of members. With normalize, a second network first learns the mean and log standard deviation of the draws and the distribution network is trained on the standardised draws.

The surrogate is then validated at n_validation fresh points against reference runs with reference_draws draws. For every served statistic the card records the surrogate’s error, the full model’s own Monte Carlo error at the operational sample count n_samples, and the equivalent sample count. The acceptance criterion is that the surrogate is no worse than the full model at n_samples draws. A statistic that misses it raises the warning AMO201 (the surrogate is still returned, and AmortizationCard.passed is False).

Queries through OperationalSurrogate.answer() fall back to the full model when the structure of the model changed, when audits marked the surrogate stale, outside the envelope, when the members disagree more than 1.2 times the 99.5 percent quantile of their disagreement on the validation points, near a decision threshold, and for a random audit_rate share of queries. Every fallback is logged with its reasons.

Parameters:

model (Any) – A ProbabilityModel or CompiledModel.
attribute (str) – Name of the scalar attribute to amortize.
inputs (Mapping[str, tuple[float, float]]) – Operational envelope, {parameter name: (low, high)} in model units. The inputs are Param leaves of the model (an operational input such as a wind speed is declared as a Param).
params (Mapping[str, float] | None) – Values (model units) of other parameters, held fixed. Parameters not given keep their initial values.
n_points (int) – Design points of the training data.
n_draws (int) – Full model draws per design point.
n_samples (int | None) – Operational sample count of a full query, the accuracy yardstick and the sample count of fallbacks. Default the compiled model’s n_samples, or 1024 when model is not compiled.
stats (Sequence[str]) – Served statistics: "mean", "quantile:<level>", "cvar:<level>" (upper tail, as cvar()), "prob:><value>" or "prob:<<value>".
head (str) – "histogram" (default), "quantile" (at fixed levels from 0.005 to 0.995 plus the levels of the served quantiles and CVaRs), "gaussian" or "mixture" (three components).
bins (int) – Bins of a histogram head.
tails (str | None) – Tails of a histogram or quantile head, default "exponential".
loss (Any) – Training loss, default "kl" for a histogram head, "pinball" for a quantile head and "crps" otherwise. Any pretrain() loss or mix.
hidden (Sequence[int]) – Widths of the hidden layers. Three layers of 128 plus the head stay within the free tier. Smaller networks train faster but rarely reach the accuracy of a few thousand full model draws.
activation (str) – Activation of the hidden layers.
members (int) – Ensemble members. Their disagreement is a fallback gate.
n_steps (int) – Training steps.
learning_rate (float) – Peak learning rate.
batch_size (int) – Mini batch size.
n_validation (int) – Validation points.
reference_draws (int | None) – Draws per validation point of the reference runs, default max(32 * n_samples, 4 * n_draws), rounded down to a multiple of n_samples and at least 2 * n_samples.
thresholds (Mapping[str, float] | None) – Decision thresholds by statistic, for example {"quantile:0.95": 120.0}. A query whose statistic lies within margin times the tolerated error of its threshold runs the full model.
margin (float) – Multiplier of the tolerated error around thresholds.
audit_rate (float) – Share of queries replayed on the full model as an audit.
disagreement (float | None) – Fixed disagreement threshold instead of the calibrated one.
check_structure (bool) – Recompute the training model’s structure hash at every guard and fall back when it changed.
normalize (bool) – Learn the location and scale of the attribute over the envelope with a second small network and let the distribution network learn the shape of the standardised attribute. This keeps the bins of a histogram head (and the knots of a quantile head) where the probability is at every point of the envelope, which matters when the spread changes a lot across it.
seed (int) – Seed of the design, the draws, the training and the audits.

Returns:

An OperationalSurrogate.

Raises:

prophys.diagnostics.ProphysError – (AMO101) For an unknown attribute, an envelope input that is not a scalar model parameter, an empty or inverted range, an unknown statistic or head, or a threshold on a statistic that is not served. (AMO102) For an attribute that is not scalar. (NN124) When design points give non finite draws.

Return type:

OperationalSurrogate

Example

op = prp.amortize(model, "loss", {"wind": (5.0, 25.0), "load": (0.5, 1.5)},
                  n_samples=4096, thresholds={"quantile:0.99": 12.0})
print(op.card)
ans = op.answer({"wind": [12.0, 30.0], "load": [1.0, 1.0]})
ans.route, ans.reasons      # the second row is outside the envelope

prophys.break_even(model, attribute, inputs, *, params=None, n_samples=None, n_points=2000, n_draws=256, n_steps=3000, hidden=(128, 128, 128), bins=64, members=3, batch_size=256, fallback_rate=0.05, queries=None, seed=0)[source]#

Estimate whether amortizing an attribute pays off, before training anything.

Times one full model query at the operational sample count in the middle of the envelope, the draws of up to 64 design points (extrapolated to n_points), a few training steps of a histogram head network (extrapolated to n_steps) and one surrogate query, then returns the number of queries after which the up front cost is recovered. As a rule of thumb amortizing pays when a full query costs 100 ms or more and thousands of queries follow. A model that answers in milliseconds is better served directly.

Parameters:

model (Any) – A ProbabilityModel or CompiledModel.
attribute (str) – The attribute to amortize.
inputs (Mapping[str, tuple[float, float]]) – Operational envelope, as for amortize().
params (Mapping[str, float] | None) – Fixed values of other parameters.
n_samples (int | None) – Operational sample count of a full query, default as for amortize().
n_points (int) – Planned design points.
n_draws (int) – Planned draws per design point.
n_steps (int) – Planned training steps.
hidden (Sequence[int]) – Planned hidden widths.
bins (int) – Planned histogram bins.
members (int) – Planned ensemble members.
batch_size (int) – Planned mini batch size.
fallback_rate (float) – Expected share of queries running the full model, audits included.
queries (float | None) – Expected number of queries, sets pays_off.
seed (int) – Seed of the timing runs.

Returns:

A BreakEvenEstimate.

Raises:

prophys.diagnostics.ProphysError – (AMO101) For an unknown attribute or an invalid envelope, as for amortize().

Return type:

BreakEvenEstimate

class prophys.engine.OperationalSurrogate(**kw)[source]#

A trained surrogate of a whole model’s predictive distribution, with guarded fallback.

Built by amortize(). query() returns the surrogate’s statistics, distribution() its predictive distribution, guard() the gate values, and answer() the guarded answer: the surrogate where every gate passes, the full model otherwise, with the reason of every fallback appended to log.

Parameters:: kw (Any)

model#: The full ProbabilityModel.

compiled#: Its compiled form used for fallbacks.

attribute#: The amortized attribute.

names#: Envelope parameter names, the column order of array inputs.

box#: Envelope bounds [d, 2] in model units.

surrogate#: The trained Surrogate.

card#: The AmortizationCard.

stats#: Served statistic names.

thresholds#: Decision thresholds by statistic.

margin#: Multiplier of the tolerated error around thresholds.

audit_rate#: Share of queries replayed on the full model as an audit.

log#: One record (a dict) per query row that ran the full model, with its inputs, reasons and the surrogate’s and full model’s values.

audits#: The records of log whose only reason is an audit.

stale#: True once at least five audits show squared errors averaging more than four times the validated ones. Every later query then runs the full model.

n_queries#: Query rows answered so far.

distribution(x)[source]#

Predictive distribution of the surrogate at x, batched over rows.

Parameters:: x (Any) – A mapping from envelope parameter to values, or an array [n, d] in the order of names.
Returns:: A closed form distribution of the ensemble mean output with cdf, quantile, mean, cvar, log_prob and crps methods on arrays, rescaled by the learned location and scale when the surrogate was trained with normalize.

query(x)[source]#

Surrogate statistics at x without any guard.

Parameters:: x (Any) – Queries, as for distribution().
Returns:: Statistic name to values [n].
Return type:: dict

guard(x, model=None)[source]#

Gate values and fallback reasons of a batch of queries, without answering them.

The audit draw is not part of the guard, answer() adds it.

Parameters:

x (Any) – Queries, as for distribution().
model (Any) – The current model, to compare its structure with the one the surrogate was trained on. Without it the training model’s structure is checked again when check_structure was set.

Returns:

A Guard.

Return type:

Guard

full(x, model=None)[source]#

Statistics of the full model at x with the operational sample count.

Parameters:

x (Any) – Queries, as for distribution().
model (Any) – A ProbabilityModel or CompiledModel to run instead of the training one. Envelope inputs and fixed values that are not parameters of it are dropped.

Returns:

Statistic name to values [n].

Return type:

answer(x, model=None)[source]#

Guarded answers: the surrogate where every gate passes, the full model otherwise.

Rows go to the full model when the structure of the model differs from the training model, when audits marked the surrogate stale, when the row lies outside the envelope, when the ensemble members disagree, when a statistic lies within margin times the tolerated error of its decision threshold, or when the row is drawn for an audit (a random audit_rate share). Each such row is appended to log with its reasons. Rows drawn only for an audit are also recorded in audits, which then decide whether the surrogate is stale.

Parameters:

x (Any) – Queries, as for distribution().
model (Any) – The current model, default the training model.

Returns:

An OperationalAnswer.

Return type:

OperationalAnswer

property fallback_rate: float#: Share of answered queries that ran the full model so far.

fallback_summary()[source]#

Number of fallbacks per reason so far (a row can count under several reasons).

Return type:: dict

break_even(queries=None, fallback_rate=None)[source]#

Break even estimate from the costs measured while amortizing.

Parameters:

queries (float | None) – Expected number of queries over the surrogate’s life.
fallback_rate (float | None) – Share of queries running the full model, default the observed rate (audits included), or the audit rate before any query.

Returns:

A BreakEvenEstimate.

Return type:

BreakEvenEstimate

export()[source]#

Export the surrogate network with the amortization card in its training card.

Returns:: A SurrogateExport.

class prophys.engine.OperationalAnswer(values, route, reasons, surrogate_values, seconds)[source]#

Answers to a batch of queries from OperationalSurrogate.answer().

Parameters:

values (dict)
route (list)
reasons (list)
surrogate_values (dict)
seconds (float)

values: dict#: Statistic name to values [n], from the surrogate or the full model per row.

route: list#: "surrogate" or "full" per row.

reasons: list#: Fallback reasons per row, empty for rows the surrogate served.

surrogate_values: dict#: What the surrogate predicted for every row, fallbacks included.

seconds: float#: Wall time of the whole batch.

class prophys.engine.Guard(in_support, disagreement, threshold_margin, reasons)[source]#

Gate values per row of a batch of queries, from OperationalSurrogate.guard().

Parameters:

in_support (ndarray)
disagreement (ndarray)
threshold_margin (ndarray)
reasons (list)

in_support: ndarray#: Whether each row lies inside the trained envelope.

disagreement: ndarray#

the mean distance between member and ensemble quantiles over half the ensemble’s 16 to 84 percent range (zero with one member).

Type:: Ensemble disagreement per row

threshold_margin: ndarray#: Smallest distance to a decision threshold in units of margin times the tolerated error. Below one means too close to decide (inf without thresholds).

reasons: list#: Fallback reasons per row, an empty tuple when the surrogate may serve it. Audits are not included.

property route: list[str]#: "surrogate" or "full" per row.

class prophys.engine.AmortizationCard(attribute, inputs, structure_hash, data_hash, n_points, n_draws, n_samples, reference_draws, n_validation, stats, surrogate_rmse, monte_carlo_rmse, n_equivalent, accepted, crps, crps_floor, coverage, pit, disagreement_threshold, data_seconds, train_seconds, full_query_seconds, surrogate_query_seconds)[source]#

Training and acceptance record of an OperationalSurrogate.

Errors compare the surrogate with reference runs of the full model that use reference_draws draws per validation point. The Monte Carlo error is that of the full model itself at n_samples draws, the sample count an operational query would use, estimated from chunks of n_samples reference draws. A statistic passes the accuracy criterion when the surrogate’s error is no larger than that Monte Carlo error, which is when its equivalent sample count n_equivalent is at least n_samples.

Parameters:

attribute (str)
inputs (dict)
structure_hash (str)
data_hash (str)
n_points (int)
n_draws (int)
n_samples (int)
reference_draws (int)
n_validation (int)
stats (list)
surrogate_rmse (dict)
monte_carlo_rmse (dict)
n_equivalent (dict)
accepted (dict)
crps (float)
crps_floor (float)
coverage (dict)
pit (list)
disagreement_threshold (float)
data_seconds (float)
train_seconds (float)
full_query_seconds (float)
surrogate_query_seconds (float)

attribute: str#: The amortized attribute.

inputs: dict#: Operational envelope, {parameter: [low, high]} in model units.

structure_hash: str#: Structure hash of the model at training time (prophys.export.structure_hash()).

data_hash: str#: SHA256 digest of the design points and training draws.

n_points: int#: Design points the surrogate was trained on.

n_draws: int#: Monte Carlo draws of the full model per design point.

n_samples: int#: Operational sample count of the full model, the accuracy yardstick.

reference_draws: int#: Draws per validation point of the reference runs.

n_validation: int#: Validation points.

stats: list#: Served statistics, such as "quantile:0.95".

surrogate_rmse: dict#: Root mean square error of the surrogate per statistic, reference noise removed.

monte_carlo_rmse: dict#: Root mean square Monte Carlo error of the full model at n_samples draws, per statistic.

n_equivalent: dict#: Equivalent full model sample count per statistic, the draws at which the full model’s Monte Carlo error equals the surrogate’s error (inf for a surrogate error of zero).

accepted: dict#: Whether the surrogate is at least as accurate as the full model at n_samples draws, per statistic.

crps: float#: Mean CRPS of the surrogate’s predictive distribution against at most 4096 reference draws per point.

crps_floor: float#: Mean CRPS of the reference distribution against itself, the lowest reachable value.

coverage: dict#: Share of reference draws inside the surrogate’s central intervals, keyed "50", "90" and "98" by nominal percent.

pit: list#: PIT histogram of the reference draws, ten bins (0.1 each when calibrated).

disagreement_threshold: float#: Ensemble disagreement above which a query falls back (inf for a single member).

data_seconds: float#: Wall time of generating the training draws.

train_seconds: float#: Wall time of training the surrogate.

full_query_seconds: float#: Measured wall time of one full model query at n_samples draws in the middle of the envelope.

surrogate_query_seconds: float#: Measured wall time of guarding one query in the middle of the envelope.

property passed: bool#: Whether every served statistic passes the accuracy criterion.

class prophys.engine.BreakEvenEstimate(full_query_seconds, surrogate_query_seconds, fallback_rate, data_seconds, train_seconds, effective_query_seconds, speedup, break_even_queries, queries, pays_off, recommendation)[source]#

Whether amortizing a model pays off, from measured or estimated costs.

The up front cost is data generation plus training. Each served query then costs the surrogate query plus fallback_rate times the full query (fallbacks and audits). Amortizing pays off after break_even_queries queries. The recommendation advises keeping the full model when a full query takes under 100 ms and the speedup is below ten.

Parameters:

full_query_seconds (float)
surrogate_query_seconds (float)
fallback_rate (float)
data_seconds (float)
train_seconds (float)
effective_query_seconds (float)
speedup (float)
break_even_queries (float)
queries (float | None)
pays_off (bool | None)
recommendation (str)

full_query_seconds: float#: Wall time of one full model query.

surrogate_query_seconds: float#: Wall time of one guarded surrogate query.

fallback_rate: float#: Assumed or observed share of queries that run the full model, audits included.

data_seconds: float#: Wall time of generating the training draws.

train_seconds: float#: Wall time of training.

effective_query_seconds: float#: Expected cost of one query through the operational surrogate.

speedup: float#: Full model query time over effective_query_seconds.

break_even_queries: float#: Queries after which the up front cost is recovered, inf if never.

queries: float | None#: Expected number of queries, if given.

pays_off: bool | None#: Whether queries exceeds break_even_queries, None without queries.

recommendation: str#: Verdict in one sentence.

prophys.engine.FALLBACK_REASONS = ('structure_changed', 'stale', 'outside_support', 'disagreement', 'near_threshold', 'audit')#

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

Results#

class prophys.engine.CalibrationResult(params, raw_params, loss_history, final_loss, n_steps, diagnostics=<factory>, surrogate_priors=<factory>)[source]#

Result of calibrate(), finetune() or fit_distribution().

Pass raw_params on to any evaluation (compiled.expectation(name, params=fit.raw_params)), to parameter_uncertainty() or to export_model(). Read params for the fitted values themselves.

Parameters:

params (dict[str, Any])
raw_params (dict[str, Any])
loss_history (list[float])
final_loss (float)
n_steps (int)
diagnostics (list)
surrogate_priors (dict)

params: dict[str, Any]#: Fitted parameters in constrained (graph visible) space.

raw_params: dict[str, Any]#: Fitted parameters in unconstrained (optimizer) space, frozen ones included at the values they were held at.

loss_history: list[float]#: Loss (negative log-likelihood plus any prior term) at every step, before that step’s update.

final_loss: float#: Last entry of loss_history (nan when no step was taken).

n_steps: int#: Number of optimizer steps taken.

diagnostics: list#: Convergence and parameter drift findings as prophys.diagnostics objects. Those that are not errors were also emitted as warnings.

surrogate_priors: dict#: Pretrained priors of every surrogate in the model, filled only when finetune() ran with surrogate_feedback.

property converged: bool#: Whether the fit converged, that is no “did not converge” (CAL201) or “diverged” (CAL202) finding.

surrogate_shift(surrogate)[source]#

Shift of a fine tuned surrogate away from its pretrained prior.

Parameters:: surrogate – A pretrained surrogate in the fitted model.
Returns:: Per leaf name, half the squared Mahalanobis distance in nats of the fitted value from the pretrained mean under the pretrained precision. The linear head is reported as a whole under the key "head". Leaves missing from raw_params are left out.
Return type:: dict[str, float]

class prophys.engine.OptimizationResult(params, loss_history, final_objective, n_steps)[source]#

Result of a design optimization with optimize().

Parameters:

params (dict[str, Any])
loss_history (list[float])
final_objective (float)
n_steps (int)

params: dict[str, Any]#: Optimized parameter values by name, in constrained space (the values the model uses).

loss_history: list[float]#: Penalized objective at every step, before that step’s update.

final_objective: float#: Last entry of loss_history (nan when no step was taken). It includes the constraint penalty, so it equals the objective only when every constraint is satisfied.

n_steps: int#: Number of optimizer steps taken.

Data input#

The readers of prophys.io are documented on Reading Data; the parsing conventions they share are listed here.

prophys.io.tabular.MISSING_TOKENS = frozenset({'', '-', '--', '?', 'missing', 'n/a', 'na', 'nan', 'none', 'null'})#: Strings a CSV field may hold to mean “not measured”. Compared case-insensitively after stripping whitespace.

prophys.io.tabular.CENSORED_PREFIXES = ('<=', '≤', '<', 'n.d.', 'nd', 'bdl')#: Prefixes marking a left-censored (“below detection limit”) record, as analytical laboratories conventionally write them. Longer prefixes come first so "<=0.5" is not read as "<" followed by "=0.5".

Plotting#

prophys.plot.distribution(compiled, attribute_name, n_samples=4096, seed=0, ax=None, bins=60)[source]#

Plot a density histogram of Monte Carlo samples of one attribute of a compiled model.

Parameters:

compiled – A CompiledModel.
attribute_name (str) – Name of the attribute to sample.
n_samples (int) – Number of samples drawn.
seed (int) – Seed of the PRNG key.
ax – Existing Axes to draw into. A new figure is created by default.
bins (int) – Number of histogram bins.

Returns:

The matplotlib Axes (use ax.figure for the parent Figure).

prophys.plot.field(field_obj, ax=None, cmap='viridis')[source]#

Plot a scalar Field as a heatmap in its own frame’s coordinates.

Parameters:

field_obj – A scalar Field. It is evaluated with no bound inputs, so its values must not depend on unbound Input or RandomVariable leaves.
ax – Existing Axes to draw into. A new figure is created by default.
cmap (str) – Matplotlib colormap name.

Returns:

The matplotlib Axes.

prophys.plot.risk_map(expr_fn, grid, ax=None, cmap='magma')[source]#

Evaluate a spatial expression at every point of a grid and plot it as a heatmap.

Parameters:

expr_fn – Callable expr_fn(PointList) -> Expr returning one value per point, typically closing over a calibrated model, e.g. lambda pts: pts.closest_distance(houses).
grid – A Grid.
ax – Existing Axes to draw into. A new figure is created by default.
cmap (str) – Matplotlib colormap name.

Returns:

The matplotlib Axes.

prophys.plot.calibration(result, ax=None)[source]#

Plot the loss history of a calibration or optimization run.

Parameters:

result – A result with a loss_history, such as a CalibrationResult or OptimizationResult.
ax – Existing Axes to draw into. A new figure is created by default.

Returns:

The matplotlib Axes.

prophys.plot.graph(expr, filename=None)[source]#

Draw the node structure of an expression graph with graphviz.

Requires the graphviz Python package.

Parameters:

expr (Any) – The expression to draw.
filename (str | None) – If given, also render the graph to this file name (in graphviz’s default output format, which sets the file suffix).

Returns:

The graphviz.Digraph.

Units#

How units are declared and checked is explained on Units & Dimensions.

class prophys.units.Unit(symbol, dimension=Dimension(length=0.0, mass=0.0, time=0.0, current=0.0, temperature=0.0, amount=0.0, luminosity=0.0), scale=1.0, offset=0.0, tag=None)[source]#

A unit of measurement: a symbol, a dimension, and a scale to SI.

scale is the multiplier taking a value in this unit to the SI base unit of its dimension: km has scale 1000, mm has scale 0.001. offset handles the two affine temperature scales, where conversion is not a pure multiplication. It is zero for every other unit.

tag distinguishes dimensionless units that must not be interchanged, such as currencies, counts and arbitrary indices. Two tagged units are compatible only when their tags match, so EUR never silently becomes USD.

Parameters:

symbol (str)
dimension (Dimension)
scale (float)
offset (float)
tag (str | None)

property is_affine: bool#: Whether conversion involves an offset as well as a scale, which makes the unit unusable in a product or a ratio (there is no meaningful “degrees Celsius per metre”).

class prophys.units.Dimension(length=0.0, mass=0.0, time=0.0, current=0.0, temperature=0.0, amount=0.0, luminosity=0.0)[source]#

Exponents of the seven SI base dimensions.

Compared by value, so m/s and km/h share a dimension while m and s do not. Exponents are rationals in principle but are stored as floats so that sqrt(m) (which appears in spectral density units) remains representable.

Parameters:

length (float)
mass (float)
time (float)
current (float)
temperature (float)
amount (float)
luminosity (float)

property is_dimensionless: bool#: True when every base dimension exponent is zero (a pure number).

exception prophys.units.UnitError[source]#: Raised when two quantities’ units are incompatible, or a unit string cannot be parsed.

prophys.units.DIMENSIONLESS#: The dimension of a pure number.

prophys.units.parse(text)[source]#

Parse a unit expression such as "m", "km/h", "kg*m/s^2".

An empty string, None, "1" and "-" all give the dimensionless unit, so an unannotated model keeps working unchanged. A Unit passes through, which lets every API accept either form.

Parameters:: text (str | Unit | None)
Return type:: Unit

prophys.units.convert(value, from_, to)[source]#

Convert value from unit from_ to unit to.

Works elementwise on anything supporting arithmetic (a Python float, a NumPy array, a JAX array or a symbolic Expr), so a converted quantity can go straight back into a model graph. Affine temperature scales are converted through kelvin.

Raises:

UnitError – When the units are incompatible or cannot be parsed.

Parameters:

value (Any)
from_ (str | Unit | None)
to (str | Unit | None)

Return type:

prophys.units.conversion_factor(from_, to)[source]#

Return the multiplier taking a value from unit from_ into unit to.

Defined only for purely multiplicative units. An affine scale (degC, degF) has no single factor, so convert() must be used for those.

Parameters:

from_ (str | Unit | None)
to (str | Unit | None)

Return type:

float

prophys.units.compatible(a, b)[source]#

Return whether two units share a dimension and a tag, and so convert into each other.

Unparseable units are reported as incompatible rather than raising.

Parameters:

a (str | Unit | None)
b (str | Unit | None)

Return type:

bool

prophys.units.check_compatible(a, b, context='')[source]#

Raise UnitError unless a and b are compatible.

context is prefixed to the message. Pass what was being combined so the error names the model construct rather than only the units.

Parameters:

a (str | Unit | None)
b (str | Unit | None)
context (str)

Return type:

None

prophys.units.register_unit(unit, *, aliases=())[source]#

Add a unit to the registry under its symbol and any aliases.

Use it for domain units this module does not ship, such as a currency, a counting unit or an index specific to a site:

register_unit(Unit("CHF", tag="currency"))

Parameters:

unit (Unit)
aliases (Iterable[str])

Return type:

None

prophys.units.describe(text)[source]#

Return a description of a unit on one line: its symbol, dimension and SI scale, plus tag and offset if any.

Used in error messages and in ModelPackage metadata.

Parameters:: text (str | Unit | None)
Return type:: str

prophys.units.same_unit(units)[source]#

Return the first of a collection of unit labels after checking that all are compatible with it.

units may be a mapping, whose values are the labels, or an iterable. An empty collection gives the dimensionless unit.

Raises:: UnitError – When any label is incompatible with the first.
Parameters:: units (Mapping[str, Any] | Iterable[Any])
Return type:: Unit

Documented elsewhere#

These APIs have their reference documentation on the narrative page that explains them, so each object has exactly one home:

Page	Covers
Structures	`prophys.structures` — points, lines, routes, rays, polygons, ellipses, polyhedra, fields, grids, meshes and networks
Transformations	`Transform` and the transformation classes of `prophys.transformations`
Distributions	The distribution families of `prophys.distributions`
Units & Dimensions	How units are parsed, converted and checked (reference above)
Observations	`Observations` — censored, interval and missing records
Reading Data	`prophys.io` — CSV, tables, wind roses, rasters, NetCDF
Training & Design Optimization	`calibrate()`, `finetune()`, `fit_distribution()`, `optimize()`, `CompiledModel`
Uncertainty & Error	Monte-Carlo error (`assess()`) and parameter uncertainty (`parameter_uncertainty()`, `propagate()`, `profile_likelihood()`)
Sensitivity Analysis	`local_sensitivity()`, `sobol_sensitivity()`
Model Validation	`information_criteria()`, `waic()`, `loo()`, `compare()`, `calibration()`, `crps()`, `posterior_predictive_check()`, `report_card()`
Posterior Inference	`prophys.inference` — NUTS posteriors over a compiled model
Export Format	`ModelPackage`, `RunManifest` and the reproducibility helpers

Learned surrogates#

class prophys.nn.Surrogate(name, inputs, output='', layers=(), *, outputs=None, hidden=None, activation='silu', monotone=None, standardize=True, members=1, seed=0)[source]#

A trainable neural network surrogate with named, unit checked inputs, usable in a model.

Calling the surrogate on expressions (s(v=wind, rho=density)) returns a graph node, so a surrogate is used inside a model like any other expression. Every trainable weight is an ordinary Param named <name>/<layer>/<tensor> (e.g. "power/L0/W"), so pretraining with pretrain(), joint fine tuning with finetune(), parameter uncertainty, sensitivity analysis and export treat the weights like any other parameter.

Free tier: up to 4 weighted layers built from Dense, Residual, Fourier, RFF, LayerNorm and the output heads, plus any Physics layers. The automatically appended output layer, every dense layer of a Residual and every head with weights count. Larger or other architectures (attention, LoRA, monotone, convex, Lipschitz layers) warn on construction (LIC101) and need a license to evaluate.

Parameters:

name (str) – Surrogate name, the prefix of every parameter name.
inputs (Sequence[str] | Mapping[str, str]) – Input names, or a mapping of input name to unit string ({"v": "m/s", "rho": "kg/m^3"}).
output (str) – Unit of the single output.
layers (Sequence[L.Layer]) – Layer sequence from prophys.nn.layers. A linear output layer of the output width is appended when the stack does not end in a head. Default: two Dense(64, activation) layers.
outputs (Sequence[str] | Mapping[str, str] | None) – Several named outputs instead of one: a list of names (sharing unit output) or a mapping of name to unit.
hidden (Sequence[int] | None) – Shortcut for layers: widths of Dense hidden layers with activation. Give either layers or hidden.
activation (str) – Activation of the hidden layers.
monotone (Mapping[str, int] | None) – Input name to +1 (increasing) or -1 (decreasing). The network must then be built from Monotone, Bounded and monotone activations.
standardize (bool) – Insert input and output standardization, whose constants pretrain() sets from the training data.
members (int) – Number of anchored ensemble members, evaluated together. Their spread is the surrogate’s epistemic uncertainty.
seed (int) – Seed of the weight initialization.

Raises:

prophys.diagnostics.ProphysError – With an NN1xx code for an invalid architecture (unknown layer, bad widths, a layer that breaks declared monotonicity or convexity, both layers and hidden given) or an unparseable unit.

Example

>>> import prophys as prp
>>> from prophys import nn
>>> power = nn.Surrogate("power", {"v": "m/s", "rho": "kg/m^3"}, "kW", hidden=(64, 64))
>>> nn.pretrain(power, {"v": v_scada, "rho": rho_scada}, p_scada)
>>> energy = prp.UncertainAttribute("energy", prp.Gaussian(power(v=wind, rho=1.225), 50.0))

param_leaves()[source]#

Return every trainable Param leaf of the network, including parameters of compiled physics layers.

Return type:: list[Param]

physics_exprs()[source]#

Return the output expressions of the physics layers that run outside the network program.

Return type:: list[Expr]

params()[source]#

Return every trainable Param leaf by name, including those of physics layers.

Return type:: dict[str, Param]

default_values()[source]#

Return the current value of every trainable leaf by name, in unconstrained space.

Return type:: dict[str, Any]

head_names()[source]#

Return the parameter names of the final linear layer (weight, then bias), or an empty list.

Return type:: list[str]

head_matrix(values)[source]#

Return the weight [K, H] and bias [K] of the final linear layer from leaf values by name.

Parameters:: values (Mapping[str, Any])
Return type:: tuple[ndarray, ndarray]

set_values(values)[source]#

Replace the initial values of named leaves (unconstrained values by name), e.g. after pretraining.

Unknown names are ignored.

Parameters:: values (Mapping[str, Any])
Return type:: None

set_constant(segment, value)[source]#

Set a frozen constant such as "in/shift" for every ensemble member.

Raises:

prophys.diagnostics.ProphysError – (NN101) If the surrogate has no such constant.

Parameters:

segment (str)
value (Any)

Return type:

None

constant(segment)[source]#

Return the value of a frozen constant such as "in/scale", one row per member.

Parameters:: segment (str)
Return type:: ndarray

output_scale()[source]#

Return the scale of the output standardization (one when there is none).

Return type:: Any

set_standardization(x_mean, x_std, y_mean=None, y_std=None)[source]#

Set the input and output standardization constants (nothing happens when standardize is off).

Inputs are mapped to (x - x_mean) / x_std and outputs back by y * y_std + y_mean. pretrain() calls this with the moments of the training data.

Parameters:

x_mean (Any)
x_std (Any)
y_mean (Any)
y_std (Any)

Return type:

None

post_step(values)[source]#

Apply the projection run after every training step, the identity for networks.

Parameters:: values (Mapping[str, Any])
Return type:: Mapping[str, Any]

schedule(step, n_steps)[source]#

Return extra evaluation environment entries for a training step, none for networks.

Parameters:

step (int)
n_steps (int)

Return type:

property n_params: int#: Number of trainable network weights, over all ensemble members.

apply(values, x, ctx=None, env=None)[source]#

Evaluate the network numerically on a batch of inputs, per ensemble member.

Parameters:

values (Mapping[str, Any]) – Unconstrained leaf values by name overriding the current ones (None or {} for the current weights).
x (Any) – Raw (unstandardized) inputs [B, in] in the declared input order and units, or [B, rows, cols] tokens for a FieldSurrogate.
ctx (Any) – Flat context tokens of a FieldSurrogate with a context.
env (Mapping[str, Any] | None) – Extra evaluation environment for physics layers.

Returns:

Raw head outputs [E, B, out] ([E, B, rows, cols] for a field output), with E the number of members.

Return type:

features(values, x, ctx=None)[source]#

Return the features [E, B, H] feeding the final linear layer, in standardized units.

Raises:

prophys.diagnostics.ProphysError – (NN301) Unless the surrogate is a single network program ending in a linear Dense head.

Parameters:

values (Mapping[str, Any])
x (Any)
ctx (Any)

Return type:

derivatives(values, x, directions=None, ctx=None)[source]#

Compute first and second directional derivatives of the outputs with respect to the inputs.

Exact Taylor mode derivatives in the native network interpreter, as used by physics informed pretraining.

Parameters:

values (Mapping[str, Any] | None) – Unconstrained leaf values by name (None for the current weights).
x (Any) – Raw inputs [B, in].
directions (Any) – Directions [K, in]. Defaults to every unit input axis.
ctx (Any) – Flat context tokens, if the surrogate takes a context.

Returns:

(first, second), each [E, B, K, out].

Raises:

prophys.diagnostics.ProphysError – (NN101) If the surrogate has a Physics layer outside the network program.

Return type:

tuple[Any, Any]

named(*args, **kwargs)[source]#

Return one point prediction node per named output, keyed by output name.

Parameters:

args (Any)
kwargs (Any)

Return type:

dict[str, Expr]

raw(*args, **kwargs)[source]#

Return the node of the full raw head output (every quantile, the mean and scale, mixture logits).

Parameters:

args (Any)
kwargs (Any)

Return type:

SurrogateCall

members_output(*args, **kwargs)[source]#

Return the node of one point prediction per ensemble member, on the last axis.

Parameters:

args (Any)
kwargs (Any)

Return type:

SurrogateCall

distribution(*args, **kwargs)[source]#

Return the predictive distribution of a probabilistic head as a prophys distribution.

A Gaussian for a GaussianHead, a Mixture for a MixtureHead, a QuantileDistribution for a QuantileHead and a BarDistribution for a HistogramHead, usable wherever a distribution is, e.g. in an UncertainAttribute. The outputs of ensemble members are averaged first (for a histogram head that is the mixture of the members’ distributions).

Raises:

prophys.diagnostics.ProphysError – (NN110) For any other head.

Parameters:

args (Any)
kwargs (Any)

epistemic(*args, **kwargs)[source]#

Return the point prediction plus a draw of the surrogate’s own (epistemic) uncertainty.

The uncertainty enters through a random variable named <name>/epistemic: a random ensemble member when there are several, otherwise a Gaussian perturbation with the standard deviation of the fitted head posterior. Monte Carlo evaluation of the model then propagates the surrogate’s uncertainty, and sensitivity analyses report its share.

Parameters:

args (Any)
kwargs (Any)

Return type:

in_support(*args, **kwargs)[source]#

Return a node that is one inside the region the surrogate was trained on and zero outside.

Parameters:

args (Any)
kwargs (Any)

Return type:

support_distance(*args, **kwargs)[source]#

Return a node with the distance to the training data relative to the support threshold.

Values above one flag extrapolation.

Parameters:

args (Any)
kwargs (Any)

Return type:

lipschitz_bound()[source]#

Return an upper bound on the maximum norm Lipschitz constant in raw input and output units.

Available when every layer is a Lipschitz layer or a tanh or relu activation, None otherwise.

Return type:: float | None

class prophys.nn.FieldSurrogate(name, rows, cols, output='', layers=(), *, context=None, standardize=True, members=1, seed=0)[source]#

A surrogate whose input is a set or field of tokens rather than named scalars.

The input is rows tokens with cols features each, with an optional context set of shape context (rows, cols). Tokens come from field_tokens(), point_tokens() or sample_set(), so the input can be an interpolated field, a point cloud or the Monte Carlo sample of an upstream quantity, which makes the surrogate a map from distributions to outputs. Call it as s(tokens) or s(tokens, context). The other methods are those of Surrogate.

Parameters:

name (str) – Surrogate name, the prefix of every parameter name.
rows (int) – Number of input tokens.
cols (int) – Features per token.
output (str) – Output unit.
layers (Sequence[L.Layer]) – Layer sequence. Defaults to Dense(64, "gelu"), PoolMean, Dense(1, "identity"). Token wise layers apply per token and a pooling layer reduces the set.
context (tuple[int, int] | None) – (rows, cols) of an optional context token set.
standardize (bool) – Insert input and output standardization.
members (int) – Number of ensemble members.
seed (int) – Seed of the weight initialization.

raw(tokens, context=None)[source]#

Return the node of the full raw head output for tokens and an optional context.

Parameters:

tokens (Any)
context (Any)

Return type:

SurrogateCall

distribution(tokens, context=None)[source]#

Return the predictive distribution of a probabilistic head for tokens.

Parameters:

tokens (Any)
context (Any)

class prophys.nn.TreeEnsemble(name, inputs, output='', *, n_trees=100, depth=4, outputs=1, learn_features=False, monotone=None, link='identity', straight_through=True, members=1, seed=0)[source]#

A differentiable tree ensemble surrogate with named, unit checked inputs.

The prediction is a base value plus the sum of n_trees complete, axis aligned trees of depth depth, optionally passed through a logistic link. Splits are softsign relaxations with a learned sharpness per node: in a model compiled with mode="opt" routing is smooth, and with mode="eval" it is hard, exactly as in gradient boosting libraries. Every threshold, sharpness, leaf value and the base are ordinary Param leaves named <name>/thresholds, <name>/log_scale, <name>/leaves and <name>/base (plus <name>/logits with learn_features), so a tree ensemble is pretrained with pretrain() and fine tuned, calibrated and analysed like any other part of the model. It supports the same graph node methods as Surrogate (call, members_output(), epistemic(), in_support(), support_distance()). Trained XGBoost and LightGBM models import with from_xgboost() and from_lightgbm().

The free tier evaluates ensembles of up to 20 trees of depth at most 4 with fixed split features (note the default n_trees=100). A larger ensemble, or one with learned split features, emits a LIC101 warning at construction and is refused at evaluation without a license.

Parameters:

name (str) – Surrogate name, the prefix of every parameter name.
inputs (Sequence[str] | Mapping[str, str]) – Input names, or a mapping from input name to unit string.
output (str) – Output unit.
n_trees (int) – Number of trees.
depth (int) – Depth of every tree, 1 to 12.
outputs (int) – Number of output columns.
learn_features (bool) – Learn the split feature of every node as a softmax over the inputs instead of fixing it at random.
monotone (Mapping[str, int] | None) – Mapping from input name to +1 (increasing) or -1 (decreasing). Enforced after every optimizer step by projecting the leaves (see post_step()). Requires fixed split features.
link (str) – "identity" or "logistic", applied to the summed margin.
straight_through (bool) – During smooth training, route hard in the forward pass and take softsign gradients, so the trained ensemble is exactly the hard ensemble it deploys as.
members (int) – Number of independently initialised ensemble members.
seed (int) – Seed of the split features and initial values.

Raises:

prophys.diagnostics.ProphysError – For an invalid depth, tree or output count, an unknown link, a bad monotone entry, or monotone with learn_features.

Example

trees = nn.TreeEnsemble("rail_break", {"temp": "degC", "age": "yr"}, "1/km/yr",
                        n_trees=20, depth=4, monotone={"age": 1})
nn.pretrain(trees, records, breaks)
rate = trees(temp=t, age=a)          # an expression node

param_leaves()[source]#

Every trainable Param leaf of the ensemble, as a list.

Return type:: list[Param]

physics_exprs()[source]#

Physics layer expressions of the surrogate, always empty for trees.

Return type:: list

params()[source]#

Trainable Param leaves by full name (<name>/leaves and so on).

Return type:: dict[str, Param]

default_values()[source]#

Current value of every leaf, by name (pretrained values once trained).

Tree leaves carry no constraining transform, so these are both the unconstrained and the constrained values.

Return type:: dict[str, Any]

set_values(values)[source]#

Replace the current values of named leaves, ignoring unknown names.

Parameters:: values (Mapping[str, Any]) – Mapping from full leaf name to new value.
Return type:: None

head_names()[source]#

Names of the head leaves, the leaf values and the base.

These are the leaves fine tuning trains with surrogate_scope="head".

Return type:: list[str]

head_matrix(values)[source]#

Leaves as a [K, T L] matrix and the base, the linear head of hard routing.

Parameters:: values (Mapping[str, Any])
Return type:: tuple[ndarray, ndarray]

property n_params: int#: Total number of trainable scalars across all leaves and members.

output_scale()[source]#: Scale of the output standardisation, always one for trees.

set_standardization(*args, **kwargs)[source]#

Accept standardisation constants and ignore them.

Trees are invariant to monotone rescaling of their inputs, so a tree ensemble is not standardised. The method exists for the common surrogate interface used by pretrain().

Return type:: None

apply(values, x, ctx=None, env=None)[source]#

Evaluate every member numerically on raw inputs.

Parameters:

values – Leaf values by name. Missing leaves take their current values.
x – Raw inputs [B, in] in declared input order.
ctx – Unused, accepted for the common surrogate interface.
env – Evaluation environment. Routing is hard unless the environment selects smooth (optimisation) mode.

Returns:

Outputs [E, B, out] for E members, after the link.

features(values, x, ctx=None)[source]#: One hot hard leaves [1, B, T L] of the first member, the features of the leaf posterior.

schedule(step, n_steps)[source]#

Temperature annealed geometrically from one to tau_opt during training.

Parameters:

step (int)
n_steps (int)

Return type:

post_step(values)[source]#

Project the leaves back onto the monotone constraints.

Applied after every optimizer step by pretrain(). It is the identity when no monotone constraints were declared.

Parameters:: values – Leaf values by name.
Returns:: The values with the leaf values projected.

monotone_violation(values=None)[source]#

Largest leaf ordering violation of the monotone constraints.

Return type:: float

data_init(x, y)[source]#

Data driven initialisation before pretraining.

Sets the thresholds at random quantiles of each node’s feature, the sharpness from the feature spread and the base at the target mean. Imported ensembles only get the sharpness.

Parameters:

x (ndarray)
y (ndarray)

Return type:

None

members_output(*args, **kwargs)[source]#

Node with one point prediction per ensemble member on the last axis.

Inputs are given as for a call of the ensemble.

epistemic(*args, **kwargs)[source]#

Point prediction plus a draw of the ensemble’s own uncertainty.

The draw enters through a random variable named <name>/epistemic: with several members it selects one member at random, otherwise it draws from the Bayesian posterior of the leaf values fitted by pretrain(). See prophys.nn.Surrogate.epistemic().

in_support(*args, **kwargs)[source]#: Node that is one inside the training support and zero outside.

support_distance(*args, **kwargs)[source]#

Node with the distance to the training data relative to the support threshold.

Values above one flag extrapolation.

classmethod from_xgboost(booster, name='xgb', inputs=None, output='')[source]#

Import a trained XGBoost booster as a pretrained tree ensemble.

Supported objectives are reg:squarederror, reg:absoluteerror, reg:pseudohubererror and binary:logistic (logistic link). Shallower branches are padded to the deepest tree’s depth, so hard routing reproduces the booster’s predictions.

Parameters:

booster – An xgboost.Booster.
name (str) – Surrogate name.
inputs (Sequence[str] | Mapping[str, str] | None) – Input names (or a mapping from name to unit) in feature order. Defaults to f0, f1, ....
output (str) – Output unit.

Returns:

A TreeEnsemble holding the booster’s trees.

Raises:

prophys.diagnostics.ProphysError – For an unsupported objective.

classmethod from_lightgbm(booster, name='lgbm', inputs=None, output='')[source]#

Import a trained LightGBM booster as a pretrained tree ensemble.

Supported objectives are regression, regression_l1, huber and binary (logistic link), with numerical splits only.

Parameters:

booster – A lightgbm.Booster.
name (str) – Surrogate name.
inputs (Sequence[str] | Mapping[str, str] | None) – Input names (or a mapping from name to unit) in feature order. Defaults to f0, f1, ....
output (str) – Output unit.

Returns:

A TreeEnsemble holding the booster’s trees.

Raises:

prophys.diagnostics.ProphysError – For an unsupported objective or a categorical split.

prophys.nn.pretrain(surrogate, inputs, targets, *, context=None, loss=None, target_kind='value', sigma_ratio=0.75, beta=0.5, val_fraction=0.2, n_steps=2000, learning_rate=0.003, batch_size=256, weight_decay=0.0, anchor_strength=1.0, residuals=(), patience=None, eval_every=25, seed=0, posterior=True, prior_probes=64)[source]#

Train a surrogate on its own data before it enters a model.

Sets the input and output standardisation from the training split, then trains with AdamW (global norm clipping at 10, cosine learning rate decay) on mini batches, evaluating the validation loss every eval_every steps and keeping the weights with the lowest one. The kept weights become the surrogate’s initial values. Afterwards it records on surrogate.pretrained everything later stages use: the training card, the support region of the training inputs (for in_support()), the Bayesian posterior of a linear head (for epistemic()) and a diagonal Gauss Newton prior for fine tuning inside the full model (prophys.engine.finetune() with surrogate_feedback).

Ensembles (members > 1) are anchored to their random initial weights (Pearce et al. 2020) with strength anchor_strength. Problems found along the way (constant inputs, too few samples, divergence, overfitting, low R2, miscalibrated intervals, a validation split outside the training support) are emitted as diagnostics warnings and listed on the result.

Parameters:

surrogate – A Surrogate, FieldSurrogate or TreeEnsemble, trained in place.
inputs (Any) – A mapping from input name to values (each optionally a (values, unit) pair, converted to the declared unit), an array [N, in] in declared input order, or tokens [N, rows, cols] for a field surrogate.
targets (Any) – Targets [N] or [N, out]. Field outputs [N, rows, cols] are flattened per row.
context (Any) – Context tokens of a field surrogate declared with a context set, [N, ctx_rows, ctx_cols] or one shared [1, ...].
loss (str | Mapping[str, float] | None) – Loss name, or a mapping from loss name to weight for a weighted mix such as {"hl_gauss": 1.0, "mse": 0.1}. Defaults to the head’s natural loss ("mse" for a linear head, "pinball", "gaussian", "mixture" or "hl_gauss" for probabilistic heads). Point losses: "mse", "mae", "bce", "crossentropy", "poisson". Probabilistic losses: "pinball" (quantile head), "gaussian" and "beta_nll" (Gaussian head, plain and beta weighted negative log likelihood), "mixture" (mixture head), "hl_gauss", "histogram_ce" and "kl" (histogram head: cross entropy against the HL Gauss smoothed target, against the hard histogram target, and the KL divergence to the target bin probabilities) and "crps" (closed form continuous ranked probability score of any probabilistic head). See prophys.nn.losses.
target_kind (str) – Form of targets. "value" (default): one value per row. "samples": [N, M] Monte Carlo draws per row, e.g. from sample(). Every loss is averaged over the draws, and histogram losses use the histogram of the draws (the Rao Blackwellised target). "probs": [N, K] probabilities of the bins of a HistogramHead with explicit borders. "cdf": [N, K + 1] (or [N, K - 1] inner) CDF values at those borders, e.g. from cdf(). Bin targets distill a distribution with no Monte Carlo noise at all.
sigma_ratio (float) – Standard deviation of the HL Gauss smoothing as a multiple of the mean width of the inner bins (0.75 after Farebrother et al. 2024).
beta (float) – Exponent of the "beta_nll" weighting (Seitzer et al. 2022). 0 is the plain likelihood, 1 weights the mean like a squared error.
val_fraction (float) – Fraction of rows held out for validation (no split when fewer than five rows are given).
n_steps (int) – Maximum number of optimizer steps.
learning_rate (float) – Peak learning rate of the cosine schedule.
batch_size (int) – Mini batch size (capped at the size of the training split).
weight_decay (float) – Decoupled weight decay on the surrogate’s own weights. Learnable physical constants of Physics layers are exempt.
anchor_strength (float) – Strength of the anchored ensemble penalty.
residuals (Sequence[PhysicsResidual]) – Physics informed loss terms, see PhysicsResidual.
patience (int | None) – Stop after this many evaluations without a new best validation loss. None trains for all n_steps.
eval_every (int) – Steps between validation evaluations.
seed (int) – Seed of the split, the mini batches and the prior probes.
posterior (bool) – Fit the Bayesian linear head posterior (single member surrogates with a linear head only).
prior_probes (int) – Random probes used to estimate the Gauss Newton prior precision. 0 skips the estimate.

Returns:

The PretrainResult training card.

Raises:

prophys.diagnostics.ProphysError – When inputs and targets disagree in length (NN127), an input is missing (NN121) or has the wrong number of columns (NN103), rows contain values that are not finite (NN124), fewer than two rows are given (NN126), the loss name is unknown (NN101) or a loss does not suit the head or the targets (NN128).

Return type:

PretrainResult

Example

power = nn.Surrogate("power", {"v": "m/s", "rho": "kg/m^3"}, "kW", hidden=(64, 64))
card = nn.pretrain(power, {"v": v_scada, "rho": rho_scada}, p_scada, n_steps=3000)
print(card)

A histogram head distilled from 256 model draws per input:

load = nn.Surrogate("load", ["v"], layers=[nn.Dense(64), nn.Dense(64), nn.HistogramHead(64)])
nn.pretrain(load, {"v": v_grid}, draws, target_kind="samples", loss="kl")

class prophys.nn.PretrainResult(surrogate, loss, n_train, n_validation, steps, best_step, train_loss, validation_loss, rmse, mae, r2, crps, coverage90, noise_std, data_hash, created, history, converged=True, diagnostics=<factory>, target_kind='value', nll=None, coverage=None, pit=None, tail_pinball=None, kl=None)[source]#

Training card returned by pretrain().

The same object is stored on the surrogate (surrogate.pretrained["card"]) and travels with it into an exported model package. Metrics are computed on the validation split in the raw units of the targets, or on the training split when there was no validation split.

Parameters:

surrogate (str)
loss (str)
n_train (int)
n_validation (int)
steps (int)
best_step (int)
train_loss (float)
validation_loss (float)
rmse (float)
mae (float)
r2 (float)
crps (float | None)
coverage90 (float | None)
noise_std (float | None)
data_hash (str)
created (str)
history (list)
converged (bool)
diagnostics (list)
target_kind (str)
nll (float | None)
coverage (dict | None)
pit (list | None)
tail_pinball (dict | None)
kl (float | None)

surrogate: str#: Name of the trained surrogate.

loss: str#

"mse", "mae", "bce", "crossentropy", "poisson", "pinball", "gaussian", "mixture", "hl_gauss", "histogram_ce", "kl", "crps" or "beta_nll". A weighted mix reads like "1*hl_gauss + 0.1*mse".

Type:: Loss that was minimized

n_train: int#: Rows in the training split.

n_validation: int#: Rows in the validation split (0 when fewer than five rows were given).

steps: int#: Optimizer steps actually taken (fewer than requested after early stopping or divergence).

best_step: int#: Step whose weights were kept, the one with the lowest validation loss.

train_loss: float#: Mini batch training loss at the last evaluation.

validation_loss: float#: Lowest validation loss, in the standardised units of the loss.

rmse: float#: Root mean square error of the point prediction.

mae: float#: Mean absolute error of the point prediction.

r2: float#: Coefficient of determination of the point prediction (nan for constant targets).

crps: float | None#: Continuous ranked probability score of a probabilistic head (Gaussian, quantile, mixture or histogram), None for point heads.

coverage90: float | None#: Fraction of targets inside the predicted 90% interval, from the head or from the Bayesian linear head posterior, None when unavailable.

noise_std: float | None#: Observation noise standard deviation estimated by the linear head posterior, in target units, None when no posterior was fitted.

data_hash: str#: SHA-256 digest of the training inputs and targets.

created: str#: UTC timestamp of the run, ISO 8601.

history: list#: [step, train_loss, validation_loss] rows, one per evaluation.

converged: bool = True#: False when training stalled or diverged (diagnostics NN201, NN202).

diagnostics: list#: Findings raised during training, as prophys.diagnostics objects.

target_kind: str = 'value'#

"value", "samples", "probs" or "cdf".

Type:: Form of the targets

nll: float | None = None#: Mean negative log predictive density of the validation targets (of every draw for sample targets), None for point heads and bin targets.

coverage: dict | None = None#: Central interval coverage by nominal level in percent, keys "50", "90" and "98", from the predictive distribution, None for point heads.

pit: list | None = None#

the fraction of validation targets in each of ten equal bins of the predictive CDF. A calibrated surrogate gives about 0.1 in every bin.

Type:: Probability integral transform histogram

tail_pinball: dict | None = None#: Pinball loss of the predicted 0.01 and 0.99 quantiles, {"0.01": .., "0.99": ..}, in target units. This is the tail accuracy that CVaR and exceedance probabilities depend on.

kl: float | None = None#: Mean KL divergence from the target bin probabilities to the predicted ones, for a histogram head, None otherwise.

class prophys.nn.PhysicsResidual(fn, points, weight=1.0, name='residual')[source]#

A physics informed residual loss term for pretrain().

Also available as prophys.nn.pretrain.Residual. At every training step the surrogate is evaluated at the collocation points, a Derivatives object is built from its outputs and exact first and second input derivatives, and weight * mean(fn(D) ** 2) is added to the data loss. Residuals need a surrogate without Physics layers that run as Python stages, since derivatives are taken through the network alone.

Example

A steady heat equation u'' = 0 enforced on 200 points:

pts = {"x": np.linspace(0.0, 1.0, 200)}
res = nn.PhysicsResidual(lambda D: D.dd("x"), pts, weight=0.1)
nn.pretrain(s, x_obs, u_obs, residuals=[res])

Parameters:

fn (Callable[[Derivatives], Any])
points (Any)
weight (float)
name (str)

fn: Callable[[Derivatives], Any]#

takes a Derivatives and returns residual values (any shape), which are driven towards zero.

Type:: Residual function

points: Any#: Collocation inputs, in the same form as the inputs of pretrain() (a mapping from input name to values, or an array [N, in]).

weight: float = 1.0#: Multiplier on the residual’s mean square in the loss.

name: str = 'residual'#: Label of the term.

class prophys.nn.Derivatives(names, x, u, d1, d2)[source]#

Outputs and input derivatives of a surrogate at collocation points.

Passed to the fn of a PhysicsResidual during pretrain(). u holds the outputs [N, out] and x the raw inputs [N, in]. Derivatives are with respect to the raw (unstandardised) inputs and are computed exactly by Taylor mode differentiation in the optimized native kernels, not by finite differences. Inputs are addressed by index or by input name. With an ensemble, the first member’s outputs and derivatives are used.

Parameters:

names (Sequence[str])
x (Any)
u (Any)
d1 (Any)
d2 (Any)

d(i)[source]#

First derivative of every output with respect to input i.

Parameters:: i – Input index or input name.
Returns:: Array [N, out] of du/dx_i.

dd(i)[source]#

Second derivative of every output with respect to input i.

Parameters:: i – Input index or input name.
Returns:: Array [N, out] of d2u/dx_i2 (no mixed derivatives).

laplacian(over=None)[source]#

Laplacian of every output, the sum of unmixed second derivatives.

Parameters:: over (Sequence[Any] | None) – Inputs (indices or names) to sum over. Defaults to every input.
Returns:: Array [N, out].

input(i)[source]#

Collocation values of input i in raw units.

Parameters:: i – Input index or input name.
Returns:: Array [N].

class prophys.nn.LinearPosterior(mean, v, s2, alpha, beta, log_evidence, n)[source]#

Gaussian posterior of a surrogate’s linear head y = [phi, 1] w.

Built by fit_linear_posterior(). pretrain() stores one on the surrogate (surrogate.pretrained["posterior"]) for single member surrogates with a linear head, where phi are the features feeding the last dense layer (or the one hot leaf indicators of a tree ensemble). It supplies the epistemic draw of epistemic(). The posterior covariance is stored in factored form: (alpha I + beta V diag(s2) V^T)^-1.

Parameters:

mean (ndarray)
v (ndarray)
s2 (ndarray)
alpha (float)
beta (float)
log_evidence (float)
n (int)

mean: ndarray#: Posterior mean of the weights, [H + 1, K] (last row the bias).

v: ndarray#: Right singular vectors of the design [phi, 1], [H + 1, r].

s2: ndarray#: Squared singular values of the design, [r].

alpha: float#: Prior precision of the weights, in standardised units.

beta: float#: Noise precision, in standardised target units.

log_evidence: float#: Log marginal likelihood at the selected alpha and beta.

n: int#: Number of training rows.

precision()[source]#

Posterior precision over the rows of mean.

Return type:: ndarray

variance(phi)[source]#

Epistemic variance of the head output [phi, 1] w for new features.

Parameters:: phi (Any) – Features [N, H].
Returns:: Variance [N] in standardised units, excluding observation noise.
Return type:: Any

property noise_std: float#: Observation noise standard deviation 1 / sqrt(beta), standardised units.

to_dict()[source]#

The posterior as a JSON serializable dict of lists and floats.

Return type:: dict

classmethod from_dict(d)[source]#

Rebuild a posterior from the output of to_dict().

Parameters:: d (dict)
Return type:: LinearPosterior

class prophys.nn.SupportRegion(lo, hi, center, precision, threshold)[source]#

Region of input space covered by a surrogate’s training data.

The intersection of a per input quantile box and a Mahalanobis ellipsoid with shrunk covariance. Built by fit_support(). pretrain() stores one on the surrogate (surrogate.pretrained["support"]) and it backs in_support() and support_distance(). For a field surrogate the rows are the per feature mean and standard deviation over the tokens.

Parameters:

lo (ndarray)
hi (ndarray)
center (ndarray)
precision (ndarray)
threshold (float)

lo: ndarray#: Lower corner of the box, [d].

hi: ndarray#: Upper corner of the box, [d].

center: ndarray#: Centre of the ellipsoid (the training mean), [d].

precision: ndarray#: Inverse of the shrunk training covariance, [d, d].

threshold: float#: Mahalanobis radius of the ellipsoid.

distance(x)[source]#

Distance of input rows relative to the support boundary.

The larger of the Mahalanobis distance divided by threshold and, outside the box, one plus the largest excess beyond the box relative to its width. At most one inside the support.

Parameters:: x (Any) – Rows [..., d].
Returns:: Distances [...].
Return type:: Any

to_dict()[source]#

The region as a JSON serializable dict of lists and floats.

Return type:: dict

classmethod from_dict(d)[source]#

Rebuild a region from the output of to_dict().

Parameters:: d (dict)
Return type:: SupportRegion

class prophys.nn.QuantileDistribution(levels, values, tails='exponential')[source]#

Predictive distribution of a quantile head: a piecewise linear quantile function with tails.

The quantile function interpolates the predicted quantiles values linearly between the levels. Below the lowest and above the highest level a tail carries the remaining probability. Its density continues the density of the adjacent segment, so there is no jump at the outer quantiles, and its shape follows tails:

"exponential" (default): an exponential tail, unbounded support. Tail quantiles, CVaR and the log density of observations beyond the predicted band are finite and principled.
"halfnormal": a half normal tail, lighter than exponential.
"linear": the linear extension of the quantile function to the levels zero and one, which bounds the support. This was the only behaviour before prophys 0.3. The log density outside the support is -1e30 as it was then.

The log density, CDF, quantile function, mean, variance, CVaR, CRPS and reparameterised samples are all closed form. This is what distribution() returns for a surrogate ending in a QuantileHead, and it can be used like any other prophys distribution.

Parameters:

levels (Sequence[float]) – Strictly increasing quantile levels inside (0, 1), at least two.
values (Any) – Predicted quantiles [..., len(levels)] (an array or an expression). A running maximum along the last axis makes them nondecreasing.
tails (str) – "exponential", "halfnormal" or "linear".

Raises:

prophys.diagnostics.ProphysError – When the levels or the tail name are invalid (NN101).

Example

q = nn.Surrogate("load", ["v"], layers=[nn.Dense(64, "silu"), nn.QuantileHead((0.05, 0.5, 0.95))])
nn.pretrain(q, v, load)
dist = q.distribution(wind)       # a QuantileDistribution
dist.cvar(0.99).evaluate()        # from the exponential tail

log_prob(x)[source]#

Log density at x, in closed form, and -1e30 outside the support of linear tails.

Parameters:: x (Any)
Return type:: Expr

class prophys.nn.BarDistribution(borders, probs, tails='exponential')[source]#

Predictive distribution of a histogram head: a bar (Riemann) distribution with tails.

Probability probs[k] is spread uniformly over the bin [borders[k], borders[k + 1]]. With tails="exponential" (default) or "halfnormal" the two outer bins are replaced by tails of the same mass anchored at the inner border, as in Muller et al. (2022), so the support is unbounded and every observation has a finite log density. The tail’s density at the border continues the density of the adjacent bin, so the tail length follows the learned probabilities rather than the (arbitrary) width of the outer bin. With two bins, half of each tail’s mass lies within one outer bin width instead, and a single bin becomes two tails about its centre. An exponential tail matched this way reproduces the mean excess of a Gaussian tail closely. A half normal tail is lighter. tails="none" keeps uniform outer bins and a bounded support.

The CDF, log density, quantile function, mean, variance, CVaR, CRPS and reparameterised samples are all closed form. This is what distribution() returns for a surrogate ending in a HistogramHead, and it is also a standalone distribution usable anywhere in a prophys model.

Parameters:

borders (Any) – Strictly increasing bin borders [K + 1] (or batched [..., K + 1]), an array or an expression.
probs (Any) – Bin probabilities [..., K] (an array or an expression), normalised to sum to one.
tails (str) – "exponential", "halfnormal" or "none".

Raises:

prophys.diagnostics.ProphysError – When the borders do not increase or the tail name is unknown (NN101).

Example

d = nn.BarDistribution(np.linspace(-3, 3, 65), probs)
d.quantile(0.99).evaluate(), d.cvar(0.99).evaluate()

Probabilistic training losses and distribution targets for surrogates.

These are the loss functions behind the probabilistic options of pretrain(), exposed so they can be checked and reused. They work on plain JAX arrays and broadcast over leading batch axes.

Histogram losses (Imani and White 2018, Farebrother et al. 2024): the cross entropy of predicted bin probabilities against a target distribution over the same bins. hl_gauss_probs() builds the HL Gauss target, a Gaussian around the observation integrated over each bin. samples_to_probs() gives the (optionally smoothed) histogram of Monte Carlo draws and cdf_to_probs() the bin probabilities from CDF values at the borders.
kl_divergence(), which equals the cross entropy minus the target’s entropy, so it is zero exactly at the target.
Closed form continuous ranked probability scores of the Gaussian (crps_gaussian()) and Gaussian mixture (crps_mixture()) heads, and the sample based crps_ensemble(), with the fair estimator of Ferro (2014).
beta_nll(), the Gaussian negative log likelihood with the gradient weighting of Seitzer et al. (2022).
energy_score() for multivariate samples.

prophys.nn.losses.hl_gauss_probs(y, borders, sigma)[source]#

HL Gauss target: a Gaussian around each observation integrated over the bins.

The mass of the Gaussian below the first inner border goes to the first bin and the mass above the last inner border to the last bin, so the target sums to one exactly and observations outside the borders still have a valid target (the outer bins carry the tails).

Parameters:

y (Any) – Observations [...].
borders (Any) – Increasing bin borders [K + 1] (or broadcastable [..., K + 1]).
sigma (Any) – Standard deviation of the smoothing Gaussian, scalar or [...]. Zero gives the hard one hot histogram target.

Returns:

Target probabilities [..., K].

Return type:

prophys.nn.losses.samples_to_probs(samples, borders, sigma=0.0)[source]#

Histogram of Monte Carlo draws over the bins, the Rao Blackwellised target.

Each draw contributes its (optionally HL Gauss smoothed) target and the contributions are averaged, so the result is the empirical distribution of the draws projected onto the bins.

Parameters:

samples (Any) – Draws [..., M].
borders (Any) – Bin borders [K + 1].
sigma (Any) – HL Gauss smoothing standard deviation. Zero bins the draws.

Returns:

Target probabilities [..., K].

Return type:

prophys.nn.losses.cdf_to_probs(cdf, bins=None)[source]#

Bin probabilities from CDF values at the bin borders.

Parameters:

cdf (Any) – CDF values at the inner borders [..., K - 1]. With bins given, values at every border [..., K + 1] are accepted too. The outer bins take all the mass beyond the inner borders, matching the tails of a bar distribution. Values are clipped to [0, 1] and made nondecreasing.
bins (int | None) – Number of bins K, needed only to recognise values at all K + 1 borders.

Returns:

Probabilities [..., K] that sum to one.

Return type:

prophys.nn.losses.cross_entropy(pred_probs, target_probs, eps=0.0)[source]#

Cross entropy -sum_k t_k log p_k over the last axis.

Parameters:

pred_probs (Any) – Predicted probabilities [..., K].
target_probs (Any) – Target probabilities [..., K].
eps (float) – Floor of the predicted probabilities inside the logarithm. Zero uses the smallest positive normal number of the dtype.

Returns:

The cross entropy [...].

Return type:

prophys.nn.losses.kl_divergence(target_probs, pred_probs, eps=0.0)[source]#

Kullback Leibler divergence KL(target || pred) over the last axis, with 0 log 0 = 0.

Its gradient in the prediction equals that of cross_entropy(). The value is the cross entropy minus the target entropy, clipped at zero, so it is nonnegative and zero exactly when the two agree.

Parameters:

target_probs (Any)
pred_probs (Any)
eps (float)

Return type:

prophys.nn.losses.crps_gaussian(mu, sigma, y)[source]#

Closed form CRPS of a Gaussian N(mu, sigma^2) at observation y (Gneiting et al. 2005).

Parameters:

mu (Any)
sigma (Any)
y (Any)

Return type:

prophys.nn.losses.crps_mixture(logits, mu, sigma, y)[source]#

Closed form CRPS of a Gaussian mixture (Grimit et al. 2006).

E|X - y| - E|X - X'| / 2 with both expectations in closed form.

Parameters:

logits (Any) – Unnormalised component log weights [..., C].
mu (Any) – Component means [..., C].
sigma (Any) – Component standard deviations [..., C].
y (Any) – Observations [...].

Return type:

prophys.nn.losses.crps_ensemble(samples, y, fair=True)[source]#

CRPS of an ensemble of draws [..., M] at observation y.

With fair=True the pairwise term is divided by M (M - 1) rather than M^2 (Ferro 2014), which removes the bias of a finite ensemble.

Parameters:

samples (Any)
y (Any)
fair (bool)

Return type:

prophys.nn.losses.energy_score(samples, y, fair=True)[source]#

Energy score of multivariate draws [..., M, D] at an observation [..., D].

E||X - y|| - E||X - X'|| / 2, the multivariate generalisation of the CRPS. fair works as for crps_ensemble().

Parameters:

samples (Any)
y (Any)
fair (bool)

Return type:

prophys.nn.losses.gaussian_nll(mu, sigma, y)[source]#

Gaussian negative log likelihood without the constant log sqrt(2 pi).

Parameters:

mu (Any)
sigma (Any)
y (Any)

Return type:

prophys.nn.losses.beta_nll(mu, sigma, y, beta=0.5)[source]#

The beta NLL of Seitzer et al. (2022): the Gaussian NLL weighted by stop_gradient(sigma^(2 beta)).

beta = 0 is the plain NLL and beta = 1 gives the mean the gradient of a mean squared error. The minimiser in expectation stays the true mean and variance for every beta.

Parameters:

mu (Any)
sigma (Any)
y (Any)
beta (float)

Return type:

prophys.nn.losses.pinball(levels, q, y)[source]#

Pinball (quantile) loss of predicted quantiles [..., L] at y, averaged over the levels.

Parameters:

levels (Any)
q (Any)
y (Any)

Return type:

prophys.nn.pfn.prior_borders(bins, features=4, context=64, queries=32, tasks=512)[source]#

Bin borders of a bar head at equally spaced quantiles of standardised prior targets.

Draws tasks tasks from the mechanism prior with a fixed seed, standardises each task’s query targets with its context statistics (as InContextRegressor does) and places bins + 1 borders at their quantiles (Muller et al. 2022). Deterministic.

Returns:

Strictly increasing borders [bins + 1].

Parameters:

bins (int)
features (int)
context (int)
queries (int)
tasks (int)

Return type:

ndarray

prophys.nn.field_tokens(field, channels=1, coords=True)[source]#

Token rows of a raster field: cell centre coordinates followed by the cell values.

Parameters:

field – A raster structure with values [..., ny, nx] (or [..., ny, nx, channels]), origin and cell size.
channels (int) – Number of value channels per cell.
coords (bool) – Whether to prepend the world coordinates x, y of each cell centre.

Returns:

An expression of shape [..., ny * nx, 2 + channels] (or [..., ny * nx, channels] without coordinates). Leading batch axes of the values are kept.

Raises:

ProphysError – NN123 when the field’s channel count differs from channels.

Return type:

prophys.nn.point_tokens(points, features=None)[source]#

Token rows of a point set: coordinates followed by optional per point features.

Parameters:

points – A point structure with coords [n, 2], or an expression or array of coordinates.
features (Any) – Optional features [n, f] or [n].

Returns:

An expression of shape [n, 2 + f] (the coordinates alone when no features are given).

Return type:

prophys.nn.geometry_tokens(geometry)[source]#

Token rows of a geometry, one row [x0, y0, x1, y1] per segment or edge.

Works for a line list, a polyline, a polygon (edges including the closing one) or a polygon list.

Returns:: An expression of shape [n_segments, 4].
Raises:: ProphysError – NN123 for any other structure.
Return type:: Expr

prophys.nn.sample_set(expr, n, seed=0)[source]#

Turn the Monte Carlo sample of an upstream quantity into a token set.

Returns a SampleSet whose value is [n, k]: n draws of expr (flattened to k columns), drawn independently of the model’s own Monte Carlo draw with fixed random numbers so the set is smooth in the parameters. Feeding it to a FieldSurrogate makes the surrogate a map from a distribution to an output.

Parameters:

expr (Any) – The expression to sample, typically depending on random variables.
n (int) – Sample size.
seed (int) – Seed of the fixed random numbers.

Return type:

SampleSet

Example

>>> tokens = nn.sample_set(load, 256)               # [256, 1]
>>> damage = nn.FieldSurrogate("damage", rows=256, cols=1)(tokens)

prophys.nn.load_state_dict(surrogate, state, mapping, strict=True)[source]#

Copy pretrained tensors from a PyTorch style state dictionary into a surrogate.

Sets state[mapping[leaf]] as the initial value of every mapped leaf of surrogate. Leaf names are given without the surrogate prefix, for example L0/attn/x/W. Tensors may be PyTorch tensors or arrays. Each is reshaped to its leaf’s shape after its size is checked.

Parameters:

surrogate – The surrogate to update, in place.
state (Mapping[str, Any]) – The state dictionary.
mapping (Mapping[str, str]) – Leaf name to state dictionary key, for example from transolver_mapping().
strict (bool) – When true, a key missing from state is an error. When false, it is skipped.

Returns:

Sorted full names of the leaves that were set.

Raises:

ProphysError – NN101 for an unknown leaf or (with strict) a missing key, NN103 for a size mismatch.

Return type:

list[str]

Example

>>> mapping = nn.transolver_mapping("L1", "blocks.0")
>>> nn.load_state_dict(s, torch.load("transolver.pt"), mapping)

prophys.nn.transolver_mapping(ours, theirs)[source]#

Map a TransolverBlock layer’s leaves to the keys of a Transolver block module.

The block module has ln_1, Attn, ln_2 and mlp children, as in the official implementation.

Parameters:

ours (str) – Leaf prefix of the layer, such as "L1".
theirs (str) – Key prefix of the block, such as "blocks.0".

Returns:

Mapping of leaf name to state dictionary key, for load_state_dict().

Return type:

dict[str, str]

prophys.nn.physics_attention_mapping(ours, theirs)[source]#

Map a SliceAttention layer’s leaves to the keys of a Transolver physics attention module.

Parameters:

ours (str) – Leaf prefix of the layer, such as "L1".
theirs (str) – Key prefix of the Physics_Attention_Irregular_Mesh module, such as "blocks.0.Attn".

Returns:

Mapping of leaf name to state dictionary key, for load_state_dict().

Return type:

dict[str, str]

prophys.nn.DPOT(name, grid, in_channels, in_timesteps, out_channels=None, out_timesteps=1, *, patch_size=8, embed_dim=64, depth=4, n_blocks=4, modes=32, mlp_ratio=1.0, out_layer_dim=32, normalize=False, act='gelu', time_agg='exp_mlp', output='', members=1, seed=0)[source]#

Build a DPOT field surrogate with the hyperparameters of the reference model.

Needs a license.

Parameters:

name (str) – Surrogate name, the prefix of every weight.
grid (tuple[int, int]) – Spatial grid (X, Y), multiples of patch_size.
in_channels (int) – Channels per input snapshot.
in_timesteps (int) – Input snapshots per sample.
out_channels (int | None) – Channels per output snapshot. Defaults to in_channels.
out_timesteps (int) – Output snapshots per sample.
patch_size (int) – Side of the square patches.
embed_dim (int) – Latent width, divisible by 8 and by n_blocks.
depth (int) – Number of Fourier neural operator blocks.
n_blocks (int) – Channel blocks of the spectral filter.
modes (int) – Fourier modes kept.
mlp_ratio (float) – Block MLP width relative to embed_dim.
out_layer_dim (int) – Width of the output decoder.
normalize (bool) – Normalise inputs per channel and map outputs back.
act (str) – Activation name.
time_agg (str) – "mlp" or "exp_mlp".
output (str) – Output unit.
members (int) – Ensemble members.
seed (int) – Seed of the weight initialisation.

Returns:

A DPOTSurrogate.

Return type:

DPOTSurrogate

Example

>>> op = nn.DPOT("flow", grid=(64, 64), in_channels=4, in_timesteps=10,
...              patch_size=8, embed_dim=64, modes=16)
>>> y = op(op.tokens(history))          # history: [..., 64, 64, 10, 4]

class prophys.nn.DPOTSurrogate(name, net, output='', *, members=1, seed=0)[source]#

A field surrogate running a DPOT network.

The surrogate is called on tokens [..., X*Y*T, C] (from tokens() or dpot_tokens() of fields [..., X, Y, T, C]) and returns [..., X*Y, out_timesteps * out_channels], which fields() reshapes to [..., X, Y, out_timesteps, out_channels]. Input and output standardisation is off. Use normalize on the network instead. Build it with DPOT() or from_dpot().

Parameters:

name (str) – Surrogate name, the prefix of every weight.
net (DPOTNet) – The network configuration.
output (str) – Output unit.
members (int) – Ensemble members.
seed (int) – Seed of the weight initialisation.

tokens(field)[source]#

Tokens [..., X*Y*T, C] of fields [..., X, Y, T, C], the surrogate’s input layout.

Parameters:: field (Any)
Return type:: Any

fields(y)[source]#

Reshape outputs [..., X*Y, out_timesteps * out_channels] to snapshots.

Returns:: Array of shape [..., X, Y, out_timesteps, out_channels].
Parameters:: y (Any)
Return type:: Any

class prophys.nn.DPOTNet(grid, in_channels, in_timesteps, out_channels, out_timesteps=1, patch_size=8, embed_dim=64, depth=4, n_blocks=4, modes=32, mlp_ratio=1.0, out_layer_dim=32, normalize=False, act='gelu', time_agg='exp_mlp')[source]#

The whole DPOT network as one surrogate layer.

Expects X * Y * in_timesteps rows of in_channels columns (the layout of dpot_tokens()) and returns X * Y rows of out_timesteps * out_channels columns. Parameter names follow the PyTorch module, for example blocks/3/filter/w1 or patch_embed/0/W, so checkpoints convert by name. The classification head of the reference is not built. Usually created through DPOT() or from_dpot() rather than directly.

Parameters:

grid (tuple[int, int])
in_channels (int)
in_timesteps (int)
out_channels (int)
out_timesteps (int)
patch_size (int)
embed_dim (int)
depth (int)
n_blocks (int)
modes (int)
mlp_ratio (float)
out_layer_dim (int)
normalize (bool)
act (str)
time_agg (str)

grid: tuple[int, int]#: Spatial grid (X, Y), both multiples of patch_size.

in_channels: int#: Channels of each input snapshot.

in_timesteps: int#: Input snapshots per sample.

out_channels: int#: Channels of each output snapshot.

out_timesteps: int = 1#: Output snapshots per sample.

patch_size: int = 8#: Side of the square patches the grid is split into.

embed_dim: int = 64#: Latent width, divisible by 8 and by n_blocks.

depth: int = 4#: Number of Fourier neural operator blocks.

n_blocks: int = 4#: Channel blocks of the block diagonal spectral filter.

modes: int = 32#: Fourier modes kept by the spectral filter.

mlp_ratio: float = 1.0#: Hidden width of each block’s MLP relative to embed_dim.

out_layer_dim: int = 32#: Width of the output decoder.

normalize: bool = False#: Whether inputs are normalised per channel and outputs mapped back.

act: str = 'gelu'#

gelu, tanh, sigmoid, relu, softplus, silu or elu.

Type:: Activation

time_agg: str = 'exp_mlp'#: Time aggregation, "mlp" or "exp_mlp".

check(where='')[source]#

Structural findings of the configuration, as a list of diagnostics (empty when valid).

Checks that the grid is a multiple of the patch size, that embed_dim splits into the group norm’s eight groups and into n_blocks, and that the activation, time aggregation and normalisation settings are supported.

Parameters:: where (str)
Return type:: list

emit(b, prefix)[source]#

Emit the DPOT operations and weights into the network being built.

Raises:

ProphysError – When check() reports an error, or NN103 when the incoming tensor does not have the token layout.

Parameters:

b (Builder)
prefix (str)

Return type:

None

prophys.nn.from_dpot(state, *, modes, out_timesteps=1, act='gelu', name='dpot', output='')[source]#

Build a DPOT surrogate carrying the weights of a released checkpoint.

Reads the configuration with dpot_config(), builds the surrogate with DPOT() and loads the weights with dpot_state().

Parameters:

state (Mapping[str, Any]) – The checkpoint’s state dictionary, for example from torch.load.
modes (int) – Fourier modes the checkpoint was trained with.
out_timesteps (int) – Output snapshots per sample.
act (str) – Activation of the checkpoint.
name (str) – Surrogate name.
output (str) – Output unit.

Returns:

A DPOTSurrogate with the checkpoint’s weights as its initial values.

Return type:

DPOTSurrogate

Example

>>> state = torch.load("model_Ti.pth", map_location="cpu")
>>> op = nn.from_dpot(state, modes=32)

prophys.nn.dpot_config(state, *, modes, out_timesteps=1, act='gelu')[source]#

Hyperparameters of a DPOT checkpoint, read from its tensor shapes.

The Fourier modes, the output time steps and the activation are not stored in a checkpoint and are given here.

Parameters:

state (Mapping[str, Any]) – The checkpoint’s state dictionary of PyTorch tensors or arrays. A "model" entry and module. prefixes are unwrapped.
modes (int) – Fourier modes the checkpoint was trained with.
out_timesteps (int) – Output snapshots the output layer is split into.
act (str) – Activation of the checkpoint.

Returns:

Keyword arguments for DPOT().

Raises:

ProphysError – NN103 when the output maps do not split into out_timesteps.

Return type:

prophys.nn.dpot_state(surrogate, state)[source]#

Leaf values of surrogate from a DPOT state dictionary.

Convolution kernels flatten in PyTorch order, the transposed output convolution is permuted to the token layout and the classification head is skipped. Every shape is checked.

Parameters:

surrogate (DPOTSurrogate) – A DPOTSurrogate of matching configuration.
state (Mapping[str, Any]) – The checkpoint’s state dictionary.

Returns:

Mapping of full leaf name to array, ready for surrogate.set_values.

Raises:

ProphysError – NN103 when a tensor’s size does not match its leaf.

Return type:

prophys.nn.dpot_tokens(field)[source]#

Flatten fields [..., X, Y, T, C] to DPOT tokens [..., X*Y*T, C].

Accepts arrays and prophys expressions. An expression yields an expression node.

Parameters:: field (Any)
Return type:: Any

class prophys.nn.InContextRegressor(features=4, context=64, queries=32, width=64, heads=4, depth=3, levels=(0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975), head='quantile', bins=128, borders=None, tails='exponential', seed=0)[source]#

Prior fitted network that predicts the posterior predictive distribution from a context set in one pass.

fit() stores observed rows as the context. predict_quantiles() then returns predictive quantiles at new inputs with one forward pass and no training, and predict_distribution() the whole predictive distribution. context and queries are fixed sizes of the network: larger context sets are subsampled evenly, smaller ones are padded by repetition, and query rows are processed in chunks of queries. Inputs and targets are standardised with the context statistics. Needs a license (prefix attention).

Two output heads are available. head="quantile" (the shipped network) predicts the quantiles at levels and is trained with the pinball loss. head="bar" predicts a BarDistribution over bins bins at quantiles of the prior’s targets with tails, trained with cross entropy (Muller et al. 2022), which is the KL divergence to the posterior predictive up to a constant. Its quantiles come from inverting the closed form CDF.

Parameters:

features (int) – Maximum number of input columns.
context (int) – Context rows per pass.
queries (int) – Query rows per pass.
width (int) – Latent width of the transformer.
heads (int) – Attention heads, which must divide width.
depth (int) – Number of transformer blocks.
levels (Sequence[float]) – Quantile levels of the quantile head, and the levels predict_quantiles() reports for either head.
head (str) – "quantile" or "bar".
bins (int) – Bins of the bar head.
borders (Sequence[float] | None) – Bin borders of the bar head in standardised target units. Defaults to prior_borders().
tails (str) – Tails of the bar head, "exponential" or "halfnormal".
seed (int) – Seed of the weight initialisation.

Example

>>> pfn = nn.InContextRegressor.pretrained()
>>> q = pfn.fit(x_obs, y_obs).predict_quantiles(x_new)   # [n, 9]

fit(x, y)[source]#

Store the observed rows as the context and return self, without any training.

Parameters:

x (Any) – Inputs [n, f] (or [n] for one column), f <= features.
y (Any) – Targets [n].

Return type:

InContextRegressor

predict_distribution(x)[source]#

Predictive distribution at query inputs, in closed form and in the units of the targets.

A bar head gives its bar distribution. A quantile head gives the piecewise linear quantile function with exponential tails (see QuantileDistribution).

Parameters:: x (Any) – Query inputs [m, f] (or [m] for one column).
Returns:: A distribution object batched over the m rows, with cdf, log_prob, quantile, mean, variance, cvar, crps and sample methods on arrays.
Raises:: ProphysError – NN101 when fit() was not called, NN103 when there are more input columns than features.

predict_quantiles(x)[source]#

Predictive quantiles at query inputs, one column per level.

Parameters:: x (Any) – Query inputs [m, f] (or [m] for one column).
Returns:: Array [m, len(levels)] in the units of the targets.
Raises:: ProphysError – NN101 when fit() was not called, NN103 when there are more input columns than features.
Return type:: ndarray

predict(x)[source]#

Predictive medians [m] at query inputs x (the level closest to 0.5).

Parameters:: x (Any)
Return type:: ndarray

save(path)[source]#

Save the configuration and weights to a NumPy .npz file at path.

Parameters:: path (str)
Return type:: None

classmethod load(path)[source]#

Load a network written by save().

Parameters:: path (str)
Return type:: InContextRegressor

classmethod pretrained()[source]#

Load the network shipped with prophys, trained on the mechanism prior.

The shipped weights are encrypted and are decrypted in memory only while a valid license is active. The network uses prefix attention, a licensed surrogate feature.

Raises:: ProphysError – NN101 when no pretrained network is installed, LIC102 without a valid license.
Return type:: InContextRegressor

prophys.nn.train_pfn(model, steps=2000, tasks=16, learning_rate=0.001, seed=0, log_every=0, loss=None, sigma_ratio=0.75)[source]#

Train an InContextRegressor on fresh synthetic tasks from the mechanism prior.

Every step draws tasks new tasks with sample_prior_task(), each with between 8 and model.context context rows and model.queries query rows, and minimises the loss on the query rows with AdamW (cosine decay, gradient clipping). A quantile head minimises the pinball loss, a bar head the cross entropy of the bin probabilities ("histogram_ce", Muller et al. 2022), its HL Gauss smoothed variant ("hl_gauss") or the closed form CRPS ("crps"). The trained weights are stored on model.

Parameters:

model (InContextRegressor) – The network to train, in place.
steps (int) – Optimiser steps.
tasks (int) – Tasks per step.
learning_rate (float) – Peak learning rate.
seed (int) – Seed of the task generator.
log_every (int) – Print the mean loss every this many steps. Zero is silent.
loss (str | None) – Loss of a bar head, "histogram_ce" (default), "hl_gauss" or "crps". A quantile head always uses the pinball loss.
sigma_ratio (float) – HL Gauss smoothing as a multiple of the mean width of the inner bins.

Returns:

The loss of every step.

Return type:

list[float]

prophys.nn.sample_prior_task(rng, n, features)[source]#

Draw one synthetic regression task from the mechanism prior.

Inputs are uniform columns or Gaussian columns of standard deviation one or 0.7, of which a random subset is active. The latent function comes from one of seven mechanism families (dose response, hazard power law, decay with distance, piecewise linear, seasonal, smooth random field, interaction), is standardised, and is observed through Gaussian, multiplicative, Student t or Poisson count noise.

Parameters:

rng (Generator) – NumPy random generator.
n (int) – Number of rows.
features (int) – Number of input columns.

Returns:

A PriorTask with x of shape [n, features] and y of shape [n].

Return type:

PriorTask

class prophys.nn.PriorTask(x, y, family)[source]#

One synthetic regression task drawn from the mechanism prior.

Returned by sample_prior_task().

Parameters:

x (ndarray)
y (ndarray)
family (str)

x: ndarray#: Inputs of shape [n, features].

y: ndarray#: Noisy targets of shape [n].

family: str#

"logistic", "power", "decay", "piecewise", "seasonal", "smooth" or "interaction".

Type:: Mechanism family of the latent function

prophys.nn.find_surrogates(model)[source]#

Find every learned surrogate used by a probabilistic model.

Walks the expression graph of every attribute’s distribution, including upstream random variables, and collects each Surrogate, FieldSurrogate and TreeEnsemble it reaches. This is how prophys.engine.finetune() finds the surrogates to anchor with surrogate_feedback.

Parameters:: model – A ProbabilityModel.
Returns:: A dict from surrogate name to surrogate.
Return type:: dict[str, Any]

class prophys.nn.SurrogateCall(surrogate, inputs, ctx=None, mode='point')[source]#

Expression node evaluating a surrogate on input expressions.

Returned by calling a surrogate (s(...), s.raw(...), s.members_output(...)) and not usually constructed directly. It behaves like any other expression: arithmetic, use as a distribution parameter, evaluation and differentiation.

Parameters:

surrogate (Surrogate)
inputs (list[Expr])
ctx (Expr | None)
mode (str)

batch_inputs(env)[source]#

Evaluate the inputs and return flat input rows [B, in], the batch shape and the flat context.

Return type:: tuple[Any, tuple, Any]

prophys.nn.fit_linear_posterior(phi, y, alphas=None, betas=None, mean=None)[source]#

Fit a Bayesian linear regression by maximising the model evidence.

Regresses y on the features phi plus a bias column with an isotropic Gaussian prior of precision alpha and Gaussian noise of precision beta, choosing both by maximising the exact log evidence over a fixed grid (MacKay 1992). The search is exact and deterministic.

With mean given, the posterior is instead the Laplace approximation at those trained weights: the noise precision is the maximum likelihood value of their residuals and only the prior precision is chosen by the evidence.

Parameters:

phi (Any) – Features [N, H].
y (Any) – Targets [N] or [N, K].
alphas (Sequence[float] | None) – Candidate prior precisions. Defaults to 49 values spaced evenly in log over 1e-6 .. 1e6.
betas (Sequence[float] | None) – Candidate noise precisions. Defaults to 51 values spaced evenly in log over 1e-2 .. 1e8. Unused when mean is given.
mean (Any) – Trained weights [H + 1, K] (bias last) to centre the posterior on.

Returns:

The LinearPosterior.

Return type:

LinearPosterior

prophys.nn.fit_support(x, coverage=0.99, shrinkage=0.1)[source]#

Fit the support region of a set of training rows.

The box spans the central coverage quantiles of each column. The ellipsoid is centred on the mean, uses the sample covariance shrunk towards its mean variance by shrinkage, and its radius is the coverage quantile of the training rows’ Mahalanobis distances.

Parameters:

x (Any) – Training rows [N, d].
coverage (float) – Fraction of training rows inside the box per column and inside the ellipsoid.
shrinkage (float) – Weight of the isotropic target in the covariance, in [0, 1].

Returns:

The SupportRegion.

Return type:

SupportRegion

class prophys.nn.SampleSet(expr, n, seed=0)[source]#

Expression node holding a fresh Monte Carlo sample [n, k] of an expression.

The sample is drawn with fixed random numbers (common random numbers, set by seed), so it is smooth in the model parameters. Every random variable expr depends on is drawn anew, independently of the outer Monte Carlo draw of the model. Create it with sample_set().

Parameters:

expr (Any) – The expression to sample.
n (int) – Sample size (rows).
seed (int) – Seed of the fixed random numbers.

Layers#

Layer catalogue of network surrogates.

Layers are declarative descriptions passed to Surrogate or FieldSurrogate. Each layer contributes its operations to one network program that runs in the optimized native kernels, and declares its named weights with their initial values. The weights of layer k of surrogate s are ordinary Param leaves named s/Lk/<tensor> (for example power/L0/W), so they train, fine tune and export like any other parameter. Weight matrices follow the PyTorch Linear layout [out, in], so state dictionaries convert without transposes (see load_state_dict()).

A layer maps a tensor of rows by cols values per sample. For a plain Surrogate there is one row, whose columns are the inputs. A FieldSurrogate carries one row per token.

Free tier: networks of up to four weighted layers built from Dense, Residual, Fourier, RFF and LayerNorm, with any activation and any output head (QuantileHead, HistogramHead, GaussianHead, MixtureHead, Bounded, Choice, Positive). LoRA adapters, Monotone, InputConvex, Lipschitz and every attention layer need a license (see License).

class prophys.nn.layers.Activation(name='silu')[source]#

Elementwise activation function applied as its own layer.

The shape is unchanged. Use it after a head or a linear layer, for example Activation("softplus") for a positive output (see Positive()).

Parameters:: name (str) – Activation name, one of ACTIVATIONS.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Bounded(low, high)[source]#

Output head squashing into the open interval (low, high) with a sigmoid.

Bounds are in output units and the shape is unchanged. Preserves monotonicity. Place it after a linear layer of the output width.

Parameters:

low (float) – Lower bound.
high (float) – Upper bound, above low.

monotone_preserving = True#: The layer keeps elementwise monotonicity of its input.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Choice(temperature=1.0)[source]#

Softmax output head giving the choice probabilities of a discrete choice model.

Normalises the previous layer’s units to probabilities that sum to one, softmax(h / temperature). The shape is unchanged.

Parameters:: temperature (float) – Softmax temperature.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.CrossAttention(heads=4)[source]#

Multi head cross attention from the rows to a context set.

Each row attends to the context rows (for example geometry tokens or a sensor set) declared with context=(rows, cols) on the FieldSurrogate. The shape is unchanged. Needs a license.

Parameters:: heads (int) – Number of attention heads, which must divide the column count.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Dense(width, activation='silu', bias=True, lora_rank=0, lora_alpha=1.0)[source]#

Fully connected layer: an affine map followed by an activation.

Maps the cols columns of each row to width columns, act(W h + b) with W of shape [width, cols]. Weights are named <prefix>/W and <prefix>/b. A Dense layer with an "identity" activation at the end of the stack is the linear head that pretrain() and fine tuning treat specially.

With lora_rank above zero the base weight and bias are frozen and a low rank update (lora_alpha / lora_rank) B A is trained instead (Hu et al. 2022), with A of shape [lora_rank, cols] and B of shape [width, lora_rank] initialised to zero. LoRA adapters need a license. A plain Dense layer is on the free tier.

Parameters:

width (int) – Output columns.
activation (str) – Activation name, one of ACTIVATIONS.
bias (bool) – Whether to add a bias vector.
lora_rank (int) – Rank of the LoRA update. Zero trains the full weight.
lora_alpha (float) – LoRA scale numerator.

Example

>>> layers = [nn.Dense(64, "silu"), nn.Dense(64, "silu"), nn.Dense(1, "identity")]

property monotone_preserving#

an unconstrained dense layer can break monotonicity.

Type:: Always False

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Fourier(features=16, scale=1.0, trainable=False, keep_input=True)[source]#

Fourier feature embedding sin(2 pi B h) and cos(2 pi B h) (Tancik et al. 2020).

B of shape [features, cols] is drawn from Normal(0, scale) and is frozen unless trainable. The output has 2 * features columns, plus the cols input columns when keep_input is true. Helps a network represent high frequency responses. On the free tier.

Parameters:

features (int) – Number of frequencies.
scale (float) – Standard deviation of the random frequencies.
trainable (bool) – Whether the frequency matrix is trained.
keep_input (bool) – Whether to pass the input through next to the features.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.GaussianHead[source]#

Mean and scale output head for a heteroscedastic Gaussian likelihood.

Outputs two columns: the mean and an unconstrained scale that is mapped through softplus to the standard deviation.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.InputConvex(widths=(32, 32), outputs=1, activation='softplus')[source]#

Input convex neural network (Amos et al. 2017), convex in every input.

Hidden layers of widths units and a linear output layer of outputs units. Weights between hidden layers are kept non negative and every layer sees the network input again, so each output is a convex function of the input. Must see the input rows (place it before any pooling or reshape) and may only be followed by softplus or relu activations. Needs a license.

Parameters:

widths (Sequence[int]) – Widths of the hidden layers.
outputs (int) – Output columns.
activation (str) – "softplus" or "relu".

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.LayerNorm[source]#

Layer normalisation over the columns of each row, with learned gain and bias.

Weights <prefix>/g (initialised to one) and <prefix>/b (zero). The shape is unchanged. On the free tier.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Lipschitz(width, activation='tanh', cap=1.0)[source]#

Dense layer with a bounded Lipschitz constant in the maximum norm.

The absolute sum of each weight row is kept at most cap, so the layer changes its output by at most cap times the largest input change. A network of Lipschitz layers and activations of slope at most one has a guaranteed bound, reported by lipschitz_bound(). Needs a license.

Parameters:

width (int) – Output columns.
activation (str) – One of tanh, relu, sigmoid, softplus, leaky_relu or identity (slope at most one).
cap (float) – Bound on the absolute row sums of the weight.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.MixtureHead(components=3)[source]#

Gaussian mixture output head for multimodal predictive distributions.

Outputs 3 * components columns: the component logits, the means and the unconstrained scales, in that order.

Parameters:: components (int) – Number of mixture components.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Monotone(width, activation='relu', output=False)[source]#

Constrained monotone layer (Runje and Shankaranarayana 2023).

Weights reaching a constrained input are made non negative through an absolute value, and the units split into convex, concave and saturated thirds of an unbounded activation, so the network stays a universal monotone approximator and can extrapolate trends. Directions come from the surrogate’s monotone mapping for the first monotone layer, later monotone layers constrain every unit. A surrogate with a monotone mapping must be built from Monotone, Bounded and monotone activations only. Needs a license.

Parameters:

width (int) – Output columns.
activation (str) – Activation of the hidden units, one of ACTIVATIONS.
output (bool) – True for a linear output layer (identity activation).

Example

>>> s = nn.Surrogate("hazard", ["age", "load"], "1/yr",
...                  layers=[nn.Monotone(32), nn.Monotone(32)],
...                  monotone={"age": 1, "load": 1})

monotone_preserving = True#: The layer keeps elementwise monotonicity of its input.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Physics(fn, name='physics', native=None)[source]#

A layer whose body is a prophys symbolic expression.

fn receives one symbolic column per incoming unit and returns an expression or a list of expressions, one per outgoing column. It may use any prophys operation and its own Param leaves, which train together with the network weights. Expressions built from elementary operations (arithmetic, powers, exp, log, sqrt, sin, cos, tan, tanh, sigmoid, softplus, erf, erfc, abs, sign, relu, minimum, maximum, clip and where on a comparison) with scalar constants and scalar parameters compile into native code, so the surrogate stays one native program that exports, runs on the GPU and has exact derivatives. Other expressions run as a prophys stage between native segments.

Parameters:

fn (Callable[[...], Any]) – Function of the incoming columns returning one expression or a list of expressions.
name (str) – Label of the layer.
native (bool | None) – True forces native compilation (and raises NN112 when the expression cannot compile), False forbids it, None decides automatically.

Example

>>> power = nn.Physics(lambda v, rho: 0.5 * rho * v**3)

outputs(n_in)[source]#

Placeholder inputs for n_in columns and the output expressions fn builds from them.

Parameters:: n_in (int)

compile_reason(n_in)[source]#

Why the expression cannot compile natively, none when it can.

Parameters:: n_in (int)
Return type:: str | None

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.PoolAttention(heads=4)[source]#

Attention pooling of the rows with a learned seed query (Lee et al. 2019).

Reduces rows by cols to one row of cols. Needs a license.

Parameters:: heads (int) – Number of attention heads.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.PoolMean[source]#

Mean over the rows, a permutation invariant summary of a set.

Reduces rows by cols to one row of cols. Has no weights.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

prophys.nn.layers.Positive()[source]#

Softplus output activation for a positive quantity.

Returns Activation("softplus"). Append it after the linear output layer.

Example

>>> layers = [nn.Dense(32), nn.Dense(1, "identity"), nn.Positive()]

Return type:: Activation

class prophys.nn.layers.PrefixAttention(heads=4, keys=1)[source]#

Self attention in which every row attends only to the first keys rows.

This is the attention mask of prior fitted networks, with the context rows first and the query rows after them (see InContextRegressor). The shape is unchanged and the column count must be divisible by heads. Needs a license.

Parameters:

heads (int) – Number of attention heads.
keys (int) – Number of leading rows every row may attend to.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.QuantileHead(levels=(0.05, 0.25, 0.5, 0.75, 0.95), tails='exponential')[source]#

Noncrossing quantile output head.

A dense layer predicts the lowest quantile and positive increments, so the len(levels) output columns are always ordered. Train it with pretrain() (pinball or CRPS loss) and use distribution() to get the predictive distribution, a QuantileDistribution.

Parameters:

levels (Sequence[float]) – Quantile levels, strictly increasing inside (0, 1).
tails (str) – Tails of the predictive distribution beyond the outer levels: "exponential" (default), "halfnormal" or "linear" (the bounded linear extension used before prophys 0.3). See QuantileDistribution.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.RFF(features=128, lengthscale=1.0)[source]#

Random Fourier features of an RBF kernel (Rahimi and Recht 2007).

Maps the input to features random cosine features with frozen random frequencies and phases and a learned log length scale (<prefix>/log_lengthscale). A linear head on top is a sparse spectrum Gaussian process, whose Laplace posterior over the head is exact. On the free tier.

Parameters:

features (int) – Number of random features.
lengthscale (float) – Initial kernel length scale.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Reshape(rows, cols)[source]#

Reinterpret each sample’s tensor as rows by cols values.

Has no weights. The element count must stay the same.

Parameters:

rows (int) – New number of rows.
cols (int) – New number of columns.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Residual(activation='silu', depth=2, width=None)[source]#

Residual block h + f(h) with depth dense layers inside.

The inner layers use activation except the last, which is linear, maps back to the input width and starts with weights scaled by 0.1, so the block begins close to the identity. The output has as many columns as the input. On the free tier, where each inner dense layer counts towards the limit of weighted layers.

Parameters:

activation (str) – Activation of the inner layers, one of ACTIVATIONS.
depth (int) – Number of dense layers inside the block.
width (int | None) – Width of the hidden inner layers. Defaults to the input width.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.SelfAttention(heads=4)[source]#

Multi head self attention over the rows of a set or field.

Every row attends to every row. The result is equivariant to row permutations and the shape is unchanged. The column count must be divisible by heads. Needs a license.

Parameters:: heads (int) – Number of attention heads.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.SetBlock(heads=4, mlp_ratio=2)[source]#

Pre norm set transformer block: self attention and a GELU MLP, each with a residual connection.

The shape is unchanged. Needs a license.

Parameters:

heads (int) – Number of attention heads, which must divide the column count.
mlp_ratio (int) – Hidden width of the MLP relative to the columns.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.SliceAttention(heads=4, slices=16, head_dim=None)[source]#

Transolver physics attention (Wu et al. 2024), linear in the number of rows.

Rows are softly assigned to slices learned states, the slice tokens attend to each other and the result is spread back to the rows. The shape is unchanged. Weight names match physics_attention_mapping(). Needs a license.

Parameters:

heads (int) – Number of attention heads.
slices (int) – Number of slice tokens.
head_dim (int | None) – Width per head. Defaults to cols // heads.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.TransolverBlock(heads=4, slices=16, mlp_ratio=2)[source]#

Pre norm Transolver block: slice attention and a GELU MLP, each with a residual connection.

h + SliceAttention(LayerNorm(h)) followed by h + MLP(LayerNorm(h)) with a hidden width of mlp_ratio times the columns. The shape is unchanged. Weight names match transolver_mapping(). Needs a license.

Parameters:

heads (int) – Number of attention heads.
slices (int) – Number of slice tokens.
mlp_ratio (int) – Hidden width of the MLP relative to the columns.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.HistogramHead(bins=64, borders='quantile', tails='exponential')[source]#

Histogram (bar distribution) output head for a full predictive distribution.

A dense layer of width bins followed by a softmax gives the probability of each bin of the target. The predictive distribution is a BarDistribution: uniform within each bin, with exponential (or half normal) tails replacing the two outer bins, so the support is unbounded and the CDF, quantiles, CVaR, CRPS and log density are closed form. Tail risk measures therefore come from learned tail probabilities rather than an extrapolation.

Train it with pretrain(). The natural losses are "hl_gauss" (default, Farebrother et al. 2024), "histogram_ce", "kl" for distribution targets and "crps". The head has one output. The program is the same dense layer and softmax as a Choice head, so it runs everywhere a surrogate runs, on the free tier too.

Parameters:

bins (int) – Number of bins K, at least one.
borders (Any) – How the K + 1 bin borders are placed: "quantile" (default) at equally spaced quantiles of the training targets, "uniform" equally spaced between the smallest and largest training target, or an explicit increasing array of K + 1 borders in the units of the target (needed for target_kind="probs" or "cdf").
tails (str) – "exponential" (default), "halfnormal" or "none" (bounded support). See BarDistribution.

Example

>>> s = nn.Surrogate("load", ["v"], layers=[nn.Dense(64, "silu"), nn.Dense(64, "silu"), nn.HistogramHead(64)])
>>> nn.pretrain(s, v, load, loss="hl_gauss")
>>> s.distribution(v=wind).cvar(0.99)

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Called by the surrogate while it assembles its network program. User code does not call it.

Parameters:

b (Builder)
prefix (str)

Return type:

None

class prophys.nn.layers.Layer[source]#

Base class of every surrogate layer.

A layer is a small declarative object. Subclasses implement emit, which appends the layer’s operations and allocates its weights on a Builder. prefix (such as "L2") is the name under which its weights are allocated.

monotone_preserving = False#: The layer keeps elementwise monotonicity of its input.

convex_preserving = False#: The layer keeps convexity of its input.

emit(b, prefix)[source]#

Emit this layer’s operations and weights into the network being built.

Parameters:

b (Builder) – The builder of the network program.
prefix (str) – Name prefix of this layer’s weights, such as "L2".

Return type:

None

class prophys.nn.layers.Builder(rows, cols, ctx_rows=0, ctx_cols=0, seed=0)[source]#

Accumulates the operations and weight tensors of one network program.

A surrogate creates one builder per ensemble member and passes it to the Layer.emit of each layer in turn. The builder tracks the current tensor shape (rows, cols), the weight segments allocated so far and flags that record whether the network still preserves monotonicity or convexity. Only needed when writing a new layer type.

Parameters:

rows (int) – Rows of the input tensor per sample (one for a point surrogate).
cols (int) – Columns of the input tensor.
ctx_rows (int) – Rows of the optional context set, zero when there is none.
ctx_cols (int) – Columns of the context set.
seed (int) – Seed of the weight initialisation.

where#: Location of the layer being emitted, for diagnostics.

record_layer: list[str]#: Layer label of every record, to map planner errors to layers.

alloc(name, shape, init, *, trainable=True, role='weight')[source]#

Allocate a named weight tensor and return its offset in the parameter vector.

Parameters:

name (str) – Tensor name, unique within the network, such as L0/W.
shape (Sequence[int]) – Tensor shape.
init (Any) – Initial value, broadcast to shape.
trainable (bool) – False for a frozen constant.
role (str) – "weight", "bias", "const", "table" or "phys".

Raises:

ProphysError – NN101 when name is already allocated.

Return type:

int

op(code, fields)[source]#

Append one operation, named by its key in the operation table, with its numeric fields.

Parameters:

code (str)
fields (Sequence[float])

Return type:

None

fail(code, message, hint='')[source]#

Raise a diagnostic error located at the layer being emitted.

Parameters:

code (str)
message (str)
hint (str)

Return type:

None

slot()[source]#

Reserve one of the eight residual slots and return its index.

Raises:: ProphysError – NN109 when all eight slots are in use.
Return type:: int

program()[source]#

The assembled network program as a flat float64 array (header followed by the operations).

Return type:: ndarray

glorot(fan_out, fan_in, act)[source]#

Random initial weight matrix [fan_out, fan_in] suited to activation act.

He normal for ReLU like activations (relu, leaky_relu, elu, silu, gelu), the SIREN style uniform range sqrt(6 / fan_in) for sin, and Glorot uniform otherwise.

Parameters:

fan_out (int)
fan_in (int)
act (str)

Return type:

ndarray

prophys.nn.layers.ACTIVATIONS#: Activation names accepted by every layer’s activation argument, mapped to their code in the network program. "linear" is an alias of "identity" and "swish" of "silu".

Explainability#

class prophys.explain.ShapValues(names, values, base, method, interactions=None)[source]#

Feature attributions of individual predictions.

Returned by every local explanation method. For each of N explained predictions with K outputs, base + values.sum(axis=-1) reproduces the prediction (exactly for TreeSHAP, exact Shapley values and exact KernelSHAP, and up to sampling or quadrature error for the other methods). Arrays stay NumPy arrays for slicing, and to_dict() converts them.

Parameters:

names (list)
values (Any)
base (Any)
method (str)
interactions (Any)

names: list#: Feature names, in the order of the last axis of values.

values: Any#: Attributions [N, K, d] of N predictions with K outputs to d features.

base: Any#

the mean prediction over the background, the prediction at the baseline (integrated gradients), or the expected value under the node covers (TreeSHAP).

Type:: Reference value [K]

method: str#: Name of the method that produced the attributions.

interactions: Any = None#: Pairwise Shapley interaction indices [N, K, d, d] when computed (shapley_values()), else None.

to_dict()[source]#

The attributions as a dict ready for JSON, arrays as nested lists.

Return type:: dict

importance()[source]#

Global importance: mean absolute attribution per feature.

Returns:: A dict from feature name to the mean of |values| over predictions and outputs.
Return type:: dict[str, float]

prophys.explain.shapley_values(model, x, background, names=None, interactions=True)[source]#

Exact interventional Shapley values and pairwise interactions, for up to 16 features.

Every one of the 2^d coalitions is evaluated: features in the coalition are taken from the explained row, the others from each background row, and the predictions are averaged over the background. This costs 2^d * N * B model evaluations.

Parameters:

model (Any) – A surrogate, tree ensemble or callable (see as_function()).
x (Any) – Rows [N, d] to explain.
background (Any) – Background rows [B, d] defining the reference distribution.
names (Sequence[str] | None) – Feature names. Defaults to the surrogate’s input names or x0, x1, ....
interactions (bool) – Also compute the pairwise Shapley interaction indices.

Returns:

ShapValues with values [N, K, d], base the mean prediction over the background, and interactions [N, K, d, d] when requested.

Raises:

prophys.diagnostics.ProphysError – When d > 16 (NN101). Use kernel_shap() instead.

Return type:

prophys.explain.kernel_shap(model, x, background, n_samples=2048, names=None, seed=0)[source]#

Approximate interventional Shapley values by KernelSHAP.

Coalitions are drawn in complementary pairs, with sizes from the Shapley kernel distribution and uniform regression weights (Covert and Lee 2021). Coalition values are interventional averages over the background, as in shapley_values(). When n_samples >= 2^d - 2 every coalition is enumerated and the result is exact.

Parameters:

model (Any) – A surrogate, tree ensemble or callable (see as_function()).
x (Any) – Rows [N, d] to explain.
background (Any) – Background rows [B, d].
n_samples (int) – Number of sampled coalitions.
names (Sequence[str] | None) – Feature names. Defaults to the surrogate’s input names or x0, x1, ....
seed (int) – Seed of the coalition sampling.

Returns:

ShapValues with values [N, K, d] and base the mean prediction over the background.

Return type:

prophys.explain.tree_shap(ensemble, x, background=None)[source]#

Exact path dependent TreeSHAP attributions of a tree ensemble.

Explains the hard routed ensemble at its current leaf values (first member). The node covers that weight the paths come from an imported booster (from_xgboost(), from_lightgbm()) or, when background is given, are counted on those rows. With a logistic link the attributions are of the margin, before the link.

Parameters:

ensemble – A TreeEnsemble.
x (Any) – Raw inputs [N, d] to explain.
background (Any) – Rows [B, d] to count node covers on.

Returns:

ShapValues with values [N, K, d] and base the expected margin under the covers.

Raises:

prophys.diagnostics.ProphysError – When no covers are available (NN302).

Return type:

Example

shap = explain.tree_shap(trees, x_test, background=x_train)
print(shap.importance())

prophys.explain.integrated_gradients(model, x, baseline=None, steps=64, names=None)[source]#

Integrated gradients attributions of a differentiable model.

Integrates the gradient along the straight path from baseline to each row with the trapezoid rule (Sundararajan et al. 2017). The attributions sum to the difference between the prediction and the baseline prediction, up to the quadrature error.

Parameters:

model (Any) – A differentiable surrogate or callable (see as_function()).
x (Any) – Rows [N, d] to explain.
baseline (Any) – Baseline row [d] or [1, d]. Defaults to zero.
steps (int) – Number of trapezoid intervals.
names (Sequence[str] | None) – Feature names. Defaults to the surrogate’s input names or x0, x1, ....

Returns:

ShapValues with values [N, K, d] and base the prediction at the baseline.

Return type:

prophys.explain.expected_gradients(model, x, background, n_samples=256, names=None, seed=0)[source]#

Expected gradients attributions of a differentiable model.

Integrated gradients averaged over n_samples baselines drawn from the background, each with one uniformly drawn point on the path (Erion et al. 2021). They approximate interventional Shapley values for differentiable models at a cost independent of the number of features.

Parameters:

model (Any) – A differentiable surrogate or callable (see as_function()).
x (Any) – Rows [N, d] to explain.
background (Any) – Background rows [B, d].
n_samples (int) – Number of (baseline, path point) draws, shared by all rows.
names (Sequence[str] | None) – Feature names. Defaults to the surrogate’s input names or x0, x1, ....
seed (int) – Seed of the draws.

Returns:

ShapValues with values [N, K, d] and base the mean prediction over the background.

Return type: