%% Machine Learning - Classification
% Table of Contents
% Ingesting: Reading the data matrix
%     Data Matrix
% Prepare: Preprocessing Data
%     Process Null Values
%         Remove rows that have a missing target values
%         Processing Missing Feature Values
%     Handle Outliers
%         Definition of an Outlier - IQR
%         Matlab functions for handling outliers
%         Customized Specifications of Outliers
% Feature Scaling
% Feature Scaling and Normalization
%     Min-Max Scaler
%     Standard Scaler
%         Matlab functions for scaling and normalization
% Analyzing Features
%     Feature Cross-Correlation Plots
% Data Visualization
% Working with full Iris DataSet
% Function for reading in datamatrix
% Function for Preprocessing
% Feature Correlations
% Correlation Plots - Independence of Features
% Multidimensional Scaling  - Visualizing in 2D
% Model Learning: Model Training
%     Split data matrix into training and testing sets
%     Training a Model
%     Estimation of training success
%     Hyperparameter Tuning
% Benchmarking Model on Previously Unseen Test Data
% Confusion Matrix - Evaluates the performance for each class
%     Visual representaiton on confusion matrix
%     Numerical Values -- Computing F1 Metric
% Receiver Operating Curve
%  Appendix
%     Functions
%         Function f1score
%
%
%% Ingesting: Reading the data matrix
% In its canonical form, the data matrix is set with samples along its rows, which is direction labelled 1, and columns as the features or variables, which is the direction labelled 2. Generally, the datamatrix has one column labelled Y which is the target variable.
%
%% Data Matrix
% The canonical form for the data matrix is to have
% samples along the rows
% features along the colums
% Individual vectors X_k = \{x_{k,1}, x_{k,2}, \ldots x_{k,p}\}is a p-dimensional vector representing the k^{th}sample.
%
% \[ \mathbf{X} = \left( \begin{array}{cccccc}
% x_{1,1} ~~&~~ x_{1,2} ~~&~~ x_{1,3} ~~&~~ \ldots ~~&~~ x_{1,p-1} ~~&~~ x_{1,p} \\
% x_{2,1} & x_{2,2} & x_{2,3} & \ldots & x_{2,p-1} & x_{2,p} \\
%  &  & \ldots &  &  & \\
% x_{i,1} & x_{i,2} & x_{i,3} & \ldots & x_{i,p-1} & x_{i,p} \\
%  &  & \ldots &  &  & \\
% x_{n,1} & x_{n,2} & x_{n,3} & \ldots & x_{n,p-1} & x_{n,p} \\
% \end{array} \right) \]
% The data matrix in the following example is in its canonical form.
% Sometimes, the ingested data has samples along its columns, and feature names as row names. In this case, the data matrix will need to be transposed to bring it into its canonical form.
% X = transpose(X)
% Or simply:
% X = X'
% Note that as such, the table that is read in, the variable dm in the example below, cannot be directly be transposed since a table stores variable names detected from the first row of the cvs file along its columns.
%

dm = readtable('Data/iris_subset.csv');

target = "species";
varnames = dm.Properties.VariableNames;
features = varnames (~strcmp(varnames, target));

X = dm{:, features}
y = dm{:, target}

%
%% Prepare: Preprocessing Data
% Remove or otherwise process missing or empty values
% Remove or otherwise process outliers
% Scale Data - Minmax Scaling or Standard Scaling
%
%
%% Process Null Values
%
%% Remove rows that have a missing target values
%

% find and remove rows where target is missing
 missing_target = ismissing(y);
 X(missing_target,:)=[];
 y(missing_target) = [];
 X, y

%

% Assign numerical values (1, 2, 3, ..) to each 
% unique non-null target value 

[labels, ~, y] = unique(y);
disp ("Target Values Coded as 1, 2, 3");
labels

%
%% Processing Missing Feature Values
% Strategy 1: One strategy to handle null values is to remove rows (samples) containing any null value. Function rmmissing can remove rows (dimension = 1), or columns (dimension = 2) that have missing values
% X = rmmissing (X, 1);	% remove rows with missing data
% X = rmmissing (X, 2); 	% remove columns with missing data
%
%
% Strategy 2: (1) Missing values are to replace the missing value with a value that is listed closest to the missing value in the data matrix. (2) Filling the missing values by creating a piecewise linear function, or (3) by fitting a piecewise cubic spline between the known points and interpolating for the missing value.
% X = fillmissing(X, "nearest") % copies value from neighbors
% X = fillmissing(X, "linear")  % linear interpolation of data
% X = fillmissing(X, "spline")  % spline interpolation of data
%

X

%

X = fillmissing(X, "linear")

%
%% Handle Outliers
%
%
%% Definition of an Outlier - IQR
% There are three standardized definitions of data outliers.
% First, the definition is based on Inter-Quartile Range or IQR which is defined as follows:
% IQR = Q3 - Q1~~~~~~\mathit{where~Q1~is~lower~quartile,~and~Q3~the~upper~quartile}
% Values less than the lower gate or greater than upper gate are considered as outliers
%
% \mathrm{Lower~Outlier~Gate} = Q1 - 1.5 \times IQR \\
% \mathrm{Upper~Outlier~Gate} =  Q3 + 1.5 \times IQR
% Second, the definition based on mean and standard deviations:
%
% \mathrm{Lower~Outlier~Gate} = \mu - 3 \cdot \sigma \\
% \mathrm{Upper~Outlier~Gate} = \mu + 3 \cdot \sigma
%
%% Matlab functions for handling outliers
% In addition to providing a function isoutlier(), Matlab also has functions for replacing outlier values with its filloutlier() function.
% The following functions use Approach Three for identifying the outliers and replacing their values with chosen data value.
% X = filloutliers(X, "center")  % median value
% X = filloutliers(X, "clip")    % clips to upper/lower threshold
% X = filloutliers(X, "nearest") % Nearest of, or
% X = filloutliers(X, "linear")  % Linear interpolation of, or
% X = filloutliers(X, "spline")  % Cubic spline interpolation
%                                % of non-outlier values
% The following function use a defined method to look for an outlier, and replace it with the chosen data value.
% In (1) outliers are defined using the "mean" method these are points more than three standard deviations from the mean, and replaced by the nearest element that is not an outlier.
% In (2) the outliers are determined by the MAD method described above and replaced by linearly interpolating over non-outlier values.
% And in (3), the outliers are determined by the inter-quartile gates defined in method(1).
% X = filloutliers (X, "nearest", "mean")         % (1)
% X = filloutliers (X, "linear", "median")        % (2)
% X = filloutliers (X, "spline", "quartiles")     % (3)
%% Customized Specifications of Outliers
% Outliers may be identified and processed in a customized manner as well.
%
%% Feature Scaling
%
%% Feature Scaling and Normalization
% After the null values and outliers are processed, the next preprocessing phase aims at normalization and scaling of the dataframe. This sets the data values on somewhat of a leveled playing field so that relative disparity in the range of values used for a specific feature does not bias any machine learning algorithms utilized downstream. In both scaling and normalization, the values of the numeric variables are transformed so that the transformed data points have specific and desirable properties.
%
%% Min-Max Scaler
% Data scaling transforms the data so that it fits within a specific scale, such as 0--100 or 0--1. Usually, 0-1 is utilized and essentially it brings the data values into a predefined range for all features, which is useful for metric-driven algorithms such as k-Nearest Neighbors and Support Vector Machines.
% Minmax Transformation takes applied to each feature, or the column of the datamatrix:
%             X_{scl} = \frac{X - X_{min}}{X_{max} - X_{min}}
%% Standard Scaler
% Stadard scaler performs a normalization of the such that the transformed data values follow a normal distribution. Normalizing the data is often desirable for a class of statistical machine learning algorithms that expect that the input data will be normally distributed, such as linear discriminant analysis and naive Bayesian. Normalization transforms data values to their corresponding Z- scores with the transformed value X_{new} defined as follows using the mean \mu_X and standard deviation \sigma_X of the observed values for feature X:
%             X_{new} = \frac{X - \mu_X}{\sigma_X}
%
%% Matlab functions for scaling and normalization
% The function normalize(X, dim, method) may instead be also used for peforming both the scaling and normalization operations discussed above. The function may be used to operate along dimension 1 or dimension 2. The following examples illustrate the use of this function along direction 1 since the canonical datamatrix is represented with columns as features.
% X = normalize (X, 1, 'zscore') % results in zscores
% X = normalize (X, 1, 'center') % results in mean = 0
% X = normalize (X, 1, 'range')  % scales to [0,1] interval
%

X  = normalize(X, 1, 'zscore')              % standard scaler

%
%% Analyzing Features
%
% One of the key goal in developing robust models is to use maintain build models using non-redudant features. It is also important to analyze how the feature values are distributed.
%
%% Feature Cross-Correlation Plots
%

figure()
p = corrplot(dm, TestR="on", DataVariables=features);

%
%% Data Visualization
% Requires data to be reduced to two dimensions. A useful technique for visualization is Multidimensional Scaling.
% MDS problem is posed as having a collection of 𝑛 objects with pairwise distances of \{d_{ij}\}, 1 ≤ 𝑖, 𝑗 ≤ 𝑛, which are then represented in the Euclidean space with points:
%         \mathbf{y_1, y_2 \ldots y_n} \in \mathrm \Re^k, such that:
%         ||\mathbf{y_i} - \mathbf{y_j}|| = d_{ij} ~~~ \forall i, j
% The idea in multidimensional scaling is to begin with a matrix of pairwise distances and then determine the location of points in a 2D or 3D space that substantially satisfies the distance criteria imposed by the pairwise matrix.
% A function for performing Classical MDS for dimensionality reduction is included in Matlab’s statistical and machine learning toolbox. The following example illustrates the use of this function cmdscale() that transforms a distance matrix into set of vectors Y.
% [Y, ev] = cmdscale (D, 2)    % transform to a 2D space
%
% First step in to develop a distance matrix. The datamatrix X \mathrm{~is~a~} N \times p  matrix corresponsing to N samples each being a p dimensional vector. A scatter matrix is constructed as X^TX is used as the input argument to the function cmdscale().
%

D = pdist(X_norm);

%

[Y, E] = cmdscale(D, 2);

%

figure()
scatter(Y(:,1), Y(:,2), 50, y, 'filled')

%
%% Working with full Iris DataSet
%
% The iris data set is a commonly utilized dataset in machine learning. The data set comprises 150 measurements, where each measurement includes four features, namely the length and width of the petals, and the length and width of the sepals of three categories of flowers. The three categories of the dataset are iris setosa, iris versicolor, and iris verginica. The data set is well balanced, as there are approximately 50 samples from each of the three categories.
%
% The outer parts of the flower which is often colored green that encloses a developing bud is called the sepal. The colored part of a flower comprises of petals.
%
%
%% Function for reading in datamatrix
% Reads in the csv file, removes rows that have an undefined value for target variable, and also replaces target values with unique integer keys.
% [X, y, features, labels] = read_datamatrix ("iris_subset.csv", "species");
%

[X, y, features, labels] = read_datamatrix("Data/iris.csv", "species");
disp("Target Labels coded to 1, 2, 3");
labels

%
%% Function for Preprocessing
% Removes all X rows, and y rows where y values are missing
% Replaces all missing X values with mean values of the column/feature
% Outliers are replaced with median value
% Apply min-max scaler
%

[X, y] = preprocess (X, y);

%
%% Feature Correlations

figure();

tc{1} = [1 0 0]; % red
tc{2} = [0 1 0]; % green
tc{3} = [0 0 1]; % blue
colors = cell2mat({tc{y}}')

%

c=1;
for i = 1:(size(X,2)-1)
    for j = (i+1):size(X,2)
        subplot(2, 3, c);
        scatter (X(:,i), X(:,j), 20, colors, 'filled');
        grid;
        title("(" + features{i} + ", " + features{j} + ")", ...
            'Interpreter', 'None');
        c = c+1;
    end
end

%
%% Correlation Plots - Independence of Features
%

figure()
p = corrplot(X,  'VarNames', {'PtL', 'PtW', 'SpL', 'SpW'});

%
%% Multidimensional Scaling  - Visualizing in 2D
%

figure();
D = pdist(X, "minkowski", 10);
[Y, E] = cmdscale(D, 2);
scatter(Y(:,1), Y(:,2), 50, colors, 'filled');
grid

%
% Perform Multidimensional scaling with cityblock distance metric.
%

figure()
D = pdist(X, "cityblock");
[Y, E] = cmdscale(D, 2);
scatter(Y(:,1), Y(:,2), 50, colors, 'filled');
grid

%% Model Learning: Model Training

% Read  and prepare the data 
[X, y] = preprocess (X, y);

%
%% Split data matrix into training and testing sets
%

rng(10);

%

% Could use a number for partitioning also
sample_cnt = size(X, 1)
part = cvpartition(sample_cnt, "HoldOut", 0.3)

%

feats = 1:4;
X_trn = X(part.training, feats);
X_tst = X(part.test, feats);
y_trn = y(part.training);
y_tst = y(part.test);

%
%% Training a Model
% k - Nearest Neighbor Model
% Based on the idea of similar items are close to each other in the feature space.
%
%

mdl = fitcknn(X_trn, y_trn, 'NumNeighbors', 3);

%
%% Estimation of training success
% In order to get a base estimation of success of training, we can compute an error rate for finding how well does the model work on classifiying the training data itself
%

[predict_trn, ~] = predict(mdl, X_trn);
cp_trn = classperf(y_trn, predict_trn);
fprintf ("Performance on Training Data Set: %5.3f", ...
                        (1-cp_trn.ErrorRate)*100);

%
% Alternatively, MATLAB's function resubLoss provides this information
%

train_error = resubLoss(mdl);
fprintf ("Performance on Training Data Set: %5.3f", ...
            (1-train_error)*100);

%
%% Hyperparameter Tuning
% Evaluate the value of parameter "k" that yeilds the best performance
%

kvals = [3, 5, 7, 9, 11, 15];
mdl_accuracy = dictionary;
for k = kvals
    mdl = fitcknn(X_trn, y_trn, 'NumNeighbors', k);
    mdl_accuracy{k} = (1-resubLoss(mdl));
end

%

mdl_accuracy = [mdl_accuracy{kvals}];

%

figure()
% plot a bar graph of accuracy vs. k
bar(categorical(kvals), mdl_accuracy)

%
% Choose the value of "k" that yields the best accuracy
%

[~, indx_max] = max(mdl_accuracy);
k_opt = kvals(indx_max);
fprintf ("Optimal value for neighborhood k is %d", k_opt);

%
% Next train a model with this optimal value of "k"
%

final_mdl = fitcknn(X_trn, y_trn, 'NumNeighbors', k_opt);
train_error = resubLoss(final_mdl);
fprintf ("Performance on Training Data Set for k = %d is %5.3f", ...
                       k_opt, (1-train_error)*100);

%
%% Benchmarking Model on Previously Unseen Test Data
%

[predict_tst, predict_scores] = predict(final_mdl, X_tst);
cp_tst = classperf(y_tst, predict_tst);
fprintf ("Performance on Testing Data Set: %5.3f", ...
                    (1-cp_tst.ErrorRate)*100);

%
%% Confusion Matrix - Evaluates the performance for each class
%
%% Visual representaiton on confusion matrix
%

confusionchart(y_tst, predict_tst)

%
%% Numerical Values -- Computing F1 Metric
%

cm_test = confusionmat (y_tst, predict_tst)
perf_score = f1score(cm_test);
perf_score.F1
fprintf("Overall F1 score is %4.3f", mean(perf_score.F1));

%
%% Receiver Operating Curve
%

figure()
rocObj = rocmetrics(y_tst, predict_scores, mdl.ClassNames);
plot(rocObj)
rocObj.AUC

%

[predict_trn, trn_scores] = predict(final_mdl, X_trn);
figure()
rocObj = rocmetrics(y_trn, trn_scores, mdl.ClassNames);
plot(rocObj)

%
%%  Appendix
%% Functions
% preprocess performs the following functions.
% Removes all X rows, and y rows where y values are missing
% Replaces all missing X values with linear interpolation
% Outliers are replaced with median value
% Apply min-max scaler
%

function [X, y] = preprocess(X, y)
    % remove rows where target is missing
    missing_target = ismissing(y);
    X(missing_target,:)=[];
    y(missing_target) = [];
    
    % linear interpolation for missing features
    X = fillmissing(X, "linear");

    % replace outliers with median
    X = filloutliers(X, "center");

    % Normalize with min-max scaler
    X = normalize(X, 1, 'range');
end

% Function reads in csvfile treating target as the target variable, and treating the remaining variables  as numeric features.
% Unique value of target variables are replaced with values of 1, 2, 3, etc.

function [X, y, features, labels] = read_datamatrix(csvfile, target)

    dm = readtable (csvfile);
    varnames = dm.Properties.VariableNames;

    features = varnames (~strcmp(varnames, target));

    X = dm{:, features};
    y = dm{:, target};
    missing_target = ismissing(y);
    X(missing_target,:)=[];
    y(missing_target) = [];

    [labels, ~, y] = unique(y);
end

%% Function f1score
%

%% Function f1score (confusion_matrix)
% Computes F1 Scores using data in confusion matrix
% confusion_matrix - N x N confusion matrix

%

function F1scores = f1score(confusion_matrix)
    tp = diag(confusion_matrix) + eps;

    tp_plus_fp = sum(confusion_matrix, 1)';
    tp_plus_fn = sum(confusion_matrix, 2);

    precision = tp ./ tp_plus_fp;
    recall = tp ./ tp_plus_fn;

    F1scores = harmmean([precision' ; recall']);
    F1scores = table(F1scores','VariableNames',{'F1'});

end

%%
%%