forked from WilkAndy/cnn_microseismic
-
Notifications
You must be signed in to change notification settings - Fork 0
/
learn.py
executable file
·342 lines (293 loc) · 16.2 KB
/
learn.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
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
#!/usr/bin/env python
# Name: learn.py
# Author: Andy Wilkins, [email protected], +61 7 3327 4497, Queensland Centre for Advanced Technologies, PO Box 883, Kenmore, Qld, 4069, Australia
# Year: 2019
# Software required: python2 software stack, listed below
# Language: python
# Program size: 21kB
import os
import sys
import numpy as np
from optparse import OptionParser, OptionValueError
from pandas import read_csv
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from keras.layers import Flatten
from keras.layers.convolutional import Conv2D
from keras.layers.convolutional import MaxPooling2D
from keras.optimizers import SGD
from keras.utils import np_utils
from keras.callbacks import EarlyStopping
from keras import backend as K
K.set_image_dim_ordering('th') # (channels, rows, cols) in X
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
# parse command line
p = OptionParser(usage="""usage: %prog [options] <data> <output_file_basename>
Attempts to learn info about <data>.
<data> is assumed to be a CSV file with a header indicated by "#"
<data> is assumed to have a line of the form
#Duration d NumChannels n
where d and n are integers that define the 'width' and 'height' of the 'images' in <data>
Here we use a CNN to attempt to distinguish a true even from a false-positive event
Eg:
./%prog -v -a -p -- cnn_data.csv out
""")
p.add_option("-v", action="store_true", dest="verbose", help="Verbose output")
p.add_option("-p", action="store_true", dest="plot", help="Plot accuracy and loss")
p.add_option("-a", action="store_true", dest="augment", help="Augment the training dataset using shifts, etc, of the data.")
p.add_option("--mask", action="store", default=0.1, type="float", help="Cut off this fraction of <data> from each end. Eg, if the 'width' of each 'image' is 100, and mask=0.1, then 10 pixels are cut from the beginning and the end of each image. The purpose of this is so that the images can be shifted to generate more data. Default=%default")
p.add_option("--min_training", action="store", default=2000, type="int", help="If the number of training images is less than this quantity, generate more by random shifts of the images. This option is irrelevant if the -a flag is not chosen. Default=%default")
p.add_option("--epochs", action="store", default=50, type="int", help="Number of epochs to use. Default=%default")
p.add_option("--permute", action="store", default=0.5, type="float", help="When augmenting the training set (with the -a flag), permute the data from original Joey channels (1,2,3), (8,9), (10, 11), with frequency given by this value, otherwise perform shifts. Eg, --permute=0.4 means 40% of the augmentations will be a random channel permute, while 60% of the augmentations will be data shifts. Default=%default")
(opts, args) = p.parse_args()
if len(args) != 2:
sys.stderr.write("Incorrect number of arguments. Run with '-h' to get help on this program\n")
sys.exit(1)
(data_file, output_file) = args
if opts.mask > 0.4:
sys.stderr.write("Currently you cannot cut more than 40% from each image.\nYou must choose mask less than 0.4\n")
sys.exit(1)
seed = 7
np.random.seed(seed)
if opts.verbose: sys.stdout.write("Obtaining Duration and NumChannels from " + data_file + "... ");
image_width = None
image_height = None
f = open(data_file, 'r')
for line in f:
if image_width != None and image_height != None:
break
if not line.strip():
continue
line = line.strip();
if not line.startswith("#"):
continue
if line.startswith("#Duration"):
line = line.split()
image_width = int(line[1])
image_height = int(line[3])
f.close()
if image_width == None or image_height == None:
sys.stderr.write("Cannot determine parameters. Exiting\n")
else:
if opts.verbose: sys.stdout.write("Done. (Width = " + str(image_width) + ", Height = " + str(image_height) + ")\n")
if opts.verbose: sys.stdout.write("Reading data file " + data_file + "... ");
data = read_csv(data_file, header=None, comment='#').values
if opts.verbose: sys.stdout.write("Done.\n")
if opts.verbose: sys.stdout.write("Extracting inputs and outputs from data... ");
num_inputs = len(data[0]) - 1
X_orig = data[:, 0:num_inputs]
Y_orig = data[:, num_inputs]
# reshape into the size of the "images"
X_orig = X_orig.reshape(X_orig.shape[0], 1, image_height, image_width)
if opts.verbose: sys.stdout.write("Done.\n")
def create_cnn_model(training_images):
if opts.verbose: sys.stdout.write(" Defining model... ");
num_training = len(training_images)
height = len(training_images[0][0])
width = len(training_images[0][0][0])
model = Sequential()
version = 10 # Other versions removed for clarity. Older versions had different Conv2D, etc, layers
if opts.verbose: sys.stdout.write("Using version " + str(version) + ".")
if version == 10:
# experimenting with different Conv2D
# trying to reduce over-fitting compared with version==7
model.add(Conv2D(20, (height, 10), input_shape = (1, height, width), activation='relu'))
model.add(Dropout(0.2))
model.add(Conv2D(10, (1, 30), input_shape = (1, height, width), activation='relu'))
model.add(Dropout(0.2))
model.add(Flatten())
model.add(Dense(int(0.1 * num_training), kernel_initializer='glorot_uniform', activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(1, kernel_initializer='glorot_uniform', activation='sigmoid'))
if opts.verbose: sys.stdout.write(" Done.\n");
if opts.verbose: sys.stdout.write(" Compiling model... ");
# stocastic gradient descent
# with initial learningrate=0.01
# and subsequent learningrate = (initial learningrate)/(1 + decay * epochs) [I think!]
# and momentum=0.9
# using nesterov
sgd = SGD(lr=0.01, momentum=0.9, decay=0.01/50, nesterov=True)
model.compile(loss='binary_crossentropy', optimizer=sgd, metrics=['accuracy'])
if opts.verbose: sys.stdout.write(" Done.\n");
return model
def augment(X, Y, num_augmentations, permute_chance):
# take X and Y, and shift a random choice of the Xs, padding with zeroes, to produce an augmented dataset
# This returns just the extra data generated, not the concatenation of the original and the extra.
extra_X = []
extra_Y = []
num_unaugmented = len(Y)
# remember that
# X[i][j][k][l] is intensity of: the i^th input 'image', the j^th colour (we only have 1 colour), the k^th row (we have 8 rows), and the l^th column (1500 columns)
img_width = len(X[0][0][0])
max_shift = int(img_width * opts.mask - 1.5)
for i in range(int(num_augmentations * (1.0 - permute_chance))):
x_to_shift = np.random.randint(0, num_unaugmented)
shift_by = np.random.randint(-max_shift, max_shift + 1)
shifted = np.roll(X[x_to_shift][0], shift_by, axis=1) # axis=1 means shift intensities of each row
shifted = np.array([shifted]) # put back 'color'
extra_X.append(shifted)
extra_Y.append(Y[x_to_shift])
for i in range(int(num_augmentations * permute_chance)):
x_to_permute = np.random.randint(0, num_unaugmented)
permute_type = np.random.ranf()
permuted = np.copy(X[x_to_permute][0])
if permute_type < 0.2:
# swap rows 0 and 1 (Joey channels 1 and 2)
permuted[[0, 1]] = permuted[[1, 0]]
if permute_type < 0.4:
# swap rows 0 and 2 (Joey channels 1 and 3)
permuted[[0, 2]] = permuted[[2, 0]]
if permute_type < 0.6:
# swap rows 1 and 2 (Joey channels 2 and 3)
permuted[[1, 2]] = permuted[[2, 1]]
elif permute_type < 0.8:
# swap rows 4 and 5 (Joey channels 8 and 9)
permuted[[4, 5]] = permuted[[5, 4]]
else:
# swap rows 6 and 7 (Joey channels 10 and 11)
permuted[[6, 7]] = permuted[[7, 6]]
extra_X.append(np.array([permuted])) # put back 'color'
extra_Y.append(Y[x_to_permute])
return (np.array(extra_X), np.array(extra_Y))
if opts.verbose: sys.stdout.write("\n\nPERFORMING A STRATIFIED k-FOLD CROSS VALIDATION WITH 10 FOLDS\n")
summary_file = open(output_file + ".txt", 'w')
summary_file.write("PERFORMING A STRATIFIED k-FOLD CROSS VALIDATION WITH 10 FOLDS\n")
summary_file.close()
kfold = StratifiedKFold(n_splits = 10, shuffle = True, random_state = seed)
fold_num = 0
fig_num = 0
for train, test in kfold.split(X_orig, Y_orig):
fold_num += 1
# split off 10% of the training data to act as validation
X_train, X_validation, Y_train, Y_validation = train_test_split(X_orig[train], Y_orig[train], test_size=0.1, random_state=seed, stratify=Y_orig[train])
indices_where_Y_train_zero = np.where(Y_train == 0)[0]
indices_where_Y_train_one = np.where(Y_train == 1)[0]
if opts.verbose: sys.stdout.write("\n\nFOLD " + str(fold_num) + ".\n " + str(len(indices_where_Y_train_one)) + " correct triggers in training (out of " + str(len(np.where(Y_orig == 1)[0])) + ") " + str(len(indices_where_Y_train_zero)) + " false-positives in training (out of " + str(len(np.where(Y_orig == 0)[0])) + ")\n")
if opts.augment:
# generate the same number of ones and zeroes in the training set
num_excess = len(indices_where_Y_train_one) - len(indices_where_Y_train_zero)
if num_excess > 0:
if opts.verbose: sys.stdout.write(" Too many ones (correct positives). Generating " + str(num_excess) + " zeroes (false positives)\n")
X_gen, Y_gen = augment(X_train[indices_where_Y_train_zero], Y_train[indices_where_Y_train_zero], num_excess, 0.0)
else:
if opts.verbose: sys.stdout.write(" Too many zeroes (false positives). Generating " + str(-num_excess) + " ones (correct positives)\n")
X_gen, Y_gen = augment(X_train[indices_where_Y_train_one], Y_train[indices_where_Y_train_one], -num_excess, 0.0)
X_train = np.concatenate((X_train, X_gen))
Y_train = np.concatenate((Y_train, Y_gen))
if opts.verbose: sys.stdout.write(" Now number of zeroes = " + str(len(np.where(Y_train == 0)[0])) + " and number of ones = " + str(len(np.where(Y_train == 1)[0])) + "\n")
# ensure there are enough examples in the training set
if len(Y_train) < opts.min_training:
if opts.verbose: sys.stdout.write(" Not enough data. Generating " + str(opts.min_training - len(Y_train)) + " extra 'images' by shifting and swapping\n")
X_gen, Y_gen = augment(X_train, Y_train, opts.min_training - len(Y_train), opts.permute)
X_train = np.concatenate((X_train, X_gen))
Y_train = np.concatenate((Y_train, Y_gen))
if opts.verbose: sys.stdout.write(" Now number of zeroes = " + str(len(np.where(Y_train == 0)[0])) + " and number of ones = " + str(len(np.where(Y_train == 1)[0])) + "\n")
# trim the training set by removing from each end
num_to_remove = int(len(X_train[0][0][0]) * opts.mask - 0.5)
if opts.verbose: sys.stdout.write(" Trimming " + str(num_to_remove) + " from each end of the training 'images'\n")
X_train = np.delete(X_train, np.s_[:num_to_remove], axis = 3)
X_train = np.delete(X_train, np.s_[-num_to_remove:], axis = 3)
new_image_width = len(X_train[0][0][0])
if opts.verbose: sys.stdout.write(" Trimming " + str(num_to_remove) + " from each end of the validation 'images'\n")
X_validation = np.delete(X_validation, np.s_[:num_to_remove], axis = 3)
X_validation = np.delete(X_validation, np.s_[-num_to_remove:], axis = 3)
if opts.verbose: sys.stdout.write(" Now the width of the images is " + str(new_image_width) + "\n")
cnn_model = create_cnn_model(X_train)
if opts.verbose: sys.stdout.write(" Training:\n")
if True:
es = EarlyStopping(monitor = 'val_loss', mode='min', patience=10, min_delta=0.0, verbose=1)
history = cnn_model.fit(X_train, Y_train, validation_data=(X_validation, Y_validation), epochs=opts.epochs, batch_size=50, verbose=1, callbacks=[es])
else:
history = cnn_model.fit(X_train, Y_train, validation_data=(X_validation, Y_validation), epochs=opts.epochs, batch_size=50, verbose=1)
if opts.plot:
fig_num += 1
plt.figure(fig_num)
plt.plot(history.history['acc'])
plt.plot(history.history['val_acc'])
plt.title('model accuracy. Fold = ' + str(fold_num))
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.savefig(output_file + "_accuracy" + str(fold_num) + ".png")
fig_num += 1
plt.figure(fig_num)
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.title('model loss. Fold = ' + str(fold_num))
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.savefig(output_file + "_loss" + str(fold_num) + ".png")
# trim the test data 'images' so they are the same size as training data
X_test = X_orig[test]
Y_test = Y_orig[test]
X_test = np.delete(X_test, np.s_[:num_to_remove], axis = 3)
X_test = np.delete(X_test, np.s_[-num_to_remove:], axis = 3)
if opts.verbose: sys.stdout.write(" On the withheld data:\n")
summary_file = open(output_file + ".txt", 'a')
summary_file.write("Fold" + str(fold_num) + "\n")
y_predict = cnn_model.predict_classes(X_test)
num_true_events = 0
num_predicted_events = 0
for i in range(len(Y_test)):
if Y_test[i] == 1:
num_true_events += 1
if y_predict[i] == 1:
num_predicted_events += 1
if opts.verbose: sys.stdout.write(" Number of true events = " + str(num_true_events) + " (total size of test data = " + str(len(Y_test)) + "). Number of predicted events = " + str(num_predicted_events) + "\n")
summary_file.write(" Number of true events = " + str(num_true_events) + " (total size of test data = " + str(len(Y_test)) + "). Number of predicted events = " + str(num_predicted_events) + "\n")
num_incorrect = 0
num_correct = 0
for i in range(len(Y_test)):
if Y_test[i] == 1:
if y_predict[i] == 1:
num_correct += 1
else:
num_incorrect += 1
if opts.verbose: sys.stdout.write(" Fraction of original true events classified as a true event " + str(100 * float(num_correct)/ float(num_correct + num_incorrect)) + "\n")
summary_file.write(" Fraction of original true events classified as a true event " + str(100 * float(num_correct)/ float(num_correct + num_incorrect)) + "\n")
num_incorrect = 0
num_correct = 0
for i in range(len(Y_test)):
if Y_test[i] == 0:
if y_predict[i] == 0:
num_correct += 1
else:
num_incorrect += 1
if opts.verbose: sys.stdout.write(" Fraction of false-positive events classified as a 'not an event' " + str(100 * float(num_correct)/ float(num_correct + num_incorrect)) + "\n")
summary_file.write(" Fraction of false-positive events classified as a 'not an event' " + str(100 * float(num_correct)/ float(num_correct + num_incorrect)) + "\n")
summary_file.close()
summary_file = open(output_file + ".txt", 'r')
tot_num_true = 0
tot_num_pred = 0
frac_true = []
frac_false = []
for line in summary_file:
if not line.strip():
continue
line = line.strip()
s = line.split("Number of true events = ")
if len(s) > 1:
tot_num_true += int(s[1].split()[0])
s = line.split("Number of predicted events = ")
if len(s) > 1:
tot_num_pred += int(s[1].split()[0])
s = line.split("Fraction of original true events classified as a true event ")
if len(s) > 1:
frac_true.append(float(s[1].split()[0]))
s = line.split("Fraction of false-positive events classified as a 'not an event' ")
if len(s) > 1:
frac_false.append(float(s[1].split()[0]))
summary_file.close()
summary_file = open(output_file + ".txt", 'a')
summary_file.write("SUMMARY\n")
summary_file.write("Total number of original true events = " + str(tot_num_true) + "\n")
summary_file.write("Total number of predicted true = " + str(tot_num_pred) + "\n")
summary_file.write("Overal fraction of original true events classified as a true event " + str(np.round(np.mean(frac_true), decimals = 1)) + " +/- " + str(np.round(np.std(frac_true), decimals = 1)) + "\n")
summary_file.write("Overal fraction of false-positive events classified as 'not an event' " + str(np.round(np.mean(frac_false), decimals = 1)) + " +/- " + str(np.round(np.std(frac_false), decimals = 1)) + "\n")
summary_file.close()
sys.exit(0)
sys.exit(0)