# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Artificial Neural Networks

from IPython.display import Image
Image(filename='Image2.png',width = "600",height = "300")

# 1. Here, each node on the left represents an input feature, the connecting lines represent
# the learned coefficients, and the node on the right represents the output, which is a
# weighted sum of the inputs.
#
# 2. This model has a lot more coefficients (also called weights) to learn: there is one
# between every input and every hidden unit (which make up the hidden layer), and
# one between every unit in the hidden layer and the output.
#
# 3. Computing a series of weighted sums is mathematically the same as computing just
# one weighted sum, so to make this model truly more powerful than a linear model,
# we need one extra trick. After computing a weighted sum for each hidden unit, a
# nonlinear function is applied to the result—usually the rectifying nonlinearity (also
# known as rectified linear unit or relu) or the tangens hyperbolicus (tanh). The result of
# this function is then used in the weighted sum that computes the output, ŷ. The two
# functions are visualized in the following Figure. The relu cuts off values below zero, while tanh
# saturates to –1 for low input values and +1 for high input values. Either nonlinear
# function allows the neural network to learn much more complicated functions than a
# linear model could:

from IPython.display import Image
Image(filename='Image3.png',width = "400",height = "300")

# For the small neural network pictured in the first figure, the full formula for computing ŷ in the case of regression would be (when using a tanh nonlinearity):   
# h[0] = tanh(w[0, 0] * x[0] + w[1, 0] * x[1] + w[2, 0] * x[2] + w[3, 0] * x[3] + b[0])   
# h[1] = tanh(w[0, 1] * x[0] + w[1, 1] * x[1] + w[2, 1] * x[2] + w[3, 1] * x[3] + b[1])   
# h[2] = tanh(w[0, 2] * x[0] + w[1, 2] * x[1] + w[2, 2] * x[2] + w[3, 2] * x[3] + b[2])   
# ŷ = v[0] * h[0] + v[1] * h[1] + v[2] * h[2] + b

# Here, **w** are the weights between the input **x** and the hidden layer **h**, and **v** are the
# weights between the hidden layer **h** and the output ŷ. The weights **v** and **w** are learned
# from data, **x** are the input features, ŷ is the computed output, and **h** are intermediate
# computations. An important parameter that needs to be set by the user is the number
# of nodes in the hidden layer. This can be as small as 10 for very small or simple datasets
# and as big as 10,000 for very complex data.

# Having large neural networks made up of many of these layers of computation is
# what inspired the term “deep learning.”

# # Data used here

# Here, we use the data **Avila** from UCI: https://archive.ics.uci.edu/ml/datasets/Avila#. The Avila data set has been extracted from 800 images of the 'Avila Bible', an XII century giant Latin copy of the Bible. The prediction task consists in associating each pattern to a copyist. In this lab, the‘avila-tr.txt'data is used. For your convenience, I have pre-processed an 'CSV' file in the Github.

# ### Data Set Information:

# CLASS DISTRIBUTION (training set)
#
# A: 4286 
#
# B: 5 
#
# C: 103 
#
# D: 352 
#
# E: 1095 
#
# F: 1961 
#
# G: 446 
#
# H: 519 
#
# I: 831 
#
# W: 44 
#
# X: 522 
#
# Y: 266

# ### Attribute Information:

# F1: intercolumnar distance 
#
# F2: upper margin 
#
# F3: lower margin 
#
# F4: exploitation 
#
# F5: row number 
#
# F6: modular ratio 
#
# F7: interlinear spacing 
#
# F8: weight 
#
# F9: peak number 
#
# F10: modular ratio/ interlinear spacing 
#
# Class: A, B, C, D, E, F, G, H, I, W, X, Y

# ### Relevant Papers:

# <NAME>, <NAME>, <NAME>, <NAME> diÂ Freca, 
# Reliable writer identification in medieval manuscripts through page layout features: The 'Avila' Bible case, Engineering Applications of Artificial Intelligence, Volume 72, 2018, pp. 99-110. 
#
# <NAME>, <NAME>, <NAME> and <NAME>, 'A Method for Scribe Distinction in Medieval Manuscripts Using Page Layout Features', Lecture Notes in Computer Science, G. Maino and G. Foresti (eds.), Springer-Verlag, vol. 6978, pp. 393-402.

# ## Import modules

# +
import pandas as pd
import matplotlib.pyplot as plt

from scipy.stats import ttest_ind

from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score, KFold
from sklearn.metrics import accuracy_score
from sklearn.datasets import make_classification
from sklearn.naive_bayes import GaussianNB
import warnings
warnings.filterwarnings("ignore")
from sklearn.neural_network import MLPClassifier
from sklearn.naive_bayes import MultinomialNB
import warnings
warnings.filterwarnings("ignore")
# -

# ## Load the data

# Load the data and show the basic information
data=pd.read_csv("files/avila-tr.csv")
print('Data size: (%.f, %.f)\n' % data.shape)
data.head()

# ## Study the Multi-layer Perceptron model and its parameters
# We will use the multi-layer percepton to understand the basic usage of artificial neural network models. Class [MLPClassifier](https://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPClassifier.html#sklearn.neural_network.MLPClassifier) implements a multi-layer perceptron (MLP) algorithm that trains using Backpropagation. MLPClassifier has many parameters to configure. The effects of some important parameters have been discussed in our lecture, e.g., activation functions, number of hidden layers, number of units in a hidden layer, gradient descent, regularization, etc. So, in this workshop we would like to empirically study these parameters. 

# ### 1. Basic usage of the model

# +
# Specify features and the target
X = data.drop(['Class'], axis = 'columns')
y = data['Class']

# Split the dataset into training data and testing data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3, random_state=42)

# Training the MLPClassifier with the default parameters (random_state=42)
clf = MLPClassifier(random_state=42)
clf.fit(X_train, y_train)

# Output the accuracy on training data and test data respectively
y_pred = clf.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print('Accuracy: %.4f \n' % accuracy)
# -

# Explore the learned MLP model
print(clf)
print('\n# of layers (including the input layer): %.f\n' % clf.n_layers_)
print('MLP structure: %.f X %.f X %.f\n' % (X.shape[1], clf.get_params()['hidden_layer_sizes'][0], clf.n_outputs_))

# ### 2. Use 10-fold cross validation to report a more robust testing performance

# Use 10-fold cross validation to validate the model
clf = MLPClassifier(random_state=42)
scores_mlp_default = cross_val_score(clf, X, y, cv=10, verbose=1)
print('Accuracy range for Multi-layer Perceptron: [%.4f, %.4f]; mean: %.4f; std: %.4f\n'
      % (scores_mlp_default.min(), scores_mlp_default.max(), scores_mlp_default.mean(), scores_mlp_default.std()))

# ### 3. Increase the number of hidden units
# Here, we draw a figure to report the testing accuray with differnt number of hidden units [10,20,30,40,50,60,70,80,90,100]. The parameter 'hidden_layer_sizes' accepts a list of numbers specifying the number of units in each hidden layer. We plot the relationship between the parameter and the accuracy score.

# For each number of hidden units, we use 10-fold cross validation to report the testing accuracy.
# N.B.: The execution time is relatively long. This is to let you experience the intenstive computaiton required by artificial neural network models.
cv_scores = []
cv_scores_std = []
hidden_unit_numbers = [[10],[20],[30],[40],[50],[60],[70],[80],[90],[100]]
for i in hidden_unit_numbers:
    clf_mlp = MLPClassifier(hidden_layer_sizes=i, random_state=42)
    scores = cross_val_score(clf_mlp, X, y, scoring='accuracy', cv=10, verbose=1)
    cv_scores.append(scores.mean())
    cv_scores_std.append(scores.std())

# Plot the relationship
plt.errorbar(hidden_unit_numbers, cv_scores, yerr=cv_scores_std, marker='x', label='Accuracy')
plt.xlabel('Size of hidden units')
plt.ylim(0.5, 0.9)
plt.ylabel('Accuracy')
plt.legend(loc='best')
plt.show()

# **It can be seen that the accuracy increases when the number of units in the hidden layer increases. The reason is the a model with a bigger number of hiden layer units has a higher complexity to capture the information in data.**

# ### Task1. Number of hidder layers
# Set the number of hidden layers as two, and each has 100 units. Please report the accuracy score and compare it with the single hidden layer case mentioend above. Note that the execution time for this task is relative long. To track the progress of the execution, you could set the verbose parameter in the cross_val_score method as 'verbose=1'.

# Try a MLP model with two hidden layers
hidden_layer_sizes=[100,100]
clf_mlp=MLPClassifier(hidden_layer_sizes=hidden_layer_sizes, random_state=42)
scores_mlp_2layers=cross_val_score(clf_mlp, X, y, scoring='accuracy', cv=10, verbose=1)



# ### 4. Choose solver for the learning process
# As discussed in our lecture, we can use gradient descent methods (standard gradient and stochastic gradient descent) to learn weights for the error minimization problem. Moreover, we can have other solvers for the optimization problem. Here, we draw a figure to report the accuray with differnt solvers. The parameter 'solver' can take values from ['lbfgs','sgd', 'adam']. 'sgd' represents stochastic gradient descent.

# Try different solvers
cv_scores = []
cv_scores_std = []
solvers = ['lbfgs', 'sgd', 'adam']
for i in solvers:
    clf_mlp = MLPClassifier(solver=i, random_state=42)
    scores = cross_val_score(clf_mlp, X, y, scoring='accuracy', cv=10)
    cv_scores.append(scores.mean())
    cv_scores_std.append(scores.std())

# Plot the relationship
plt.bar(solvers, cv_scores, yerr=cv_scores_std, label='Accuracy')
plt.xlabel('Solvers')
plt.ylim([0.5, 1])
plt.ylabel('Accuracy')
plt.legend(loc='best')
plt.show()

# **It can be seen that sgd doesn't perform as well as the other two solvers**

# ### Task2. Activation functions
# It can be seen above that the default activation function is 'ReLU'. As we discussed in our lecture, we could have different types of activation functions. Please try the possible functions 'identity’, ‘logistic’, ‘tanh’, and ‘relu’ provided by the API. Draw a figure to report the accuray with differnt activation functions.

# Try different activation functions
cv_scores = []
cv_scores_std = []
functions = ['identity','logistic', 'tanh', 'relu']
for function in functions:
    clf_mlp = MLPClassifier(random_state=42, activation=function)
    scores = cross_val_score(clf_mlp, X, y, scoring='accuracy', cv=10)
    cv_scores.append(scores.mean())
    cv_scores_std.append(scores.std())

# Plot the relationship
plt.bar(functions, cv_scores, yerr=cv_scores_std, label='Accuracy')
plt.xlabel('Activation functions')
plt.ylim([0.5, 1])
plt.ylabel('Accuracy')
plt.legend(loc='best')
plt.show()

# ### 5. Use different values of alpha
# The parameter alpha is the L2 penalty (regularization term) to overcome the overfitting issue. It balance the error caused the data and that by the model structure (number of weights). 
# Here, please draw a figure to report the accuray with differnt values of alpha [0.0001,0.001,0.01, 0.1,1].

# Try different regulaization parameters
cv_scores = []
cv_scores_std = []
alphas = [0.0001,0.001,0.01, 0.1,1]
for i in alphas:
    clf_mlp = MLPClassifier(alpha=i, random_state=42)
    scores = cross_val_score(clf_mlp, X, y, scoring='accuracy', cv=10)
    cv_scores.append(scores.mean())
    cv_scores_std.append(scores.std())

# Plot the relationship
plt.errorbar(alphas, cv_scores, yerr=cv_scores_std, marker='x', label='Accuracy')
plt.xscale('log')
plt.xlabel('alpha')
plt.ylim([0.6, 0.85])
plt.ylabel('Accuracy')
plt.legend(loc='best')
plt.show()

# **It can be seen that when alpha increases, the prediction accuracy drops, showing that the overfitting issue of the model on this dataset is not a big problem.**

# ### Task3. Increase the number of iterations
# This is related to the early stopping technique we mentioned in the lecture. Please explore when is good enough to stop the iteration of weight updating. Please draw a figure to report the accuray with differnt number of iterations [200,400,600,800,1000]. The parameter 'max_iter' can specify this setting. Note that the execution time for this task is relative long. To track the progress of the execution, you could set the verbose parameter in the cross_val_score method as 'verbose=1'.

# Try different number of iterations. Plot the relationship between the performance and the number of iterations.
cv_scores = []
cv_scores_std = []
iterations = [200, 400, 600, 800, 1000]
for i in iterations:
    clf_mlp = MLPClassifier(max_iter=i, random_state=42)
    scores = cross_val_score(clf_mlp, X, y, scoring='accuracy', cv=10)
    cv_scores.append(scores.mean())
    cv_scores_std.append(scores.std())

plt.errorbar(iterations, cv_scores, yerr=cv_scores_std, marker='x', label='Accuracy')
plt.xscale('log')
plt.xlabel('alpha')
plt.ylim([0.2, 1])
plt.ylabel('Accuracy')
plt.legend(loc='best')
plt.show()

# ## Comparing with other classification models
# ### Task4. Compare the results with Naive Bayes (GaussianNB) and K-Nearest Neighbors (K=1)
# Compare their accuracy scores and use t test to show if their perofrmance has siginficantly different with significance level 0.05.

# Compare with the two models we used before
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB
clf_knn=KNeighborsClassifier(n_neighbors=1)
clf_gnb=GaussianNB()
scores_knn=cross_val_score(clf_knn, X, y, scoring='accuracy', cv=KFold(n_splits=10, shuffle=True, random_state=42))
scores_gnb=cross_val_score(clf_gnb, X, y, scoring='accuracy', cv=KFold(n_splits=10, shuffle=True, random_state=42))
print(scores_knn.mean(), scores_knn.std())
print(scores_gnb.mean(), scores_gnb.std())

from scipy.stats import ttest_ind
t, p=ttest_ind(scores_knn, cv_scores)
print(t,p) 
t, p=ttest_ind(scores_gnb, cv_scores)
print(t,p) 



StarcoderdataJupyterScriptsDedupFiltered