import numpy as np import tensorflow as tf import random import os from keras import backend as K from numpy import genfromtxt import matplotlib.pyplot as plt import pandas as pd from keras.regularizers import l2 from keras.utils import to_categorical from keras.optimizers import SGD, Nadam, Adam from keras.layers import Input, Dense, Activation, Dropout, BatchNormalization, Flatten, merge from keras.callbacks import ModelCheckpoint, EarlyStopping from keras.models import Model, model_from_json from sklearn.metrics import roc_auc_score from sklearn.metrics import precision_recall_fscore_support from sklearn.model_selection import StratifiedShuffleSplit from sklearn.svm import SVC from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import accuracy_score from sklearn.metrics import roc_auc_score from sklearn.svm import SVC from sklearn.model_selection import GridSearchCV from keras.wrappers.scikit_learn import KerasClassifier from sklearn.model_selection import StratifiedKFold from keras.constraints import maxnorm from sklearn.model_selection import ParameterGrid from sklearn.model_selection import train_test_split from keras.layers import concatenate from sklearn.metrics import roc_auc_score,accuracy_score,roc_curve,auc import argparse from sklearn.metrics import roc_auc_score,auc,roc_curve,f1_score,accuracy_score from lifelines.utils import concordance_index def parse_args(): parser = argparse.ArgumentParser() parser.add_argument('--seed',default=155,type=int,help='random seed') parser.add_argument('--i',default=0,type=int,help='num_of_ensemble') return parser.parse_args() def create_model(neuron, activation, l2term, lr, lr_decay, rate, norm, is_a, is_c, i): neuron_a = 20 activation_a = 'relu' l2term_a = 0.0001 rate_a = 0.05 norm_a = 1 input_array = Input(shape=(20,)) hidden1_array_ = Dense(neuron_a,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_a), activation=activation_a, kernel_constraint=maxnorm(norm_a),trainable=is_a)(input_array) hidden1_array = Dropout(rate_a)(hidden1_array_) hidden2_array_ = Dense(neuron_a,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_a), activation=activation_a, kernel_constraint=maxnorm(norm_a),trainable=is_a)(hidden1_array) hidden2_array = Dropout(rate_a)(hidden2_array_) hidden3_array_ = Dense(neuron_a,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_a), activation=activation_a, kernel_constraint=maxnorm(norm_a),trainable=is_a)(hidden2_array) hidden3_array = Dropout(rate_a)(hidden3_array_) Output_array = Dense(1, activation='sigmoid',trainable=is_a)(hidden3_array) model_array = Model(inputs=input_array, outputs=Output_array) if not is_a: model_array.load_weights('../../model/bimodal_ens/breast/bimodal_array_' + str(i) + '.hdf5') else: return model_array neuron_c = 5 activation_c = 'tanh' l2term_c = 0.0001 rate_c = 0.05 norm_c = 1 input_clinical = Input(shape=(10,)) hidden1_clinical_ = Dense(neuron_c,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_c), activation=activation_c, kernel_constraint=maxnorm(norm_c),trainable=is_c)(input_clinical) hidden1_clinical = Dropout(rate_c)(hidden1_clinical_) hidden2_clinical_ = Dense(neuron_c,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_c), activation=activation_c, kernel_constraint=maxnorm(norm_c),trainable=is_c)(hidden1_clinical) hidden2_clinical = Dropout(rate_c)(hidden2_clinical_) Output_clinical = Dense(1, activation='sigmoid')(hidden2_clinical) model_clinical = Model(inputs=input_clinical, outputs=Output_clinical) if not is_c: model_clinical.load_weights('../../model/bimodal_ens/breast/bimodal_clinical_' + str(i) + '.hdf5') else: return model_clinical merge_vector = concatenate([hidden3_array,hidden2_clinical],axis=-1) # merge_vector=merge([hidden3_array,hidden2_clinical],mode='concat',concat_axis=1) hidden1_merge_ = Dense( neuron, kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term), activation=activation, kernel_constraint=maxnorm(norm))(merge_vector) hidden1_merge = Dropout(rate)(hidden1_merge_) Output_merge = Dense(1, activation='sigmoid')(hidden1_merge) model_merge = Model(inputs=[input_array, input_clinical],outputs=Output_merge) return model_merge def create_model_2(neuron, activation, l2term, lr, lr_decay, rate, norm, is_a, is_c, i): neuron_a = 20 activation_a = 'relu' l2term_a = 0.0001 rate_a = 0.05 norm_a = 1 input_array = Input(shape=(20,)) hidden1_array_ = Dense(neuron_a,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_a), activation=activation_a, kernel_constraint=maxnorm(norm_a),trainable=is_a)(input_array) hidden1_array = Dropout(rate_a)(hidden1_array_) hidden2_array_ = Dense(neuron_a,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_a), activation=activation_a, kernel_constraint=maxnorm(norm_a),trainable=is_a)(hidden1_array) hidden2_array = Dropout(rate_a)(hidden2_array_) hidden3_array_ = Dense(neuron_a,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_a), activation=activation_a, kernel_constraint=maxnorm(norm_a),trainable=is_a)(hidden2_array) hidden3_array = Dropout(rate_a)(hidden3_array_) Output_array = Dense(1, activation='sigmoid',trainable=is_a)(hidden3_array) model_array = Model(inputs=input_array, outputs=Output_array) if not is_a: model_array.load_weights('../../model/bimodal_ens/breast/bimodal_array_' + str(i) + '.hdf5') else: return model_array neuron_c = 5 activation_c = 'tanh' l2term_c = 0.0001 rate_c = 0.05 norm_c = 1 input_clinical = Input(shape=(10,)) hidden1_clinical_ = Dense(neuron_c,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_c), activation=activation_c, kernel_constraint=maxnorm(norm_c),trainable=is_c)(input_clinical) hidden1_clinical = Dropout(rate_c)(hidden1_clinical_) hidden2_clinical_ = Dense(neuron_c,kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term_c), activation=activation_c, kernel_constraint=maxnorm(norm_c),trainable=is_c)(hidden1_clinical) hidden2_clinical = Dropout(rate_c)(hidden2_clinical_) Output_clinical = Dense(1, activation='sigmoid')(hidden2_clinical) model_clinical = Model(inputs=input_clinical, outputs=Output_clinical) if not is_c: model_clinical.load_weights('../../model/bimodal_ens/breast/bimodal_clinical_' + str(i) + '.hdf5') else: return model_clinical merge_vector = concatenate([hidden3_array,hidden2_clinical],axis=-1) # merge_vector=merge([hidden3_array,hidden2_clinical],mode='concat',concat_axis=1) hidden1_merge_ = Dense( neuron, kernel_initializer='glorot_normal', kernel_regularizer=l2(l2term), activation=activation, kernel_constraint=maxnorm(norm))(merge_vector) hidden1_merge = Dropout(rate)(hidden1_merge_) Output_merge = Dense(1, activation='sigmoid')(hidden1_merge) model_merge = Model(inputs=[input_array, input_clinical],outputs=Output_merge) model_merge.load_weights('../../model/bimodal_ens/breast/bimodal_merge_' + str(i) + '.hdf5') return model_merge def main(args): # random seed settings os.environ['PYTHONHASHSEED'] = '\'' + str(args.seed) + '\'' np.random.seed(int(args.seed)) random.seed(int(args.seed)) session_conf = tf.ConfigProto(intra_op_parallelism_threads=1, inter_op_parallelism_threads=1) tf.set_random_seed(int(args.seed)) sess = tf.Session(graph=tf.get_default_graph(), config=session_conf) K.set_session(sess) # load data & normalization data = np.load('../../data/breast/metabric.npz',allow_pickle=True) x_train = data['X_rest'] y_train = data['y_rest'] x_test = data['X_test'] y_test = data['y_test'] num_clinical = 10 MIN = x_train.min(axis=0) MAX = x_train.max(axis=0) for col in [0,2,6]: x_train[:,col] = (x_train[:,col] - MIN[col]) / (MAX[col] - MIN[col]) x_test[:,col] = (x_test[:,col] - MIN[col]) / (MAX[col] - MIN[col]) for col in range(30): if col < num_clinical: continue x_train[:,col] = (x_train[:,col] - MIN[col]) / (MAX[col] - MIN[col]) x_test[:,col] = (x_test[:,col] - MIN[col]) / (MAX[col] - MIN[col]) i = int(args.i) model_a = create_model( neuron=10, activation='tanh', l2term=0.0001, lr=0.00002, lr_decay=0.004, rate=0.05, norm=1, is_a=True, is_c=False, i=i) model_a.compile(loss='binary_crossentropy', optimizer = Nadam(lr=0.0002,schedule_decay=0.004)) model_a.fit(x_train[:,num_clinical:],y_train,batch_size=10,epochs=200,verbose=1) model_a.save('../../model/bimodal_ens/breast/bimodal_array_' + str(i) + '.hdf5') model_a.load_weights('../../model/bimodal_ens/breast/bimodal_array_' + str(i) + '.hdf5') logits = model_a.predict(x_test[:,num_clinical:]) fpr,tpr,thr = roc_curve(y_test.astype(int),logits,pos_label=1) print(auc(fpr,tpr)) model_c = create_model( neuron=10, activation='tanh', l2term=0.0001, lr=0.00002, lr_decay=0.004, rate=0.05, norm=1, is_a=False, is_c=True, i=i) model_c.compile(loss='binary_crossentropy', optimizer = Nadam(lr=0.002,schedule_decay=0.004)) model_c.fit(x_train[:,:num_clinical],y_train,batch_size=10,epochs=800,verbose=1) model_c.save('../../model/bimodal_ens/breast/bimodal_clinical_' + str(i) + '.hdf5') model_c.load_weights('../../model/bimodal_ens/breast/bimodal_clinical_' + str(i) + '.hdf5') logits = model_c.predict(x_test[:,:num_clinical]) fpr,tpr,thr = roc_curve(y_test.astype(int),logits,pos_label=1) print(auc(fpr,tpr)) model_b = create_model( neuron=10, activation='tanh', l2term=0.0001, lr=0.00002, lr_decay=0.004, rate=0.05, norm=1, is_a=False, is_c=False, i=i) model_b.compile(loss='binary_crossentropy', optimizer = Nadam(lr=0.00002,schedule_decay=0.004)) model_b.fit([x_train[:,num_clinical:],x_train[:,:num_clinical]],y_train,batch_size=10,epochs=100,verbose=1) model_b.save('../../model/bimodal_ens/breast/bimodal_merge_' + str(i) + '.hdf5') model_b.load_weights('../../model/bimodal_ens/breast/bimodal_merge_' + str(i) + '.hdf5') model_b2 = create_model_2( neuron=10, activation='tanh', l2term=0.0001, lr=0.00002, lr_decay=0.004, rate=0.05, norm=1, is_a=False,is_c=False, i=i) model_b2.compile(loss='binary_crossentropy', optimizer = Nadam(lr=0.00002,schedule_decay=0.004)) model_b2.fit([x_train[:,num_clinical:],x_train[:,:num_clinical]],y_train,batch_size=10,epochs=100,verbose=1) model_b2.save('../../model/bimodal_ens/breast/bimodal_merge_2_' + str(i) + '.hdf5') model_b2.load_weights('../../model/bimodal_ens/breast/bimodal_merge_2_' + str(i) + '.hdf5') logits = model_b2.predict([x_test[:,num_clinical:],x_test[:,:num_clinical]]) fpr,tpr,thr = roc_curve(y_test.astype(int),logits,pos_label=1) print(auc(fpr,tpr)) # best threshold # valid_num = int(np.shape(x_train)[0] * 0.25) # x_valid = x_train[-valid_num:,:] # y_valid = y_train[-valid_num:] # y_pred_valid = model_b2.predict([x_valid[:,num_clinical:],x_valid[:,:num_clinical]]) # fpr,tpr,thr = roc_curve(y_valid,y_pred_valid) # thr_best = thr[np.argmax(np.subtract(tpr,fpr))] # print(thr_best) if __name__ == "__main__": args = parse_args() main(args)