GitHub

See also the examples directory for more involved examples.

Basic example with CHV method

A simple example of jump process is shown. We look at the following process of switching dynamics where

\[X(t) = (x_c(t), x_d(t)) \in\mathbb R\times\lbrace-1,1\rbrace.\]

In between jumps, $x_c$ evolves according to

\[\dot x_c(t) = 10x_c(t),\quad\text{ if } x_d(t)\text{ is even}.\]

\[\dot x_c(t) = -3x_c(t)^2,\quad \text{ otherwise}.\]

We first need to load the library. 

using Sundials, Random
using PiecewiseDeterministicMarkovProcesses
const PDMP = PiecewiseDeterministicMarkovProcesses
PiecewiseDeterministicMarkovProcesses

We then define a function that encodes the dynamics in between jumps. We need to provide the vector field of the ODE. Hence, we define a function that, given the continuous state $x_c$ and the discrete state $x_d$ at time $t$, returns the vector field. In addition some parameters can be passed with the variable parms.

function F!(ẋ, xc, xd, parms, t)
	if mod(xd[1],2)==0
		ẋ[1] = 10xc[1]
	else
		ẋ[1] = -3xc[1]^2
	end
end
F! (generic function with 1 method)

Let's consider a stochastic process with following transitions:

TransitionRateReaction numberJump
$x_d\to x_d+[1,0]$$k(x_c)$1[1]
$x_d\to x_d+[0,1]$$parms$2[1]

We implement these jumps using a 2x1 matrix nu of Integers, such that the jumps on each discrete component xd are given by nu * xd. Hence, we have nu = [1 0;0 -1].

The rates of these reactions are encoded in the following function.

k(x) = 1 + x

function R!(rate, xc, xd, parms, t, issum::Bool)
	# rate function
	if issum == false
	# in the case, one is required to mutate the vector `rate`
		rate[1] = k(xc[1])
		rate[2] = parms[1]
		return 0.
	else
	# in this case, one is required to return the sum of the rates
		return k(xc[1]) + parms[1]
	end
end

# initial conditions for the continuous/discrete variables
xc0 = [1.0]
xd0 = [0, 0]

# matrix of jumps for the discrete variables, analogous to chemical reactions
nu = [1 0 ; 0 -1]

# parameters
parms = [50.]
1-element Vector{Float64}:
 50.0

We define a problem type by giving the characteristics of the process F, R, Delta, nu, the initial conditions, and the timespan to solve over:

Random.seed!(8) # to get the same result as this simulation!
problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (0.0, 10.0))
PDMPProblem{Float64, Int64, Vector{Float64}, Vector{Int64}, PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(Main.F!), VariableRate{typeof(Main.R!)}, PiecewiseDeterministicMarkovProcesses.Jump{Int64, Matrix{Int64}, typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)}, Vector{Float64}, Vector{Int64}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}}, Random.TaskLocalRNG}([0.0, 10.0], PiecewiseDeterministicMarkovProcesses.PDMPJumpTime{Float64, Int64}(0.0, 0.0, 0, 0.0, [0.0, 0.0], false, 0), [0.0], [1.0;;], [0; 0;;], Float64[], PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(Main.F!), VariableRate{typeof(Main.R!)}, PiecewiseDeterministicMarkovProcesses.Jump{Int64, Matrix{Int64}, typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)}, Vector{Float64}, Vector{Int64}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}}(Main.F!, VariableRate{typeof(Main.R!)}(Main.R!), PiecewiseDeterministicMarkovProcesses.Jump{Int64, Matrix{Int64}, typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)}([1 0; 0 -1], PiecewiseDeterministicMarkovProcesses.Delta_dummy), [1.0], [0, 0], [1.0], [0, 0], PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}([0.0, 0.0], [0.0, 0.0], Any[], true), [50.0]), Random.TaskLocalRNG())

After defining the problem, you solve it using solve.

sol =  PDMP.solve(problem, CHV(CVODE_BDF()))
PDMPResult{Float64, RecursiveArrayTools.VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, RecursiveArrayTools.VectorOfArray{Int64, 2, Vector{Vector{Int64}}}}([0.0, 0.016584640908818088, 0.021832926104404297, 0.045989251418264696, 0.05839841766949395, 0.07404882595197447, 0.07636215473113918, 0.07647559650345033, 0.10249771152070208, 0.11475192733986933  …  9.842673638516434, 9.845242498578735, 9.858079614209279, 9.866281266880032, 9.88435958387877, 9.946876145348343, 9.947493781307273, 9.964080347237013, 9.976804951803679, 10.0], [1.0 1.1803928764763592 … 0.4751111085748859 0.4600304173185194], [0 1 … 27 27; 0 0 … -510 -510], [6.9076201220579e-310], (false, true), 539, 0)

In this case, we chose to sample pb with the CHV algorithm where the flow in between jumps is integrated with the solver CVODE_BDF() from DifferentialEquations.jl.

We can then plot the solution as follows:

# plotting
using Plots
Plots.plot(sol.time,sol.xc[1,:],label="xc")
Example block output

Basic example with the rejection method

The previous method is useful when the total rate function varies a lot. In the case where the total rate is mostly constant in between jumps, the rejection method is more appropriate.

The rejection method assumes some a priori knowledge of the process one wants to simulate. In particular, the user must be able to provide a bound on the total rate. More precisely, the user must provide a constant bound in between the jumps. To use this method, R_tcp! must return sum(rate), bound_rejection. Note that this means that in between jumps, one has:

sum(rate)(t) <= bound_rejection

function R2!(rate, xc, xd, parms, t, issum::Bool)
	# rate function
	bound_rejection = 1. + parms[1] + 15  # bound on the total rate
	if issum == false
	# in the case, one is required to mutate the vector `rate`
		rate[1] = k(xc[1])
		rate[2] = parms[1]
		return 0., bound_rejection
	else
	# in this case, one is required to return the sum of the rates
		return k(xc[1]) + parms[1], bound_rejection
	end
end
R2! (generic function with 1 method)

We can now simulate this process as follows

Random.seed!(8) # to get the same result as this simulation!
problem = PDMP.PDMPProblem(F!, R2!, nu, xc0, xd0, parms, (0.0, 1.0))
sol =  PDMP.solve(problem, Rejection(CVODE_BDF()))
PDMPResult{Float64, RecursiveArrayTools.VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, RecursiveArrayTools.VectorOfArray{Int64, 2, Vector{Vector{Int64}}}}([0.0, 0.017237222770056303, 0.03631074954381572, 0.039849508374552645, 0.039938841973049324, 0.04518341084742461, 0.05938144897549995, 0.09845524187549184, 0.11659246172483147, 0.1508858215902958  …  0.8464246768800252, 0.8509412543037606, 0.8650341632184299, 0.9043182659618885, 0.9127339458823895, 0.9186168012457484, 0.9561986633910544, 0.9634157569710974, 0.9861912733467151, 1.0], [1.0 1.1881209662553882 … 0.36669262595664703 0.3612633432297831], [0 0 … 1 1; 0 -1 … -49 -49], [6.90762700107025e-310, 52.18812096625539, 52.43779249159507, 52.4895840399464, 52.49091548389138, 52.57119478398121, 52.810888000103105, 53.676624431308355, 54.2089043572537, 55.52159157397741  …  51.44614850519544, 51.43331677237729, 51.43078755191963, 51.42308188903923, 51.40298831475079, 51.39892961554119, 51.39614062333924, 51.37920412025487, 51.3761161445769, 51.366692625956645], (false, true), 51, 24)

In this case, we chose to sample pb with the Rejection algorithm where the flow in between jumps is integrated with the solver CVODE_BDF() from DifferentialEquations.jl.

We can then plot the solution as follows:

# plotting
using Plots
Plots.plot(sol.time,sol.xc[1,:],label="xc")
Example block output