-
Notifications
You must be signed in to change notification settings - Fork 0
/
scenarios_utils.py
104 lines (88 loc) · 4.45 KB
/
scenarios_utils.py
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
import scipy.interpolate
import itertools
import datetime
import numpy as np
import pandas as pd
def pick_scenario(setup, scn_id):
if setup.scale == 'Regions':
scenarios_specs = {
'vaccpermonthM': [3, 15, 150], # ax.set_ylim(0.05, 0.4)
# 'vacctotalM': [2, 5, 10, 15, 20],
'newdoseperweek': [125000, 250000, 479700, 1e6, 1.5e6, 2e6],
'epicourse': ['U', 'L'] # 'U'
}
elif setup.scale == 'Provinces':
if setup.nnodes == 107:
scenarios_specs = {
'vaccpermonthM': [3, 15, 150], # ax.set_ylim(0.05, 0.4)
# 'vacctotalM': [2, 5, 10, 15, 20],
'newdoseperweek': [125000, 250000, 479700, 1e6, 1.5e6, 2e6],
'epicourse': ['U', 'L'] # 'U'
}
elif setup.nnodes == 10:
scenarios_specs = {
'newdoseperweek': [125000*10, 250000*10],
'vaccpermonthM': [14*10, 15*10],
'epicourse': ['U'] # 'U', 'L'
}
# Compute all permutatios
keys, values = zip(*scenarios_specs.items())
permuted_specs = [dict(zip(keys, v)) for v in itertools.product(*values)]
specs_df = pd.DataFrame.from_dict(permuted_specs)
if setup.nnodes == 107:
specs_df = specs_df[((specs_df['vaccpermonthM'] == 15.0) | (specs_df['newdoseperweek'] == 479700.0))].reset_index(drop=True) # Filter out useless scenarios
# scn_spec = permuted_specs[scn_id]
scn_spec = specs_df.loc[scn_id]
tot_pop = setup.pop_node.sum()
scenario = {'name': f"{scn_spec['epicourse']}-r{int(scn_spec['vaccpermonthM'])}-t{int(scn_spec['newdoseperweek'])}-id{scn_id}",
'newdoseperweek': scn_spec['newdoseperweek'],
'rate_fomula': f"({scn_spec['vaccpermonthM'] * 1e6 / tot_pop / 30}*pop_nd)"
}
# Build beta scenarios:
if scn_spec['epicourse'] == 'C':
scenario['beta_mult'] = np.ones((setup.nnodes, setup.ndays))
elif scn_spec['epicourse'] == 'U':
course = scipy.interpolate.interp1d([0, 10, 40, 80, 100, 1000], [1.4, 1.2, .9,.8, 1.2, .75], kind='linear')
course = scipy.interpolate.interp1d([0, 10, 40, 80, 100, 1000],
[1.8 * .97, 1.5 * .97, .7 * .97, 1.2 * .97, 1.2 * .97, .75 * .97],
kind='cubic')
course = course(np.arange(0, setup.ndays))
scenario['beta_mult'] = np.ones((setup.nnodes, setup.ndays)) * course
elif scn_spec['epicourse'] == 'L':
course = scipy.interpolate.interp1d([0, 10, 40, 80, 100, 1000], [.8, .4, 1.2, .7, .75, .75], kind='linear')
course = scipy.interpolate.interp1d([0, 10, 40, 80, 100, 1000], [1.8, .9, .8, .68, .75, .75], kind='cubic')
course = course(np.arange(0, setup.ndays))
scenario['beta_mult'] = np.ones((setup.nnodes, setup.ndays)) * course
return scenario
def build_scenario(setup, scenario, strategy = None):
M = setup.nnodes
N = setup.ndays - 1
control_initial = np.zeros((M, N))
maxvaccrate_regional = np.zeros((M, N))
unvac_nd = np.copy(setup.pop_node)
delivery_national = np.zeros(N)
if strategy is None:
strategy = setup.pop_node
stockpile = 0
for k in range(N):
if (setup.start_date + datetime.timedelta(days=k)).weekday() == 0: # if monday
delivery_national[k] = scenario['newdoseperweek']
else:
delivery_national[k] = 0
stockpile += delivery_national[k]
for nd in range(M):
pop_nd = setup.pop_node[nd]
maxvaccrate_regional[nd, k] = eval(scenario['rate_fomula'])
to_allocate = stockpile * strategy[nd] / strategy.sum()
to_allocate = min(to_allocate, maxvaccrate_regional[nd, k]*.9, unvac_nd[nd]*.9)
control_initial[nd, k] = to_allocate
stockpile -= to_allocate
unvac_nd[nd] -= to_allocate
stockpile_national_constraint = np.cumsum(delivery_national)
for k in range(N):
if (setup.start_date + datetime.timedelta(days=k)).weekday() != 6: # if NOt sunday
stockpile_national_constraint[k] = np.inf
if stockpile_national_constraint[-1] == np.inf:
# np.nanmax(stockpile_national_constraint[stockpile_national_constraint != np.inf]) + scenario['newdoseperweek']
stockpile_national_constraint[-1] = np.cumsum(delivery_national).max()
return maxvaccrate_regional, delivery_national, stockpile_national_constraint, control_initial