-
Notifications
You must be signed in to change notification settings - Fork 7
/
chv.jl
121 lines (91 loc) · 3.96 KB
/
chv.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
### WARNING This is an old ODE solver which is not based on an iterator implementation. We keep it until LSODA has an iterator implementation
include("chvdiffeq.jl")
function solve(problem::PDMPProblem, algo::CHV{Tode}; verbose::Bool = false, ind_save_d = -1:1, ind_save_c = -1:1, n_jumps = Inf64, reltol = 1e-7, abstol = 1e-9, save_positions = (false, true), save_rate = false, finalizer = finalize_dummy) where {Tode <: Symbol}
verbose && println("#"^30)
ode = algo.ode
@assert ode in [:cvode, :lsoda, :adams, :bdf]
verbose && printstyled(color=:red,"--> Start CHV method\n")
# initialise the problem. If I call twice this solve function, it should give the same result...
init!(problem)
# we declare the characteristics for convenience
caract = problem.caract
ratecache = caract.ratecache
ti, tf = problem.tspan
n_jumps += 1 # to hold initial vector
nsteps = 1 # index for the current jump number
xc0 = caract.xc0
xd0 = caract.xd0
# Set up initial simulation time
t = ti
X_extended = similar(xc0, length(xc0) + 1)
for ii in eachindex(xc0)
X_extended[ii] = xc0[ii]
end
X_extended[end] = ti
#useful to use the same array, as it can be used in CHV(ode)
Xd = caract.xd
if ind_save_c[1] == -1
ind_save_c = 1:length(xc0)
end
if ind_save_d[1] == -1
ind_save_d = 1:length(xd0)
end
xc_hist = VectorOfArray([copy(xc0)[ind_save_c]])
xd_hist = VectorOfArray([copy(xd0)[ind_save_d]])
res_ode = zeros(2, length(X_extended))
nsteps += 1
# define the ODE flow, this leads to big memory saving
if ode == :cvode || ode == :bdf
Flow = (X0_,Xd_,Δt,r_) -> Sundials.cvode((tt,x,xdot) -> algo(xdot, x, caract, tt), X0_, [0., Δt], abstol = abstol, reltol = reltol, integrator = :BDF)
elseif ode==:adams
Flow = (X0_,Xd_,Δt,r_) -> Sundials.cvode((tt,x,xdot) -> algo(xdot, x, caract, tt), X0_, [0., Δt], abstol = abstol, reltol = reltol, integrator = :Adams)
elseif ode==:lsoda
Flow = (X0_,Xd_,Δt,r_) -> LSODA.lsoda((tt,x,xdot,data) -> algo(xdot, x, caract, tt), X0_, [0., Δt], abstol = abstol, reltol = reltol)
end
# we use the first time interval from the one generated by the constructor PDMPProblem
δt = problem.simjptimes.tstop_extended
# Main loop
while (t < tf) && (nsteps < n_jumps)
verbose && println("--> t = ", t," - δt = ", δt, ",nstep = ", nsteps)
res_ode .= Flow(X_extended, Xd, δt, ratecache.rate)
verbose && println("--> ode solve is done!")
# this holds the new state of the continuous component
@inbounds for ii in eachindex(X_extended)
X_extended[ii] = res_ode[end, ii]
end
# this is the next jump time
t = res_ode[end, end]
caract.R(ratecache.rate, X_extended, Xd, caract.parms, t, false)
# jump time:
if (t < tf) && nsteps < n_jumps
# Update event
ev = pfsample(ratecache.rate)
# we perform the jump
affect!(caract.pdmpjump, ev, X_extended, Xd, caract.parms, t)
verbose && println("--> Which reaction? => ", ev)
verbose && println("--> xd = ", Xd)
# save state, post-jump
pushTime!(problem, t)
push!(xc_hist, X_extended[ind_save_c])
push!(xd_hist, Xd[ind_save_d])
save_rate && push!(problem.rate_hist, caract.R(ratecache.rate, X_extended, Xd, caract.parms, t, true)[1])
finalizer(ratecache.rate, caract.xc, caract.xd, caract.parms, t)
δt = - log(rand())
else
if ode in [:cvode, :bdf, :adams]
res_ode_last = Sundials.cvode((tt, x, xdot)->caract.F(xdot,x,Xd,caract.parms,tt), xc_hist[end], [problem.time[end], tf], abstol = 1e-9, reltol = 1e-7)
else#if ode==:lsoda
res_ode_last = LSODA.lsoda((tt, x, xdot, data)->caract.F(xdot,x,Xd,caract.parms,tt), xc_hist[end], [problem.time[end], tf], abstol = 1e-9, reltol = 1e-7)
end
t = tf
# save state
pushTime!(problem, tf)
push!(xc_hist, res_ode_last[end,ind_save_c])
push!(xd_hist, Xd[ind_save_d])
end
nsteps += 1
end
verbose && println("--> Done")
verbose && println("--> xc = ", xd_hist[:,1:nsteps-1])
return PDMPResult(problem.time, xc_hist, xd_hist, problem.rate_hist, save_positions, length(problem.time), 0)
end