ContinuationForOptimisation/figures.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
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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
using Plots, LaTeXStrings
using BenchmarkTools, BenchmarkPlots, StatsPlots
using CSV, DataFrames

# Simulation timings
t = BenchmarkTools.load("results/simTimings/data.json")[1]
l = @layout [a b]

plot(t["Small"]["ODE"]["ODE - Standard"], label = nothing, st=:box, dpi=300, size=(397,250), linecolor=:match,
markerstrokewidth=0, title="Small Perturbation", whisker_range=0)
plot!(t["Small"]["ODE"]["ODE - Tracking"], label = nothing, st=:box, linecolor=:match,
    markerstrokewidth=0, whisker_range=0)
plot!(t["Small"]["Cont"]["Cont - Shooting"], st=:box, xaxis=nothing, linecolor=:match,
markerstrokewidth=0, ylabel="Time (s)", yformatter=x->x/1e9, ylim=(0, 0.75e9), whisker_range=0)
xaxis!(xticks=([1, 2, 3], ["Standard", "Tracking", "Continuation"]), xminorticks=false)
plotA = yaxis!(minorgrid=true)
plot(t["Large"]["ODE"]["ODE - Standard"], st=:box, dpi=300, size=(397,300), title="Large Perturbation",
linecolor=:match, markerstrokewidth=0, whisker_range=0)
plot!(t["Large"]["ODE"]["ODE - Tracking"], st=:box, linecolor=:match,
    markerstrokewidth=0, whisker_range=0)
plot!(t["Large"]["Cont"]["Cont - Shooting"], st=:box, legend=nothing, xaxis=nothing, linecolor=:match, markerstrokewidth=0,
ylabel="", yformatter=x->x/1e9, ylim=(0, 0.75e9), whisker_range=0)
xaxis!(xticks=([1, 2, 3], ["Standard", "Tracking", "Continuation"]), xminorticks=false)
plotB = yaxis!(minorgrid=true)

plot(plotA, plotB, layout=l, size=(539,250), dpi=300, margins=2Plots.mm)
annotate!(0, 0.8e9, text("A", 12, :black), subplot=1)
annotate!(0, 0.8e9, text("B", 12, :black), subplot=2)
savefig("results/simTimings/simTimings.pdf")

# MCMC
file_types = ["results/mcmc/cont_", "results/mcmc/trackingODE_", "results/mcmc/standardODE_"]
plots = []
numSamples = 40000
for file_type in file_types
    global plots
    # Read the chain data from the CSV file
    data = CSV.read(file_type*"chain.csv", DataFrame)
    accepts = data.Accept
    lls = data.ℓ
    chain = Matrix([data.gNa data.gK data.gL data.σ])

    # Plot results
    plot_params = (linewidth=2., dpi=300)
    paramNames = ["θ₁" "θ₂" "θ₃" "σ"]

    # Plot acceptance rate
    plot([mean(accepts[max(i-499,1):i]) for i in 1:numSamples], xlabel="Iteration",
    ylim=(0,1), label="Acceptance Rate", xlim = (0,numSamples), xticks=([0,20000,40000],["0","2×10⁵","4×10⁵"]); plot_params...)
    vline!([numSamples*0.25+0.5], label="Burn In", color=:red, linewidth=1.5, linestyle=:dot)
    plots = [plots... vline!([numSamples*0.1+0.5], label="Adaption", color=:green, linewidth=1.5, linestyle=:dot, legend=nothing)]

    # Plot log likelihood
    plot(lls, xlabel="Iteration", xlim=(0,numSamples),
    label="ℓ", xticks=([0,20000,40000],["0","2×10⁵","4×10⁵"]); plot_params...)
    vline!([numSamples*0.25+0.5], label="Burn In", color=:red, linewidth=1.5, linestyle=:dot)
    plots = [plots... vline!([numSamples*0.1+0.5], label="Adaption", color=:green, linewidth=1.5, linestyle=:dot, legend=nothing)]

    # Remove burn in stage to get posterior distribution
    burnIn = round(Int, numSamples*0.25)
    posterior = chain[burnIn+1:end, :]

    # Plot parameter convergence
    pTrueWithNoise = [1.0, 1.0, 1.0, 2.0]
    order = [4, 3, 1, 2]
    plot(chain[:,order]./pTrueWithNoise[order]', label="", xticks=([0,20000,40000],["0","2×10⁵","4×10⁵"]),
    xlabel="Iteration", xlim=(0,numSamples);
    plot_params...)
    # Hodge podge of lines in the right order for the legend
    for i in 1:4
        plot!([-1], [-1], label=paramNames[i], color=findfirst(i.==order); plot_params...)
    end
    vline!([numSamples*0.25+0.5], label="Burn In", color=:red, linewidth=1.5, linestyle=:dot)
    ylims!(0.7,1.2)
    plots = [plots... vline!([numSamples*0.1+0.5], label="Adaption", color=:green, linewidth=1.5, linestyle=:dot, legend=nothing)]

    # Plot posterior histograms
    p = corrplot(posterior, label=paramNames, size=(539,500), xrot=90, fillcolor=:thermal)
    for i in 1:4
        for j in 1:4
            if i != j
                scatter!(p, [pTrueWithNoise[i]], [pTrueWithNoise[j]], subplot=(j-1)*4+i, label="", color=:red, marker=:x)
            end
            if i == j
                vline!(p, [pTrueWithNoise[i]], subplot=(j-1)*4+i, label="", color=:red)
            end
        end
    end
    for (i, subplot) in enumerate([4,8,12])
        yaxis!(p, yformatter=x -> x, subplot=subplot, ymirror=true, yaxisposition=:right, ylabel=paramNames[i])
    end
    for subplot in [2, 3, 7]
        yaxis!(p, subplot=subplot, ymirror=true, yaxisposition=:right)
    end
    for subplot in [1,6,11,16]
        yaxis!(p, yformatter=x -> x, subplot=subplot)
    end
    xaxis!(p, xticks=([0.92, 0.96, 1.00, 1.04], ["0.92", "0.96", "1.00", "1.04"]), subplot=15)
    xaxis!(p, xticks=([1.8, 2.0, 2.2], ["1.8", "2.0", "2.2"]), subplot=16)
    yaxis!(p, ylabel="", subplot=1)

    for col in 3:4
        for row in 1:col-1
            if col == 3
                xaxis!(p, xticks=([0.92, 0.96, 1.00, 1.04], []), subplot=(row-1)*4+col)
            else
                xaxis!(p, xticks=([1.8, 2.0, 2.2], []), subplot=(row-1)*4+col)
            end
        end
    end

    savefig(file_type*"posterior.pdf")
end

l = @layout [a b c]

function plotter(plots, title)
    plot(plots..., layout=l, size=(539,250), dpi=300, link=:all, bottom_margin=2Plots.mm, right_margin=3Plots.mm, yformatter=:none, title=["A" "B" "C"], titlelocation=:left)
    yaxis!(yformatter=x->x, ylabel=title, subplot=1)
    return plot!(legend=:bottomright, subplot=3)
end

plot_ = plotter(plots[1:3:end], "Acceptance Rate")
plot!(legend=:topright, subplot=3)
savefig("results/mcmc/acceptanceRate.pdf")
plot_ = plotter(plots[2:3:end], "Log Likelihood")
savefig("results/mcmc/logLikelihood.pdf")
plot_ = plotter(plots[3:3:end], "Normalized Parameters")
savefig("results/mcmc/convergence.pdf")

# Cipa - Simulation timings
tStandard = BenchmarkTools.load("results/cipa/simTimings/standard.json")[1]
tTracking = BenchmarkTools.load("results/cipa/simTimings/tracking.json")[1]
tCont = BenchmarkTools.load("results/cipa/simTimings/continuation.json")[1]

plot(tStandard, st=:box, dpi=300, size=(397, 300), linecolor=:match,
    markerstrokewidth=0, title="CiPA Convergence Times", label=nothing, whisker_range=0)
plot!(tTracking, st=:box, linecolor=:match, markerstrokewidth=0, label = nothing, whisker_range=0)
plot!(tCont, st=:box, legend=:bottomleft, xaxis=nothing, linecolor=:match,
    markerstrokewidth=0, ylabel="Time (s)", yformatter=x -> x / 1e9, label = nothing, ylim=(0, 0.75e12), whisker_range=0)
yaxis!(minorgrid=true)
xaxis!(xticks=([1,2,3], ["Standard", "Tracking", "Continuation"]), xminorticks=false)
savefig("results/cipa/simTimings/simTimings.pdf")