Pychastic 0.1 reference

SDE problem setting

class pychastic.sde_problem.SDEProblem(a, b, x0, tmax, exact_solution=None)

Stores a vector stochastic differential equation of the form: \(d\mathbf{x} = \mathbf{a(x)} dt + \mathbf{B(x)} d\mathbf{w}\), where a and B are vector of drifts and matrix of noise coefficients.

a vector and B matrix should be formatted such that

x = x + a(x)*dt + np.dot(B(x),dw)

is a valid code fragment (note B has as many columns as there are noise terms).

Parameters

  • a (callable) – Function describing drift term of the equation, should return np.array of length dimension.

  • b (callable) – Function describing noise term of the equation, should return np.array of size (dimension,noiseterms).

  • dimension (int) – Dimension of the space in which \(d\mathbf{x}\) sits.

  • noiseterms (int) – Dimension of the space in which \(d\mathbf{w}\) sits.

  • x0 (np.array) – Initial value of the stochastic process.

  • tmax (float) – Time at which integration should stop.

Example

>>> import numpy as np
>>> problem = pychastic.sde_problem.VectorSDEProblem(lambda x: np.array([1,1]), lambda x: np.array([[1,0.5],[0.5,1]]), 2, 2, np.array([1.5,0.5]), 1)

SDE solvers

class pychastic.sde_solver.SDESolver(scheme='euler', dt=0.01, dt_adapting_factor=10, min_dt=None, max_dt=None, error_terms=1, target_mse_density=0.01, adaptive=False)

Produces realisations of stochastic process to solve method. Controls numerical integration features via attributes.

Parameters

  • scheme ({'euler', 'milstein', 'wagner_platen'}, default: 'euler') – Type of scheme used for integration.

  • dt (float) – Step size in fixed-step integration.

solve(problem, seed=0, chunk_size=1, chunks_per_randomization=None, progress_bar=True)

Solves SDE problem given by problem. Integration parameters are controlled by attribues of SDESolver object.

Parameters

  • problem (SDEProblem) – (Vector) SDE problem to be solved.

  • seed (int, optional) – value of seed for PRNG.

  • chunk_size (int or None, optional) – Make steps in solver in chunks of chunk_size steps. If chunk_size = n then value at every nth step is returned. If None then maximal size of chunk is used only final value is returned.

  • chunks_per_randomization (int or None, optional) – Sample wiener trajectories once per chunks_per_randomization chunks. If chunks_per_randomization = n then PRNG runs at every nth chunk. If None, PRNG runs once, at beginning of simulation. Smaller values lead to less memory usage. Larger values increase speed.

  • progress_bar (True, False) – Display tqdm style progress bar during computation.

Returns

Under following keys in returned dict you’ll find:

  • time_values – (steps,) jnp.array containing timestamps coresponding to each trajectory.

  • solution_values – (steps, problem_dimension) jnp.array containing values of integrated SDE.

  • wiener_values – (steps, noise_dimension) jnp.array containing values of Wiener processes driving the SDE.

Return type

dict

Example

>>> import pychastic
>>> import jax.numpy as jnp
>>> solver = pychastic.sde_solver.SDESolver()
>>> problem = pychastic.sde_problem.SDEProblem(
... lambda x: jnp.array([1/(2*x[0]),0]),       # [1/2r,0]
... lambda x: jnp.array([
...    [jnp.cos(x[1]),jnp.sin(x[1])],           # cos phi,      sin phi
...    [-jnp.sin(x[1])/x[0],jnp.cos(x[1])/x[0]] # -sin phi / r, cos phi / r
... ]),
... x0 = jnp.array([1.0,0.0]), # r=1.0, phi=0.0
... tmax=0.02
... )
>>> solution = solver.solve(problem)
>>> (r,phi) = (solution["solution_values"][-1,:])
>>> compare = {"integrated":(r*jnp.array([jnp.cos(phi),jnp.sin(phi)])).T,"exact":solution["wiener_values"][-1,:]+problem.x0}
>>> print(compare["integrated"],compare["exact"])
solve_many(problem: pychastic.sde_problem.SDEProblem, n_trajectories=1, step_post_processing=None, seed=0, chunk_size=1, chunks_per_randomization=None, progress_bar=True)

Solves SDE problem given by problem. Integration parameters are controlled by attribues of SDESolver object.

Parameters

  • problem (SDEProblem) – (Vector) SDE problem to be solved.

  • n_trajectories (int, optional) – Number of sample paths to generate

  • step_post_processing (callable) – (Advanced) Function of with call signature f(x) returning canonical coordinates of x. Usefull when simulating process on a manifold with does not have covering map from \(\mathbb{R}^n\) such as \(SO(3)\). Post processing function has to jit with jax. To deal with branch cuts and such refer to jax.lax.cond.

  • seed (int, optional) – value of seed for PRNG.

  • chunk_size (int or None, optional) – Make steps in solver in chunks of chunk_size steps. If chunk_size = n then value at every nth step is returned. If None then maximal size of chunk is used only final value is returned.

  • chunks_per_randomization (int or None, optional) – Sample wiener trajectories once per chunks_per_randomization chunks. If chunks_per_randomization = n then PRNG runs at every nth chunk. If None, PRNG runs once, at beginning of simulation. Smaller values lead to less memory usage. Larger values increase speed.

  • progress_bar (True, False) – Display tqdm style progress bar during computation.

Returns

Under following keys in returned dict you’ll find:

  • time_values – (n_trajectories, steps) jnp.array containing timestamps coresponding to each trajectory.

  • solution_values – (n_trajectories, steps, problem_dimension) jnp.array containing values of integrated SDE.

  • wiener_values – (n_trajectories, steps, noise_dimension) jnp.array containing values of Wiener processes driving the SDE.

Return type

dict

Example

>>> import pychastic
>>> import jax.numpy as jnp
>>> solver = pychastic.sde_solver.SDESolver()
>>> problem = pychastic.sde_problem.SDEProblem(
... lambda x: jnp.array([1/(2*x[0]),0]),       # [1/2r,0]
... lambda x: jnp.array([
...    [jnp.cos(x[1]),jnp.sin(x[1])],           # cos phi,      sin phi
...    [-jnp.sin(x[1])/x[0],jnp.cos(x[1])/x[0]] # -sin phi / r, cos phi / r
... ]),
... x0 = jnp.array([1.0,0.0]), # r=1.0, phi=0.0
... tmax=0.02
... )
>>> solution = solver.solve_many(problem,n_trajectories = 100)
>>> solution_principal = solver.solve_many(problem,step_post_processing = lambda x : jnp.fmod(x,2*jnp.pi),n_trajectories = 100)
>>> (r,phi) = (solution["solution_values"][:,-1,:]).transpose()
>>> compare = {"integrated":(r*jnp.array([jnp.cos(phi),jnp.sin(phi)])).T,"exact":solution["wiener_values"][:,-1,:]+problem.x0}
>>> print(compare["integrated"][0],compare["exact"][0])

Expected values of principle integrals

pychastic.wiener_integral_moments.E(alpha) → numpy.polynomial.polynomial.Polynomial

Calculate expected value of mixed winer integral E(I_alpha). Entries in multindex alpha coresspond to:

  • 0 - time dimension

  • positive integer - wiener dimension

Parameters

alpha (array of integers) –

Returns

polynomial in time

Return type

numpy.polynomial.Polynomial

Example

Compute expected value of integral of Wiener process with respect to itself.

>>> import pychastic.wiener_integral_moments
>>> pychastic.wiener_integral_moments.E([1,1]) # returns 0
pychastic.wiener_integral_moments.E2(alpha, beta=None) → numpy.polynomial.polynomial.Polynomial

Calculate expected value of a product of two mixed winer integrals E(I_alpha*I_beta) or E(I_alpha^2) if beta=None. Entries in multindices alpha and beta coresspond to:

  • 0 - time dimension

  • positive integer - wiener dimension

Parameters

  • alpha (array of integers) – multindex of the first integral

  • beta (array of integers or None) – multindex of the second integral

Returns

polynomial in time

Return type

numpy.polynomial.Polynomial

Example

Compute variance of integral of Wiener process with respect to itself on interval [0,1].

>>> import pychastic.wiener_integral_moments
>>> mean = pychastic.wiener_integral_moments.E([1,1])(1)
>>> mean_square = pychastic.wiener_integral_moments.E2([1,1])(1)
>>> print(f'Variance is {mean_square - mean**2}') # prints 0.5

previous

Examples