-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathTransport.py
140 lines (86 loc) · 3.49 KB
/
Transport.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
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
#!/usr/bin/env python
# coding: utf-8
# In[1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pulp
import glob
from matplotlib import rcParams
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Hiragino Maru Gothic Pro', 'Yu Gothic', 'Meirio', 'Takao', 'IPAexGothic', 'IPAPGothic', 'VL PGothic', 'Noto Sans CJK JP']
# In[8]:
filenames = glob.glob('data_severe/x_*')
filenames.sort()
filenames
date_list = [filename.split('_')[-1].split('.')[0] for filename in filenames]
f_date = max(date_list)
print("forcasted data ={0}".format(f_date))
# In[9]:
def calc_transport(gamma,x_type,forecast_date):
# データの入力
if x_type == 'upper':
df_x = pd.read_csv('data_severe/x0975_{0}.csv'.format(forecast_date),index_col=0 )
if x_type == 'mean':
df_x = pd.read_csv('data_severe/x_{0}.csv'.format(forecast_date),index_col=0 )
if x_type == 'bottom':
df_x = pd.read_csv('data_severe/x0025_{0}.csv'.format(forecast_date),index_col=0 )
x = df_x.values
df_severe_beds = pd.read_csv('data_Koro/severe_beds.csv',index_col=0)
dirnames = (df_severe_beds['japan_prefecture_code']+df_severe_beds['都道府県名']).values
names = df_severe_beds['都道府県名'].values
weeks = df_severe_beds.columns[2:].values
new_week = max(weeks)
M = df_severe_beds[new_week].values
N = x.shape[1]
T = x.shape[0]
U = np.zeros((T,N,N))
L = np.kron(np.ones((1,N)),np.eye(N)) - np.kron(np.eye(N),np.ones((1,N)))
df_w = pd.read_csv('data_Kokudo/w_distance.csv',index_col=0)
W= df_w.values
k = 0
y = np.zeros(x.shape)
uv = np.zeros((T,N*N))
y[0] = x[0].copy()
w_pulp = W.T.reshape(-1)
status = np.zeros(T-1)
gammas = gamma * np.ones(T)
##################Start of OPTIMIZATION#######################################################
for k in range(T-1):
for i_calc in range(100):
A_pulp = L.copy()
b_pulp =np.trunc(gammas[k]*M) - (y[k] + x[k+1]-x[k])
# 数理モデル
m = pulp.LpProblem()
# 変数
uv_N = N*N
# x_pulp = [pulp.LpVariable('x%d'%i, lowBound=0) for i in range(uv_N)]
x_pulp = [pulp.LpVariable('x%d'%i, lowBound=0,cat='Integer') for i in range(uv_N)]
# 目的関数
m += pulp.lpDot(w_pulp,x_pulp)
# 拘束条件(Ax<=b)
for i in range(A_pulp.shape[0]):
m += pulp.lpDot(A_pulp[i],x_pulp) <= b_pulp[i]
#Build the solverModel for your preferred
solver = pulp.PULP_CBC_CMD(maxSeconds = 30)
status[k] = m.solve(solver)
print(k,gammas[k])
if status[k]==1:
gammas[k+1] = gammas[k]
break;
else:
gammas[k] = gammas[k] + 0.01
for i in range(len(x_pulp)):
uv[k,i] = pulp.value(x_pulp[i])
y[k+1] = y[k] + x[k+1] - x[k] + L.dot(uv[k])
################## End of OPTIMIZATION #######################################################
np.save('data_transport/u_{0}_{1:03}_{2}.npy'.format(x_type,int(gamma*100),forecast_date), uv)
np.save('data_transport/gammas_{0}_{1:03}_{2}.npy'.format(x_type,int(gamma*100),forecast_date), gammas)
# In[10]:
gamma_list = [0.8,1.0]
# 使う条件の設定
x_type_list = ['upper','mean','bottom']
for gamma in gamma_list:
for x_type in x_type_list:
calc_transport(gamma,x_type,f_date)
# In[ ]: