Added initial version of the Pilot 2 Benchmark 2 keras code.

bvanessen · bvanessen · commit f8a964cd10fd · 2017-03-01T07:55:13.000-08:00
diff --git a/P2B2/__init__.py b/P2B2/__init__.py
@@ -0,0 +1,4 @@
+print 'import candle_rnn'
+import os,sys
+HOME=os.environ['HOME']
+sys.path.append('%s/Work/Python/Git_Folder/caffe-tools/Newkeras/keras'%HOME)
diff --git a/P2B2/__main__.py b/P2B2/__main__.py
@@ -0,0 +1,150 @@
+import theano
+import numpy as np
+import scipy as sp
+import pickle
+import sys,os
+import glob
+import optparse
+import matplotlib
+matplotlib.use('TKAgg')
+import pylab as py
+py.ion()
+HOME=os.environ['HOME']
+def parse_list(option, opt, value, parser):
+  setattr(parser.values, option.dest, value.split(','))
+
+if __name__=="__main__":
+### Hyperparameters and model save path
+	parser=optparse.OptionParser()
+	parser.add_option("--train", action="store_true",dest="train_bool",default=False,help="Invoke training")
+	parser.add_option("--learning-rate",help="learning rate",dest="learning_rate",type=float,default=0.1)
+	parser.add_option("--noise-factor",help="noise",dest="noise_factor",type=float,default=0.0)
+	parser.add_option("--cool",action="store_true",dest="cool",default=False,help="Cool Learning Rate")
+	parser.add_option("--epochs",help="epochs",dest="epochs",type=int,default=1)
+	parser.add_option("--batch-size",help="batch size",dest="batch_size",type=int,default=1)
+	parser.add_option("--look-back",help="look back time window",dest="look_back",type=int,default=1)
+	parser.add_option("--home-dir",help="Home Directory",dest="home_dir",type=str,default='/Users/talathi1/Work/Python/Git_Folder/caffe-tools/keras')
+	parser.add_option("--save-dir",help="Save Directory",dest="save_path",type=str,default=None)
+	parser.add_option("--model-file",help="Trained Model Pickle File",dest="weight_path",type=str,default=None)
+	parser.add_option("--memo",help="Memo",dest="base_memo",type=str,default=None)
+	parser.add_option("--seed", action="store_true",dest="seed",default=False,help="Random Seed")
+	(opts,args)=parser.parse_args()
+
+	## set the seed
+	if opts.seed:
+		np.random.seed(7)
+	else:
+		np.random.seed(np.random.randint(10000))
+
+	## Set paths
+	if not os.path.isdir(opts.home_dir):
+		print ('Keras home directory not set')
+		sys.exit(0)
+	sys.path.append('home_dir')
+	
+	import candle_helper_functions as hf
+	reload(hf)
+	maps=hf.autoencoder_preprocess()
+	
+	## Import keras modules
+	from keras.optimizers import SGD,RMSprop,Adam
+	from keras.datasets import mnist
+	from keras.callbacks import LearningRateScheduler,ModelCheckpoint
+	from keras.regularizers import l2,WeightRegularizer
+	from keras import callbacks
+	from keras.layers.advanced_activations import ELU
+	from keras.preprocessing.image import ImageDataGenerator
+
+	batch_size = opts.batch_size
+##### Read Data ########
+	print ('Reading Data...')
+	data_file='%s/Research/DeepLearning/ECP CANDLE/Benchmarks/Benchmarks.git/P2B2/sim-numpy.npy'%HOME ### can code to read at the terminal
+	print 'Data File: %s' %data_file
+	print 'Data Format: [Num Samples, Num Molecules, Num Atoms, Position]'
+	
+	X=np.load(data_file) ### Data is: Samples, Molecules, Atoms, x-pos,y-pos,z-pos
+	## Take center of mass for atoms:
+	X_A=X.mean(axis=2) ## Data is: Samples, Molecules, x-pos,y-pos,z-pos
+	data=X_A[0:-1,:,2] ## only consider z-dimension
+	X_train,y_train=hf.create_dataset(data,opts.look_back,look_forward=1) ## convert data to a sequence 
+	temporal_dim=X_train.shape[1]
+	input_dim=X_train.shape[2]
+	subset_sample_weight=np.ones((X_train.shape[0],1))
+	sample_weight=np.zeros((X_train.shape[0],opts.look_back))
+	sample_weight[:,0:1]=subset_sample_weight
+
+	print('X_train type and shape:', X_train.dtype, X_train.shape)
+	print('X_train.min():', X_train.min())
+	print('X_train.max():', X_train.max())
+
+### Define Model, Solver and Compile ##########
+	print ('Define the model and compile')
+	opt=Adam(lr=opts.learning_rate)
+
+	print ('using mlp network')
+	model_type='mlp'
+	hidden_layers=[512,256,128,64,32,16]
+	recurrent_layers=[16,16,16]
+	## Model is a Autoencoder-RNN network
+	model=hf.rnn_dense_auto(weights_path=None,T=temporal_dim,D=input_dim,nonlinearity='relu',hidden_layers=hidden_layers,recurrent_layers=recurrent_layers)
+		
+	memo='%s_%s_%0.5f'%(opts.base_memo,model_type,opts.learning_rate)
+
+	print 'Autoencoder Regression problem'
+	model.compile(optimizer=opt, loss='mean_squared_error',sample_weight_mode="temporal")
+	model.summary()	## print model summary in details
+	#sys.exit(0)
+
+#### Print Compact Model Configuration ###########
+	# num_layers=len(model.layers)
+	# print '*'*10,'Model Configuration','*'*10
+	# for i  in range(len(model.layers)):	
+	# 	print i,': ',model.layers[i].name, ':', model.layers[i].output_shape[:]
+
+### Set up for Training and Validation
+	total_epochs = opts.epochs
+	initial_lrate=opts.learning_rate
+	if opts.cool:
+		drop=0.5
+	else:
+		drop=1.0
+	
+	epochs_drop=1+int(np.floor(total_epochs/3))
+		
+	def step_decay(epoch):
+		global initial_lrate,epochs_drop,drop
+		lrate = initial_lrate * np.power(drop, np.floor((1+epoch)/epochs_drop))
+		return lrate
+	lr_scheduler = LearningRateScheduler(step_decay)
+
+#### Train the Model
+	if opts.train_bool:
+		history = callbacks.History()
+		if opts.save_path !=None:
+			model_file='%s/%s.hdf5'%(opts.save_path,memo)
+			checkpointer=ModelCheckpoint(filepath=model_file, verbose=1)
+			callbacks=[history,lr_scheduler,checkpointer]
+		else:
+			callbacks=[history,lr_scheduler]
+		model.fit(X_train, y_train, batch_size=batch_size,shuffle=False,nb_epoch=total_epochs,callbacks=callbacks,verbose=1,sample_weight=sample_weight)
+		
+		loss_data={'train': history.history['loss']}
+		if opts.save_path!=None:
+			loss_file='%s/%s.pkl'%(opts.save_path,memo)
+			o=open(loss_file,'wb')
+			pickle.dump(loss_data,o)
+			o.close()
+
+		## Generate model forecast figure  
+		x=X_train[0:1]
+		xmod=x.reshape(x.shape[1],x.shape[2])
+		yf=hf.generate_timedistributed_forecast(model,x,X_train.shape[0]+opts.look_back)
+		yt=yt=y_train[:,0,:]
+		ytn=np.vstack([xmod,yt])
+		py.figure();py.plot(ytn.mean(axis=1))
+		py.hold('on');py.plot(yf.mean(axis=1))
+
+
+
+	
+	
diff --git a/P2B2/candle_helper_functions.py b/P2B2/candle_helper_functions.py
@@ -0,0 +1,234 @@
+from __future__ import absolute_import
+import theano
+import matplotlib
+if 'MACOSX' in matplotlib.get_backend().upper():
+  matplotlib.use('TKAgg')
+import pylab as py
+py.ion() ## Turn on plot visualization
+
+import gzip,pickle
+import numpy as np
+from PIL import Image
+import cv2
+import keras.backend as K
+K.set_image_dim_ordering('th')
+from keras.layers import Input, merge, TimeDistributed,LSTM,GRU,RepeatVector
+from keras.models import Sequential,Model
+from keras.layers.core import Flatten, Dense, Dropout, Activation, Reshape
+from keras.initializations import normal, identity, he_normal,glorot_normal,glorot_uniform,he_uniform
+from keras.layers.normalization import BatchNormalization
+import threading
+
+
+############# Define Data Generators ################
+class ImageNoiseDataGenerator(object):
+    '''Generate minibatches with
+    realtime data augmentation.
+    '''
+    def __init__(self,corruption_level=0.5):
+
+        self.__dict__.update(locals())
+        self.p=corruption_level
+        self.lock = threading.Lock()
+
+    def _flow_index(self, N, batch_size=32, shuffle=False, seed=None):
+        b = 0
+        total_b = 0
+        while 1:
+            if b == 0:
+                if seed is not None:
+                    np.random.seed(seed + total_b)
+
+                if shuffle:
+                    index_array = np.random.permutation(N)
+                else:
+                    index_array = np.arange(N)
+
+            current_index = (b * batch_size) % N
+            if N >= current_index + batch_size:
+                current_batch_size = batch_size
+            else:
+                current_batch_size = N - current_index
+
+            if current_batch_size == batch_size:
+                b += 1
+            else:
+                b = 0
+            total_b += 1
+            yield index_array[current_index: current_index + current_batch_size], current_index, current_batch_size
+
+    def flow(self, X, y, batch_size=32, shuffle=False, seed=None):
+        assert len(X) == len(y)
+        self.X = X
+        self.y = y
+        self.flow_generator = self._flow_index(X.shape[0], batch_size, shuffle, seed)
+        return self
+
+    def __iter__(self):
+        # needed if we want to do something like for x,y in data_gen.flow(...):
+        return self
+
+    def next(self):
+        # for python 2.x
+        # Keep under lock only the mechainsem which advance the indexing of each batch
+        # see # http://anandology.com/blog/using-iterators-and-generators/
+        with self.lock:
+            index_array, current_index, current_batch_size = next(self.flow_generator)
+        # The transformation of images is not under thread lock so it can be done in parallel
+        bX = np.zeros(tuple([current_batch_size] + list(self.X.shape)[1:]))
+        for i, j in enumerate(index_array):
+            x = self.X[j]
+            x = self.insertnoise(x,corruption_level=self.p)
+            bX[i] = x
+        bY = self.y[index_array]
+        return bX, bY
+
+    def __next__(self):
+        # for python 3.x
+        return self.next()
+
+    def insertnoise(self,x,corruption_level=0.5):
+        return np.random.binomial(1,1-corruption_level,x.shape)*x
+
+
+### Generate RNN compatible dataset
+def create_dataset(dataset, look_back=1,look_forward=1):
+    ## input is np.array of dim T,D
+    #output is np.array X: N,look_back,D and Y: N,D
+    # where N=T-look_back-1
+    assert(look_back>=look_forward)
+    dataX, dataY = [], []
+    for i in range(len(dataset)-look_back-(look_forward-1)-1):
+        a = dataset[i:(i+look_back), :]
+        dataX.append(a)
+        dataY.append(dataset[i + look_back:i+look_back+look_forward, :])
+    dataX=np.array(dataX)
+    dataY=np.array(dataY)    
+    if look_back-look_forward>0:
+        dataY_mod=np.zeros((dataY.shape[0],dataX.shape[1],dataY.shape[2]))
+        dataY_mod[:,0:dataY.shape[1],:]=dataY
+    else:
+        dataY_mod=dataY
+    return dataX, dataY_mod
+
+def generate_timedistributed_forecast(model,x,prediction_length=10):
+    ## to be used when rnn is used for sequence to sequence mapping
+    N,T,D=x.shape
+    x_data=x[0,:,:].copy()
+    x_revise=x.copy()
+    for i in range(prediction_length):
+        y_pred=model.predict(x_revise[0:1,:,:],batch_size=1)
+        yf=y_pred[:,0,:]
+        #print 'prediction:',yf
+        x_data=np.vstack([x_data,yf])
+        
+        #x=x.reshape(T,D) ## assume N=1 ... i.e. one sample
+        #print 'data before prediction:\n',x
+        x_revise[0,0:T-1,:]=x_revise[0,1:T,:]
+        x_revise[0,T-1:T,:]=yf
+        #x=x.reshape(1,T,D)
+        #print 'data after appending prediction\n',x    
+    return x_data
+
+##### Define Neural Network Models ###################
+def simple_test_rnn(T=1,D=1):
+    input_shape=(T,D)
+    input_img = Input(shape=input_shape)
+    encoder=TimeDistributed(Dense(20,activation='relu'))(input_img)
+    rnn=LSTM(10,activation='elu',return_sequences=True, stateful=False)(encoder)
+    decoder=TimeDistributed(Dense(20,activation='relu'))(rnn)
+    model=Model(input=input_img,output=decoder)
+    return model
+
+
+def dense_auto(weights_path=None,input_shape=(784,),hidden_layers=None,nonlinearity='relu'):
+    input_img = Input(shape=input_shape)
+    
+    if hidden_layers!=None:
+        if type(hidden_layers)!=list:
+            hidden_layers=list(hidden_layers)
+        for i,l in enumerate(hidden_layers):
+            if i==0: 
+                encoded=Dense(l,activation=nonlinearity)(input_img)
+            else:
+                encoded=Dense(l,activation=nonlinearity)(encoded)
+
+        for i,l in reversed(list(enumerate(hidden_layers))):
+            if i <len(hidden_layers)-1:
+                if i==len(hidden_layers)-2:
+                    decoded=Dense(l,activation=nonlinearity)(encoded)
+                else:
+                    decoded=Dense(l,activation=nonlinearity)(decoded)
+        decoded=Dense(input_shape[0])(decoded)
+    else:
+        decoded=Dense(input_shape[0])(input_img)
+
+    model=Model(input=input_img,output=decoded)
+    
+    if weights_path:
+        print('Loading Model')
+        model.load_weights(weights_path)
+    return model
+
+def rnn_dense_auto(weights_path=None,T=1,D=1,nonlinearity='relu',hidden_layers=None,recurrent_layers=None):
+    input_shape=(T,D)
+    input_img = Input(shape=input_shape)
+
+    if hidden_layers!=None:
+        if type(hidden_layers)!=list:
+            hidden_layers=list(hidden_layers)
+        for i,l in enumerate(hidden_layers):
+            if i==0: 
+                encoded=TimeDistributed(Dense(l,activation=nonlinearity))(input_img)
+            else:
+                encoded=TimeDistributed(Dense(l,activation=nonlinearity))(encoded)
+
+        for i,l in enumerate(recurrent_layers):
+            if i==0:
+                rnn=LSTM(l,return_sequences=True, stateful=False)(encoded)
+            else:
+                rnn=LSTM(l,return_sequences=True, stateful=False)(rnn)
+
+        for i,l in reversed(list(enumerate(hidden_layers))):
+            if i <len(hidden_layers)-1:
+                if i==len(hidden_layers)-2:
+                    decoded=TimeDistributed(Dense(l,activation=nonlinearity))(rnn)
+                else:
+                    decoded=TimeDistributed(Dense(l,activation=nonlinearity))(decoded)
+        decoded=TimeDistributed(Dense(D))(decoded)
+    else:
+        decoded=TimeDistributed(Dense(D))(input_img)
+    
+    model=Model(input=input_img,output=decoded)
+
+    if weights_path:
+        print('Loading Model')
+        model.load_weights(weights_path)
+    return model
+
+
+class autoencoder_preprocess():
+    def __init__(self,img_size=(784,),noise_factor=0.):
+        self.noise=noise_factor
+        self.img_size=img_size
+        self.lock = threading.Lock()
+
+    def add_noise(self,X_train):
+        ## Add noise to input data
+        np.random.seed(100)
+        ind=np.where(X_train==0)
+        rn=self.noise*np.random.rand(np.shape(ind)[1])
+        X_train[ind]=rn
+        return X_train
+    
+    def renormalize(self,X_train,mu,sigma):
+        X_train=(X_train-mu)/sigma
+        X_train = X_train.astype("float32")
+        return X_train
+
+def get_activations(model, layer, X_batch):
+    get_activations = K.function([model.layers[0].input, K.learning_phase()], model.layers[layer].output)
+    activations = get_activations([X_batch,0])
+    return activations
+
+
