<jupyter_start><jupyter_text>Pre-processing, dimensionalization, visualization and image sorting<jupyter_code>import numpy as np
import pandas as pd
import random
import cv2
import matplotlib.pylab as plt
import os
from os import listdir, path
from PIL import Image as PImage # converting images into arrays
def load_images(path, width=100, height=100, inter=cv2.INTER_LANCZOS4) :
    """Read all images from a directory, converts them to all to the same dimensions,
    and return a DataFrame where each row corresponds to an image.
    
    Parameters
    ----------
    path (str): path (abolute or relative in the current working directory) of the directory of the images 
    width (int): (default=100)
    height (int): (default=100)
    inter : flag that takes one of the following cv2 interpolation methods/algorithm for image resize:
            cv2.INTER_LANCZOS4
            cv2.INTER_NEAREST (default)
            cv2.INTER_LINEAR
            cv2.INTER_AREA
            cv2.INTER_CUBIC

    Returns
    ------
    df (pandas DataFrame): 
        rows : images
        columns: 1st) integers representing the season each photo was taken
                      corresponding to the inital letter of each image file
                      (0:winter:W, 1:fall:F, 2:spring:S)
                rest) RGB color data for each pixel
    """
    path_abs = os.getcwd()+path
    images = os.listdir(path_abs)
    
    paths=[str(path_abs)+'/'+img for img in images]
    seasons=[]
    for filename in images:
        seasons.append(filename[0])
        dic = {'W':0, 'F':1, 'S':2}
        seasons=[dic.get(n, n) for n in seasons]
        
    dim = (width, height)
    x = [] # images as arrays
    for i,patth in enumerate(paths):
        # Read and resize image
        img_size_full = cv2.imread(patth)
        img_size_new = cv2.resize(img_size_full, dim, interpolation=inter)
        img_flat = img_size_new.reshape(-1,width*height*3).tolist()[0]
        obs_full = np.append(seasons[i], img_flat)
        x.append(obs_full)
        
    df = pd.DataFrame(x)
    
    names=['label']
    for i in range(1,(width+1)):
        for j in range(1,(height+1)):
            for c in ('B', 'G', 'R'):
                names.append(str(i)+','+str(j)+':'+c)
    df.columns = names
    df.index = [str(path)+'/'+img for img in os.listdir(path_abs)]
    return(df)
df=load_images('/images')
df
def pca_image_space_visualization(x, width=16., height=8.):
    """Read a labeled dataframe, calculate the two first principal components of the data
    and return a plot representing the 2-dimensional data.
    
    Parameters
    ----------
    x (pandas DataFrame) : dataframe i) labeled with label in column 0 
                                    ii) in the form of output of function load_images
                                        representing flattened image data
    width (float) : width of the plot (default=16.)
    height (float) : height of the plot (default=8.)
                                        
    Returns
    ------
    plot : in which images are visualized in the two-dimensional space resulting 
           from the projection of the data in the first two principal components
    """    
    from sklearn.preprocessing import StandardScaler
    # Standardizing the features
    features = x.iloc[:,1:].values
    features_st = StandardScaler().fit_transform(features)
    
    from sklearn.decomposition import PCA
    # initalize and apply PCA transformation
    pca = PCA(n_components=2)
    principalComponents = pca.fit_transform(features_st)
    # create new DataFrame with PCA values
    principalDf = pd.DataFrame(data = principalComponents
                               , columns = ['principal component 1', 'principal component 2'])
    paths = principalDf.index = list(x.index)
    finalDf = pd.concat([principalDf, x.iloc[:,0]], axis = 1)
    
    pca_image_space_visualization.explained_variance_ratio = pca.explained_variance_ratio_
                         
    import matplotlib.pyplot as plt
    from matplotlib.offsetbox import OffsetImage, AnnotationBbox
    pca1=list(finalDf.iloc[:,0].values)
    pca2=list(finalDf.iloc[:,1].values)
    img_paths=[os.getcwd()+patth for patth in paths]
    fig, ax = plt.subplots(figsize = (width, height))
    ax.scatter(pca1, pca2)
    ax.axis('off')
    for pca1, pca2, img_path in zip(pca1, pca2, img_paths):
        ab = AnnotationBbox(OffsetImage(plt.imread(img_path), zoom=0.25), (pca1, pca2))
        ax.add_artist(ab)
pca_image_space_visualization(df)<jupyter_output><empty_output><jupyter_text>In the above plot, the images located closer together are positively related to each other in some way, while images located far apart across the PCs are negatively related (always in respect of the linearly uncorrelated variables PC1 and PC2).But what about generalization of these or other conclusions from this plot for the original image space?<jupyter_code>print('The above plot represents image data that explain only the',
      "{:.2%}".format(pca_image_space_visualization.explained_variance_ratio.sum()),
      'of the variance/information of the initial data.')<jupyter_output>The above plot represents image data that explain only the 50.94% of the variance/information of the initial data.<jupyter_text>In that sense, we can say that the 2-D graph, using just PC1 and PC2, would be a 50.94% approximation or 50.94% accurate representation of the original image space.Let's insist a bit more:<jupyter_code>print('We can see that the x axis —corresponding to the first  principal component— represents the',
      "{:.2%}".format(pca_image_space_visualization.explained_variance_ratio[0]),
      ', \nwhile the y axis —corresponding to the second principal component— represents the',
      "{:.2%}".format(pca_image_space_visualization.explained_variance_ratio[1]),
      '\nof the variance(information) of the initial data.')<jupyter_output>We can see that the x axis —corresponding to the first  principal component— represents the 42.08% , 
while the y axis —corresponding to the second principal component— represents the 8.86% 
of the variance(information) of the initial data.<jupyter_text>At this point, we should keep in mind two facts and a graph:- The PCs are ranked (as first, second and so on) based on the percentage of variance/information of the original data explained by them. Thus the first principal component retains always more data information than the second one.- The PCs are linearly uncorrelated; the second principal component does not explain variance that has retained earlier by the first one.<jupyter_code>#variance0=[]

#for i in range(1,30):
#    pca = PCA(n_components=i)
#    pca.fit_transform(X)
#    variance0.append(pca.explained_variance_ratio_[i-1])
    
plt.grid(b=True, which='major', color='#666666', linestyle='-')
# Show the minor grid lines with very faint and almost transparent grey lines
plt.minorticks_on()
plt.grid(b=True, which='minor', color='#999999', linestyle='-', alpha=0.2)
plt.title('Scree plot (data variance explained by each PC)')
plt.xlabel('Sequence of PC')
plt.ylabel('Explained Variance')
plt.plot(range(1,30), variance0)<jupyter_output><empty_output><jupyter_text>Given, all the more, that in our case the first principal component explains/retains 3.75 times more information than the second one, then it is expected that is along this axis we perceive a much more semantically meaningful variance of its values, or else, it's logical it's along this axis that we have a separable clustering of winter images (obviously PC1 tends to be interpreted as the amount of white colour) even if fall and spring remain inseparable in a unified cluster. In that sense, it may be more correct and useful to have a graphical representation where the length of the axes is not set randomly but based justly on this ratio of retained variance between of the first two principal components<jupyter_code>pc1_exlained_variance = pca_image_space_visualization.explained_variance_ratio[0]
pc2_exlained_variance = pca_image_space_visualization.explained_variance_ratio[1]
r=pc2_exlained_variance/pc1_exlained_variance
pca_image_space_visualization(df, 16, 16*r)<jupyter_output><empty_output><jupyter_text>1-nn classification with the large data We will first develop a 1-nn classifier, so, in the first place we have to standarize our data<jupyter_code>from sklearn.preprocessing import StandardScaler
# standardize the features
X = df.iloc[:,1:].values
X = StandardScaler().fit_transform(X)
y= df.iloc[:,0].values
# import the knn classifier library and the cross-validation library
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score

# perform 5-fold cross-validation 1-nn classification and calculate the average accuracy
knn = KNeighborsClassifier(n_neighbors=1)
accuracy_1nn = cross_val_score(knn, X, y, cv=5, scoring='accuracy').mean()
print('for an 5-fold cross-validaiton 1-nn classification the average accuracy is:\n', accuracy_1nn)<jupyter_output>for an 5-fold cross-validaiton 1-nn classification the average accuracy is:
 0.7666666666666667<jupyter_text>At this point, we have to be careful because accuracy with a binary classifier is measured as:\begin{equation*}Acc =\frac{\text{Number of correct predictions}}{\text{Total number of predictions}} =\frac{TP + TN}{P + N} =\frac{TP + TN}{TP + FP + TN +FP}\end{equation*}Where P = Positives, N = Negatives, TP = True Positives, TN = True Negatives, FP = False Positives, and FN = False Negatives.But that is not our case, since we have to deal with a problem of a multiclass/3-class classification problem.The **multiclass classification accuracy**  is calculated as the average accuracy per class:\begin{equation*}Acc =\frac{\sum_{i=1}^k \frac{tp_i+tn_i}{tp_i+tn_i+fp_i+fn_i}}{k},\quad\quad \text{for $k:$ total number of classes}.\end{equation*} --- 1-nn classification with the low-dimensional data (using PCA decomposition)  \begin{equation*}|| X - W[t] C [t] || _F ^ 2 / || X || _F ^ 2 \end{equation*}\begin{equation*}ε (|| Χ - W [t] C [t] || _F ^ 2 - || Χ - W [t-1] C [t-1] || _F ^ 2) / || X || _F ^ 2 <ε \text{ with ε = 0.01 or 0.001 or 0.0001)}\end{equation*} At this point we should note that we will choose the number of components so as to retain at least 95% of the data variance:<jupyter_code># import PCA decomposition algorithm from the appropriate scikit-learn module
from sklearn.decomposition import PCA

# the data (X,y) are already standarized so we proceed directly to the decomposition
pca = PCA(n_components=.95)
principal_components = pca.fit_transform(X)
print('The minimun number of principal components that explain at least 95% of data variance is:\n',
      pca.n_components_,
      ',\nand indeed, the retained data variance (information) is:\n',
      "{:.2%}".format(pca.explained_variance_ratio_.sum()),
     '.\nSo, we have here a percentage significantly higher compared to the ',
      "{:.2%}".format(pca.explained_variance_ratio_[0:2].sum()),
      ' of the two prncipal components we saw above, as well as a significant reduction',
      ' in the number of features/estimators, from 30,000 to',
      pca.n_components_)<jupyter_output>The minimun number of principal components that explain at least 95% of data variance is:
 22 ,
and indeed, the retained data variance (information) is:
 95.59% .
So, we have here a percentage significantly higher compared to the  50.94%  of the two prncipal components we saw above, as well as a significant reduction  in the number of features/estimators, from 30,000 to 22<jupyter_text>With these 22-dimensional data we proceed again with the development of a 1-nn classifier.<jupyter_code># perform 5-fold cross-validation 1-nn classification and calculate the average accuracy
knn_pca = KNeighborsClassifier(n_neighbors=1)
accuracy_1nn_pca = cross_val_score(knn_pca, principal_components, y, cv=5, scoring='accuracy').mean()
print('for a 5-fold cross-validaiton 1-nn classification (PCA-95%variance) the average accuracy is:\n',
      accuracy_1nn_pca)<jupyter_output>for a 5-fold cross-validaiton 1-nn classification (PCA-95%variance) the average accuracy is:
 0.7333333333333333<jupyter_text>We thus observe a reduction in the average classification accuracy, that is not only expected, but also very similar to the reduction in the variance of the data with which we trained the *1-nn_pca* classifier. --- SVM classification with the large data This time we will try the four different kernel functions so as to choose the best one base on their classification accuracy<jupyter_code>from sklearn import svm
svm_rbf = svm.SVC(kernel='rbf', random_state=0, gamma='auto')
svm_linear = svm.SVC(kernel='linear', random_state=0, gamma='auto')
svm_poly = svm.SVC(kernel='poly', random_state=0, gamma='auto')
svm_sigmoid = svm.SVC(kernel='sigmoid', random_state=0, gamma='auto')
accuracy_svm_rbf = cross_val_score(svm_rbf, X, y, cv=5, scoring='accuracy').mean()
accuracy_svm_linear = cross_val_score(svm_linear, X, y, cv=5, scoring='accuracy').mean()
accuracy_svm_poly = cross_val_score(svm_poly, X, y, cv=5, scoring='accuracy').mean()
accuracy_svm_sigmoid = cross_val_score(svm_sigmoid, X, y, cv=5, scoring='accuracy').mean()
print('accuracy_svm_rbf    : ',
      "{:.2f}".format(accuracy_svm_rbf),
      '\naccuracy_svm_linear : ',
      "{:.2f}".format(accuracy_svm_linear),
      '\naccuracy_svm_poly   : ',
      "{:.2f}".format(accuracy_svm_poly),
      '\naccuracy_svm_sigmoid: ',
      "{:.2f}".format(accuracy_svm_sigmoid))<jupyter_output>accuracy_svm_rbf    :  0.80 
accuracy_svm_linear :  0.80 
accuracy_svm_poly   :  0.67 
accuracy_svm_sigmoid:  0.80<jupyter_text>--- SVM classification with the low-dimensional data (using PCA decomposition)  And we will do the same but this time with the 22-dimensional data from the PCA decomposition:<jupyter_code>accuracy_svm_pca_rbf = cross_val_score(svm_rbf, principal_components, y, cv=5, scoring='accuracy').mean()
accuracy_svm_pca_linear = cross_val_score(svm_linear, principal_components, y, cv=5, scoring='accuracy').mean()
accuracy_svm_pca_poly = cross_val_score(svm_poly, principal_components, y, cv=5, scoring='accuracy').mean()
accuracy_svm_pca_sigmoid = cross_val_score(svm_sigmoid, principal_components, y, cv=5, scoring='accuracy').mean()
print('accuracy_svm_pca_rbf    : ',
      "{:.2f}".format(accuracy_svm_pca_rbf),
      '\naccuracy_svm_pca_linear : ',
      "{:.2f}".format(accuracy_svm_pca_linear),
      '\naccuracy_svm_pca_poly   : ',
      "{:.2f}".format(accuracy_svm_pca_poly),
      '\naccuracy_svm_pca_sigmoid: ',
      "{:.2f}".format(accuracy_svm_pca_sigmoid))<jupyter_output>accuracy_svm_pca_rbf    :  0.43 
accuracy_svm_pca_linear :  0.80 
accuracy_svm_pca_poly   :  0.60 
accuracy_svm_pca_sigmoid:  0.80<jupyter_text>Further Investigation and DiscussionFrom the accuracy scores of the eight SVM classifiers we end up with two remarks:- considering altogether, both approaches andd all the kernel functions, the linear and sigmoid kernel functions seem to have better accuracy scores than the radial basis and the polynomial one,- regarding the accuracy score, the classification using linear and sigmoid kernel functions doesn't seem to be affected by whether we rely on all image data information or only on 95.6% of it.The second remark makes us wonder if there is a limit of retained data variance and number of principal components below which the performance of the classifiers (in terms of at least of accuracy) is lower than that of the classifiers trained with all the information and if there is, in general, some pattern (linear or not) that characterises this relation  "number of PC/explained data variance ~ accuracy of the classifier". After all, although we chose the number of principal components with a solid and widely accepted reasoning (selecting those that hold a minimum of 95% of the information), it would not be bad to be able to take advantage of even more of the basic advantages of PCA approach (especially the reduction in training time, in noise and overfitting, as well as from feature extraction aspect) without losing in terms of prediction performance.<jupyter_code>variance1=[]
variance2=[]

for i in range(1,23):
    pca = PCA(n_components=i)
    pca.fit_transform(X)
    variance1.append(pca.explained_variance_ratio_.sum())
    variance2.append(pca.explained_variance_ratio_.sum()/i)
plt.figure(figsize=(12,5))

plt.subplot(1, 2, 1)
# Show the major grid lines with dark grey lines
plt.grid(b=True, which='major', color='#666666', linestyle='-')
# Show the minor grid lines with very faint and almost transparent grey lines
plt.minorticks_on()
plt.grid(b=True, which='minor', color='#999999', linestyle='-', alpha=0.2)
plt.axhline(y=0.95, color='r', linestyle='-')
plt.text(4.5, 0.96, '95% cut-off threshold', color = 'red', fontsize=11)
plt.title('Total Variance retained vs. number of PCs')
plt.xlabel('Number of PCs')
plt.ylabel('Total Explained Variance')
plt.plot(range(1,23), variance1)


plt.subplot(1, 2, 2)
# Show the major grid lines with dark grey lines
plt.grid(b=True, which='major', color='#666666', linestyle='-')
# Show the minor grid lines with very faint and almost transparent grey lines
plt.minorticks_on()
plt.grid(b=True, which='minor', color='#999999', linestyle='-', alpha=0.2)
plt.axhline(y=.9559/22, color='r', linestyle='-')
plt.text(.5, 0.051, '95% cut-off threshold', color = 'red', fontsize=11)
plt.title('Average Explained Variance per PC vs. number of PCs')
plt.xlabel('Number of PCs')
plt.ylabel('Average Explained Variance per PC')
plt.plot(range(1,23), variance2)<jupyter_output><empty_output><jupyter_text>To this end, the two above graphs:- the cumulative sum of data variance (eigenvalues) for different numbers of principal components (eigenvectors),- the average (per PC) cumulative sum of data variance of the above eigenvalues against the number of PCs (dividing each sum of eigenvalues by the number of corresponding PCs (eigenvectors) before plotting it against the PCs' number) provide a good indication of when we hit the point of diminishing returns (i.e., little variance is gained by retaining additional eigenvalues), and that is around the number 5 of PCs.After all, given the size of our data set, we are able to calculate (always with a 5-fold cross-validation) the accuracy of the SVM classifier for the four different kernel functions for all different numbers of principal components.<jupyter_code>svm_rbf = svm.SVC(kernel='rbf', random_state=0, gamma='auto')
accuracies_svm_rbf =  []
svm_linear = svm.SVC(kernel='linear', random_state=0, gamma='auto')
accuracies_svm_linear = []
svm_poly = svm.SVC(kernel='poly', random_state=0, gamma='auto')
accuracies_svm_poly = []
svm_sigmoid = svm.SVC(kernel='sigmoid', random_state=0, gamma='auto')
accuracies_svm_sigmoid = []

for i in range(1,23):
    pcs = PCA(n_components=i).fit_transform(X)
    accuracies_svm_rbf.append(cross_val_score(svm_rbf, pcs, y, cv=5,scoring='accuracy').mean())
    accuracies_svm_linear.append(cross_val_score(svm_linear, pcs, y, cv=5, scoring='accuracy').mean())
    accuracies_svm_poly.append(cross_val_score(svm_poly, pcs, y, cv=5, scoring='accuracy').mean())
    accuracies_svm_sigmoid.append(cross_val_score(svm_sigmoid, pcs, y, cv=5, scoring='accuracy').mean())
line1, = plt.plot(range(1,23), accuracies_svm_rbf, label="rbf", color='blue')
line2, = plt.plot(range(1,23), accuracies_svm_poly, label="poly", color='red')
line3, = plt.plot(range(1,23), accuracies_svm_linear, label="linear", color='green')
line4, = plt.plot(range(1,23), accuracies_svm_sigmoid, label="sigmoid", color='orange')
# Create a legend for the first line.
first_legend = plt.legend(bbox_to_anchor=(1, 1), title="Kernel functions")
# Add the legend manually to the current Axes.
ax = plt.gca().add_artist(first_legend)
# Show the major grid lines with dark grey lines
plt.grid(b=True, which='major', color='#666666', linestyle='-')
# Show the minor grid lines with very faint and almost transparent grey lines
plt.minorticks_on()
plt.grid(b=True, which='minor', color='#999999', linestyle='-', alpha=0.2)
plt.title('Classification Accuracy per number of PCs')
plt.xlabel('Number of PCs')
plt.ylabel('Classification Accuracy')
plt.show()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup