-
Notifications
You must be signed in to change notification settings - Fork 2
/
TF_MLP_GridSearch.py
executable file
·215 lines (174 loc) · 7.78 KB
/
TF_MLP_GridSearch.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
"""
PURPOSE:
Fully connected MLP Neural Network Regression implemented in TensorFlow (TF)
INPUTS:
REQUIRED:
-x File with genotype information
-y File with values you want to predict
-label Name of column in y with the value you want to predict (i.e. trait of interest)
-save Name to include in RESULTS file (i.e. what dataset are you running)
-cv File with CV folds specified
OPTIONAL:
-arc Desired NN architecture as comma separated layer sizes (i.e. 100,50 or 200,200,50)
-act What activation function to use (sigmoid (default), relu, elu)
-epochs Number of epochs to train on (default = 1000)
-lr Learning rate (default = 0.01)
-beta Regularization parameter (default = 0.01)
-JobID Which cv fold from the cv file do you want to run?
OUTPUTS:
-RESULTS Summary of results from the run located in the dir where the script was called.
Results will be appended to this file as they complete. Use -save to give
a run a unique identifier.
EXAMPLE ON HPCC:
Log on to development node with GPUs:
$ ssh dev-intel16-k80
Load linuxbrew, modules required by TF, & activate the TF python environment
$ source /opt/software/tensorflow/1.1.0/load_tf
Run example MLP (files in /mnt/home/azodichr/GitHub/TF-GenomicSelection/):
$ python TF_MLP_GridSearch.py -x geno.csv -y pheno.csv -label Yld_Env1 -cv CVFs.csv -save wheat -arc 100,50,20
"""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import sys, os
import numpy as np
import tensorflow as tf
import pandas as pd
from datetime import datetime
import timeit
start_time = timeit.default_timer()
timestamp = datetime.now().strftime('%Y-%m-%d %H:%M:%S')
# FUNCTIONS
def multilayer_perceptron(x, weights, biases, layer_number, activation_function):
layer = x
for l in range(1,layer_number+1):
weight_name = 'h' + str(l)
bias_name = 'b' + str(l)
layer = tf.add(tf.matmul(layer, weights[weight_name]), biases[bias_name])
if activation_function.lower() == 'sigmoid':
layer = tf.nn.sigmoid(layer)
elif activation_function.lower() == 'relu':
layer = tf.nn.relu(layer)
elif activation_function.lower() == 'elu':
layer = tf.nn.elu(layer)
else:
print("Given activation function is not supported")
quit()
out_layer = tf.matmul(layer, weights['out']) + biases['out']
return out_layer
#### Set default values #####
activation_function = 'sigmoid'
training_epochs = 1000
arc = 100,50
learning_rate = 0.01
beta = 0.01 # regularization parameter
SAVE = 'test'
for i in range (1,len(sys.argv),2):
if sys.argv[i] == "-x":
X_file = sys.argv[i+1]
if sys.argv[i] == "-y":
Y_file = sys.argv[i+1]
if sys.argv[i] == "-cv":
CVs = sys.argv[i+1]
if sys.argv[i] == "-JobID":
JobID = int(sys.argv[i+1])
if sys.argv[i] == "-label":
LABEL = sys.argv[i+1]
if sys.argv[i] == "-save":
SAVE = sys.argv[i+1]
if sys.argv[i] == "-act":
activation_function = sys.argv[i+1]
if sys.argv[i] == "-epochs":
training_epochs = int(sys.argv[i+1])
if sys.argv[i] == "-lr":
learning_rate = float(sys.argv[i+1])
if sys.argv[i] == "-beta":
beta = float(sys.argv[i+1])
if sys.argv[i] == "-arc": # Desired layer sizes comma separated (i.e. 100,50,20)
arc = sys.argv[i+1]
# Read in the desired architecture
arc = arc.strip().split(',')
archit = []
for a in arc:
archit.append(int(a))
layer_number = len(archit)
# Read in geno and pheno and remove non target phenotypes
x = pd.read_csv(X_file, sep=',', index_col = 0)
y = pd.read_csv(Y_file, sep=',', index_col = 0)
y = y[[LABEL]]
yhat = np.zeros(shape = y.shape)
final_training_error = []
final_testing_error = []
final_accuracy = []
cv_folds = pd.read_csv(CVs, sep=',', index_col=0)
for c_fold in range(1,3):
print('Starting cv set: ' + str(c_fold))
#cv = cv_folds['cv_' + str(JobID)]
cv = cv_folds['cv_' + str(c_fold)]
num_cvs = np.ptp(cv) + 1 # Range of values in cv (PeakToPeak)
training_error = []
for i in range(1,num_cvs+1):
print("Predicting cv fold %i" % i)
X_train = x[cv != i]
X_test = x[cv == i]
y_train = y[cv != i]
y_test = y[cv == i]
n_input = X_train.shape[1]
n_samples = X_train.shape[0]
n_classes = y_train.shape[1]
# TF Graph Input
nn_x = tf.placeholder(tf.float32, [None, n_input])
nn_y = tf.placeholder(tf.float32, [None, n_classes])
# Store layers weight & bias (default: mean=0, sd = 1)
weights = {}
biases = {}
weights['h1'] = tf.Variable(tf.random_normal([n_input, archit[0]]))
biases['b1'] = tf.Variable(tf.random_normal([archit[0]]))
for l in range(1,layer_number):
w_name = 'h' + str(l+1)
b_name = 'b' + str(l+1)
weights[w_name] = tf.Variable(tf.random_normal([archit[l-1], archit[l]]))
biases[b_name] = tf.Variable(tf.random_normal([archit[l]]))
weights['out'] = tf.Variable(tf.random_normal([archit[-1], n_classes]))
biases['out'] = tf.Variable(tf.random_normal([n_classes]))
# Construct model
pred = multilayer_perceptron(nn_x, weights, biases, layer_number, activation_function)
# Define loss and optimizer
loss = tf.reduce_mean(tf.square(pred - nn_y)) # Mean squared error
regularizer = tf.nn.l2_loss(weights['h1']) + tf.nn.l2_loss(weights['h2'])
loss = tf.reduce_mean(loss + beta * regularizer)
optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate).minimize(loss)
# Launch the graph
sess = tf.Session()
init = tf.global_variables_initializer()
sess.run(init)
for epoch in range(training_epochs):
sess.run(optimizer, feed_dict = {nn_x:X_train, nn_y:y_train})
c = sess.run(loss,feed_dict = {nn_x:X_train, nn_y:y_train})
if (epoch+1) % 250 == 0:
print("Epoch:", '%04d' % (epoch+1), "Cost=", "{:.9f}".format(c))
if epoch+1 == training_epochs:
training_error.append(c)
print('Final mse for training cv_%i: %.5f' % (i, c))
# Predict test set and add to yhat output
y_pred = sess.run(pred, feed_dict={nn_x: X_test})
yhat[cv == i] = y_pred
testing_mse = np.mean((np.array(y)[:,0] - yhat[:,0])**2)
cor = np.corrcoef(np.array(y)[:,0], yhat[:,0])
stop_time = timeit.default_timer()
final_training_error.append(np.mean(training_error))
final_testing_error.append(testing_mse)
final_accuracy.append(cor[0,1])
print('###################\nRESULTS\n###################\n')
print('Training error (MSE +/- stdev): %0.5f (%0.5f)' % (np.mean(final_training_error), np.std(final_training_error)))
print('Testing error (MSE +/- stdev): %0.5f (%0.5f)' % (np.mean(final_testing_error), np.std(final_testing_error)))
print('Accuracy (correlation coef +/- stdev): %.5f (%0.5f)' % (np.mean(final_accuracy), np.std(final_accuracy)))
print('\nRun time: %s' % str(stop_time - start_time))
if not os.path.isfile('RESULTS.txt'):
out2 = open('RESULTS.txt', 'a')
out2.write('DateTime\tDFs\tDFy\tTrait\tCV_Fold\tNumHidLay\tArchit\tActFun\tEpochs\tLearnRate\tBeta\tTrainError\tTrainErrorSTD\tTestError\tTestErrorSTD\tAccuracy\tAccuracySTD\n')
out2 = open('RESULTS.txt', 'a')
out2.write('%s\t%s\t%s\t%s\tall\t%i\t%s\t%s\t%i\t%0.5f\t%0.5f\t%0.5f\t%0.5f\t%0.5f\t%0.5f\t%0.5f\t%0.5f\n' % (
timestamp, X_file, Y_file, LABEL, layer_number, arc, activation_function, training_epochs, learning_rate, beta,
np.mean(final_training_error), np.std(final_training_error), np.mean(final_testing_error),
np.std(final_testing_error), np.mean(final_accuracy), np.std(final_accuracy)))