ContinuationForOptimisation/simulationTimings-cipa.jl at main · CardiacModelling/ContinuationForOptimisation · GitHub

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
using DifferentialEquations, DiffEqCallbacks
using BenchmarkTools
using Random

include("cipa.jl")
using .Cipa

prob = Cipa.prob

const debug = false
const nParameters = 100 # How many parameter vectors to use for the benchmark

# Setup
STATE::Vector{Float64} = zeros(size(ic))
function affect!(integrator)
    error = sum(abs.(STATE .- integrator.u))
    debug && println(error)
    if error < 1e-6
        terminate!(integrator)
    end
    STATE .= integrator.u
end
cb = PeriodicCallback(affect!, Cipa.pulse_period, save_positions=(false, false))

# Parameters
Random.seed!(0)
params = [0.85 .+ rand(14) * 0.3 for _ in 1:nParameters]
@show params

# ODE Convergence - Standard
println("Standard Approach")
for i in eachindex(params)
    println("Standard for parameter vector $i")
    prob_de = prob
    param_map!(prob_de, params[i])
    sol = DifferentialEquations.solve(prob_de, Tsit5(); callback=cb, Cipa.solversettings(save=false, maxt=2e6)...)
    @show sol.t
    @show sol.u[end]
end

flush(stdout)

println("Running standard approach benchmark")
Random.seed!(0)
b = @benchmarkable DifferentialEquations.solve(prob, $Tsit5(); callback=$cb, $Cipa.solversettings(save=false, maxt=2e6)...) setup = (param_map!(prob, 0.85 .+ rand(14)*0.3))
t = run(b, seconds=3600*15, samples=nParameters, evals=1)
if length(t.times) != nParameters
    @show t
    throw("Not all parameters were run. $(length(t.times)) != $nParameters")
end
BenchmarkTools.save("results/cipa/simTimings/standard.json", t)

# ODE Convergence - Tracking
println("Tracking Approach")
for i in eachindex(params)
    println("Tracking for parameter vector $i")
    prob_de = remake(prob, u0=ic_conv)
    param_map!(prob_de, params[i])
    sol = DifferentialEquations.solve(prob_de, Tsit5(); callback=cb, Cipa.solversettings(save=false, maxt=2e6)...)
    @show sol.t
    @show sol.u[end]
end

flush(stdout)

println("Running tracking approach benchmark")
Random.seed!(0)
prob_de = remake(prob, u0=Cipa.ic_conv)
b = @benchmarkable DifferentialEquations.solve(prob_de, $Tsit5(); callback=$cb, $Cipa.solversettings(save=false, maxt=2e6)...) setup = (param_map!(prob_de, 0.85 .+ rand(14)*0.3))
t = run(b, seconds=3600*15, samples=nParameters, evals=1)
if length(t.times) != nParameters
    @show t
    throw("Not all parameters were run. $(length(t.times)) != $nParameters")
end
BenchmarkTools.save("results/cipa/simTimings/tracking.json", t)

for i in eachindex(params)
    println("Continuation Approach")
    println("Continuation for parameter vector $i")
    lc = Cipa.continuation(Cipa.ic_conv, ones(14),
    params[i], debug)
    @show lc
end

flush(stdout)

println("Running continuation approach benchmark")
Random.seed!(0)
b = @benchmarkable Cipa.continuation($Cipa.ic_conv, $ones(14),
    p, false) setup = (p = 0.85 .+ rand(14)*0.3)
t = run(b, seconds=3600*15, samples=nParameters, evals=1)
if length(t.times) != nParameters
    @show t
    throw("Not all parameters were run. $(length(t.times)) != $nParameters")
end
BenchmarkTools.save("results/cipa/simTimings/continuation.json", t)