#import theano
from theano import function, config, sandbox, shared 
import theano.tensor as T
import numpy as np
import scipy
import time

"""
Explanation: I accomplished the above by running this command at the command prompt:   
THEANO_FLAGS='mode=FAST_RUN,device=gpu,floatX=float32' jupyter notebook
End of explanation
"""


print( theano.config.device )
print( theano.config.lib.cnmem)  # cf. http://deeplearning.net/software/theano/library/config.html
print( theano.config.print_active_device)# Print active device at when the GPU device is initialized.

import os, sys
os.getcwd()
os.listdir( os.getcwd() ) 

%run gpu_test.py THEANO_FLAGS='mode=FAST_RUN,device=gpu,floatX=float32,lib.cnmem=0.85' # note lib.cnmem option for CnMem

"""
Explanation: More theano setup in jupyter notebook boilerplate
End of explanation
"""


# Load the diabetes dataset
diabetes = sklearn.datasets.load_diabetes()

diabetes_X = diabetes.data
diabetes_Y = diabetes.target

#diabetes_X1 = diabetes_X[:,np.newaxis,2]
diabetes_X1 = diabetes_X[:,np.newaxis, 2].astype(theano.config.floatX)
#diabetes_Y  = diabetes_Y.reshape( diabetes_Y.shape[0], 1)
diabetes_Y = diabetes_Y.astype(theano.config.floatX)

"""
Explanation: sample data boilerplate
End of explanation
"""


m_lin = diabetes_X1.shape[0]

"""
Explanation: Linear regression
cf. Linear Regression In Theano
1_linear_regression.py from github Newmu/Theano-Tutorials
Train on $m$ number of input data points
End of explanation
"""


#x1 = T.vector('x1')  # X1, input data, with only 1 feature, i.e. X \in \mathbb{R}^N, d=1 
#ylin = T.vector('ylin') # target variable for linear regression, so that Y \in \mathbb{R}

x1 = T.scalar('x1')  # X1, input data, with only 1 feature, i.e. X \in \mathbb{R}^N, d=1 
ylin = T.scalar('ylin') # target variable for linear regression, so that Y \in \mathbb{R}

"""
Explanation: input, output variables $x$, $y$ for Theano
End of explanation
"""


thet0_init_val = np.random.randn()
thet1_init_val = np.random.randn()

thet0 = theano.shared( value=thet0_init_val, name='thet0', borrow=True)  # \theta^0
thet1 = theano.shared( thet1_init_val, name='thet1', borrow=True)   # \theta^1


"""
Explanation: Parameters (for a linear slope)
$$  
(\theta^0, \theta^1) \in \mathbb{R}^2 
$$
End of explanation
"""


#h_thet = T.dot( thet1, x1) + thet0
# whereas, Newmu uses
h_thet = thet1 * x1 + thet0

"""
Explanation: hypothesis function $h_{\theta}$
$$ 
h_{\theta}(x) = \theta_1 x + \theta_0
$$
End of explanation
"""


# roshansanthosh uses 
#Jthet = T.sum( T.pow(h_thet-ylin,2))/(2*m_lin)

# whereas, Newmu uses
# Jthet = T.mean( T.sqr( thet_1*x1 + thet_0 - ylin ))

Jthet = T.mean( T.pow( h_thet-ylin,2))/2
#Jthet = sandbox.cuda.basic_ops.gpu_from_host( T.mean( 
#        sandbox.cuda.basic_ops.gpu_from_host( T.pow( h_thet-ylin,2))))/2

"""
Explanation: Cost function $J(\theta)$
End of explanation
"""


grad_thet0 = T.grad(Jthet, thet0)
grad_thet1 = T.grad(Jthet, thet1)


# so-called "learning rate"
gamma = 0.01

"""
Explanation: $$
\text{grad}{\theta}J(\theta) = ( \text{grad}{\theta^0} J , \text{grad}_{\theta^1} J ) 
$$
End of explanation
"""


train_lin = theano.function(inputs = [x1,ylin], outputs=Jthet, 
                        updates=[[thet1,thet1-gamma*grad_thet1],[thet0,thet0-gamma*grad_thet0]])



test_lin = theano.function([x1],h_thet)

#X1_lin_in = shared( diabetes_X1 ,'float32')
#Y_lin_out = shared( diabetes_Y, 'float32')

training_steps = 1000 # 10000

sh_diabetes_X1 = shared( diabetes_X1 , borrow=True)
sh_diabetes_Y  = shared( diabetes_Y, borrow=True)

"""
for i in range(training_steps):
    for x,y in zip( diabetes_X1, diabetes_Y):
        Jthet_val = train_lin( x, y )
        """

for i in range(training_steps):
#    for x,y in zip( sh_diabetes_X1, sh_diabetes_Y) :
#        Jthet_val = train_lin( x,y)
    Jthet_val = train_lin( sh_diabetes_X1, sh_diabetes_Y)

print(Jthet_val)

print( thet0.get_value() ); print( thet1.get_value() )



test_lin_out = np.array( [ test_lin( x ) for x in diabetes_X1 ] ) 

plt.plot(diabetes_X1,diabetes_Y,'ro')
plt.plot(diabetes_X1,test_lin_out)

if any([x.op.__class__.__name__ in ['GpuGemm','GpuGemv'] for x in train_lin.maker.fgraph.toposort()]):
    print("Used the gpu")
else:
    print(train_lin.maker.fgraph.toposort())

if np.any([isinstance(x.op,T.Elemwise) for x in train_lin.maker.fgraph.toposort()]):
    print("Used the cpu")

"""
Explanation: Note that "updates (iterable over pairs (shared_variable, new_expression) List, tuple or dict.) – expressions for new SharedVariable values" cf. Theano doc
End of explanation
"""


A = T.matrix('A')
B = T.matrix('B')
#matadd = function([A,B], A+B)
#matadd = function([A,B],sandbox.cuda.basic_ops.gpu_from_host(A+B) )
# Note: we are just defining the expressions, nothing is evaluated here!  
C = sandbox.cuda.basic_ops.gpu_from_host(A+B)
matadd = function([A,B], C)

#A = T.dmatrix('A')
#B = T.dmatrix('B')

A = T.matrix('A')
B = T.matrix('B')

C_out = A + B
matadd_CPU = function([A,B], C_out)

A_eg = shared( np.array([[8,6,9],[10,1,10]]), 'float32')
B_eg = shared( np.array([[3,10,2],[6,1,-1]]), 'float32')


A_eg_CPU = np.array([[8,6,9],[10,1,10]])
B_eg_CPU = np.array([[3,10,2],[6,1,-1]])

print(A_eg_CPU)
print( type( A_eg_CPU ))
print( A_eg_CPU.shape)
print( B_eg_CPU.shape)

print( matadd.maker.fgraph.toposort() )

print( matadd_CPU.maker.fgraph.toposort() )

matadd( A_eg, B_eg)

"""
Explanation: Linear Algebra and theano
cf. Week 1, Linear Algebra Review, Coursera, Machine Learning with Ng
I'll take this opportunity to provide a dictionary between the syntax of linear algebra math and numpy.   
Essentially, what I did was take Coursera's Week 1, Linear Algebra Review and then translated the math into theano, and in particular, running theano on the GPU.
Other reference that I used was 
https://simplyml.com/linear-algebra-shootout-numpy-vs-theano-vs-tensorflow-2/
Linear Algebra Shootout: NumPy vs. Theano vs. TensorFlow by Charanpal Dhanjal - 14/07/16 
Matrix addition
cf. Coursera, Intro. to Machine Learning, Linear Algebra Review, Addition and Scalar Multiplication
End of explanation
"""


np.random.randn( *A_eg_CPU.shape )

C_out = theano.shared( np.random.randn( *A_eg_CPU.shape).astype('float32') )

C_out.type()

#A_in = shared( A_eg_CPU, "float32")
#A_in = shared( A_eg_CPU, "float32")

A_in = shared( A_eg_CPU.astype("float32"), "float32")
B_in = shared( B_eg_CPU.astype("float32"), "float32")
#C_out_GPU = A_in + B_in
C_out_GPU = sandbox.cuda.basic_ops.gpu_from_host(A_in+B_in)


matadd_GPU = theano.function( [], C_out_GPU)

C_out_GPU_result = matadd_GPU()

C_out_GPU_result

"""
Explanation: The way to do it, to "force" on the GPU, is like this (cf. Speeding up your Neural Network with Theano and the GPU - Wild ML):
End of explanation
"""


print( matadd_GPU.maker.fgraph.toposort() )

#if np.any([isinstance(C_out_GPU.op, tensor.Elemwise ) and 
if np.any([isinstance( C_out_GPU.op, T.Elemwise ) and 
           ('Gpu' not in type( C_out_GPU.op).__name__) for x in matadd_GPU.maker.fgraph.toposort()]) :
    print('Used the cpu')
else:
    print('Used the gpu')

matadd_CPU( A_eg_CPU.astype("float32"), B_eg_CPU.astype("float32") )

type(A_eg)

print( type( numpy.asarray(rng.rand(2000)) ) )
numpy.asarray(rng.rand(2000)).shape

"""
Explanation: Notice how DIFFERENT this setup or syntax is: we have to set up tensor or matrix shared variables A_n, B_in, which are then used to define the theano function, theano.function.  "By using shared variables we ensure that they are present in the GPU memory".  cf. Linear Algebra Shootout: NumPy vs. Theano vs. TensorFlow
End of explanation
"""


A = T.matrix('A')
B = T.matrix('B')

C = sandbox.cuda.basic_ops.gpu_from_host( A + B ) # vs. 
# C = A + B  # this will result in an output array on the host, as opposed to CudaNdarray on device
matadd = function([A,B], C)

print( matadd.maker.fgraph.toposort() )

matadd( A_eg_CPU.astype("float32"), B_eg_CPU.astype("float32") )

"""
Explanation: Bottom Line: there are 2 ways of doing linear algebra on the GPU

symbolic computation with the usual arguments    

$$
A + B = C \in \text{Mat}_{\mathbb{R}}(M,N) 
$$   
$ \forall \, A, B \in \text{Mat}_{\mathbb{R}}(M,N)$
End of explanation
"""


A_in = shared( A_eg_CPU.astype("float32"), "float32")  # initialize with the input values, A_eg_CPU, anyway
B_in = shared( B_eg_CPU.astype("float32"), "float32")  # initialize with the input values B_eg_CPU, anyway

# C_out = A_in + B_in # this version will output to the host as a numpy.ndarray
# indeed, reading the graph,
"""
[GpuElemwise{add,no_inplace}(float32, float32), HostFromGpu(GpuElemwise{add,no_inplace}.0)]
"""
# this version immediately below, in 1 line, will result in a CudaNdarray on device
C_out = sandbox.cuda.basic_ops.gpu_from_host(A_in+B_in)

matadd_GPU = theano.function( [], C_out)

print( matadd_GPU.maker.fgraph.toposort() )

C_out_result = matadd_GPU()

C_out_result

"""
Explanation: with shared variables
End of explanation
"""


A_2 = np.array( [[4,5],[1,7] ])

a = T.scalar('a')

F = sandbox.cuda.basic_ops.gpu_from_host( a*A )
scalarmul = theano.function([a,A],F)

print( scalarmul.maker.fgraph.toposort() )

scalarmul( np.float32( 2.), A_2.astype("float32"))

"""
Explanation: Scalar Multiplication (on the GPU)
cf. Scalar Multiplication of Linear Algebra Review, coursera, Machine Learning Intro by Ng
End of explanation
"""


scalarmul( np.float32(2.), matadd( A_eg_CPU.astype("float32"), B_eg_CPU.astype("float32") ) )

u = T.vector('u')
v = T.vector('v')

w = sandbox.cuda.basic_ops.gpu_from_host( u + v)
vecadd = theano.function( [u,v],w)

t = sandbox.cuda.basic_ops.gpu_from_host( a * u)
scalarmul_vec = theano.function([a,u], t)


print(vecadd.maker.fgraph.toposort())  
print(scalarmul_vec.maker.fgraph.toposort())  


u_eg = np.array( [4,6,7], dtype="float32")
v_eg = np.array( [2,1,0], dtype="float32")

print( u_eg.shape)

scalarmul_vec( np.float32(0.5), u_eg )

vecadd( scalarmul_vec( np.float32(0.5), u_eg ) , scalarmul_vec( np.float32(-3.), v_eg )  )

"""
Explanation: Composition; Confirming that you can do composition of scalar multiplication on a matrix (or ring) addition
Being able to do composition is very important in math
End of explanation
"""


B_out = sandbox.cuda.basic_ops.gpu_from_host( T.dot(A,v))
AVmul = theano.function([A,v], B_out)
print(AVmul.maker.fgraph.toposort())

AVmul( np.array([[1,0,3],[2,1,5],[3,1,2]]).astype("float32"), np.array([1,6,2]).astype("float32"))

AVmul( np.array([[1,0,0],[0,1,0],[0,0,1]]).astype("float32"), np.array([1,6,2]).astype("float32"))

"""
Explanation: This was the computer equivalent to mathematical expression:  
$$ 
\left[ \begin{matrix} 4 \ 6 \ 7 \end{matrix} \right] /2 - 3 * \left[ \begin{matrix} 2 \ 1 \ 0 \end{matrix} \right]
$$   
sAxy or A-V multiplication or so-called "Gemv", or Matrix Multiplication on a vector, or linear transformation on a R-module, or vector space
i.e.   
$$
Av = B 
$$
End of explanation
"""


C_f = sandbox.cuda.basic_ops.gpu_from_host( T.dot(A,B)) 
matmul = theano.function([A,B], C_f)
print( matmul.maker.fgraph.toposort())

matmul( np.array( [[1,3],[2,4],[0,5]]  ).astype("float32"), np.array([[1,0],[2,3]]).astype("float32")  )

"""
Explanation: AB or Gemm or Matrix Multiplication, i.e. Ring multiplication
i.e. 
$$
A*B = C
$$
End of explanation
"""


Ainverse = sandbox.cuda.basic_ops.gpu_from_host( T.inv(A))
Ainv = theano.function([A], Ainverse)
print(Ainv.maker.fgraph.toposort())

Atranspose = sandbox.cuda.basic_ops.gpu_from_host( A.T)
AT = theano.function([A],Atranspose)
print(AT.maker.fgraph.toposort())

"""
Explanation: Inverse and Transpose
cf. Inverse and Transpose
End of explanation
"""


linregdata = pd.read_csv('./coursera_Ng/machine-learning-ex1/ex1/ex1data1.txt', header=None)

X_linreg_training = linregdata.as_matrix([0])  # pandas.DataFrame.as_matrix convert frame to its numpy-array representation
y_linreg_training = linregdata.as_matrix([1])
m_linreg_training = len(y_linreg_training)  # number of training examples  
print( X_linreg_training.shape, type(X_linreg_training)) 
print( y_linreg_training.shape, type(y_linreg_training)) 
print m_linreg_training

"""
Explanation: Summation, sum, mean, scan
Linear Regression (again), via Coursera's Machine Learning Intro by Ng, Programming Exercise 1 for Week 2
Boilerplate, load sample data
End of explanation
"""


# theta_linreg = T.vector('theta_linreg')
d = X_linreg_training.shape[1] # d = features

# Declare Theano symbolic variables
X = T.matrix('x')
y = T.vector('y')

"""
Explanation: Try representing $\theta$, parameters or "weights", of size $|\theta|$ which should be equal to the number of features $n$ (or $d$).
End of explanation
"""


#X_linreg_training = X_linreg_training.reshape( m_linreg_training,1)
#y_linreg_training = y_linreg_training.reshape( m_linreg_training,1)

# Instead, the training data X and test data values y are going to be represented by Theano symbolic variable above
#X_linreg = theano.shared(X_linreg_training.astype("float32"),"float32")
#y_linreg = theano.shared(y_linreg_training.astype("float32"),"float32")

#theta_0 = np.zeros( ( d+1,1)); print(theta_0)
theta_0 = np.zeros(  d+1); print(theta_0)

theta = theano.shared( theta_0.astype("float32"), "theta")

alpha = np.float32(0.01) # learning rate gamma or alpha

# Construct Theano "expression graph"

predicted_vals = sandbox.cuda.basic_ops.gpu_from_host( T.dot(X,theta) )  # h_{\theta}
m = np.float32( y_linreg_training.shape[0] ) 
J_theta = sandbox.cuda.basic_ops.gpu_from_host( 
    T.dot( (T.dot(X,theta) - y).T, T.dot(X,theta) - y)  * np.float32( 0.5 ) * np.float32( 1./ m )   
    ) # cost function



update_theta = sandbox.cuda.basic_ops.gpu_from_host( 
        theta - alpha * T.grad( J_theta, theta) )


gradientDescent = theano.function( 
                            inputs=[X,y],
                            outputs=[predicted_vals,J_theta],  
                            updates=[(theta, update_theta)], 
                            name = "gradientDescent")


print( gradientDescent.maker.fgraph.toposort() )

num_iters = 1500
J_History = []

"""
Explanation: Preprocess training data (due to numpy's treatment of arrays) (note, this is not needed, if you use pandas to choose which column(s) you want to make into a numpy array)
End of explanation
"""


input_X_linreg = np.hstack( ( np.ones((m_linreg_training,1)), X_linreg_training ) ).astype("float32")

y_linreg_training_processed = y_linreg_training.reshape( m_linreg_training,).astype("float32")

J_History = [0 for iter in range(num_iters)]
for iter in range(num_iters):
    predicted_vals_out, J_out = \
        gradientDescent(input_X_linreg.astype("float32"), y_linreg_training_processed.astype("float32") ) 
    J_History[iter] = J_out

Deg = (np.random.randn(40,10).astype("float32"), np.random.randint(size=40,low=0,high=2).astype("float32") )

Deg[0].shape

Deg[1].shape

theta.get_value()

dir( J_History[0] )

J_History[-5].gpudata

plt.plot( [ele.gpudata for ele in J_History])

"""
Explanation: Preprocess X to include intercepts
End of explanation
"""


input_X_linreg.shape

# GPU NOTE: Conversion to float32 to store them on the GPU!
X = theano.shared( input_X_linreg.astype('float32'), name='X' )
y = theano.shared( y_linreg_training.astype('float32'), name='y')

# GPU NOTE: Conversion to float32 to store them on the GPU!  
theta = theano.shared( np.vstack(theta_0).astype("float32"), name='theta')

# Construct Theano "expression graph"

predicted_vals = sandbox.cuda.basic_ops.gpu_from_host( 
    T.dot(X,theta) )  # h_{\theta}
m = np.float32( y_linreg_training.shape[0] )
# cost function J_theta, J_{\theta}
J_theta = sandbox.cuda.basic_ops.gpu_from_host( 
    (
        T.dot( (T.dot(X,theta) - y).T, T.dot(X,theta) - y) * np.float32(0.5) * np.float32( 1./m)  
    ).reshape([]) )   # cost function  # reshape is to force "broadcast" into 0-dim. scalar for cost function
    


update_theta = sandbox.cuda.basic_ops.gpu_from_host( 
        theta - alpha * T.grad( J_theta, theta) )

# Note that we removed the input values because we will always use the same shared variable
# GPU Note: Removed the input values to avoid copying data to the GPU.
gradientDescent = theano.function( 
                            inputs=[],
#                            outputs=[predicted_vals,J_theta],  
                            updates=[(theta, update_theta)], 
                            name = "gradientDescent")


print( gradientDescent.maker.fgraph.toposort() )

#J_History = [0 for iter in range(num_iters)]
for iter in range(num_iters):
    gradientDescent( ) 
    

print( np.vstack( theta_0).shape )
print( y_linreg_training.shape )

theta.get_value()

# Profiling
print( theano.config.profile )  # Do the vm/cvm linkers profile the execution time of Theano functions?
print( theano.config.profile_memory ) # Do the vm/cvm linkers profile the memory usage of Theano functions? It only works when profile=True.

theano.printing.debugprint(gradientDescent)

#print( gradientDescent.profile.print_summary() )
dir( gradientDescent.profile)

"""
Explanation: Denny Britz's way:  
http://www.wildml.com/2015/09/speeding-up-your-neural-network-with-theano-and-the-gpu/
Speeding up your Neural Network with Theano and the GPU
and his jupyter notebook
https://github.com/dennybritz/nn-theano/blob/master/nn-theano-gpu.ipynb
 nn-theano/nn-theano-gpu.ipynb
End of explanation
"""


import sys
import os

#sys.path.append( os.getcwd() + '/ML')
sys.path.append( os.getcwd() + '/ML' )

from linreg_gradDes import LinearReg, LinearReg_loaded
#from ML import LinearReg, LinearReg_loaded

"""
Explanation: Testing the Linear Regression with (Batch) Gradient Descent classes in ./ML/
End of explanation
"""


linregdata1 = pd.read_csv('./coursera_Ng/machine-learning-ex1/ex1/ex1data1.txt', header=None)
linregdata1.as_matrix([0]).shape
linregdata1.as_matrix([1]).shape

features = linregdata1.as_matrix([0]).shape[1]
numberoftraining = linregdata1.as_matrix([0]).shape[0]
LinReg_housing = LinearReg( features, numberoftraining , 0.01)

Xin   = LinReg_housing.preprocess_X( linregdata1.as_matrix([0]))
ytest = linregdata1.as_matrix([1]).flatten()

%time LinReg_housing.build_model( Xin, ytest )

LinRegloaded_housing = LinearReg_loaded( linregdata1.as_matrix([0]), linregdata1.as_matrix([1]), 
                                        features, numberoftraining )

%time LinRegloaded_housing.build_model()

print( LinReg_housing.gradientDescent.maker.fgraph.toposort() )
print( LinRegloaded_housing.gradientDescent.maker.fgraph.toposort() )


"""
Explanation: Boilerplate for sample input data
End of explanation
"""


def featureNormalize(X):
    """
    FEATURENORMALIZE Normalizes the features in X  
    FEATURENORMALIZE(X) returns a normalized version of X where  
    the mean value of each feature is 0 and the standard deviation  
    is 1.  This is often a good preprocessing step to do when 
    working with learning algorithms.
    
    """
    # You need to set these values correctly  
    X_norm = (X-X.mean(axis=0))/X.std(axis=0)
    mu = X.mean(axis=0)
    sigma = X.std(axis=0)
    
    return [X_norm, mu, sigma]

linregdata2 = pd.read_csv('./coursera_Ng/machine-learning-ex1/ex1/ex1data2.txt', header=None)


features = linregdata2.as_matrix().shape[1] - 1
numberoftraining = linregdata2.as_matrix().shape[0]
Xdat  = linregdata2.as_matrix( range(features) )
ytest = linregdata2.as_matrix( [features])

[Xnorm, mus,sigmas] = featureNormalize(Xdat)

LinReg_housing2 = LinearReg( features, numberoftraining, 0.01)
processed_X = LinReg_housing2.preprocess_X( Xnorm )

%time LinReg_housing2.build_model( processed_X, ytest.flatten(), 400)

LinRegloaded_housing2 = LinearReg_loaded( Xnorm, ytest, 
                                        features, numberoftraining )

%time LinRegloaded_housing2.build_model(  400)

"""
Explanation: Other (sample) datasets
Consider feature normalization
End of explanation
"""


# Load the diabetes dataset
diabetes = sklearn.datasets.load_diabetes()

diabetes_X = diabetes.data
diabetes_Y = diabetes.target

#diabetes_X1 = diabetes_X[:,np.newaxis,2]
diabetes_X1 = diabetes_X[:,np.newaxis, 2].astype(theano.config.floatX)
#diabetes_Y  = diabetes_Y.reshape( diabetes_Y.shape[0], 1)
diabetes_Y = np.vstack( diabetes_Y.astype(theano.config.floatX) )

features1 = 1 
numberoftraining = diabetes_Y.shape[0]

LinReg_diabetes = LinearReg( features1, numberoftraining, 0.01)


processed_X = LinReg_diabetes.preprocess_X( diabetes_X1 )

%time LinReg_diabetes.build_model( processed_X, diabetes_Y.flatten(), 10000)

LinRegloaded_diabetes = LinearReg_loaded( diabetes_X1, diabetes_Y, 
                                        features1, numberoftraining )

%time LinRegloaded_diabetes.build_model(  10000)

"""
Explanation: Diabetes data from sklearn, sci-kit learn
End of explanation
"""


features = diabetes_X.shape[1]


LinReg_diabetes = LinearReg( features, numberoftraining, 0.01)
processed_X = LinReg_diabetes.preprocess_X( diabetes_X )

%time LinReg_diabetes.build_model( processed_X, diabetes_Y.flatten(), 10000)

LinRegloaded_diabetes = LinearReg_loaded( diabetes_X, diabetes_Y, 
                                        features, numberoftraining )

%time LinRegloaded_diabetes.build_model(  10000)

"""
Explanation: Multiple number of features case:
End of explanation
"""


data_ex1data2 = pd.read_csv('./coursera_Ng/machine-learning-ex1/ex1/ex1data2.txt', header=None)
X_ex1data2 = data_ex1data2.iloc[:,0:2]
y_ex1data2 = data_ex1data2.iloc[:,2]
m_ex1data2 = y_ex1data2.shape[0]
X_ex1data2=X_ex1data2.values.astype(np.float32)
y_ex1data2=y_ex1data2.values.reshape((m_ex1data2,1)).astype(np.float32)
print(type(X_ex1data2))
print(type(y_ex1data2))
print(X_ex1data2.shape)
print(y_ex1data2.shape)
print(m_ex1data2)
print(X_ex1data2[:5])
print(y_ex1data2[:5])

((X_ex1data2[:,1] - X_ex1data2[:,1].mean())/( X_ex1data2[:,1].std()) ).std()

# feature Normalize
#X_ex1data2_norm = sklearn.preprocessing.Normalizer.transform(X_ex1data2 )
X_ex1data2_norm = (X_ex1data2  - np.mean(X_ex1data2, axis=0)) / np.std(X_ex1data2, axis=0)
print(X_ex1data2_norm[:,0].mean())
print(X_ex1data2_norm[:,0].std())
print(X_ex1data2_norm[:,1].mean())
print(X_ex1data2_norm[:,1].std())

# X_ex1data2_norm[:5];

X=T.matrix(dtype=theano.config.floatX)
y=T.matrix(dtype=theano.config.floatX)

Theta=theano.shared(np.zeros((2,1)).astype(theano.config.floatX))
b = theano.shared(np.zeros(1).astype(theano.config.floatX))

print(b.get_value().shape)

yhat = T.dot( X, Theta) + b

# L2 norm
J = np.cast[theano.config.floatX](0.5)*T.mean( T.sqr( yhat-y))

alpha=0.01  # learning rate
# sandbox.cuda.basic_ops.gpu_from_host
updateThetab = [ Theta-np.float32(alpha)*T.grad(J,Theta), b-np.float32(alpha)*T.grad(J,b)]
gradientDescent_step = theano.function(inputs=[X,y], 
                                          outputs=J,
                                          updates = zip([Theta,b],updateThetab) )


num_iters =400
JList=[]
for iter in range(num_iters):
    err = gradientDescent_step(X_ex1data2_norm,y_ex1data2)
    JList.append(err)

# Final mode:
print(Theta.get_value())
print(b.get_value())

# JList[-10:]
plt.plot(JList)
plt.show()

"""
Explanation: ex2 Linear Regression, on d=2 features
End of explanation
"""


os.getcwd()

os.listdir( './coursera_Ng/machine-learning-ex3/' )

os.listdir( './coursera_Ng/machine-learning-ex3/ex3' )

# Load saved matrices from file 
multiclscls_data = scipy.io.loadmat('./coursera_Ng/machine-learning-ex3/ex3/ex3data1.mat')

"""
Explanation: Multi-class Classification
cf. ex3, Programming Exercise 3: Multi-class Classification and Neural Networks, Machine Learning
1 Multi-class Classification
End of explanation
"""


import sys
import os

os.getcwd()

#sys.path.append( os.getcwd() + '/ML')
sys.path.append( os.getcwd() + '/ML' )

from gradDes import LogReg

# Test case for Cost function J_{\theta} with regularization

theta_t = np.vstack( np.array( [-2, -1, 1, 2]) )
X_t = np.array( [i/10. for i in range(1,16)]).reshape((3,5)).T
#X_t = np.hstack( ( np.ones((5,1)), X_t) ) # no need to preprocess the input data X with column of 1's
y_t = np.vstack( np.array( [1,0,1,0,1]))


MulClsCls_digits = LogReg( X_t, y_t, 3,5,0.01, 3.  )

MulClsCls_digits.calculate_cost()

MulClsCls_digits.z.get_value()

print( MulClsCls_digits.X.get_value() )
MulClsCls_digits.y.get_value()

calc_z_test = theano.function([], MulClsCls_digits.z)

calc_z_test()

MulClsCls_digits.theta.set_value( theta_t.astype('float32') )

calc_z_test()

MulClsCls_digits.calculate_cost()

print( 1/(1+np.exp( np.dot( -np.hstack( ( np.ones((5,1)), X_t) ), theta_t) ) )  )
h_test = 1/(1+np.exp( np.dot( -np.hstack( ( np.ones((5,1)), X_t) ), theta_t) ) ) 
print( np.dot( (h_test - y_t).T, h_test- y_t) * 0.5/5 ) # non-regularized J_theta cost term
np.dot( theta_t[1:].T, theta_t[1:]) * 3 / (2.* 5)


MulClsCls_digits.predict()

MulClsCls_digit

theano.config.floatX

"""
Explanation: import the classes from ML
End of explanation
"""


os.getcwd()

os.listdir( './coursera_Ng/machine-learning-ex3/' )

os.listdir( './coursera_Ng/machine-learning-ex3/ex3/' )

"""
Explanation: Neural Networks
Model representation
cf. 2 Neural Networks, 2.1 Model representation, ex3.pdf
End of explanation
"""


# Load saved matrices from file 
nn3_data = scipy.io.loadmat('./coursera_Ng/machine-learning-ex3/ex3/ex3weights.mat')

print( nn3_data.keys() )
print( type( nn3_data['Theta1']) )
print( type( nn3_data['Theta2']) )
print( nn3_data['Theta1'].shape )
print( nn3_data['Theta2'].shape )

Theta1[0]

"""
Explanation: $ \Theta_1, \Theta_2 $
End of explanation
"""


%load_ext tikzmagic

"""
Explanation: Feedforward
End of explanation
"""


np.random.seed(0)
s_l = 400 # (layer) size of layer l, i.e. number of nodes, units in layer l
s_lp1 = 25
al = theano.shared( np.random.randn(s_l+1,1).astype('float32'), name="al")
#alp1 = theano.shared( np.random.randn(s_lp1,1).astype('float32'), name="al")
#Thetal = theano.shared( np.random.randn( s_lp1,s_l+1).astype('float32') , name="Thetal")

# Feedforward, forward propagation
#z = T.dot( Thetal, al)
#g = T.nnet.sigmoid( z)


s_l = 25
s_lp1 = 10




rng = np.random.RandomState(99)
Theta_values = np.asarray( rng.uniform( 
    low=-np.sqrt( 6. / (s_l+ s_lp1)), 
    high=np.sqrt( 6./(s_l + s_lp1)), size=(s_lp1,s_l+1)), dtype=theano.config.floatX )
print( Theta_values.shape )
print( Theta_values.dtype )
#Theta_values *= np.float32(4)
Theta_values *= 4.

print( Theta_values.dtype)
Theta_values.shape

np.float32( 4)

"""
Explanation: $$
\begin{tikzpicture}
  \matrix (m) [matrix of math nodes, row sep=3em, column sep=4em, minimum width=2em]
  {
    \mathbb{R}^{s_l}  &  \mathbb{R}^{ s_l +1 }   & \mathbb{R}^{s_{l+1} } & \mathbb{R}^{s_{l+1} }  \
a^{(l)} & (a_0^{(l)} = 1, a^{(l)} ) & z^{(l+1)} & g(z^{(l+1)}) = a^{(l+1)} \
  };
  \path[->]
  (m-1-1) edge node [above] {$a_0^{(l)}=1$} (m-1-2)
  (m-1-2) edge node [above] {$\Theta^{(l)}$} (m-1-3)
  (m-1-3) edge node [above] {$g$} (m-1-4) 
  ;
  \path[|->]
  (m-2-1) edge node [above] {$a_0^{(l)}=1$} (m-2-2)
  (m-2-2) edge node [above] {$\Theta^{(l)}$} (m-2-3)
  (m-2-3) edge node [above] {$g$} (m-2-4) 
  ;
\end{tikzpicture}
$$
End of explanation
"""


%env

os.getcwd()

print( sys.path )

#sys.path.append( os.getcwd() + '/ML')
sys.path.append( '../DeepLearningTutorials/code/' )

#from logistic_sgd import LogisticRegression, load_data, sgd_optimization_mnist, predict
import logistic_sgd 

MNIST_MTLdat = logistic_sgd.load_data("../DeepLearningTutorials/data/mnist.pkl.gz") # list of training data

print(len(MNIST_MTLdat))
print(type(MNIST_MTLdat))
for ele in MNIST_MTLdat: print type(ele), len(ele) # test_set_x, test_set_y, valid_set_x, valid_set_y, train_set_x, 

print( MNIST_MTLdat[0][0].get_value().shape)
print( type(MNIST_MTLdat[0][1]))
print( MNIST_MTLdat[0][1].get_scalar_constant_value )

print( type( MNIST_MTLdat[1][1] ) )
MNIST_MTLdat[1][1].shape

dir(MNIST_MTLdat[0][1]) ;

import gzip
import six.moves.cPickle as pickle
with gzip.open("../DeepLearningTutorials/data/mnist.pkl.gz", 'rb') as f:
    try:
        train_set, valid_set, test_set = pickle.load(f, encoding='latin1')
    except:
        train_set, valid_set, test_set = pickle.load(f)

print( type( train_set[0] ))
print( train_set[0].shape )
print( type( train_set[1]))
print( train_set[1].shape )
print( type( valid_set[0] ))
print( valid_set[0].shape )
print( type( valid_set[1]))
print( valid_set[1].shape )
print( type( test_set[0] ))
print( test_set[0].shape )
print( type( test_set[1]))
print( test_set[1].shape )


X = train_set[0].T

pd.DataFrame(X.T).describe()

28*28

X_i = theano.shared( X.astype("float32"))

m = X_i.get_value().shape[1]

a1 = T.stack( [ theano.shared( np.ones((1,m)).astype("float32") ) , X_i ] , axis=1 )

print( type(a1) )
#print( a1.get_scalar_constant_value() )
dir(a1)
a1.get_parents()

a1.ndim

a1_0 = theano.shared( np.ones((1,m)).astype("float32"),name='a1_0')


a1 = T.stack( [a1_0,X_i], axis=0)

d = X_i.get_value().shape[0]
s_2 = d/2
rng1 = np.random.RandomState(1234)
Theta1_values = np.asarray( rng1.uniform( low=-np.sqrt(6./(d+s_2)),high=np.sqrt(6./(d+s_2)),size=(s_2,d+1)),
                           dtype=theano.config.floatX)
Theta1 = theano.shared(value=Theta1_values, name="Theta",borrow=True)

#rng1.uniform( low=-np.sqrt(6./(d+s_2)),high=np.sqrt(6./(d+s_2)),size=(s_2,d+1))
z1 = T.dot( Theta1, a1)
a2 = T.tanh(z1)

passthru1 = theano.function( [], a2)

print(d)
passthru1()

print(X.shape)
X_i = theano.shared( X.astype("float32"))
#m = X_i.get_value().shape[1]
m = X.shape[1]
print(m)
a1_0 = theano.shared( np.ones((1,m)).astype("float32"),name='a1_0')
print(a1_0.get_value().shape)
a1 = T.stack( [a1_0,X_i], axis=0)
addintercept = theano.function([],a1)

addintercept()

d = X_i.get_value().shape[0]
print(d)
s_2 = d/2
print(s_2)
rng1 = np.random.RandomState(1234)
Theta1_values = np.asarray( rng1.uniform( low=-np.sqrt(6./(d+s_2)),high=np.sqrt(6./(d+s_2)),size=(s_2,d)),
                           dtype=theano.config.floatX)
Theta1 = theano.shared(value=Theta1_values, name="Theta1",borrow=True)
b_values = np.vstack( np.zeros(s_2) ).astype(theano.config.floatX)
b1 = theano.shared(value=b_values, name='b1',borrow=True)
a1_values=np.array( np.zeros( (d,m)), dtype=theano.config.floatX)
a1 = theano.shared(value=a1_values, name='a1', borrow=True)
lin_z2 = T.dot( Theta1, a1) + T.tile(b1,(1,m))
#lin_z2 = T.dot( Theta1, a1)

test_mult = theano.function([],lin_z2)

print( type(b_values))
b_values.dtype


test_mult()

print( b1.get_value().shape )
T.tile( b1, (0,m))

"""
Explanation: From Deep Learning Tutorials of LISA lab of University of Montreal; logistic_sgd.py, mlp.py
End of explanation
"""


import sys
import os

#sys.path.append( os.getcwd() + '/ML')
sys.path.append( os.getcwd() + '/ML' )

from NN import Layer, cost_functional, cost_functional_noreg, gradientDescent_step


"""
Explanation: NN.py, load NN.py for Layer class for Neural Net for Multiple Layers
End of explanation
"""


# Load Training Data
print("Loading and Visualizing Data ... \n")
ex4data1 = scipy.io.loadmat('./coursera_Ng/machine-learning-ex4/ex4/ex4data1.mat')

ex4data1.keys()

print( ex4data1['X'].shape )
print( ex4data1['y'].shape )

test_rng = np.random.RandomState(1234)
#Theta1 = Layer( test_rng, 1, 400,25, 5000)

#help(Theta1.al.set_value); # Beginning with Theano 0.3.1, set_value will work in-place on the GPU, if ... source on CPU
Theta1.al.set_value( ex4data1['X'].T.astype(theano.config.floatX))

Theta1.alp1

print( type( Theta1.alp1 ) )
Theta2 = Layer( test_rng, 2, 25,10,5000, al=Theta1.alp1 )

Theta2.alp1

predicted = theano.function([],sandbox.cuda.basic_ops.gpu_from_host( Theta2.alp1 ) )

predicted().shape

print( ex4data1['y'].shape )
pd.DataFrame( ex4data1['y']).describe()

# recall that whereas the original labels (in the variable y) were 1, 2, ..., 10, for the purpose of training a 
# neural network, we need to recode the labels as vectors containing only values 0 or 1
K=10
m = ex4data1['y'].shape[0]
y_prob = [np.zeros(K) for row in ex4data1['y']]  # list of 5000 numpy arrays of size dims. (10,)
for i in range( m):
        y_prob[i][ ex4data1['y'][i]-1] = 1
y_prob = np.array(y_prob).T.astype(theano.config.floatX)  # size dims. (K,m)
print(y_prob.shape)

print( type(y_prob) )
type( np.asarray( y_prob, dtype=theano.config.floatX) )

help( T.nlinalg.trace )

y_sh_var = theano.shared( np.asarray( y_prob,dtype=theano.config.floatX),name='y')

h_test = Theta2.alp1
J = sandbox.cuda.basic_ops.gpu_from_host(
        (-T.nlinalg.trace( T.dot( T.log( h_test ), y_sh_var.T)) - T.nlinalg.trace( 
        T.dot( T.log( np.float32(1.)-h_test),(np.float32(1.)- y_sh_var.T ) )))/np.float32(m)
    )

print(type(J))
test_cost_func = theano.function([],J)

test_cost_func()

J_test_build = sandbox.cuda.basic_ops.gpu_from_host( -T.nlinalg.trace( T.dot( T.log(h_test),y_sh_var.T) ) )
test_cost_build_func = theano.function([], J_test_build)

test_cost_build_func()

"""
Explanation: Boilerplate sample data, from Coursera's Machine Learning Introduction
End of explanation
"""


Theta_testvals = scipy.io.loadmat('./coursera_Ng/machine-learning-ex4/ex4/ex4weights.mat')

print( Theta_testvals.keys() )
print( Theta_testvals['Theta1'].shape )
print( Theta_testvals['Theta2'].shape )
Theta1_testval = Theta_testvals['Theta1'][:,1:]
b1_testval = Theta_testvals['Theta1'][:,0:1]
print( Theta1_testval.shape )
print( b1_testval.shape )
Theta2_testval = Theta_testvals['Theta2'][:,1:]
b2_testval = Theta_testvals['Theta2'][:,0:1]
print( Theta2_testval.shape )
print( b2_testval.shape )

Theta1 = Layer( test_rng, 1, 400,25, 5000, activation=T.nnet.sigmoid)


Theta1.Theta.set_value( Theta1_testval.astype("float32"))
Theta1.b.set_value( b1_testval.astype('float32') )
Theta1.al.set_value( ex4data1['X'].T.astype('float32'))

"""
Explanation: Sanity check using ex4.m, Exercise 4 or Programming Exercise 4 from Coursera's Machine Learning Introduction by Ng
End of explanation
"""


Theta2 = Layer( test_rng, 2, 25,10,5000, al=Theta1.alp1 , activation=T.nnet.sigmoid)

Theta2.Theta.set_value( Theta2_testval.astype('float32'))
Theta2.b.set_value( b2_testval.astype('float32'))

h_test = Theta2.alp1
J = sandbox.cuda.basic_ops.gpu_from_host(
    T.mean( T.sum( 
        - y_sh_var * T.log( h_test ) - ( np.float32( 1) - y_sh_var) * T.log( np.float32(1) - h_test), axis =0), axis=0)
  )
#J = sandbox.cuda.basic_ops.gpu_from_host( 
#    T.log(h_test) * y_sh_var
#    )

test_cost_func = theano.function([],J)

test_cost_func()

print(type( y_sh_var) )
print( y_sh_var.get_value().shape )
print( type( h_test ))

checklayer2 = theano.function([], sandbox.cuda.basic_ops.gpu_from_host(Theta1.alp1))

checklayer2() 

testreg = theano.function([], T.sum( Theta1.Theta * Theta1.Theta ) )

testreg()

range(1,3)

Thetas_lst = [ Theta1.Theta, Theta2.Theta ]

T.sum( [ T.sum( theta*theta) for theta in Thetas_lst] )

cost_func_test = cost_functional(3, 1, y_prob, Theta2.alp1, [Theta1.Theta, Theta2.Theta])

cost_test = theano.function([], cost_func_test)

cost_test() # (this value should be about 0.383770)

grad_test = T.grad( cost_func_test,[Theta1.Theta, Theta2.Theta])

grad_test_test = theano.function([], grad_test)

print( type(grad_test_test() ) )
print( len( grad_test_test() ))
print( type(grad_test_test()[0] ))
print( grad_test_test()[0].shape )
print( grad_test_test()[1].shape )

print( range(6))
print( list( "Ernest") )
zip( range(6), list("Ernest"))
print( type(grad_test))

print( grad_test_test.maker.fgraph.toposort() )

0.01 * grad_test

test_update = [(Theta,sandbox.cuda.basic_ops.gpu_from_host( Theta - np.float32(0.01)*T.grad(cost_func_test, Theta)+0.0001*Theta ) ) for Theta in [Theta1.Theta, Theta2.Theta] ]

test_gradDes_step = theano.function( inputs=[], updates= test_update )

test_gradDes_step()

print( Theta1.Theta.get_value() )
print( Theta2.Theta.get_value() )

test_gradDes_step()

print( Theta1.Theta.get_value() )
print( Theta2.Theta.get_value() )

gradDes_test_res = gradientDescent_step(cost_func_test, [Theta1.Theta, Theta2.Theta], 0.01, 0.00001 )

print( type(gradDes_test_res) )
gradDes_step_test = gradDes_test_res[1]

gradDes_step_test()

print( Theta1.Theta.get_value() )
print( Theta2.Theta.get_value() )

gradDes_step_test()

print( Theta1.Theta.get_value() )
print( Theta2.Theta.get_value() )

y_prob.shape

ex4data1['y'].shape

pd.DataFrame( ex4data1['y']).describe()

print( Theta2.alp1.shape )
print( Theta2.alp1.shape.ndim )
# Theta2.alp1.shape.get_scalar_constant_value()
predicted_logreg = theano.function([],Theta2.alp1)

pd.DataFrame( predicted_logreg().T ).describe()

pd.DataFrame(predicted_logreg().T).describe().iloc[1:-1,:].plot()

print( np.argmax( predicted_logreg(), axis=0).shape )
np.vstack( np.argmax( predicted_logreg(),axis=0) ).shape

pd.DataFrame( np.vstack( np.argmax(predicted_logreg(),axis=0)) + 1).describe()

res = np.float32( ( np.vstack( np.argmax( predicted_logreg(),axis=0)) + 1 ) == ex4data1['y'] )
pd.DataFrame(res).describe()

range(1,3)

predicted_logreg().shape

print(y_prob.shape); print( np.argmax( y_prob,axis=0 ).shape)

"""
Explanation: For $\Theta^{(2)}$, the key to connecting $\Theta^{(2)}$ with $\Theta^{(1)}$ is to set the argument in class Layer with al=Theta1.alp1,
End of explanation
"""


sys.path.append( os.getcwd() + '/ML' )

from NN import Layer, cost_functional, cost_functional_noreg, gradientDescent_step, MLP

# Load Training Data
print("Loading and Visualizing Data ... \n")
ex4data1 = scipy.io.loadmat('./coursera_Ng/machine-learning-ex4/ex4/ex4data1.mat')

# recall that whereas the original labels (in the variable y) were 1, 2, ..., 10, for the purpose of training a 
# neural network, we need to recode the labels as vectors containing only values 0 or 1
K=10
m = ex4data1['y'].shape[0]
y_prob = [np.zeros(K) for row in ex4data1['y']]  # list of 5000 numpy arrays of size dims. (10,)
for i in range( m):
        y_prob[i][ ex4data1['y'][i]-1] = 1
y_prob = np.array(y_prob).T.astype(theano.config.floatX)  # size dims. (K,m)

print(ex4data1['X'].T.shape)
print(y_prob.shape)

digitsMLP = MLP(3,[400,25,10], 5000, ex4data1['X'].T, y_prob, T.nnet.sigmoid, 1., 0.1, 0.0000)

digitsMLP.train_model(100000)

digitsMLP.accuracy_log_reg()

print( digitsMLP.Thetas[0].Theta.get_value() )
digitsMLP.Thetas[1].Theta.get_value()


digitsMLP.predicted_vals_logreg()

testL1a2 = theano.function([], digitsMLP.Thetas[0].alp1 )
print( testL1a2() )
testL2a2 = theano.function([], digitsMLP.Thetas[1].al )
print( testL2a2() )


[1,2,3,4,5] + [8,1,5]

print( digitsMLP.y.shape )
y_cls_test = np.vstack( np.argmax( digitsMLP.y, axis=0) )
print( y_cls_test.shape )
pd.DataFrame( y_cls_test ).describe()

pred_y_cls_test = np.vstack( np.argmax( digitsMLP.predicted_vals_logreg() , axis=0))
print( pred_y_cls_test.shape )
pd.DataFrame( pred_y_cls_test ).describe()

np.mean( pred_y_cls_test == y_cls_test )

"""
Explanation: Summary for Neural Net with Multiple Layers for logistic regression (but can be extended to linear regression)

Load boilerplate training data:
End of explanation
"""


K=10
m = len(train_set[1])
y_train_prob = [np.zeros(K) for row in train_set[1]]  # list of 5000 numpy arrays of size dims. (10,)
for i in range( m):
        y_train_prob[i][ train_set[1][i]] = 1
y_train_prob = np.array(y_train_prob).T.astype(theano.config.floatX)  # size dims. (K,m)
print( y_train_prob.shape )

print( pd.DataFrame( y_train_prob).describe() )

m,d= train_set[0].shape
MNIST_MTL = MLP(3,[d,25,10], m, train_set[0].T, y_train_prob, T.nnet.sigmoid, 1., 0.1, 0.00001)

MNIST_MTL.accuracy_log_reg()

print( MNIST_MTL.Thetas[0].Theta.get_value() )
MNIST_MTL.Thetas[1].Theta.get_value()

MNIST_MTL.predicted_vals_logreg()

MNIST_MTL.train_model(100000)

MNIST_MTL.accuracy_log_reg()

print( MNIST_MTL.Thetas[0].Theta.get_value() )
MNIST_MTL.Thetas[1].Theta.get_value()

MNIST_MTL.predicted_vals_logreg()

"""
Explanation: Testing on MNIST, from University of Montreal, Deep Learning Tutorial, data
End of explanation
"""


import cPickle

save_file = open('./saved_models/MNIST_MTL_log_reg','wb')

for Thet in MNIST_MTL.Thetas:
    cPickle.dump( Thet.Theta.get_value(borrow=True), save_file,-1) # the -1 is for HIGHEST priority
    cPickle.dump( Thet.b.get_value(borrow=True), save_file,-1)

save_file.close()

MNIST_MTL.Thetas[0].al.set_value( valid_set[0].T.astype(theano.config.floatX) )

K=10
m = len(valid_set[1])
y_valid_prob = [np.zeros(K) for row in valid_set[1]]  # list of 5000 numpy arrays of size dims. (10,)
for i in range( m):
        y_valid_prob[i][ valid_set[1][i]] = 1
y_valid_prob = np.array(y_valid_prob).T.astype(theano.config.floatX)  # size dims. (K,m)
print( y_valid_prob.shape )

MNIST_MTL.y = y_valid_prob

MNIST_MTL.predicted_vals_logreg()

theano.function([], MNIST_MTL.Thetas[0].alp1)()

Layer1 = MNIST_MTL.Thetas[0]
Layer2 = MNIST_MTL.Thetas[1]
m = valid_set[0].shape[0]
print(m)

a2 = T.nnet.sigmoid( T.dot( Layer1.Theta, Layer1.al) + T.tile( Layer1.b, (1,m)) )
a3 = T.nnet.sigmoid( T.dot( Layer2.Theta, a2) + T.tile( Layer2.b, (1,m)) )
valid_pred = theano.function([], a3)()
print( valid_pred.shape)

pd.DataFrame( valid_pred.T).describe()

np.mean( np.vstack( np.argmax( valid_pred,axis=0)) == np.vstack( valid_set[1] ) )

X_in = T.matrix()

X_in.set_value( valid_set[0].T.astype(theano.config.floatX))

a2_giv = T.nnet.sigmoid( T.dot( Layer1.Theta, X_in) + T.tile(Layer1.b, (1,m)))
a3_giv = T.nnet.sigmoid( T.dot( Layer2.Theta, a2_giv) + T.tile( Layer2.b, (1,m)) )
valid_pred_givens = theano.function([], outputs=a3_giv, givens={ X_in: valid_set[0].T.astype(theano.config.floatX)} )

print( valid_pred_givens().shape )
pd.DataFrame( valid_pred_givens().T).describe()

np.mean( np.vstack( np.argmax( valid_pred_givens(),axis=0)) == np.vstack( valid_set[1] ) )

test_pred_givens = theano.function([], outputs=a3_giv, givens={ X_in: test_set[0].T.astype(theano.config.floatX)} )

np.mean( np.vstack( np.argmax( test_pred_givens(),axis=0)) == np.vstack( test_set[1] ) )

range(1,3)

range(3)

range(1,3-1)

"""
Explanation: Save the mode; cf. Getting Started, DeepLearning 0.1 documentation, Loading and Saving Models
End of explanation
"""


gls_data = pd.read_csv( "./kaggle/glass.csv")

gls_data.describe()

gls_data.get_values().shape

X_gls = gls_data.get_values()[:,:-1]
print(X_gls.shape)
y_gls = gls_data.get_values()[:,-1]
print(y_gls.shape)
print( y_gls[:10])
X_gls_train = gls_data.get_values()[:-14,:-1]
print(X_gls_train.shape)
y_gls_train = gls_data.get_values()[:-14,-1]
print(y_gls_train.shape)

K=7
m = len(y_gls_train)
y_gls_train_prob = [np.zeros(K) for row in y_gls_train]  # list of 5000 numpy arrays of size dims. (10,)
for i in range( m):
        y_gls_train_prob[i][ y_gls_train[i]-1] = 1
y_gls_train_prob = np.array(y_gls_train_prob).T.astype(theano.config.floatX)  # size dims. (K,m)
print( y_gls_train_prob.shape )

gls_MLP = MLP( 3, [9,8,7],200, X_gls_train.T, y_gls_train_prob, T.nnet.sigmoid, 0.01,0.05,0.0001 )

gls_MLP.accuracy_log_reg()

gls_MLP.train_model(10000)

gls_MLP.accuracy_log_reg()

gls_MLP.predicted_vals_logreg()

gls_MLP.train_model(10000)
gls_MLP.accuracy_log_reg()

ga

X_gls_test = gls_data.get_values()[-14:,:-1]
print( X_gls_test.shape )
y_gls_test = gls_data.get_values()[-14:,-1]
print( y_gls_test.shape)

gls_predict_on_test = gls_MLP.predict_on( 14, X_gls_test.T )

np.mean( np.vstack( np.argmax( gls_predict_on_test(), axis=0) ) == (y_gls_test-1) )

y_gls_test

np.vstack( np.argmax( gls_predict_on_test(), axis=0))

X_sym = T.matrix()

rng = np.random.RandomState(1234)
Thetab1 = Layer( rng, 1, 4,3,2, al = X_sym, activation=T.nnet.sigmoid)


Thetab1.alp1
Thetab1.Theta.get_value().shape

Thetab2 = Layer( rng, 2, 3,2,2, al=Thetab1.alp1, activation=T.nnet.sigmoid)


Thetab2.al = Thetab1.alp1

X_sym.shape[0]

T.tile( Thetab1.b, (1, X_sym.shape[0]))

test12comp = theano.function( [], outputs=Thetab2.alp1, givens={ X_sym : X42test} )

X42test = np.array([1,2,3,4,5,6,7,8]).reshape((4,2)).astype(theano.config.floatX)

test12comp()

X43test = np.array(range(1,13)).reshape((4,3)).astype(theano.config.floatX)

X43test

test43comp = theano.function( [], outputs=Thetab2.alp1, givens={ X_sym : X43test} )

test43comp()

print( type(Thetab1.al ))

lin_zlp1 = T.dot(Thetab1.Theta, Thetab1.al)+T.tile( Thetab1.b, (1,Thetab1.al.shape[1]))
a1p1 = Thetab1.g( lin_zlp1 )

Thetab1.al = X_sym

Thetab2.al = a1p1

lin_z2p1 = T.dot(Thetab2.Theta, Thetab2.al)+T.tile( Thetab2.b, (1, Thetab2.al.shape[1]))
a2p1 = Thetab2.g( lin_z2p1 )

test_gen_conn = theano.function([], outputs=a2p1, givens={ Thetab1.al : X42test })

test_gen_conn()

test_gen_conn = theano.function([], outputs=a2p1, givens={ Thetab1.al : X43test })

test_gen_conn()

"""
Explanation: cf. Glass Classification
End of explanation
"""


test_gen_conn = theano.function([], outputs=sandbox.cuda.basic_ops.gpu_from_host(a2p1), givens={ Thetab1.al : X42test })

test_gen_conn()

test_gen_conn = theano.function([], outputs=sandbox.cuda.basic_ops.gpu_from_host(a2p1), givens={ Thetab1.al : X43test })

test_gen_conn()

"""
Explanation: GPU test
End of explanation
"""


sys.path.append( os.getcwd() + '/ML' )

from NN import MLP

# Load Training Data
print("Loading and Visualizing Data ... \n")
ex4data1 = scipy.io.loadmat('./coursera_Ng/machine-learning-ex4/ex4/ex4data1.mat')

# recall that whereas the original labels (in the variable y) were 1, 2, ..., 10, for the purpose of training a 
# neural network, we need to recode the labels as vectors containing only values 0 or 1
K=10
m = ex4data1['y'].shape[0]
y_prob = [np.zeros(K) for row in ex4data1['y']]  # list of 5000 numpy arrays of size dims. (10,)
for i in range( m):
        y_prob[i][ ex4data1['y'][i]-1] = 1
y_prob = np.array(y_prob).T.astype(theano.config.floatX)  # size dims. (K,m)

print(ex4data1['X'].T.shape)
print(y_prob.shape)

digitsMLP = MLP( 3, [400,25,10], ex4data1['X'].T, y_prob, T.nnet.sigmoid, 1.)

digitsMLP.build_update(ex4data1['X'].T, y_prob, 0.01, 0.00001)

digitsMLP.predicted_vals_logreg()

digitsMLP.accuracy_logreg( ex4data1['X'].T, y_prob)

digitsMLP.train_model(10000)

digitsMLP.accuracy_logreg( ex4data1['X'].T, y_prob)

digitsMLP.train_model(50000)

digitsMLP.accuracy_logreg( ex4data1['X'].T, y_prob)

"""
Explanation: Summary for Neural Net with Multiple Layers for logistic regression (but can be extended to linear regression)
End of explanation
"""


import gzip
import six.moves.cPickle as pickle
with gzip.open("../DeepLearningTutorials/data/mnist.pkl.gz", 'rb') as f:
    try:
        train_set, valid_set, test_set = pickle.load(f, encoding='latin1')
    except:
        train_set, valid_set, test_set = pickle.load(f)

K=10
m = len(train_set[1])
y_train_prob = [np.zeros(K) for row in train_set[1]]  # list of 5000 numpy arrays of size dims. (10,)
for i in range( m):
        y_train_prob[i][ train_set[1][i]] = 1
y_train_prob = np.array(y_train_prob).T.astype(theano.config.floatX)  # size dims. (K,m)
print( y_train_prob.shape )

MNIST_MLP = MLP( 3,[784,49,10], train_set[0].T, y_train_prob, T.nnet.sigmoid, 1.)

MNIST_MLP.build_update( train_set[0].T, y_train_prob, 0.01, 0.0001)

MNIST_MLP.accuracy_logreg( train_set[0].T, y_train_prob)

MNIST_MLP.train_model(50000)

MNIST_MLP.accuracy_logreg( train_set[0].T, y_train_prob)

%time MNIST_MLP.train_model(100000)

MNIST_MLP.accuracy_logreg( train_set[0].T,y_train_prob)

m = len(valid_set[1])
y_valid_prob = [np.zeros(K) for row in valid_set[1]]  # list of 5000 numpy arrays of size dims. (10,)
for i in range( m):
        y_valid_prob[i][ valid_set[1][i]] = 1
y_valid_prob = np.array(y_valid_prob).T.astype(theano.config.floatX)  # size dims. (K,m)
print( y_valid_prob.shape )

m = len(test_set[1])
y_test_prob = [np.zeros(K) for row in test_set[1]]  # list of 5000 numpy arrays of size dims. (10,)
for i in range( m):
        y_test_prob[i][ test_set[1][i]] = 1
y_test_prob = np.array(y_test_prob).T.astype(theano.config.floatX)  # size dims. (K,m)
print( y_test_prob.shape )

MNIST_MLP.accuracy_logreg( valid_set[0].T,y_valid_prob)

MNIST_MLP.accuracy_logreg( test_set[0].T,y_test_prob)

MNIST_d = train_set[0].T.shape[0]
print(MNIST_d)
MNIST_MLP = MLP( 3,[MNIST_d,25,10], train_set[0].T, y_train_prob, T.nnet.sigmoid, 1.)
MNIST_MLP.build_update( train_set[0].T, y_train_prob, 0.1, 0.00001)

MNIST_MLP.accuracy_logreg( train_set[0].T, y_train_prob)

MNIST_MLP.train_model(150000)

MNIST_MLP.accuracy_logreg( train_set[0].T, y_train_prob)

MNIST_MLP.accuracy_logreg( valid_set[0].T, y_valid_prob)

MNIST_MLP.accuracy_logreg( test_set[0].T, y_test_prob)

"""
Explanation: Testing on University of Montreal LISA lab MNIST data
End of explanation
"""


