parameter_optimization.jl 4.42 KB
Newer Older
Peter Jentsch's avatar
Peter Jentsch committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
using KissABC
using CurveFit
using Loess


@views function exp_growth_rate_estimation(incident_cases)
    diff(l) = [l[i] - l[i-1] for i=2:length(l)]
    l = length(incident_cases)
    loess_interp = loess(1:l,incident_cases)
    interp_incident_cases = Loess.predict(loess_interp, collect(1.0:l))
    second_deriative =diff(diff(interp_incident_cases))
    inflection_pt = any(second_deriative.<=-1) ? findfirst(<=(-1),second_deriative) : l  #completely heuristic
    ic,exp_growth_rate = exp_fit(1:inflection_pt, incident_cases[1:inflection_pt]) #fit y = a exp(b x), b is exp growth rate
    return ic,exp_growth_rate
end


function vaccination_data()
Peter Jentsch's avatar
Peter Jentsch committed
19
20
    ymo_total_vaccinated = [0.279,0.38,0.7]
    total_vaccinated = 0.33
Peter Jentsch's avatar
Peter Jentsch committed
21
22
23
24

end

function fit_epi_parameters(p_tuple, target_growth_rate)
Peter Jentsch's avatar
Peter Jentsch committed
25
26
27
28
29
30
31
32
    p_tuple_without_vac = merge(p_tuple,
        (
            sim_length = 30,
            immunization_begin_day = 0,
            infection_introduction_day = 1,
            immunizing = false,
        )
    )
Peter Jentsch's avatar
Peter Jentsch committed
33
34
    samples = 1
    priors = Factored(Uniform(0.0001,0.005))    
Peter Jentsch's avatar
Peter Jentsch committed
35

Peter Jentsch's avatar
Peter Jentsch committed
36
    function cost(p)
Peter Jentsch's avatar
Peter Jentsch committed
37
        new_params = merge(p_tuple_without_vac, (base_transmission_probability = p[1],))
Peter Jentsch's avatar
Peter Jentsch committed
38
39
40
41
42
        out = mean_solve(samples, new_params ,DebugRecorder)
        incident_cases = mean.(out.daily_cases)
        _,exp_growth_rate = exp_growth_rate_estimation(incident_cases)
        return (target_growth_rate - exp_growth_rate)^2
    end
Peter Jentsch's avatar
Peter Jentsch committed
43

Peter Jentsch's avatar
Peter Jentsch committed
44
45
46
47
48
49
50
    out = ABCDE(priors,cost,0.01; verbose=true, nparticles=100,generations=20,  parallel = true) #this one has better NaN handling
    return out
end


function fit_behavioural_parameters(p_tuple, target_growth_rate)
    samples = 1
51
52
    p_names = (:π_base,:ρ_y,:ρ_m,:ρ_o)
    priors = Factored(Uniform(-2.0,0.0),Uniform(0.0,1.0),Uniform(0.0,1.0),Uniform(0.0,1.0))
Peter Jentsch's avatar
Peter Jentsch committed
53
54
55
    #simulation begins in august
    #30 days for opinion dynamics to stabilize, then immunization begins in september,
    #infection is introduced at the end of november 
56
    sim_length = 120
Peter Jentsch's avatar
Peter Jentsch committed
57
58
59
60
    p_tuple_adjust = merge(p_tuple,
        (
            sim_length = sim_length,
            I_0_fraction = 0.000,
61
            immunization_begin_day =60, 
62
            infection_introduction_day = 91,
Peter Jentsch's avatar
Peter Jentsch committed
63
64
65
66
            immunizing = true,
        )
    )
    target_cumulative_vac_proportion = 0.33
67
    vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
Peter Jentsch's avatar
Peter Jentsch committed
68
    ymo_vac = @SVector [0.255,0.278,0.602]
Peter Jentsch's avatar
Peter Jentsch committed
69
    function cost(p)
Peter Jentsch's avatar
Peter Jentsch committed
70
71
72
        new_params = merge(p_tuple_adjust, NamedTuple{p_names}(ntuple(i -> p[i],length(p_names))))
        out = DebugRecorder(0,sim_length)
        model = abm(new_params,out)
73

Peter Jentsch's avatar
Peter Jentsch committed
74
        target_cumulative_vaccinations = target_cumulative_vac_proportion*model.nodes
75
76
77
78
79
80
        target_ymo_vac = ymo_vac .* sum(vaccination_data[1:4]) .* target_cumulative_vaccinations
        ymo_vaccination_ts = out.new_ymo_immunization
        total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
        # display(target_ymo_vac)
        # display(total_preinfection_vaccination)
        return sum((total_preinfection_vaccination .- target_ymo_vac).^2)    
Peter Jentsch's avatar
Peter Jentsch committed
81
    end
Peter Jentsch's avatar
Peter Jentsch committed
82
#smc(priors,cost; verbose = true, nparticles = 100, parallel = true)#
83
    out = ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=100,  parallel = true) #this one has better NaN handling
Peter Jentsch's avatar
Peter Jentsch committed
84
85
86
87
88
    return out
end


function plot_behavioural_fit(particles,p_tuple)
89
    p_names = (:π_base, :ρ_y,:ρ_m,:ρ_o)
Peter Jentsch's avatar
Peter Jentsch committed
90
91
92
    sim_length = 210
    samples = 1
    p_tuple_adjust = merge(p_tuple,
93
    (
Peter Jentsch's avatar
Peter Jentsch committed
94
            sim_length = sim_length,
95
96
97
            I_0_fraction = 0.000,
            immunization_begin_day =60, 
            infection_introduction_day = 91,
Peter Jentsch's avatar
Peter Jentsch committed
98
99
100
101
102
103
104
            immunizing = true,
        )
    )
    p = map(e -> mode(e.particles),particles.P)
    display(p)
    new_params = merge(p_tuple_adjust, NamedTuple{p_names}(ntuple(i -> p[i],length(p_names))))
    out = mean_solve(samples, new_params ,DebugRecorder)
105
106
107
108
109
110
    target_cumulative_vac_proportion = 0.33
    vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
    ymo_vac = @SVector [0.255,0.278,0.602]
    ymo_vaccination_ts = mean.(out.new_ymo_immunization)
    total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
    display(total_preinfection_vaccination)
Peter Jentsch's avatar
Peter Jentsch committed
111
112
    p = plot_model(nothing,[nothing],[out],new_params.infection_introduction_day,new_params.immunization_begin_day)
    savefig(p,"behaviour_fit.pdf")
Peter Jentsch's avatar
Peter Jentsch committed
113
    return out
114
115
116
117
end

#once we have vaccine parameters, with seasonal, look at cumulative infections
#