test_stochasticIntegrators

Tests the Euler-Maruyama stochastic integrator (svStochasticIntegratorMayurama).

This module covers the deterministic (zero-noise) limit, single-realization agreement with a Python Euler-Maruyama reference, and Monte-Carlo weak-convergence of the mean and variance for the Ornstein-Uhlenbeck process. The higher-order weak Runge-Kutta methods (W2Ito1, W2Ito2, and the Roessler families) are verified separately by exact numerical equivalence to committed reference trajectories (test_stochasticIntegratorsJulia.py and test_stochasticIntegratorsPaper.py).

The Example 1 Monte Carlo validation runs in routine CI with its statistical assertion enabled. It replays a fixed ensemble of NumPy-generated Wiener increments, independent of the C++ standard library, and has no automatic retries.

class test_stochasticIntegrators.ComplexOrnsteinUhlenbeckSystem(x0: float = 0.1, y0: float = -0.1, theta_x: float = 0.1, theta_y: float = 0.073, sigma_x1: float = 0.015, sigma_x2: float = 0.011, sigma_y1: float = 0, sigma_y2: float = 0.029)

Bases: object

A process defined by

dx = -theta*x*dt + sigma_x1*dW_1 + sigma_x2*dW_2 dy = -theta*y*dt + sigma_y1*dW_1 + sigma_y2*dW_2

f(t: float, x: ndarray[tuple[Any, ...], dtype[float64]])

Drift function for the complex OU process.

property g

Return list of diffusion functions.

g1(t: float, x: ndarray[tuple[Any, ...], dtype[float64]])

First diffusion function for the complex OU process.

g2(t: float, x: ndarray[tuple[Any, ...], dtype[float64]])

Second diffusion function for the complex OU process.

class test_stochasticIntegrators.Example1System(y1_0: float = 1, y2_0: float = 1)

Bases: object

Example 1 dynamical system in:

Tang, X., Xiao, A. Efficient weak second-order stochastic Runge-Kutta methods for Itô stochastic differential equations. Bit Numer Math 57, 241-260 (2017). https://doi.org/10.1007/s10543-016-0618-9

f(t: float, x: ndarray[tuple[Any, ...], dtype[float64]])

Drift function for Example 1 system.

property g

Return list of diffusion functions.

g1(t: float, x: ndarray[tuple[Any, ...], dtype[float64]])

First diffusion function for Example 1 system.

g2(t: float, x: ndarray[tuple[Any, ...], dtype[float64]])

Second diffusion function for Example 1 system.

class test_stochasticIntegrators.ExponentialSystem(x0: float = 1)

Bases: object

A simple deterministic system with one state: dx/dt = x*t.

f(t: float, x: ndarray[tuple[Any, ...], dtype[float64]])

Drift function for the exponential system.

class test_stochasticIntegrators.OrnsteinUhlenbeckSystem(x0: float = 0.1, mu: float = 0, theta: float = 0.1, sigma: float = 0.01)

Bases: object

A process defined by

dx = theta*(mu - x)*dt + sigma*dW

f(t: float, x: ndarray[tuple[Any, ...], dtype[float64]])

Drift function for the OU process.

property g

Return list of diffusion functions.

g1(t: float, x: ndarray[tuple[Any, ...], dtype[float64]])

Diffusion function for the OU process.

mean(t: float)

E[x(t)]

var(t: float)

Var(x(t))

test_stochasticIntegrators.estimateErrorAndEmpiricalVariance(computeTrajectory: Callable[[], ndarray[tuple[Any, ...], dtype[float64]]], G: Callable[[ndarray[tuple[Any, ...], dtype[float64]]], float], estimateGOnTrajectory: float, M1: int, M2: int)

Computes the error and empirical variance according to the equations described in Section 4 of Tang & Xiao.

Parameters:

• computeTrajectory (Callable[[], npt.NDArray[np.float64]]) – A function that,

• called (when)

• from (realizes one simulation by forward propagating the system)

• paper (In the) – y_N.

• G (Callable[[npt.NDArray[np.float64]], float]) – A differentiable function

• number. (that takes in a state vector and outputs a single)

• estimateGOnTrajectory (float) – The estimate of the application of the function

• system. (G on the random variable that is the last state of the propagated)

• paper – E[G(y(t_N))]

• M1 (int) – Number of batches.

• M2 (int) – Number of trajectories per batch.

Returns:

The integrator error (\hat{\mu} in the paper), and the empirical variance (\hat{\sigma}_\mu^2 in the paper).

Return type:

tuple[float, float]

test_stochasticIntegrators.eulerMayuramaIntegrate(f: Callable[[float, ndarray[tuple[Any, ...], dtype[float64]]], ndarray[tuple[Any, ...], dtype[float64]]], g_list: List[Callable[[float, ndarray[tuple[Any, ...], dtype[float64]]], ndarray[tuple[Any, ...], dtype[float64]]]], x0: ndarray[tuple[Any, ...], dtype[float64]], dt: float, tf: float, rng_seed: int)

Euler-Mayurama integrator for the vector SDE:

dX = f(t,X) dt + sum_k g_list[k](t,X) dW_k

Parameters:

• f – Drift function.

• g_list – List of diffusion functions.

• x0 – Initial state.

• dt – Time step.

• tf – Final time.

• rng_seed – Random seed.

Returns:

Array of state trajectories, including time as the first column.

test_stochasticIntegrators.getBasiliskSim(method: Literal['EulerMayurama'], dt: float, x0: ndarray[tuple[Any, ...], dtype[float64]], f: Callable[[float, ndarray[tuple[Any, ...], dtype[float64]]], ndarray[tuple[Any, ...], dtype[float64]]], g: List[Callable[[float, ndarray[tuple[Any, ...], dtype[float64]]], ndarray[tuple[Any, ...], dtype[float64]]]], seed: int | None)

Set up and return a Basilisk simulation for a given SDE and integrator method.

Parameters:

• method – Integration method (only “EulerMayurama”).

• dt – Time step.

• x0 – Initial state.

• f – Drift function.

• g – List of diffusion functions.

• seed – RNG seed (or None for random).

Returns:

Tuple of (scSim, stateModel, integratorObject, stateLogger).

test_stochasticIntegrators.test_deterministic(method: Literal['EulerMayurama'], plot: bool = False)

Test deterministic integration (no diffusion) for Euler-Maruyama. Compares Basilisk and Python implementations against the analytical solution for the exponential system dx/dt = x*t.

Parameters:

• method – Integration method (only “EulerMayurama”).

• plot – If True, plot the relative error.

test_stochasticIntegrators.test_deterministicLimitAllIntegrators(integratorClassName: str)

Every stochastic integrator must reduce to a correct ODE solver when there is no noise. This runs each integrator on a ZERO-NOISE (m==0) system and checks it against the analytic solution.

This check is deliberately oracle-INDEPENDENT: unlike the Julia/paper equivalence tests (which replay prescribed increments through the same recurrence the reference used), the exponential system has a known closed-form solution, so a wrong drift tableau (alpha / A0 / c0) or a broken deterministic step is caught here even if the committed reference happened to share the same mistake.

It also guards the m==0 code path: an integrator that reads a per-source noise entry (dW/dZ) outside its k < m loop reads an empty Eigen vector when there are no noise sources, which this test exercises directly.

The system is dx/dt = -x (with g = []), exact solution x(t) = x0 * exp(-t). A decaying linear drift is used (rather than the growing dx/dt = x*t) so that a mis-scaled drift weight produces a clearly bounded, easily-diagnosed error.

test_stochasticIntegrators.test_example1(method: Literal['EulerMayurama'], plot: bool = False)

Test Example 1 system from Tang & Xiao (2017) for all integrator methods. Compares Basilisk and Python implementations for a single realization.

Parameters:

• method – Integration method.

• plot – If True, plot the state trajectories.

test_stochasticIntegrators.test_example1_second_moment_equations()

Verify the moment equations underlying Example 1’s analytic reference.

Ito’s formula gives the instantaneous rates of the squared components. These must satisfy the closed moment equations whose solution, from unit initial conditions, is used by the Monte Carlo validation below.

test_stochasticIntegrators.test_ou(method: Literal['EulerMayurama'], plot: bool = False)

Test Ornstein-Uhlenbeck process integration for all integrator methods. Compares Basilisk and Python implementations for a single realization.

Parameters:

• method – Integration method.

• plot – If True, plot the state trajectories.

test_stochasticIntegrators.test_ouComplex(method: Literal['EulerMayurama'], plot: bool = False)

Test integration of a two-dimensional coupled Ornstein-Uhlenbeck process. Compares Basilisk and Python implementations for a single realization.

Parameters:

• method – Integration method.

• plot – If True, plot the state trajectories.

test_stochasticIntegrators.test_validateExample1(method: Literal['EulerMayurama'])

Check Example 1’s empirical second moment with known discretization bias.

Run 1,000 trajectories in ten batches over the five-second interval used by the manual validation. Prescribed Wiener increments from NumPy’s frozen RandomState algorithm make the ensemble repeatable across C++ standard libraries, up to floating-point roundoff. The assertion is always enabled and failures are not retried.

Correct the error against the continuous analytic solution by the exact Euler-Maruyama discretization bias. The remaining sampling error must be within two estimated standard errors of the grand mean. The estimator returns the variance between batch means, so divide it by the number of batches before taking its square root. This bound is a regression criterion for this fixed ensemble, not a confidence guarantee for newly drawn samples.

Parameters:

method – Integration method.

test_stochasticIntegrators.test_validateOu(method: Literal['EulerMayurama'], figureOfMerit: Literal['mean', 'variance'], tf: float = 0.1)

Validate the weak accuracy of the integrators for the Ornstein-Uhlenbeck process. Compares the empirical mean or variance of the final state to the analytical value using multiple Monte Carlo batches.

Parameters:

• method – Integration method.

• figureOfMerit – “mean” or “variance” to validate.