@@ -6,24 +6,16 @@ weave_options:
66---
77
88```julia
9- using JumpProcesses, Graphs, Plots, Statistics, BenchmarkTools
9+ using JumpProcesses, Graphs, OrdinaryDiffEq, Plots, Statistics, BenchmarkTools
10+ fmt = :png
1011```
1112
1213# Model and example solutions
1314
1415```julia
1516function reset_history!(h; start_time = nothing)
16- if start_time === nothing
17- start_time = -Inf
18- end
1917 @inbounds for i = 1:length(h)
20- hi = h[i]
21- ix = 0
22- while ((ix + 1) <= length(hi)) && hi[ix+1] <= start_time
23- ix += 1
24- end
25- end
26- h[i] = ix == 0 ? eltype(h)[] : hi[1:ix]
18+ h[i] = eltype(h)[]
2719 end
2820 nothing
2921end
@@ -109,15 +101,15 @@ for (i, method) in enumerate(methods)
109101 else
110102 sol = solve(jump_prob, Tsit5())
111103 end
112- push!(plots, plot(sol.t, sol[1:V, :]', title=shortlabels[i], legend=false, format=fmt))
104+ push!(plots, plot(sol.t, sol[1:V, :]', title=shortlabels[i], legend=false, format=fmt))
113105end
114106plot(plots..., layout=(1, 2), format=fmt)
115107```
116108
117109# Benchmarking performance of the methods
118110
119111```julia
120- Vs = append!([1], 5:5:100 )
112+ Vs = append!([1], 5:5:95 )
121113Gs = [erdos_renyi(V, 0.2) for V in Vs]
122114nsims = 50
123115benchmarks = Vector{Vector{BenchmarkTools.Trial}}()
@@ -128,22 +120,21 @@ _u = nothing
128120_h = nothing
129121_jump_prob = nothing
130122_stepper = nothing
131- _trial = nothing
132123for method in methods
133124 println("Method: $(method)")
134125 push!(benchmarks, Vector{BenchmarkTools.Trial}())
135126 _bs = benchmarks[end]
136127 for (i, V) in enumerate(Vs)
137128 trial = try
138- _G = Gs[i]
139- _g = [neighbors(_G, i) for i in 1:nv(_G)]
140- _u = [0. for i in 1:nv(_G)]
141- _h = [eltype(tspan)[] for _ in 1:nv(_G)]
142- _jump_prob = hawkes_problem(p, method; u=u , tspan=tspan, g=g , h=h )
129+ global _G = Gs[i]
130+ global _g = [neighbors(_G, i) for i in 1:nv(_G)]
131+ global _u = [0. for i in 1:nv(_G)]
132+ global _h = [eltype(tspan)[] for _ in 1:nv(_G)]
133+ global _jump_prob = hawkes_problem(p, method; u=_u , tspan=tspan, g=_g , h=_h )
143134 if typeof(method) <: QueueMethod
144- _stepper = SSAStepper()
135+ global _stepper = SSAStepper()
145136 else
146- _stepper = Tsit5()
137+ global _stepper = Tsit5()
147138 end
148139 @benchmark(
149140 solve(_jump_prob, _stepper),
@@ -155,25 +146,22 @@ for method in methods
155146 catch
156147 BenchmarkTools.Trial(BenchmarkTools.Parameters(samples=nsims, evals=1, seconds=60))
157148 end
158- push!(_bs, _trial )
149+ push!(_bs, trial )
159150 if (V == 1 || V % 20 == 0) && (length(_bs[end]) > 0)
160151 println("\tV = $V, b = $(_bs[end])")
161152 end
162153 end
163- time = Vector{Vector{Float64}}(undef, nsims)
164- run_benchmark!(time, jump_prob, stepper)
165- push!(benchmarks, time)
166154end
167155```
168156
169157```julia
170158plot(yscale=:log10, xlabel="V", ylabel="Time (ns)", legend_position=:outertopright)
171159for (i, method) in enumerate(methods)
172160 _bs, _Vs = [], []
173- for (i , b) in enumerate(benchmarks)
161+ for (j , b) in enumerate(benchmarks[i] )
174162 if length(b) > 0
175163 push!(_bs, median(b.times))
176- push!(_Vs, Vs[i ])
164+ push!(_Vs, Vs[j ])
177165 end
178166 end
179167 plot!(_Vs, _bs, label="$method")
0 commit comments