"""
Helper functions for FNN with ESP
=================================================================
Author: <NAME> <<EMAIL>>

This is a pure python implementation of a Feed Forward Deep network with ESP
activation. See main text for details.

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
#Math Libraries
import numpy as np

#Visualization libraries
import matplotlib.pyplot as plt

#metrics
import sklearn
import sklearn.datasets
from sklearn.metrics import r2_score

np.random.seed(1)

#=====================
# DNN Initialization
#=====================

#Sigmoid Function, aka the Fermi-Dirac distribution
def sigmoid(beta, h):
    s= 1/(1+ np.exp(-beta*h))
    return s

#Parameter initialization
def initialize_parameters(layers):

    '''
    Initialise the parameters of the model:

    L-- number of layers in the network (excluding the ouput)
    W-- weight matrix, dim: (l, l-1) initialized to a small number drawn from a standard normal distribution
        mean 0 and std 1.
        Different initializations:  *np.sqrt(2/layers[l-1]), sum_j omega_{i,j}=1
        The second initialization makes sure that the sum of the weights is 1 for each column

    b-- bias vector, dim: (l,1)
    beta-- inverse "temperature". initialized according to randn, but beta>0. We Initialise beta small, i.e.
           high temperature. Note that each unit has its own beta as it attains only local equilibrium.
           Another possible initialization of beta is to 1 for each unit.
    '''
    np.random.seed(1)

    parameters = {}
    L = len(layers)            # number of layers in the network

    for l in range(1, L):

        parameters['w' + str(l)] = np.random.randn(layers[l], layers[l-1])
        norm= np.sum(parameters['w' + str(l)],1).reshape(layers[l],1)
        parameters['w' + str(l)] = parameters['w' + str(l)]/norm

        parameters['b' + str(l)] = np.zeros((layers[l], 1))
        parameters['beta'+str(l)]= np.random.randn(layers[l],1)*0.1

    return parameters


#Post activation function
def act(h, beta, activation):
    """
    Activation functions:

    esp  --  expected signal propagation
    relu  -- zero noise limit of esp
    sigma -- sigmoid/Fermi-Dirac distribution

    """
    if activation == 'esp':
        A = h*sigmoid(beta, h)

    elif activation== 'relu':
        A= np.maximum(0,h)

    elif activation== 'sigmoid':
        A= sigmoid(beta, h)

    return A


#--------Forward propagation----------------------------------------------------------------

def FW_prop(X, parameters, activation):

    """
    Arguments:
    X-- placeholder of the input data.
    parameters-- dictionary of parameters, layer by layer, in the network.
    activations-- list of activation functions to apply to the pre-activation outputs

    Evaluates:
    A_prev --activation of the previous layer, used in the fwd pass
    cache_linear["Z"+str(l)]-- dictionary of pre-activation outputs
    cache_act["A"+str(l)]-- dictionary of post-activation outputs

    Returns:
    caches-- array containing all the post and pre- activation values, layer by layer

    """

    cache_linear={} #dictionary, cache of the linear outputs
    cache_act={} #dictionary, cache of activations
    caches = [] #list of output caches

    L= len(parameters)//3+1 # number of layers (this is counted in a different way than later)

    a_prev= X #
    caches.append((X,0) ) #Add the input as the 0th layer. This is used for backpropagation

    for l in range(1,L):

        cache_linear["h"+str(l)] = np.dot(parameters["w"+str(l)], a_prev)+ parameters["b"+str(l)]
        cache_act["a"+str(l)] = act(cache_linear["h"+str(l)], parameters['beta'+str(l)], activation[l-1])
        a_prev= cache_act["a"+str(l)]

        caches.append((cache_act["a"+str(l)], cache_linear["h"+str(l)]))  #Returns all the caches in the FWD pass

    return  caches


#---------------Objective function-----------------------------------------------------------

def obj(AL, Y, activation):

    """
    Arguments:
    AL -- activation of the output layer
    Y --  ground truth.
    activation -- activation function to apply. The last value discriminates between a classificatin (sigmoid) append
    a regression (esp ot relu) taks.

    Returns:
    cost -- cost function
    """

    L= len(activation) #number of layers
    m = Y.shape[1] #number of training examples

    last = activation[L-1]

    if last == 'sigmoid': #use cross entropy loss function
          cost = np.sum(-(np.multiply(Y,np.log(AL))+np.multiply((1-Y),np.log(1-AL)) ))/m

    elif last == 'esp' or last == 'relu': #use minimum squared error (L2 loss)
          cost = np.sum((AL-Y)**2)/(2*m)

    cost = np.squeeze(cost) #to make sure loss is a scalar.

    return cost

#reorder the parameters 
def layer_par(parameters):

    L= int(len(parameters)/3)+1

    W = [parameters['w'+str(l)] for l in range(1,L)]
    Beta = [parameters['beta'+str(l)] for l in range(1,L)]

    return W, Beta

#-------------------------One step of backpropagation-------------------------
def Back_one_step(da, dt, Wl , Betal, hl, a_prev, actl):

    '''
    Evaluate one step of backpropagation:

    Arguments:
    da -- recursive parameters to evaluate dW and db
    dt -- recursive parameters to evaluate dbeta

    Wl    -- weight of the (l-1)th unit
    Betal -- (inverse) noise value of the lth unit
    hl--  pre-activation value of the lth unit
    a_prev--  post-activation value of the previous unit

    Evaluate:
    g -- derivative of the esp at layer l
    f -- derivative of the sigmoid at layer l
    dc -- recursive parameter to compute dW and db for 'esp'
    ds -- recursive parameter to compute dbeta for esp or everything for 'sigmoid'
    da_prev -- value of dA feeding the next layer
    dt_prev -- value of dt feeding the next layer


    Return:
    grads-- dictionary containing the gradients

    '''

    #Initialise variables
    s = sigmoid(Betal, hl) #Activation value of the lth unit.


    if actl == 'esp':
        #Building Blocks
        g = s+hl*Betal*s*(1-s) #derivative of esp w.r.t. w ot b
        f = (hl**2)*s*(1-s) #derivative of esp w.r.t. beta
        dc = da*g

    elif actl == 'sigmoid':
        #Building Blocks
        g = hl*Betal*s*(1-s) #derivative of sigmoid w.r.t. w ot b
        f = hl*s*(1-s) #derivative of sigmoid w.r.t. beta
        dc = da*s*(1-s)

    # gradients of the parameters
    dWl = np.dot(dc,a_prev.T)
    dbl = np.sum(dc, axis=1, keepdims=True)
    dBetal= np.sum(dt*f, axis=1, keepdims=True)

    #Next
    da_prev = np.dot(Wl.T, dc) #update dWl, dbl for the next BW step
    dt_prev = np.dot(Wl.T, (dt*g)) #update dbeta for the next BW step

    assert (da_prev.shape == a_prev.shape)
    assert (dWl.shape == Wl.shape)
    assert (dBetal.shape == Betal.shape)

    return da_prev, dt_prev, dWl, dbl, dBetal

#Back prop through the all network

def BKW_prop(Y, caches, parameters, activation):

    """
    Implement the backward propagation

    Arguments:
    Y -- ground truth
    caches -- list of caches of pre/post activation values of each unit for every layer [1,L]
    parameters -- list of parameters in the model for every layer

    Returns:
    grads -- A dictionary with the gradients
             grads["dA" + str(l)] = ...
             grads["dt" + str(l)] = ...
             grads["dW" + str(l)] = ...
             grads["db" + str(l)] = ...
             grads["dbeta" + str(l)] = ...
    """

    grads = {}
    a_prev =[]

    L= len(parameters)//3 # number of layers (integer)

    W, Beta = layer_par(parameters) #obtain W and beta for each layer

    m  = Y.shape[1] # number of training examples

    AL = caches[L][0] #post-activation value of the last unit

    actL= activation[L-1]

    assert(Y.shape == AL.shape) #make sure that the dimenstions of the final layer and labels coincides

    # Initializing the backpropagation
    if actL == 'sigmoid':
       delta = - ( np.divide(Y, AL) - np.divide(1 - Y, 1 - AL) )/m

    elif actL == 'esp':
        delta = (AL - Y)/m

    # Initialise the output layer

    grads["da" + str(L)], grads["dt" + str(L)], grads["dw" + str(L)], grads["db" + str(L)], grads["dbeta" + str(L)] = Back_one_step(delta, delta, W[L-1] , Beta[L-1], caches[L][1], caches[L-1][0],actL )

    for l in reversed(range(L-1)):

        da_tmp, dt_tmp ,dWl_tmp, dbl_tmp, dbetal_tmp = Back_one_step( grads["da" + str(l+2)], grads["dt" + str(l+2)], W[l] , Beta[l], caches[l+1][1], caches[l][0] , activation[l])

        grads["da" + str(l + 1)] = da_tmp
        grads["dt" + str(l + 1)] = dt_tmp
        grads["dw" + str(l + 1)] = dWl_tmp
        grads["db" + str(l + 1)] = dbl_tmp
        grads["dbeta" + str(l + 1)] = dbetal_tmp

    return grads

#=================================
#Update rule for the parameters
#=================================
#This is a simple gradient descent update rule. More advanced update rules can be implemented.


def update_parameters(parameters, grads, learning_rate):
    """
    Update parameters using pure gradient descent

    Arguments:
    parameters --  dictionary containing W, b, beta for each layer
    grads --  dictionary containing the gradients.

    Returns:
    parameters -- python dictionary containing the updated parameters
                  parameters["W" + str(l)] = ...
                  parameters["b" + str(l)] = ...
                  parameters["beta" + str(l)] = ...
    """

    L = len(parameters)//3 # number of layers (integer)

    # Update paramters
    for l in range(L):
        parameters['w'+str(l+1)]= parameters['w'+str(l+1)] - learning_rate* grads["dw" + str(l+1)]

        parameters['b'+str(l+1)]= parameters['b'+str(l+1)] - learning_rate* grads["db" + str(l+1)]
        parameters['beta'+str(l+1)] = parameters['beta'+str(l+1)] - learning_rate* grads["dbeta" + str(l+1)]

    return parameters


#=======================
#Main Program
#=======================


def Run_DNN(x, y, layers, activation, learning_rate, num_iterations, print_cost=False):

    """
    Run the DNN to find the optimal set of paramters

    Arguments:
    X -- data, iput marix
    Y -- true "label" vector
    layers -- list containing the input size and each layer size
    learning_rate -- learning rate of the gradient descent update rule
    num_iterations -- number of iterations of the optimization loop
    print_cost -- if True, it prints the cost every 100 steps

    Returns:
    parameters -- parameters learnt by the model. They can then be used in the predictor to find the ouput
    """
    #Initialise parameters
    parameters=0

    costs = []                         # store the value of the objective function for each epoch
    gradients= []                      # store the value of the gradients for each epoch

    #Take the transpose of the inputs to abide to our convention
    X = x.T
    Y = y.T

    f , m = X.shape # f: number of features, m: number of training examples

    #Initialise network
    network = np.append(f, layers)
    L= len(network)

    # Parameters initialization.
    parameters= initialize_parameters(network)

    # Loop (gradient stochastic gradient descent)
    for i in range(0, num_iterations):

        # Forward propagation
        caches = FW_prop(X, parameters, activation)

        # Compute cost.
        cost = obj(caches[L-1][0] , Y, activation)

        # Backward propagation. It outputs the value of the gradients for the current epoch. This can be used to study the dynamics
        #of the optimization algo
        grads = BKW_prop(Y, caches, parameters, activation)

        # Update parameters.
        parameters = update_parameters(parameters, grads, learning_rate)

        # Print the cost every 1000 training example
        if print_cost and i % 1000 == 0:
            print ("Cost after iteration %i: %f" %(i, cost))

        #store the value of cost and grads every 1000 epochs (chanche this value)
        if print_cost and i % 1000 == 0:
            costs.append(cost)
            gradients.append(grads)

    # plot the cost
    plt.plot(np.squeeze(costs))
    plt.ylabel('cost')
    plt.xlabel('iterations (per tens)')
    plt.title("Learning rate =" + str(learning_rate))
    plt.show()

    return parameters, gradients, costs


#======================
#Prediction functions
#======================

#Prediction and metrics for Regression Models

def predict_reg(x, y, parameters, activation):

    """
    Arguments:
    X-- feature matrix from training/test set
    Y-- ground truth, label matrix   from training/test set
    activation-- array contaiing the type of activation function of each layer
    parameters-- dictionary containing the trained parameters

    Return:
    Yp-- predicted label matrix
    accuracy-- prediction accuracy
    """

    #Take transpose of input
    X = x.T
    Y = y.T

    L = len(parameters)//3
    m = X.shape[1]

    (n_y, m) = y.shape #shape of the ground truth matrix

    caches = FW_prop(X, parameters, activation)

    #An= np.transpose(cache_A["A"+str(L-1)]) #last activation
    Yp = caches[L][0] # predicted output

    accuracy= r2_score(y, Yp.T)

    return  Yp.T , accuracy

#Prediction and metrics for Classification Models

def predict_class(x, y, parameters, activation):
    """
    This function is used to predict the results of a  n-layer neural network.

    Arguments:
    X -- data set of examples you would like to label
    parameters -- parameters of the trained model

    Returns:
    p -- predictions for the given dataset X
    """
    #Take transpose of input
    X = x.T
    Y = y.T

    L = len(parameters)//3
    m = X.shape[1]
    Yp = np.zeros((1,m), dtype = np.int)

    # Forward propagation
    caches = FW_prop(X, parameters, activation)

    p = caches[L][0] # predicted probability of 0/1

    # convert probs to 0/1 predictions. If p>1/2 Yp=1 otherwise 0. Note that a dynamic threshold coulb be implemented
    for i in range(0, m):
        if p[0,i] > 0.5:
            Yp[0,i] = 1
        else:
            Yp[0,i] = 0

    #evaluate accuracy for single label prediction
    accuracy = np.mean((Yp[0,:] == Y[0,:]) )

    return p, Yp, accuracy

#General Prediction function

def predict(x, y, parameters, activation):

    '''
    Prediction function. It uses either predict_reg or predict_class according to wether the ouput activation
    is esp (relu) or sigmoid

    Returns:
    Predicted value and accuracy. Other metrics can be used than the one proposed here.

    '''

    L = len(parameters)//3

    actL = activation[L-1]

    if actL == 'sigmoid':
        p, Yp, accuracy = predict_class(x, y, parameters, activation)
        model = [p, Yp, accuracy]

    elif actL == 'esp':

        Yp, accuracy = predict_reg(x, y, parameters, activation)
        model = [Yp, accuracy]

    print('accuracy=', accuracy)

    return model

"""
Helper functions for FNN with ESP
=================================================================
Author: <NAME> <<EMAIL>>

This is a pure python implementation of a Feed Forward Deep network with ESP
activation. See main text for details.

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
#Math Libraries
import numpy as np

#Visualization libraries
import matplotlib.pyplot as plt

#metrics
import sklearn
import sklearn.datasets
from sklearn.metrics import r2_score

np.random.seed(1)

#=====================
# DNN Initialization
#=====================

#Sigmoid Function, aka the Fermi-Dirac distribution
def sigmoid(beta, h):
    s= 1/(1+ np.exp(-beta*h))
    return s

#Parameter initialization
def initialize_parameters(layers):

    '''
    Initialise the parameters of the model:

    L-- number of layers in the network (excluding the ouput)
    W-- weight matrix, dim: (l, l-1) initialized to a small number drawn from a standard normal distribution
        mean 0 and std 1.
        Different initializations:  *np.sqrt(2/layers[l-1]), sum_j omega_{i,j}=1
        The second initialization makes sure that the sum of the weights is 1 for each column

    b-- bias vector, dim: (l,1)
    beta-- inverse "temperature". initialized according to randn, but beta>0. We Initialise beta small, i.e.
           high temperature. Note that each unit has its own beta as it attains only local equilibrium.
           Another possible initialization of beta is to 1 for each unit.
    '''
    np.random.seed(1)

    parameters = {}
    L = len(layers)            # number of layers in the network

    for l in range(1, L):

        parameters['w' + str(l)] = np.random.randn(layers[l], layers[l-1])
        norm= np.sum(parameters['w' + str(l)],1).reshape(layers[l],1)
        parameters['w' + str(l)] = parameters['w' + str(l)]/norm

        parameters['b' + str(l)] = np.zeros((layers[l], 1))
        parameters['beta'+str(l)]= np.random.randn(layers[l],1)*0.1

    return parameters


#Post activation function
def act(h, beta, activation):
    """
    Activation functions:

    esp  --  expected signal propagation
    relu  -- zero noise limit of esp
    sigma -- sigmoid/Fermi-Dirac distribution

    """
    if activation == 'esp':
        A = h*sigmoid(beta, h)

    elif activation== 'relu':
        A= np.maximum(0,h)

    elif activation== 'sigmoid':
        A= sigmoid(beta, h)

    return A


#--------Forward propagation----------------------------------------------------------------

def FW_prop(X, parameters, activation):

    """
    Arguments:
    X-- placeholder of the input data.
    parameters-- dictionary of parameters, layer by layer, in the network.
    activations-- list of activation functions to apply to the pre-activation outputs

    Evaluates:
    A_prev --activation of the previous layer, used in the fwd pass
    cache_linear["Z"+str(l)]-- dictionary of pre-activation outputs
    cache_act["A"+str(l)]-- dictionary of post-activation outputs

    Returns:
    caches-- array containing all the post and pre- activation values, layer by layer

    """

    cache_linear={} #dictionary, cache of the linear outputs
    cache_act={} #dictionary, cache of activations
    caches = [] #list of output caches

    L= len(parameters)//3+1 # number of layers (this is counted in a different way than later)

    a_prev= X #
    caches.append((X,0) ) #Add the input as the 0th layer. This is used for backpropagation

    for l in range(1,L):

        cache_linear["h"+str(l)] = np.dot(parameters["w"+str(l)], a_prev)+ parameters["b"+str(l)]
        cache_act["a"+str(l)] = act(cache_linear["h"+str(l)], parameters['beta'+str(l)], activation[l-1])
        a_prev= cache_act["a"+str(l)]

        caches.append((cache_act["a"+str(l)], cache_linear["h"+str(l)]))  #Returns all the caches in the FWD pass

    return  caches


#---------------Objective function-----------------------------------------------------------

def obj(AL, Y, activation):

    """
    Arguments:
    AL -- activation of the output layer
    Y --  ground truth.
    activation -- activation function to apply. The last value discriminates between a classificatin (sigmoid) append
    a regression (esp ot relu) taks.

    Returns:
    cost -- cost function
    """

    L= len(activation) #number of layers
    m = Y.shape[1] #number of training examples

    last = activation[L-1]

    if last == 'sigmoid': #use cross entropy loss function
          cost = np.sum(-(np.multiply(Y,np.log(AL))+np.multiply((1-Y),np.log(1-AL)) ))/m

    elif last == 'esp' or last == 'relu': #use minimum squared error (L2 loss)
          cost = np.sum((AL-Y)**2)/(2*m)

    cost = np.squeeze(cost) #to make sure loss is a scalar.

    return cost

#reorder the parameters 
def layer_par(parameters):

    L= int(len(parameters)/3)+1

    W = [parameters['w'+str(l)] for l in range(1,L)]
    Beta = [parameters['beta'+str(l)] for l in range(1,L)]

    return W, Beta

#-------------------------One step of backpropagation-------------------------
def Back_one_step(da, dt, Wl , Betal, hl, a_prev, actl):

    '''
    Evaluate one step of backpropagation:

    Arguments:
    da -- recursive parameters to evaluate dW and db
    dt -- recursive parameters to evaluate dbeta

    Wl    -- weight of the (l-1)th unit
    Betal -- (inverse) noise value of the lth unit
    hl--  pre-activation value of the lth unit
    a_prev--  post-activation value of the previous unit

    Evaluate:
    g -- derivative of the esp at layer l
    f -- derivative of the sigmoid at layer l
    dc -- recursive parameter to compute dW and db for 'esp'
    ds -- recursive parameter to compute dbeta for esp or everything for 'sigmoid'
    da_prev -- value of dA feeding the next layer
    dt_prev -- value of dt feeding the next layer


    Return:
    grads-- dictionary containing the gradients

    '''

    #Initialise variables
    s = sigmoid(Betal, hl) #Activation value of the lth unit.


    if actl == 'esp':
        #Building Blocks
        g = s+hl*Betal*s*(1-s) #derivative of esp w.r.t. w ot b
        f = (hl**2)*s*(1-s) #derivative of esp w.r.t. beta
        dc = da*g

    elif actl == 'sigmoid':
        #Building Blocks
        g = hl*Betal*s*(1-s) #derivative of sigmoid w.r.t. w ot b
        f = hl*s*(1-s) #derivative of sigmoid w.r.t. beta
        dc = da*s*(1-s)

    # gradients of the parameters
    dWl = np.dot(dc,a_prev.T)
    dbl = np.sum(dc, axis=1, keepdims=True)
    dBetal= np.sum(dt*f, axis=1, keepdims=True)

    #Next
    da_prev = np.dot(Wl.T, dc) #update dWl, dbl for the next BW step
    dt_prev = np.dot(Wl.T, (dt*g)) #update dbeta for the next BW step

    assert (da_prev.shape == a_prev.shape)
    assert (dWl.shape == Wl.shape)
    assert (dBetal.shape == Betal.shape)

    return da_prev, dt_prev, dWl, dbl, dBetal

#Back prop through the all network

def BKW_prop(Y, caches, parameters, activation):

    """
    Implement the backward propagation

    Arguments:
    Y -- ground truth
    caches -- list of caches of pre/post activation values of each unit for every layer [1,L]
    parameters -- list of parameters in the model for every layer

    Returns:
    grads -- A dictionary with the gradients
             grads["dA" + str(l)] = ...
             grads["dt" + str(l)] = ...
             grads["dW" + str(l)] = ...
             grads["db" + str(l)] = ...
             grads["dbeta" + str(l)] = ...
    """

    grads = {}
    a_prev =[]

    L= len(parameters)//3 # number of layers (integer)

    W, Beta = layer_par(parameters) #obtain W and beta for each layer

    m  = Y.shape[1] # number of training examples

    AL = caches[L][0] #post-activation value of the last unit

    actL= activation[L-1]

    assert(Y.shape == AL.shape) #make sure that the dimenstions of the final layer and labels coincides

    # Initializing the backpropagation
    if actL == 'sigmoid':
       delta = - ( np.divide(Y, AL) - np.divide(1 - Y, 1 - AL) )/m

    elif actL == 'esp':
        delta = (AL - Y)/m

    # Initialise the output layer

    grads["da" + str(L)], grads["dt" + str(L)], grads["dw" + str(L)], grads["db" + str(L)], grads["dbeta" + str(L)] = Back_one_step(delta, delta, W[L-1] , Beta[L-1], caches[L][1], caches[L-1][0],actL )

    for l in reversed(range(L-1)):

        da_tmp, dt_tmp ,dWl_tmp, dbl_tmp, dbetal_tmp = Back_one_step( grads["da" + str(l+2)], grads["dt" + str(l+2)], W[l] , Beta[l], caches[l+1][1], caches[l][0] , activation[l])

        grads["da" + str(l + 1)] = da_tmp
        grads["dt" + str(l + 1)] = dt_tmp
        grads["dw" + str(l + 1)] = dWl_tmp
        grads["db" + str(l + 1)] = dbl_tmp
        grads["dbeta" + str(l + 1)] = dbetal_tmp

    return grads

#=================================
#Update rule for the parameters
#=================================
#This is a simple gradient descent update rule. More advanced update rules can be implemented.


def update_parameters(parameters, grads, learning_rate):
    """
    Update parameters using pure gradient descent

    Arguments:
    parameters --  dictionary containing W, b, beta for each layer
    grads --  dictionary containing the gradients.

    Returns:
    parameters -- python dictionary containing the updated parameters
                  parameters["W" + str(l)] = ...
                  parameters["b" + str(l)] = ...
                  parameters["beta" + str(l)] = ...
    """

    L = len(parameters)//3 # number of layers (integer)

    # Update paramters
    for l in range(L):
        parameters['w'+str(l+1)]= parameters['w'+str(l+1)] - learning_rate* grads["dw" + str(l+1)]

        parameters['b'+str(l+1)]= parameters['b'+str(l+1)] - learning_rate* grads["db" + str(l+1)]
        parameters['beta'+str(l+1)] = parameters['beta'+str(l+1)] - learning_rate* grads["dbeta" + str(l+1)]

    return parameters


#=======================
#Main Program
#=======================


def Run_DNN(x, y, layers, activation, learning_rate, num_iterations, print_cost=False):

    """
    Run the DNN to find the optimal set of paramters

    Arguments:
    X -- data, iput marix
    Y -- true "label" vector
    layers -- list containing the input size and each layer size
    learning_rate -- learning rate of the gradient descent update rule
    num_iterations -- number of iterations of the optimization loop
    print_cost -- if True, it prints the cost every 100 steps

    Returns:
    parameters -- parameters learnt by the model. They can then be used in the predictor to find the ouput
    """
    #Initialise parameters
    parameters=0

    costs = []                         # store the value of the objective function for each epoch
    gradients= []                      # store the value of the gradients for each epoch

    #Take the transpose of the inputs to abide to our convention
    X = x.T
    Y = y.T

    f , m = X.shape # f: number of features, m: number of training examples

    #Initialise network
    network = np.append(f, layers)
    L= len(network)

    # Parameters initialization.
    parameters= initialize_parameters(network)

    # Loop (gradient stochastic gradient descent)
    for i in range(0, num_iterations):

        # Forward propagation
        caches = FW_prop(X, parameters, activation)

        # Compute cost.
        cost = obj(caches[L-1][0] , Y, activation)

        # Backward propagation. It outputs the value of the gradients for the current epoch. This can be used to study the dynamics
        #of the optimization algo
        grads = BKW_prop(Y, caches, parameters, activation)

        # Update parameters.
        parameters = update_parameters(parameters, grads, learning_rate)

        # Print the cost every 1000 training example
        if print_cost and i % 1000 == 0:
            print ("Cost after iteration %i: %f" %(i, cost))

        #store the value of cost and grads every 1000 epochs (chanche this value)
        if print_cost and i % 1000 == 0:
            costs.append(cost)
            gradients.append(grads)

    # plot the cost
    plt.plot(np.squeeze(costs))
    plt.ylabel('cost')
    plt.xlabel('iterations (per tens)')
    plt.title("Learning rate =" + str(learning_rate))
    plt.show()

    return parameters, gradients, costs


#======================
#Prediction functions
#======================

#Prediction and metrics for Regression Models

def predict_reg(x, y, parameters, activation):

    """
    Arguments:
    X-- feature matrix from training/test set
    Y-- ground truth, label matrix   from training/test set
    activation-- array contaiing the type of activation function of each layer
    parameters-- dictionary containing the trained parameters

    Return:
    Yp-- predicted label matrix
    accuracy-- prediction accuracy
    """

    #Take transpose of input
    X = x.T
    Y = y.T

    L = len(parameters)//3
    m = X.shape[1]

    (n_y, m) = y.shape #shape of the ground truth matrix

    caches = FW_prop(X, parameters, activation)

    #An= np.transpose(cache_A["A"+str(L-1)]) #last activation
    Yp = caches[L][0] # predicted output

    accuracy= r2_score(y, Yp.T)

    return  Yp.T , accuracy

#Prediction and metrics for Classification Models

def predict_class(x, y, parameters, activation):
    """
    This function is used to predict the results of a  n-layer neural network.

    Arguments:
    X -- data set of examples you would like to label
    parameters -- parameters of the trained model

    Returns:
    p -- predictions for the given dataset X
    """
    #Take transpose of input
    X = x.T
    Y = y.T

    L = len(parameters)//3
    m = X.shape[1]
    Yp = np.zeros((1,m), dtype = np.int)

    # Forward propagation
    caches = FW_prop(X, parameters, activation)

    p = caches[L][0] # predicted probability of 0/1

    # convert probs to 0/1 predictions. If p>1/2 Yp=1 otherwise 0. Note that a dynamic threshold coulb be implemented
    for i in range(0, m):
        if p[0,i] > 0.5:
            Yp[0,i] = 1
        else:
            Yp[0,i] = 0

    #evaluate accuracy for single label prediction
    accuracy = np.mean((Yp[0,:] == Y[0,:]) )

    return p, Yp, accuracy

#General Prediction function

def predict(x, y, parameters, activation):

    '''
    Prediction function. It uses either predict_reg or predict_class according to wether the ouput activation
    is esp (relu) or sigmoid

    Returns:
    Predicted value and accuracy. Other metrics can be used than the one proposed here.

    '''

    L = len(parameters)//3

    actL = activation[L-1]

    if actL == 'sigmoid':
        p, Yp, accuracy = predict_class(x, y, parameters, activation)
        model = [p, Yp, accuracy]

    elif actL == 'esp':

        Yp, accuracy = predict_reg(x, y, parameters, activation)
        model = [Yp, accuracy]

    print('accuracy=', accuracy)

    return model

Helper functions for FNN with ESP ================================================================= Author: <NAME> <<EMAIL>>  This is a pure python implementation of a Feed Forward Deep network with ESP activation. See main text for details.  # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at # #     http://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. #Math Libraries #Visualization libraries #metrics #===================== # DNN Initialization #===================== #Sigmoid Function, aka the Fermi-Dirac distribution #Parameter initialization Initialise the parameters of the model:      L-- number of layers in the network (excluding the ouput)     W-- weight matrix, dim: (l, l-1) initialized to a small number drawn from a standard normal distribution         mean 0 and std 1.         Different initializations:  *np.sqrt(2/layers[l-1]), sum_j omega_{i,j}=1         The second initialization makes sure that the sum of the weights is 1 for each column      b-- bias vector, dim: (l,1)     beta-- inverse "temperature". initialized according to randn, but beta>0. We Initialise beta small, i.e.            high temperature. Note that each unit has its own beta as it attains only local equilibrium.            Another possible initialization of beta is to 1 for each unit. # number of layers in the network #Post activation function Activation functions:      esp  --  expected signal propagation     relu  -- zero noise limit of esp     sigma -- sigmoid/Fermi-Dirac distribution #--------Forward propagation---------------------------------------------------------------- Arguments:     X-- placeholder of the input data.     parameters-- dictionary of parameters, layer by layer, in the network.     activations-- list of activation functions to apply to the pre-activation outputs      Evaluates:     A_prev --activation of the previous layer, used in the fwd pass     cache_linear["Z"+str(l)]-- dictionary of pre-activation outputs     cache_act["A"+str(l)]-- dictionary of post-activation outputs      Returns:     caches-- array containing all the post and pre- activation values, layer by layer #dictionary, cache of the linear outputs #dictionary, cache of activations #list of output caches # number of layers (this is counted in a different way than later) # #Add the input as the 0th layer. This is used for backpropagation #Returns all the caches in the FWD pass #---------------Objective function----------------------------------------------------------- Arguments:     AL -- activation of the output layer     Y --  ground truth.     activation -- activation function to apply. The last value discriminates between a classificatin (sigmoid) append     a regression (esp ot relu) taks.      Returns:     cost -- cost function #number of layers #number of training examples #use cross entropy loss function #use minimum squared error (L2 loss) #to make sure loss is a scalar. #reorder the parameters #-------------------------One step of backpropagation------------------------- Evaluate one step of backpropagation:      Arguments:     da -- recursive parameters to evaluate dW and db     dt -- recursive parameters to evaluate dbeta      Wl    -- weight of the (l-1)th unit     Betal -- (inverse) noise value of the lth unit     hl--  pre-activation value of the lth unit     a_prev--  post-activation value of the previous unit      Evaluate:     g -- derivative of the esp at layer l     f -- derivative of the sigmoid at layer l     dc -- recursive parameter to compute dW and db for 'esp'     ds -- recursive parameter to compute dbeta for esp or everything for 'sigmoid'     da_prev -- value of dA feeding the next layer     dt_prev -- value of dt feeding the next layer       Return:     grads-- dictionary containing the gradients #Initialise variables #Activation value of the lth unit. #Building Blocks #derivative of esp w.r.t. w ot b #derivative of esp w.r.t. beta #Building Blocks #derivative of sigmoid w.r.t. w ot b #derivative of sigmoid w.r.t. beta # gradients of the parameters #Next #update dWl, dbl for the next BW step #update dbeta for the next BW step #Back prop through the all network Implement the backward propagation      Arguments:     Y -- ground truth     caches -- list of caches of pre/post activation values of each unit for every layer [1,L]     parameters -- list of parameters in the model for every layer      Returns:     grads -- A dictionary with the gradients              grads["dA" + str(l)] = ...              grads["dt" + str(l)] = ...              grads["dW" + str(l)] = ...              grads["db" + str(l)] = ...              grads["dbeta" + str(l)] = ... # number of layers (integer) #obtain W and beta for each layer # number of training examples #post-activation value of the last unit #make sure that the dimenstions of the final layer and labels coincides # Initializing the backpropagation # Initialise the output layer #================================= #Update rule for the parameters #================================= #This is a simple gradient descent update rule. More advanced update rules can be implemented. Update parameters using pure gradient descent      Arguments:     parameters --  dictionary containing W, b, beta for each layer     grads --  dictionary containing the gradients.      Returns:     parameters -- python dictionary containing the updated parameters                   parameters["W" + str(l)] = ...                   parameters["b" + str(l)] = ...                   parameters["beta" + str(l)] = ... # number of layers (integer) # Update paramters #======================= #Main Program #======================= Run the DNN to find the optimal set of paramters      Arguments:     X -- data, iput marix     Y -- true "label" vector     layers -- list containing the input size and each layer size     learning_rate -- learning rate of the gradient descent update rule     num_iterations -- number of iterations of the optimization loop     print_cost -- if True, it prints the cost every 100 steps      Returns:     parameters -- parameters learnt by the model. They can then be used in the predictor to find the ouput #Initialise parameters # store the value of the objective function for each epoch # store the value of the gradients for each epoch #Take the transpose of the inputs to abide to our convention # f: number of features, m: number of training examples #Initialise network # Parameters initialization. # Loop (gradient stochastic gradient descent) # Forward propagation # Compute cost. # Backward propagation. It outputs the value of the gradients for the current epoch. This can be used to study the dynamics #of the optimization algo # Update parameters. # Print the cost every 1000 training example #store the value of cost and grads every 1000 epochs (chanche this value) # plot the cost #====================== #Prediction functions #====================== #Prediction and metrics for Regression Models Arguments:     X-- feature matrix from training/test set     Y-- ground truth, label matrix   from training/test set     activation-- array contaiing the type of activation function of each layer     parameters-- dictionary containing the trained parameters      Return:     Yp-- predicted label matrix     accuracy-- prediction accuracy #Take transpose of input #shape of the ground truth matrix #An= np.transpose(cache_A["A"+str(L-1)]) #last activation # predicted output #Prediction and metrics for Classification Models This function is used to predict the results of a  n-layer neural network.      Arguments:     X -- data set of examples you would like to label     parameters -- parameters of the trained model      Returns:     p -- predictions for the given dataset X #Take transpose of input # Forward propagation # predicted probability of 0/1 # convert probs to 0/1 predictions. If p>1/2 Yp=1 otherwise 0. Note that a dynamic threshold coulb be implemented #evaluate accuracy for single label prediction #General Prediction function Prediction function. It uses either predict_reg or predict_class according to wether the ouput activation     is esp (relu) or sigmoid      Returns:     Predicted value and accuracy. Other metrics can be used than the one proposed here.