-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathactive_learning_2Dexample.py
169 lines (128 loc) · 4.95 KB
/
active_learning_2Dexample.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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 12 14:17:28 2019
@author: fsc
"""
import numpy as np
import matplotlib.pyplot as plt
from pyDOE import lhs
from models_para_tf import Eikonal2DnetCV2RPF
import entropy_estimators as ee
np.random.seed(1234)
def plot_ensemble(T_star, CV_star, X_train, Y_train, filename = None):
plt.set_cmap('jet_r')
scale = np.linspace(0,0.75, 16)
scaleCV = np.linspace(0.9,1.5, 16)
plt.close('all')
fig = plt.figure(1)
fig.set_size_inches((10,15))
plt.subplot(321)
plt.contourf(X_m, Y_m, T_star.mean(1).reshape(X_m.shape), scale)
plt.colorbar()
plt.scatter(X_train, Y_train)
plt.xlim([0,1])
plt.ylim([0,1])
plt.title('predicted activation times [ms]')
plt.subplot(322)
plt.contourf(X_m, Y_m, (T_star.mean(1) - T[:,0]).reshape(X_m.shape))
plt.colorbar()
plt.scatter(X_train, Y_train)
plt.xlim([0,1])
plt.ylim([0,1])
plt.title('activation time error [ms]')
plt.subplot(323)
plt.contourf(X_m, Y_m, CV_star.mean(1).reshape(X_m.shape))
plt.colorbar()
plt.scatter(X_train, Y_train)
plt.xlim([0,1])
plt.ylim([0,1])
plt.title('predicted conduction velocity [mm/ms]')
plt.subplot(324)
plt.contourf(X_m, Y_m, (CV[:,0] - CV_star.mean(1)).reshape(X_m.shape))
plt.colorbar()
plt.scatter(X_train, Y_train)
plt.xlim([0,1])
plt.ylim([0,1])
plt.title('conduction velocity error [mm/ms]')
plt.set_cmap('jet')
plt.subplot(325)
plt.contourf(X_m, Y_m, T_star.std(1).reshape(X_m.shape))
plt.colorbar()
plt.scatter(X_train, Y_train)
plt.xlim([0,1])
plt.ylim([0,1])
plt.title('activation times std [ms]')
ent = np.apply_along_axis(lambda x: ee.entropy(x[:,None]), 1, T_star)
plt.subplot(326)
plt.contourf(X_m, Y_m, ent.reshape(X_m.shape))
plt.colorbar()
plt.scatter(X_train, Y_train)
plt.xlim([0,1])
plt.ylim([0,1])
plt.title('activation times entropy [-]')
plt.tight_layout()
if filename is not None:
plt.savefig(filename)
if __name__ == "__main__":
def exact(X, Y):
return np.minimum(np.sqrt(X**2 + Y**2), 0.7*np.sqrt((X - 1)**2 + (Y - 1)**2))
def CVexact(X, Y):
mask = np.less_equal(np.sqrt(X**2 + Y**2), 0.7*np.sqrt((X - 1)**2 + (Y - 1)**2))
return mask*1.0 + ~mask*1.0/0.7
# create the data
N_grid = 50
x = y = np.linspace(0,1,N_grid)[:,None]
# start training with a low number of samples
N_train = 10
X_m, Y_m = np.meshgrid(x,y)
X = X_m.flatten()[:,None]
Y = Y_m.flatten()[:,None]
T = exact(X,Y)
CV = CVexact(X,Y)
X_train_all = lhs(2, N_train)
X_train = X_train_all[:,:1]
Y_train = X_train_all[:,1:]
T_train = exact(X_train, Y_train)
# network parameters
X_pde = X
Y_pde = Y
layers = [2,20,20,20,20,20,1]
CVlayers = [2,5,5,5,5,1]
Batch = 30 # this is the number of networks to train in parallel
CVmax = 1.5
model = Eikonal2DnetCV2RPF(X_pde, Y_pde, X_train, Y_train, T_train,
layers, CVlayers, Batch, C = CVmax, alpha = 1e-7, alphaL2 = 1e-9)
#start training the model with the initial dataset
model.train_Adam_minibatch(20000 + 5000*(N_train - 10), size = 50)
T_star, CV_star = model.predict(X,Y)
plot_ensemble(T_star, CV_star, X_train, Y_train, filename = 'results/AL_NNpara_0.pdf')
N_AL = 40 # number of samples to acquire with active learning
# store how the error is evolving
errorsT = [np.sqrt(np.average((T_star.mean(1) - T[:,0])**2))]
errorsCV = [np.average(np.abs(CV_star.mean(1) - CV[:,0]))]
T_stars = [T_star.mean(1)]
CV_stars = [CV_star.mean(1)]
print('RMSE:',errorsT[-1])
print('RMSE CV:',errorsCV[-1])
# list of available candidates for sample during active learning
available = list(range(X.shape[0]))
# start active learning
for i in range(N_AL):
# compute the entropy for the available candidates
ent = np.apply_along_axis(lambda x: ee.entropy(x[:,None]), 1, T_star[available])
idx_next = ent.argmax() # find the point of maximum entropy
# add it to the dataset
x_next, y_next = X[available][idx_next], Y[available][idx_next]
T_next = exact(x_next, y_next)
model.add_point(x_next, y_next, T_next)
available.remove(available[idx_next])
# and continue training
model.train_Adam_minibatch(5000, size = 96)
# predict and evaluate the error
T_star, CV_star = model.predict(X,Y)
plot_ensemble(T_star, CV_star, model.x_e, model.y_e, filename = 'results/AL_NNpara_%i.pdf' % (i+1))
errorsT.append(np.sqrt(np.average((T_star.mean(1) - T[:,0])**2)))
errorsCV.append(np.average(np.abs(CV_star.mean(1) - CV[:,0])))
print(i,'RMSE:',errorsT[-1])
print(i,'RMSE CV:',errorsCV[-1])