%% Mass Spectromtery Data Analysis - Mouse Pacreatic
%
% Table of Contents
% Read spectrograph file
% Plot of 3 randomly selected spectras
% Envelopes of the spectrograms
% Create Predictors and Response Variables
% Perform Preprocessing
% Partition the data into training and test sets
% SVM Classifier
% SVM Classifier - Training
% SVM Classifier - Testing
% Dim Reduction with PCA
% Training with Principal Components
% Testing with Principal Components
% Classification Learner App
% Functions
%     Function read_mass_spec
%     Function preprocess
%     Function f1score
%% Read spectrograph file

rng("default")
DATA = ".";
control_data = DATA + "/Data/Pancreas/Control";
cancer_data = DATA + "/Data/Pancreas/Cancerous";

%

normal = read_mass_spec(control_data, "CSV");
cancer = read_mass_spec (cancer_data, "CSV");

%

% Get sample counts
c_sample_cnt = size(cancer, 2);
n_sample_cnt = size(normal, 2);

% Collect m/z and intensities of cancer/normal
c_m_z = [cancer.M_Z]';
c_intensity = [cancer.Intensity]';
n_m_z = [normal.M_Z]';
n_intensity = [normal.Intensity]';

%
%% Plot of 3 randomly selected spectras

% m_z stays constant across spectrographs
m_z = n_m_z(1,:);

xAxisLabel = 'Mass/Charge (M/Z)';       % x label for plots
yAxisLabel = 'Ion Intensity';           % y label for plots
c_sample_indx = randi(c_sample_cnt, 3, 1);
n_sample_indx = randi(n_sample_cnt, 3, 1);

%

figure;

n_index = randi(n_sample_cnt)
plot(m_z, n_intensity(randi(n_index),:), 'b');
axis([0 10000 0 5000])
grid
xlabel('M/Z')
ylabel('Intensity')
set(gca, 'FontSize', 14); % Applies to current axes

%

figure; hold on;

hC = plot(m_z, c_intensity(c_sample_indx,:), 'r');
hN = plot(m_z, n_intensity(n_sample_indx,:), 'b');
set(gca, 'FontSize', 14); % Applies to current axes

intensity = [c_intensity; n_intensity];

%


xlabel(xAxisLabel); ylabel(yAxisLabel);
axis([min(m_z) max(m_z) min(min(intensity)) max(max(intensity))/10])
legend([hC(1),hN(1)],{'Pre-Pancreatic Cancer', 'Control'})
title('Sample Spectrograms')
set(gca, 'FontSize', 14); % Applies to current axes
hold off;

%
%% Envelopes of the spectrograms

% mean, min, and max for normal samples
hold on
mean_N = mean(n_intensity);
max_N = max(n_intensity); 
min_N = min(n_intensity); 

% mean, min, and max for cancer samples
mean_C = mean(c_intensity);  
max_C = max(c_intensity); 
min_C = min(c_intensity); 

figure; hold on;
hC = plot(m_z,mean_C,'r--');
hN = plot(m_z,mean_N,'b--');

gC = plot(m_z,[max_C' min_C'],'r-');
gN = plot(m_z,[max_N' min_N'],'b-');
xlabel(xAxisLabel); ylabel(yAxisLabel);
axis([min(m_z) max(m_z) min(min(intensity)) max(max(intensity))/10])
legend([hN, hC, gN(1), gC(1)],{'Control Group Avg.','Pre-Pancreas Cancer Group Avg',...
                            'Control Envelope','Pre-Pancreas Cancer Envelope'},...
                            'Location', 'NorthWest')
title('Group Average and Group Envelopes')
set(gca, 'FontSize', 14); % Applies to current axes
hold off;

%
%% Create Predictors and Response Variables

% Response Output variables
y = zeros(size(c_intensity, 1) + size(n_intensity, 1), 1);
% Normal samples -- output = 0, Abnormal -- output = 1
y (1:end) = 0;
y (1:size(c_intensity, 1)) = 1;
% Predictor Matrix: Samples x Features 
X = [c_intensity; n_intensity];

%
%% Perform Preprocessing
%

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

%
%% Partition the data into training and test sets

per_eval = 0.20;          % training size for cross-validation
rng('default'); 
dataset_size = numel(y);
cv  = cvpartition(dataset_size, 'holdout', per_eval)
X_trn = X(cv.training, :);
y_trn = boolean(y(cv.training));
X_tst = X(cv.test, :);
y_tst = boolean(y(cv.test));

%
%% SVM Classifier
%
%% SVM Classifier - Training

Mdl = fitcsvm(X_trn, y_trn);
trainError = resubLoss(Mdl);
sprintf('Training Error: %4.3f', trainError)

%
%% SVM Classifier - Testing

[y_pred, scores] = predict(Mdl, X_tst);
cm = confusionchart(y_tst, y_pred);
cm.FontSize = 20
classperf(y_tst, y_pred)
cmat = confusionmat(y_tst, y_pred);
f1score(cmat)

%
%% Dim Reduction with PCA

% Retain 99% variance

%

[coeff,score,~,~, explained,mu] = pca(X);

%

ev = min(find(cumsum(explained)>99))
Xpc = score(:, 1:ev);

%
%% Training with Principal Components

Xpc_trn = Xpc(cv.training, :);
Xpc_tst = Xpc(cv.test, :);
Mdl_pc = fitcsvm(Xpc_trn, y_trn);
trainError_pc = resubLoss(Mdl_pc);
sprintf('Training Error: %4.3f', trainError_pc)

%
%% Testing with Principal Components

[y_pred_pc, scores_pc] = predict(Mdl_pc, Xpc_tst);
cm = confusionchart(y_tst, y_pred_pc);
cm.FontSize = 20
classperf(y_tst, y_pred_pc)
cmat_pc = confusionmat(y_tst, y_pred_pc);
f1score(cmat_pc)

%
%% Classification Learner App
%

classificationLearner

%
%% Functions
%% Function read_mass_spec
% Function reads mass spectra from a folder
%

function [Y, cnt] = read_mass_spec (folderlist, ext)
    cnt = 0;
    for i = 1:length(folderlist)
        folder = folderlist{i};
        info = dir(folder);
        all_files = {info.name};
        data_files = cellfun(@(x) contains(upper(x), upper(ext)), all_files);
        files = all_files(data_files);
    
        % Read all spectrographs -- skip the first (header) row
        warning('OFF', 'MATLAB:table:ModifiedAndSavedVarnames');
        for j = 1:numel(files)
            filename = folder + "/" + files{j};
            spec = readtable(filename);
            cnt = cnt + 1;
            Y(cnt).M_Z = spec.M_Z;
            Y(cnt).Intensity = spec.Intensity;
        end
    end
end

%
%% Function 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 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