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], RecursiveArrayTools.VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.0]]), RecursiveArrayTools.VectorOfArray{Int64, 2, Vector{Vector{Int64}}}([[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, 0.0, 0.0, 0.0, 0.0], Any[]), [50.0]), Random.TaskLocalRNG())

After defining the problem, you solve it using solve.

sol =  PDMP.solve(problem, CHV(CVODE_BDF()))
PiecewiseDeterministicMarkovProcesses.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], RecursiveArrayTools.VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.0], [1.1803928764763592], [1.1588551769471922], [1.0690727773986204], [1.0281525468228008], [0.9808052592926326], [0.9741744649029329], [0.9738516366722972], [0.905045040549314], [0.8759014941508346]  …  [0.5874152795862426], [0.5847681769519445], [0.5718889124236182], [0.5639533407669047], [0.5472158155769341], [0.4962820869071921], [0.4958261871265907], [0.4838873931864888], [0.4751111085748859], [0.4600304173185194]]), RecursiveArrayTools.VectorOfArray{Int64, 2, Vector{Vector{Int64}}}([[0, 0], [1, 0], [1, -1], [1, -2], [1, -3], [1, -4], [1, -5], [1, -6], [1, -7], [1, -8]  …  [27, -502], [27, -503], [27, -504], [27, -505], [27, -506], [27, -507], [27, -508], [27, -509], [27, -510], [27, -510]]), [6.907861429626e-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()))
PiecewiseDeterministicMarkovProcesses.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], RecursiveArrayTools.VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.0], [1.1881209662553882], [1.4377924915950695], [1.4895840399464029], [1.4909154838913776], [1.5711947839812146], [1.8108880001031022], [2.676624431308355], [3.2089043572536973], [4.5215915739774095]  …  [0.43331677237729055], [0.4307875519196267], [0.42308188903923305], [0.40298831475079305], [0.3989296155411953], [0.3961406233392419], [0.37920412025486727], [0.37611614457690024], [0.36669262595664703], [0.3612633432297831]]), RecursiveArrayTools.VectorOfArray{Int64, 2, Vector{Vector{Int64}}}([[0, 0], [0, -1], [0, -2], [0, -3], [0, -4], [0, -5], [0, -6], [0, -7], [0, -8], [1, -8]  …  [1, -41], [1, -42], [1, -43], [1, -44], [1, -45], [1, -46], [1, -47], [1, -48], [1, -49], [1, -49]]), [6.907822294336e-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