<jupyter_start><jupyter_text>Neural Network in pure python<jupyter_code># import libraries
import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Import training and test set<jupyter_code>Z_training = np.loadtxt(fname = 'training_set.txt', delimiter=',')
Z_training.shape
Z_test = np.loadtxt(fname = 'test_set.txt', delimiter=',')
Z_test.shape
x_training, y_training = Z_training.T
x_test, y_test = Z_test.T<jupyter_output><empty_output><jupyter_text>Training and test set visualization<jupyter_code>plt.plot(x_training,y_training, 'b.', label = 'training set')
plt.plot(x_test,y_test, 'r.', label = 'test set')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>There is an interval [-3,-2] where x data are missing; this is likely to be the interval where the real test will be done. Also x > 2 could be a good interval to be tested. We should see in the end that the predictions of the final trained model in those intervals aren't crazy.  Define other activation functions<jupyter_code>from scipy.special import expit
# expit(x) = 1/(1+exp(-x))
def sigmoid(x):
    return expit(x)

# 1st derivative
#sigmoid_der = lambda x: sigmoid(x) * (1 - sigmoid(x))

def sigmoid_der(x):
    return sigmoid(x) * (1 - sigmoid(x))
# vectorized implementation

def ReLU(x):
    x = x.reshape((len(x),1))
    zeros = np.zeros(x.shape)
    z = np.stack((x,zeros), axis=1)
    return np.max(z, axis=1)

def ReLU_der(x):
    z = np.zeros(x.shape)
    z[x>0] = 1
    return z
# vectorized implementation

def LeakyReLU(x):
    x = x.reshape((len(x),1))
    ones = x*0.1
    z = np.stack((x,ones), axis=1)
    return np.max(z, axis=1)

def LeakyReLU_der(x):
    z = np.ones(x.shape)
    z[x<0] = 0.1
    return z<jupyter_output><empty_output><jupyter_text>Show learned function<jupyter_code>def show_forecast_vs_groundtruth(NN, x_test, y_test):
    x_axis = np.linspace(-4,4,4000)
    plt.plot(x_test, [NN.forward(x) for x in x_test], 'b.', label = 'predictions')
    plt.plot(x_test, y_test, '.r', label = 'ground truth')
    plt.plot(x_axis, [NN.forward(x) for x in x_axis], 'g--', lw = 0.7, label =  'learned function' )
    plt.legend()
    plt.show()<jupyter_output><empty_output><jupyter_text>Define NN class with more methodsThe idea is to incorporate all the parameters of the training inside the NN, so that then we can implement an automatic grid search with cross validation.Skelethon  of the class:```pythonclass Network():        def __init__(self, N_neurons, act_func, act_der_func, lr, n_epochs, \                 en_decay=True, lr_final=1e-4, early_stopping = True, tol = 1e-3, \                 en_penalty=False, penalty=1e-4, en_grad_clipping=True, grad_treshold=10)        def load_weights(self, WBh1, WBh2, WBo)        def save_weights(self)        def forward(self, x, additional_out=False)        def update(self, x, label)        def train(self, x_train, y_train, x_val, y_val, train_log=False, verbose=False)        def evalutate_mean_loss(self, x_test, y_test)            def plot_weights(self)```<jupyter_code>#%% Network class

class Network():
    
    def __init__(self, N_neurons, act_func, act_der_func, lr, n_epochs, \
                 en_decay=True, lr_final=1e-4, early_stopping = True, tol = 1e-3, \
                 en_penalty=False, penalty=1e-4, en_grad_clipping=True, grad_treshold=10):
        
        Ni, Nh1, Nh2, No = N_neurons
        ### WEIGHT INITIALIZATION (Xavier)
        # Initialize hidden weights and biases (layer 1)
        Wh1 = (np.random.rand(Nh1, Ni) - 0.5) * np.sqrt(12 / (Nh1 + Ni))
        Bh1 = np.zeros([Nh1, 1])
        self.WBh1 = np.concatenate([Wh1, Bh1], 1) # Weight matrix including biases
        # Initialize hidden weights and biases (layer 2)
        Wh2 = (np.random.rand(Nh2, Nh1) - 0.5) * np.sqrt(12 / (Nh2 + Nh1))
        Bh2 = np.zeros([Nh2, 1])
        self.WBh2 = np.concatenate([Wh2, Bh2], 1) # Weight matrix including biases
        # Initialize output weights and biases
        Wo = (np.random.rand(No, Nh2) - 0.5) * np.sqrt(12 / (No + Nh2))
        Bo = np.zeros([No, 1])
        self.WBo = np.concatenate([Wo, Bo], 1) # Weight matrix including biases
        
        ### ACTIVATION FUNCTION
        self.act = act_func
        self.act_der = act_der_func
        
        # TRAINING PARAMETERS
        self.lr = lr
        self.n_epochs = n_epochs
        self.en_decay = en_decay
        self.tolerance = tol
        if en_decay == True:
            self.lr_decay = (lr_final / lr)**(1 / n_epochs)
        else:
            self.lr_decay = None
        
        if early_stopping == True:
            self.early_stopping = True
        else:
            self.early_stopping = False
        
        if (en_penalty == 'L1'):
            self.penalty = penalty
            self.penalty_type = 'L1'
        elif (en_penalty == 'L2'):
            self.penalty = penalty
            self.penalty_type = 'L2'
        else:
            self.penalty = None
            self.penalty_type = None
            
        if en_grad_clipping == True:
            self.grad_treshold = grad_treshold
        else:
            self.grad_treshold = None
    
    def load_weights(self, WBh1, WBh2, WBo):
        self.WBh1 = WBh1
        self.WBh2 = WBh2
        self.WBo = WBo
        
    def save_weights(self):
        np.save("WBh1", self.WBh1)
        np.save("WBh2", self.WBh2)
        np.save("WBo", self.WBo)
        
    def forward(self, x, additional_out=False):
        
        # Convert to numpy array
        x = np.array(x)
        
        ### Hidden layer 1
        # Add bias term
        X = np.append(x, 1)
        # Forward pass (linear)
        H1 = np.matmul(self.WBh1, X)
        # Activation function
        Z1 = self.act(H1)
        
        ### Hidden layer 2
        # Add bias term
        Z1 = np.append(Z1, 1)
        # Forward pass (linear)
        H2 = np.matmul(self.WBh2, Z1)
        # Activation function
        Z2 = self.act(H2)
        
        ### Output layer
        # Add bias term
        Z2 = np.append(Z2, 1)
        # Forward pass (linear)
        Y = np.matmul(self.WBo, Z2)
        # NO activation function
        
        if additional_out:
            return Y.squeeze(), Z2
        
        return Y.squeeze()
        
    def update(self, x, label):
        
        # Convert to numpy array
        X = np.array(x)
        
        ### Hidden layer 1
        # Add bias term
        X = np.append(X, 1)
        # Forward pass (linear)
        H1 = np.matmul(self.WBh1, X)
        # Activation function
        Z1 = self.act(H1)
        
        ### Hidden layer 2
        # Add bias term
        Z1 = np.append(Z1, 1)
        # Forward pass (linear)
        H2 = np.matmul(self.WBh2, Z1)
        # Activation function
        Z2 = self.act(H2)
        
        ### Output layer
        # Add bias term
        Z2 = np.append(Z2, 1)
        # Forward pass (linear)
        Y = np.matmul(self.WBo, Z2)
        # NO activation function
        
        # Evaluate the derivative terms
        D1 = Y - label
        D2 = Z2
        D3 = self.WBo[:,:-1]
        D4 = self.act_der(H2)
        D5 = Z1
        D6 = self.WBh2[:,:-1]
        D7 = self.act_der(H1)
        D8 = X
        
        # Layer Error
        Eo = D1
        Eh2 = np.matmul(Eo, D3) * D4
        Eh1 = np.matmul(Eh2, D6) * D7
         
        # Derivative for weight matrices
        dWBo = np.matmul(Eo.reshape(-1,1), D2.reshape(1,-1))
        dWBh2 = np.matmul(Eh2.reshape(-1,1), D5.reshape(1,-1))
        dWBh1 = np.matmul(Eh1.reshape(-1,1), D8.reshape(1,-1))
        
        if self.grad_treshold != None:
            # compute the norm of the whole gradient
            grad_norm = np.linalg.norm(np.concatenate([dWBo.flatten(),dWBh2.flatten(),dWBh1.flatten()]))
            # if the norm exceeds the threshold, we rescale component wise so that the norm is clipped to the treshold
            if grad_norm > self.grad_treshold :
                dWBo = dWBo/grad_norm*self.grad_treshold 
                dWBh1 = dWBh1/grad_norm*self.grad_treshold 
                dWBh2 = dWBh2/grad_norm*self.grad_treshold 
      
        if self.penalty_type == 'L1':
            
            # Compute penalty 
            WBh1_sgn = np.ones(self.WBh1.shape)
            WBh1_sgn[self.WBh1<0] = -1

            WBh2_sgn = np.ones(self.WBh2.shape)
            WBh2_sgn[self.WBh2<0] = -1

            WBo_sgn = np.ones(self.WBo.shape)
            WBo_sgn[self.WBo<0] = -1
        
            # Update the weights
            self.WBh1 -= self.lr * (dWBh1 + self.penalty * WBh1_sgn)
            self.WBh2 -= self.lr * (dWBh2 + self.penalty * WBh2_sgn)
            self.WBo -= self.lr * (dWBo + self.penalty * WBo_sgn)
            
            # Compute the loss
            loss = (Y - label)**2/2 + \
                   self.penalty * (np.abs(self.WBh1).sum() + np.abs(self.WBh2).sum() + np.abs(self.WBo).sum())
            
        elif self.penalty_type == 'L2':
        
            # Update the weights
            self.WBh1 -= self.lr * (dWBh1 + self.penalty * self.WBh1)
            self.WBh2 -= self.lr * (dWBh2 + self.penalty * self.WBh2)
            self.WBo -= self.lr * (dWBo + self.penalty * self.WBo)
            
            # Compute the loss
            loss = (Y - label)**2/2 + \
                   self.penalty *( (self.WBh1**2).sum() + (self.WBh2**2).sum() + (self.WBo**2).sum() )
        
        else: # no penalty case
            
            # Update the weights
            self.WBh1 -= self.lr * dWBh1
            self.WBh2 -= self.lr * dWBh2
            self.WBo -= self.lr * dWBo
            
            # Compute penalty error
            loss = (Y - label)**2/2     
            
        return loss
    
    def train(self, x_train, y_train, x_val, y_val, train_log=False, verbose=False):
    
        if (train_log == True) or (self.early_stopping == True):
            train_loss_log = []
            val_loss_log = []
        if (self.early_stopping == True):
            last_check = 40
        for num_ep in range(self.n_epochs):
            # Learning rate decay
            if self.en_decay:
                self.lr *= self.lr_decay
            # Train single epoch (sample by sample, no batch for now)
            train_loss_vec = [self.update(x, y) for x, y in zip(x_train, y_train)]
            avg_train_loss = np.mean(train_loss_vec)
            # Validation network
            y_val_est = np.array([self.forward(x) for x in x_val])
            avg_val_loss = np.mean((y_val_est - y_val)**2/2) # just mean squared error
            # Log
            train_loss_log.append(avg_train_loss)
            val_loss_log.append(avg_val_loss)
            
            if verbose == True:
                print('Epoch %d - lr: %.5f - Train loss: %.4f - Test loss: %.4f' % \
                      (num_ep + 1, self.lr, avg_train_loss, avg_val_loss))
            
            if self.early_stopping and (num_ep > 100):
                # if the minimum loss of the last 20 epochs is greater than the mean of the previous 20 epochs
                # in a range of tolerance specified, then either stop if the learning is adaptive
                # or adapte the learning rate untill it is smaller or equal to 10^-4
                if np.mean(val_loss_log[-40:-20]) - self.tolerance < np.min(val_loss_log[-20:]):
                    if (not self.en_decay) and (self.lr >= 1e-4): #only case in which we do not stop
                        if num_ep - last_check > 20: #wait at least 20 epochs between adjustment of learning rates
                            if verbose:
                                print("Training is not improving. Reducing by 10 the learning rate.")
                            self.lr = self.lr/10 
                            last_check = num_ep
                        else:
                            continue
                    else:
                        if (not self.en_decay) and num_ep - last_check < 20:
                            continue
                        else:
                            if verbose:
                                print("Training interrupted by early stopping.")
                            break 
                else:
                    continue
            
        if train_log == True:
            return train_loss_log, val_loss_log
    
    
    def evaluate_mean_loss(self, x_test, y_test, save=False):
        y_predicted = np.array([self.forward(x) for x in x_test])
        average_loss = np.mean((y_predicted - y_test)**2/2)
        if save:
            np.savetxt('predictions.txt', y_predicted)
        return average_loss
        
    def plot_weights(self):
    
        fig, axs = plt.subplots(3, 1, figsize=(12,9))
        axs[0].hist(self.WBh1.flatten(), 20)
        axs[0].set_title('First hidden layer weights')
        axs[1].hist(self.WBh2.flatten(), 50)
        axs[1].set_title('Second hidden layer weights')
        axs[2].hist(self.WBo.flatten(), 20)
        axs[2].set_title('Output layer weights')
        #plt.legend()
        #plt.grid()
        plt.tight_layout()
        plt.show()<jupyter_output><empty_output><jupyter_text>L1 penalty termIn the initialization if the flag `en_penalty` is set to `'L1'`, then the L1 penalty is used.New loss function:$$ L = \frac{1}{2}(Y-\hat{Y})^2 + \lambda \left(\sum_{i,j}|W_{ij}^{(o)}| + \sum_{i,j}|W_{ij}^{(h1)}| + \sum_{i,j}|W_{ij}^{(h2)}| \right)$$Thus there is an additional term in the backpropagation. For example:$$\frac{\partial L}{\partial W_{ij}^{(o)}} = (Y - \hat{Y}) + \lambda \cdot sgn(W_{ij}^{(o)}) $$  L2 penalty termIn the initialization if the flag `en_penalty` is set to `'L2'`, then the L2 penalty is used.New loss function:$$ L = \frac{1}{2}(Y-\hat{Y})^2 + \lambda \left(\sum_{i,j}(W_{ij}^{(o)})^2 + \sum_{i,j}(W_{ij}^{(h1)})^2 + \sum_{i,j}(W_{ij}^{(h2)})^2 \right)$$Thus there is an additional term in the backpropagation. For example:$$\frac{\partial L}{\partial W_{ij}^{(o)}} = (Y - \hat{Y}) + 2\lambda \cdot W_{ij}^{(o)}$$  Gradient clippingTo avoid exploding gradients, if a gradient is greater in absolute value than a certain threshold `grad_tresh`, clip it to the treshold.Gradient clipping has a severe issue, that is the fact that it is altering the learning rule. This has the positive effect to avoid exploding gradients, but it doesn't hold anymore the convexity property of the loss function, hence for some values of the gradient treshold the training loss oscillates a lot.<jupyter_code># example
grad = np.array([-1,-1000,2000,1,2,3])
print("Initial grad: \n", grad)
# set a treshold
grad_tresh = 10

grad_norm = np.linalg.norm(grad)

if grad_norm > grad_tresh:
    grad_clipped = grad/grad_norm*grad_tresh
print("Clipped grad: \n",grad_clipped)<jupyter_output>Initial grad: 
 [   -1 -1000  2000     1     2     3]
Clipped grad: 
 [-4.47212925e-03 -4.47212925e+00  8.94425849e+00  4.47212925e-03
  8.94425849e-03  1.34163877e-02]<jupyter_text>Initialize new NN class passing parameters<jupyter_code># just to see if it works
x_train = x_training
y_train = y_training
x_val = x_test
y_val = y_test
# single choice of parameters
parameters = dict(N_neurons = [1,50,20,1], 
                  act_func = LeakyReLU,
                  act_der_func = LeakyReLU_der, 
                  lr = 0.005,
                  n_epochs = 2000,
                  en_decay = True,
                  lr_final = 1e-4,
                  early_stopping = True,
                  tol = 1e-3,
                  en_penalty = 'L1',
                  penalty = 0.0005,
                  en_grad_clipping = True,
                  grad_treshold = 10)
# pass a dictionary as keyword=argument list to the __init__ function
net = Network(**parameters)
train_loss_log, val_loss_log = net.train(x_train, y_train, x_val, y_val, train_log=True, verbose=True)
plt.plot(np.arange(len(train_loss_log)), train_loss_log)
plt.plot(np.arange(len(val_loss_log[10:])), val_loss_log[10:])
plt.xlabel('Number of epochs', fontsize = 15)
plt.ylabel('MSE validation loss', fontsize = 15)
plt.title('Early stopping demonstration', fontsize = 16)
plt.show()
show_forecast_vs_groundtruth(net,x_train, y_train)<jupyter_output><empty_output><jupyter_text>Some manual testing to cut down options without grid search and to tune the early stopping algorithm.Also we can see that the architecture is not influencing that much the performances, whereas I was using a penalty too high.```python starting configurationparameters = dict(N_neurons = [1,50,5,1],                   act_func = LeakyReLU,                  act_der_func = LeakyReLU_der,                   lr = 0.005,                  n_epochs = 2000,                  en_decay = True,                  lr_final = 1e-4,                  early_stopping = True,                  tol = 1e-4,                  en_penalty = 'L2',                  penalty = 0.05,                  en_grad_clipping = True,                  grad_treshold = 10)``` Changes in number of neuronsN_neurons = [1,50,5,1]Epoch 1000 - lr: 0.00071 - Train loss: 2.4792 - Test loss: 1.6337 Epoch 2000 - lr: 0.00010 - Train loss: 2.3872 - Test loss: 1.5329N_neurons = [1,50,10,1]Epoch 1000 - lr: 0.00071 - Train loss: 2.4803 - Test loss: 1.6321 Epoch 2000 - lr: 0.00010 - Train loss: 2.3891 - Test loss: 1.5353N_neurons = [1,50,20,1]Epoch 1000 - lr: 0.00071 - Train loss: 2.4802 - Test loss: 1.6322 Epoch 2000 - lr: 0.00010 - Train loss: 2.3889 - Test loss: 1.5351N_neurons = [1,20,20,1]Epoch 1000 - lr: 0.00071 - Train loss: 2.4802 - Test loss: 1.6330 Epoch 2000 - lr: 0.00010 - Train loss: 2.3886 - Test loss: 1.5347N_neurons = [1,100,20,1]Epoch 1000 - lr: 0.00071 - Train loss: 2.4794 - Test loss: 1.6338 Epoch 2000 - lr: 0.00010 - Train loss: 2.3876 - Test loss: 1.5333 Changes in the penalty termN_neurons = [1,100,20,1] en_penalty = 'L2' penalty = 0.1                   Epoch 1000 - lr: 0.00071 - Train loss: 3.7074 - Test loss: 1.8263 Epoch 2000 - lr: 0.00010 - Train loss: 3.5992 - Test loss: 1.6757N_neurons = [1,100,20,1] en_penalty = 'L1' penalty = 0.01 Epoch 1000 - lr: 0.00071 - Train loss: 1.1644 - Test loss: 1.2625 Epoch 1440 - lr: 0.00030 - Train loss: 1.1430 - Test loss: 1.2519 Training interrupted by early stopping.N_neurons = [1,100,20,1] en_penalty = 'L2' penalty = 0.01 Epoch 239 - lr: 0.00313 - Train loss: 1.3674 - Test loss: 1.5084Training interrupted by early stopping. Finally found a good range for the penalty termN_neurons = [1,100,20,1] en_penalty = 'L2' penalty = 0.001 Epoch 1000 - lr: 0.00071 - Train loss: 0.7452 - Test loss: 0.7281 Epoch 2000 - lr: 0.00010 - Train loss: 0.6805 - Test loss: 0.6667```pythonparameters = dict(N_neurons = [1,100,20,1],                   act_func = sigmoid,                  act_der_func = sigmoid_der,                   lr = 0.005,                  n_epochs = 2000,                  en_decay = True,                  lr_final = 1e-4,                  early_stopping = True,                  tol = 1e-4,                  en_penalty = 'L1',                  penalty = 0.0005,                  en_grad_clipping = True,                  grad_treshold = 10)```Epoch 1000 - lr: 0.00071 - Train loss: 0.9419 - Test loss: 0.7879 Epoch 1612 - lr: 0.00021 - Train loss: 0.9186 - Test loss: 0.7788 Training interrupted by early stopping.  K-fold cross validationIdea:- split training data in K subsets- use K-1 to train and 1 to validate (like a test set)- cycle over all K subset and store the validation losses- rank the hyperparameters used as the mean validation loss<jupyter_code>def get_K_folds(x, y, K=5):
    """
    Splits training set x with associated labels y in K (default 5) subsets.
    
    Parameters
    ----------
    x : numpy array, training set 
    y : numpy array, training labels
    K : int, number of folds (optional)
    
    Returns
    -------
    x_subsets : list of K numpy arrays
        Contains the K training subsets
    y_subsets : list of K numpy arrays
        Contains the labels of the K subsets
    """
    m = x.shape[0]
    permutation = np.random.permutation(m)
    x_shuffled = x[permutation]
    y_shuffled = y[permutation]
    
    subset_length = int(m/K)
    x_subsets = []
    y_subsets = []
    for i in range(subset_length-1):
        x_subsets.append(x_shuffled[i*subset_length: (i+1)*subset_length])
        y_subsets.append(y_shuffled[i*subset_length: (i+1)*subset_length])
    
    x_subsets.append(x_shuffled[(K-1)*subset_length:])
    y_subsets.append(x_shuffled[(K-1)*subset_length:])
    
    return x_subsets, y_subsets
x_subsets, y_subsets = get_K_folds(x_training, y_training)
def evaluate_HP(model, params, x_training, y_training, K):
    """
    Train and evaluates a model given some parameters using a K-fold cross validation.
    
    Parameters
    ----------
    model      : class of the predictor
        Must be instantiated as model(**params)
        Must have a "train" method like:
        net.train(x_train, y_train, x_val, y_val, train_log=True, verbose=False)
    params     : dictionary {'model_HP':value}
        hyper-parameters of the model
    x_training : numpy array, training set 
    y_training : training labels
    K          : int, number of folds
    
    """
    Kfolds_val_losses = []
    epochs_length = []
    x_subsets, y_subsets = get_K_folds(x_training, y_training, K)
    
            
    for i in range(K):
        print("Computing fold %d out of %d..."%(i+1,K))
        
        # instantiate new network
        net = model(**params)
        
        # validation sets for this fold
        x_val = x_subsets[i]
        y_val = y_subsets[i]
        
        # training sets
        x_train = np.array([])
        y_train = np.array([])
        for j in range(K):
            if j != i:
                if len(x_train) == 0:
                    x_train = x_subsets[j]
                    y_train = y_subsets[j]
                else:
                    x_train = np.concatenate((x_train, x_subsets[j]), axis=0)
                    y_train = np.concatenate((y_train, y_subsets[j]), axis=0)
        
        
        _, val_loss_log = net.train(x_train, y_train, x_val, y_val, train_log=True, verbose=False)
        # either choose the last element or the smaller (could be a good choice if early stopping is implemented)
        Kfolds_val_losses.append(val_loss_log[-1])
        epochs_length.append(len(val_loss_log)) #record the real number of training epochs

    HP_score = np.array(Kfolds_val_losses).mean()
    HP_deviation = np.array(Kfolds_val_losses).std()
    mean_epochs = int(np.mean(epochs_length))
    return HP_score, HP_deviation, mean_epochs
# single choice of parameters
parameters = dict(N_neurons = [1,20,20,1], 
                  act_func = LeakyReLU,
                  act_der_func = LeakyReLU_der, 
                  lr = 0.01,
                  n_epochs = 50,
                  en_decay = True,
                  lr_final = 0.001,
                  en_penalty = 'L1',
                  penalty = 1e-2,
                  early_stopping = True)
score, dev, mean_epochs = evaluate_HP(Network, parameters,x_training, y_training, K=5)
print(score, dev)
def print_parameters(params):
    print("Parameters: ")
    print('='*75)
    for key in params:
        if (key == 'lr') or (key == 'tol'):
            print(key, '\t \t', params[key])
        else:
            print(key, '\t', params[key])
    print('='*75)

def print_HP_score(params,score,dev):
    print_parameters(params)
    print("Loss: %.4f +/- %.4f"%(score,dev))
    

print_HP_score(parameters, score, dev)<jupyter_output>Parameters: 
===========================================================================
N_neurons 	 [1, 20, 20, 1]
act_func 	 <function LeakyReLU at 0x7febe5b414d0>
act_der_func 	 <function LeakyReLU_der at 0x7febe5b415f0>
lr 	 	 0.01
n_epochs 	 50
en_decay 	 True
lr_final 	 0.001
en_penalty 	 L1
penalty 	 0.01
early_stopping 	 True
===========================================================================
Loss: 1.3229 +/- 0.4566<jupyter_text>Grid Search Cross Validation<jupyter_code>def GridSearchCV(model, x_training, y_training, params_comb, K_folds=5):
    """
    Parameters
    ----------
    model      : class of the predictor
        Must be instantiated as model(**params)
        Must have a "train" method like:
        net.train(x_train, y_train, x_val, y_val, train_log, verbose)
    params_comb: list of dictionaries {'model_HP':value}
        all combinations of hyper-parameters to be evaluated
    x_training : numpy array, training set 
    y_training : training labels
    K_folds    : int, number of folds
    
    Returns
    -------
    scoring    : list of all the mean cross-validation losses 
    deviations : list of all the std deviations of the cross-validation losses 
    mean_epochs: list of the mean number of epochs of training for each 
                 configuration of HPs
    """
    import time
    scoring = []
    deviations = []
    mean_epochs = []
    combinations = len(params_comb)
    init_time = time.time()
    for i in range(combinations):
        start = time.time()
        params = params_comb[i]
        HP_score, HP_deviation, epochs = evaluate_HP(model, params, x_training, y_training, K=K_folds)
        scoring.append(HP_score)
        deviations.append(HP_deviation)
        mean_epochs.append(epochs)
        
        print('\nParameters configuration %d out of %d'%(i+1,len(params_comb)))
        print_HP_score(params,HP_score,HP_deviation)
        finish = time.time()
        print('Time for evaluation: %.1f s'%(finish-start))
        time_to_finish = (combinations-i-1)*(finish-init_time)/(i+1)
        
        if time_to_finish > 3600:
            print('Estimated time to finish : %.2f h'%(time_to_finish/3600), '\n')
        elif time_to_finish < 3600 and time_to_finish > 60:
            print('Estimated time to finish : %.2f min'%(time_to_finish/60), '\n')
        else:
            print('Estimated time to finish : %.2f s'%(time_to_finish), '\n')
    return scoring, deviations, mean_epochs<jupyter_output><empty_output><jupyter_text>Now we need to explore some portions of the parameter space and then rank the parameters configurations in order of scoring.Possible/reasonable options for all parameters:- 3 activation functions (with respective derivatives);- 1 learning rate: 0.005 -> 0.0001;- 4 possible numbers of neurons in the first hidden layer, e.g. 10, 20, 50, 100;- 3 possible numbers of neurons in the second hidden layer, e.g. 10, 20, 50 (but I just want the configurations in which Nh2 <= Nh1, which makes 9 total configs);- 4 combinations with penalties: L1 with 0.001, L1 with 0.0005, L2 with 0.001, L2 with 0.0005;Early stopping will be always on. This will enable us to use a great number of epochs (2000) without overfitting.Total number of combinations: $3*9*4 = 108$  I am not sure if exists some library to further automate the enumeration of all the possible parameters, but not having found one, I built my own implementation of that. The goal is to have a list of dictionaries, each of them containing all the keywords needed to specify the model initialization.<jupyter_code>import itertools
# we start making list of dictionaries of all the possibilities for every group of keywords

act_func_dict = [{'act_func':sigmoid,'act_der_func':sigmoid_der},
                  {'act_func':ReLU,'act_der_func':ReLU_der},
                  {'act_func':LeakyReLU,'act_der_func':LeakyReLU_der}]

lr_dict = [{'lr': 0.005, 'en_decay':True,'lr_final':0.0001}]

n_epochs_dict = [{'n_epochs':2000, 'early_stopping':True, 'tol':1e-3}]

regularization_dict = [{'en_penalty':'L1', 'penalty':0.001},
                       {'en_penalty':'L1', 'penalty':0.0005},
                       {'en_penalty':'L2', 'penalty':0.001},
                       {'en_penalty':'L2', 'penalty':0.0005}]

N_neurons_dict =[{'N_neurons': [1, 10, 10, 1]},
                 {'N_neurons': [1, 20, 10, 1]},
                 {'N_neurons': [1, 20, 20, 1]},
                 {'N_neurons': [1, 50, 10, 1]},
                 {'N_neurons': [1, 50, 20, 1]},
                 {'N_neurons': [1, 50, 50, 1]},
                 {'N_neurons': [1, 100, 10, 1]},
                 {'N_neurons': [1, 100, 20, 1]},
                 {'N_neurons': [1, 100, 50, 1]}]<jupyter_output><empty_output><jupyter_text>The important fact here is that some key-values pairs must be tried only together with other pairs, thus we must treat them as a single case. It's difficult and tedious to explain the code below here, I will just say that it works and there is the possibility to read all the combinations formatted.<jupyter_code>def combine_dictionaries(*dictionaries):
    # *var is used to pass an arbitrary number of positional variables 
    global_dict = list(itertools.product(*dictionaries))
    params_list = []
    for i in range(len(global_dict )):
        merged_dict = {}
        for d in global_dict [i]:
            merged_dict = dict(merged_dict, **d)
        params_list.append(merged_dict)
    return params_list
params_list = combine_dictionaries(N_neurons_dict, act_func_dict, lr_dict, n_epochs_dict, regularization_dict)
print("Number of combinations: ", len(params_list))

flag = True #set to True to see all combinations
if flag == True:
    for params in params_list:
        print()
        print_parameters(params)<jupyter_output>Number of combinations:  108

Parameters: 
===========================================================================
N_neurons 	 [1, 10, 10, 1]
act_func 	 <function sigmoid at 0x7feadaa23710>
act_der_func 	 <function sigmoid_der at 0x7feadaa28560>
lr 	 	 0.005
en_decay 	 True
lr_final 	 0.0001
n_epochs 	 2000
early_stopping 	 True
tol 	 	 0.001
en_penalty 	 L1
penalty 	 0.001
===========================================================================

Parameters: 
===========================================================================
N_neurons 	 [1, 10, 10, 1]
act_func 	 <function sigmoid at 0x7feadaa23710>
act_der_func 	 <function sigmoid_der at 0x7feadaa28560>
lr 	 	 0.005
en_decay 	 True
lr_final 	 0.0001
n_epochs 	 2000
early_stopping 	 True
tol 	 	 0.001
en_penalty 	 L1
penalty 	 0.0005
===========================================================================

Parameters: 
===========================================================================
N_neurons 	 [1, 10, 10, [...]<jupyter_text>Using grid search to study single parametersBefore selecting the final parameters I want to shrink their numbers by testing them keeping all other variables fixed. This is useful to make heuristics more evident.In particular I want to see what changes changing architecture and changing the penalty weight.<jupyter_code># we start making list of dictionaries of all the possibilities for every group of keywords

act_func_dict = [{'act_func':sigmoid,'act_der_func':sigmoid_der},
                  {'act_func':ReLU,'act_der_func':ReLU_der},
                  {'act_func':LeakyReLU,'act_der_func':LeakyReLU_der}]
lr_dict = [{'lr': 0.005, 'en_decay':True,'lr_final':0.0001}]

n_epochs_dict = [{'n_epochs':1000, 'early_stopping':True, 'tol':1e-3}]

regularization_dict = [{'en_penalty':'L2', 'penalty':0.1},
                       {'en_penalty':'L2', 'penalty':0.05},
                       {'en_penalty':'L2', 'penalty':0.01},
                       {'en_penalty':'L2', 'penalty':0.005},
                       {'en_penalty':'L2', 'penalty':0.001},
                       {'en_penalty':'L2', 'penalty':0.0005}]

N_neurons_dict =[{'N_neurons': [1, 10, 10, 1]},
                 {'N_neurons': [1, 20, 10, 1]},
                 {'N_neurons': [1, 20, 20, 1]},
                 {'N_neurons': [1, 50, 10, 1]},
                 {'N_neurons': [1, 50, 20, 1]},
                 {'N_neurons': [1, 50, 50, 1]},
                 {'N_neurons': [1, 100, 10, 1]},
                 {'N_neurons': [1, 100, 20, 1]},
                 {'N_neurons': [1, 100, 50, 1]}]<jupyter_output><empty_output><jupyter_text>Regularization term<jupyter_code>regularizations_list = combine_dictionaries(N_neurons_dict[2:3], act_func_dict[:1], lr_dict, n_epochs_dict, regularization_dict)
print("Number of combinations: ", len(regularizations_list))

flag = True #set to True to see all combinations
if flag == True:
    for params in regularizations_list:
        print()
        print_parameters(params)
%%time
reg_scoring, reg_deviations, reg_mean_epochs = GridSearchCV(Network, x_training, y_training, regularizations_list, K_folds=5)
penalties = np.array([0.1,0.05,0.01,0.005,0.001,0.0005])
plt.plot(penalties,reg_scoring,'-o', label = 'Sigmoid \nNh1=20 \nNh2=20')
plt.xscale('log')
plt.xlabel('Penalty term in logscale', fontsize=15)
plt.ylabel('MSE validation loss', fontsize=15)
plt.title('Effect of regularization on accuracy', fontsize=16)
plt.legend(fontsize=13)
plt.show()<jupyter_output><empty_output><jupyter_text>Number of neurons<jupyter_code>neurons_list = combine_dictionaries(N_neurons_dict, act_func_dict[:1], lr_dict, n_epochs_dict, regularization_dict[-1:])
print("Number of combinations: ", len(neurons_list))

flag = True #set to True to see all combinations
if flag == True:
    for params in neurons_list:
        print()
        print_parameters(params)
%%time
neurons_scoring, neurons_deviations, neurons_mean_epochs = GridSearchCV(Network, x_training, y_training, neurons_list, K_folds=5)
#np.save('neurons_scoring', neurons_scoring) # to make sure of not losing the results of the grid search
neurons_scoring = np.load('neurons_scoring.npy')
#! pip install seaborn
import seaborn as sns
score_matrix = np.zeros((3,4))
index_dict = {10:0, 20:1, 50:2, 100:3}
for i in range(len(N_neurons_dict)):
    Nh1 = N_neurons_dict[i]['N_neurons'][1]
    Nh2 = N_neurons_dict[i]['N_neurons'][2]
    index1 = index_dict[Nh1]
    index2 = index_dict[Nh2]
    score_matrix[index2,index1] = neurons_scoring[i]
help(sns.heatmap)
score_matrix
help(ax.set_ylim)
help(sns.heatmap)
plt.figure(figsize = (10,7))

ax = sns.heatmap(score_matrix, annot = True, annot_kws={'fontsize':18},
                 cbar_kws = {'label' : 'MSE validation loss'})

ax.set_ylim(bottom=3, top=0)

x_index = ['%d'%n for n in [10,20,50,100]]
y_index = ['%d'%n for n in [10,20,50]]

fontdict=dict(fontsize=18)

ax.set_xticklabels(x_index, rotation = 0, **fontdict)
ax.set_xlabel('Number of neurons Nh1', fontsize = 18)

ax.set_yticklabels(y_index, rotation = 0, **fontdict)
ax.set_ylabel('Number of neurons Nh2', fontsize = 18)


ax.figure.axes[-1].yaxis.label.set_size(18)
plt.title("Accuracy dependence on the architecture", fontsize = 20)
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Final grid search<jupyter_code># we start making list of dictionaries of all the possibilities for every group of keywords

act_func_dict = [{'act_func':sigmoid,'act_der_func':sigmoid_der},
                  {'act_func':ReLU,'act_der_func':ReLU_der},
                  {'act_func':LeakyReLU,'act_der_func':LeakyReLU_der}]

lr_dict = [{'lr': 0.005, 'en_decay':True,'lr_final':0.0001}]

n_epochs_dict = [{'n_epochs':2000, 'early_stopping':True, 'tol':1e-3}]

regularization_dict = [{'en_penalty':'L1', 'penalty':0.001},
                       {'en_penalty':'L1', 'penalty':0.0005},
                       {'en_penalty':'L2', 'penalty':0.001},
                       {'en_penalty':'L2', 'penalty':0.0005}]

N_neurons_dict =[{'N_neurons': [1, 10, 10, 1]},
                 {'N_neurons': [1, 20, 10, 1]},
                 {'N_neurons': [1, 20, 20, 1]},
                 {'N_neurons': [1, 50, 10, 1]},
                 {'N_neurons': [1, 50, 20, 1]},
                 {'N_neurons': [1, 50, 50, 1]},
                 {'N_neurons': [1, 100, 10, 1]},
                 {'N_neurons': [1, 100, 20, 1]},
                 {'N_neurons': [1, 100, 50, 1]}]
params_list = combine_dictionaries(N_neurons_dict, act_func_dict, lr_dict, n_epochs_dict, regularization_dict)
print("Number of combinations: ", len(params_list))

flag = True #set to True to see all combinations
if flag == True:
    for params in params_list:
        print()
        print_parameters(params)
%%time
scoring, deviations, mean_epochs = GridSearchCV(Network, x_training, y_training, params_list, K_folds=5)<jupyter_output>Computing fold 1 out of 5...
Computing fold 2 out of 5...
Computing fold 3 out of 5...
Computing fold 4 out of 5...
Computing fold 5 out of 5...

Parameters configuration 1 out of 108
Parameters: 
===========================================================================
N_neurons 	 [1, 10, 10, 1]
act_func 	 <function sigmoid at 0x7febf1acbb90>
act_der_func 	 <function sigmoid_der at 0x7febf1ad14d0>
lr 	 	 0.005
en_decay 	 True
lr_final 	 0.0001
n_epochs 	 2000
early_stopping 	 True
tol 	 	 0.001
en_penalty 	 L1
penalty 	 0.001
===========================================================================
Loss: 0.9533 +/- 0.2689
Time for evaluation: 138.1 s
Estimated time to finish : 4.10 h 

Computing fold 1 out of 5...
Computing fold 2 out of 5...
Computing fold 3 out of 5...
Computing fold 4 out of 5...
Computing fold 5 out of 5...

Parameters configuration 2 out of 108
Parameters: 
===========================================================================
N_neurons 	 [1, 10, 10, 1]
[...]<jupyter_text>For 16 combinations CPU times: user 34min 25s, sys: 1.38 s, total: 34min 27s Wall time: 34min 30sFor 108 combinations Expected time 3h 53 min<jupyter_code>np.save('HP_scoring', scoring) # to make sure of not losing the results of the grid search
np.save('HP_deviations', deviations) 
np.save('mean_epochs', mean_epochs) 
scoring = np.load('HP_scoring.npy')
deviations = np.load('HP_deviations.npy')
mean_epochs = np.load('mean_epochs.npy')
# if there are some trainings in which the gradient exploded, take care of them
nan_mask = np.isnan(scoring)
print('Worst finite result: ', np.max(scoring[~nan_mask]))
nan_substitute = np.max(scoring[~nan_mask])
scoring[nan_mask] = nan_substitute
# my metric to choose the best model
composite_loss = scoring + 2*deviations
best_params = params_list[np.argmin(composite_loss)]
best_n_epochs = mean_epochs[np.argmin(composite_loss)]
print_HP_score(score=scoring[np.argmin(composite_loss)], 
               dev = deviations[np.argmin(composite_loss)], 
               params=best_params)
# if the early stopping was used, in the final training don't use it 
#but train for the mean number of epochs used in the cross validation
if 'early_stopping' in best_params.keys():
    best_params['early_stopping'] = False
    best_params['n_epochs'] = best_n_epochs
print_HP_score(score=scoring[np.argmin(composite_loss)], 
               dev = deviations[np.argmin(composite_loss)], 
               params=best_params)<jupyter_output>Parameters: 
===========================================================================
N_neurons 	 [1, 20, 10, 1]
act_func 	 <function ReLU at 0x7febe5b41290>
act_der_func 	 <function ReLU_der at 0x7febe5b413b0>
lr 	 	 0.005
en_decay 	 True
lr_final 	 0.0001
n_epochs 	 599
early_stopping 	 False
tol 	 	 0.001
en_penalty 	 L2
penalty 	 0.001
===========================================================================
Loss: 0.8550 +/- 0.0770<jupyter_text>Final training of the best model<jupyter_code>best_model = Network(**best_params)
train_loss_log, test_loss_log = best_model.train(x_training, y_training, x_test, y_test, train_log=True, verbose=True)
plt.plot(np.arange(len(train_loss_log)), train_loss_log, label ='training loss')
plt.plot(np.arange(len(test_loss_log)), test_loss_log, label ='test loss')
plt.legend(fontsize=13)
plt.xlabel('Number of epochs', fontsize = 16)
plt.ylabel('MSE loss', fontsize = 16)
plt.title('Final model training')
plt.show()
best_model.plot_weights()
def show_forecast_vs_groundtruth(NN, x_test, y_test):
    x_axis = np.linspace(-4,4,4000)
    plt.plot(x_test, [NN.forward(x) for x in x_test], 'b.', label = 'predictions')
    plt.plot(x_test, y_test, '.r', label = 'ground truth')
    plt.plot(x_axis, [NN.forward(x) for x in x_axis], 'g--', lw = 0.7, label =  'learned function' )
    plt.xlabel('x axis', fontsize=16)
    plt.ylabel('y axis', fontsize=16)
    plt.legend(fontsize=13)
    plt.show()
show_forecast_vs_groundtruth(best_model,x_test, y_test)
show_forecast_vs_groundtruth(best_model,x_train, y_train)
import pickle
# Save parameters
with open("dict.txt", "wb") as file: 
    pickle.dump(best_params, file)
# Save weights
best_model.save_weights()<jupyter_output><empty_output><jupyter_text>Other finalists<jupyter_code>import pandas as pd
params_list = np.array(params_list)
params_list = params_list.reshape((-1,1))
scoring = scoring.reshape((-1,1))
deviations = deviations.reshape((-1,1))
composite_loss = scoring + 2*deviations
data = np.stack((params_list,scoring,deviations,composite_loss), axis=1).reshape((-1,4))
df = pd.DataFrame(data=data, columns=['parameters','loss','dev','score'])
df_best = df.sort_values(by='score', ascending=True)[:5] #best 5 functions
best_configs = df_best.to_numpy()
for i in range(len(best_configs)):
    print_HP_score(params=best_configs[i,0], score=best_configs[i,1], dev=best_configs[i,2])
    print()<jupyter_output>Parameters: 
===========================================================================
N_neurons 	 [1, 20, 10, 1]
act_func 	 <function ReLU at 0x7febe5b41290>
act_der_func 	 <function ReLU_der at 0x7febe5b413b0>
lr 	 	 0.005
en_decay 	 True
lr_final 	 0.0001
n_epochs 	 599
early_stopping 	 False
tol 	 	 0.001
en_penalty 	 L2
penalty 	 0.001
===========================================================================
Loss: 0.8550 +/- 0.0770

Parameters: 
===========================================================================
N_neurons 	 [1, 100, 50, 1]
act_func 	 <function ReLU at 0x7febe5b41290>
act_der_func 	 <function ReLU_der at 0x7febe5b413b0>
lr 	 	 0.005
en_decay 	 True
lr_final 	 0.0001
n_epochs 	 2000
early_stopping 	 True
tol 	 	 0.001
en_penalty 	 L2
penalty 	 0.001
===========================================================================
Loss: 0.8259 +/- 0.1020

Parameters: 
===========================================================================
N_neurons 	 [1, 50[...]
{'fraction_non_alphanumeric': 0.14489277317310661, 'fraction_numerical': 0.04536062887405894, 'mean_word_length': 2.7112638815441565, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 95, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '515787', 'license': 'MIT', 'n_tokens_mistral': 15807, 'n_tokens_neox': 14018, 'n_words': 3635, 'path': 'Nicola_Dainese_Ex1/Source code/Nicola_Dainese_Lab2.ipynb', 'repo_name': 'nicoladainese96/NeuralNetworks_DeepLearning', 'size': '42107'}