%% TD2 - Deep (Feedforward) Neural Network on the Breast Cancer Dataset (MATLAB)
% -------------------------------------------------------------------------  This script shows the complete workflow:
% 1) Load data
% 2) Put data in the format expected by trainNetwork
% 3) Split into training and validation sets 4
% ) Define a deep fully-connected network (a.k.a. a "multilayer perceptron")
% 5) Train the network
% 6) Evaluate performance with a confusion matrix
% IMPORTANT : MATLAB has TWO "neural network" worlds: - Classic Neural Network Toolbox functions (patternnet, feedforwardnet, ...) often expect: features in rows, samples in columns. -
%
%  Deep Learning Toolbox (trainNetwork + layers) typically expects: samples/observations in rows, features in columns for tabular data.
% The built-in dataset cancer_dataset uses the classic orientation: cancerInputs is 9 x 699 (9 features, 699 samples) cancerTargets is 2 x 699 (2 classes, one-hot encoding)
% For trainNetwork with tabular data, we will transpose cancerInputs to get 699 x 9. -------------------------------------------------------------------------

clear; clc; close all;

% For teaching and debugging, it is useful to make results reproducible.
% This fixes the random seed used for:
%   - the random train/validation split
%   - the random weight initialization
% (If you want different results each run, comment this out.)
rng(0);

%
%% 1) Load the dataset
% You can load the .mat file if it is on the MATLAB path. If "load cancer_dataset.mat" fails on your computer, try: load cancer_dataset (without .mat) or (on many MATLAB versions): [cancerInputs, cancerTargets] = cancer_dataset;

load cancer_dataset.mat

% Sanity check: look at the size of the loaded variables
fprintf('Original cancerInputs size  : %d x %d (features x samples)\n', size(cancerInputs,1), size(cancerInputs,2));
fprintf('Original cancerTargets size : %d x %d (classes  x samples)\n\n', size(cancerTargets,1), size(cancerTargets,2));

%
%% 2) Put predictors X in the format expected by trainNetwork
% For TABULAR data, trainNetwork expects: N x P (N observations, P features) So we transpose cancerInputs.

X = cancerInputs.';          % 699 x 9

% OPTIONAL: feature scaling / normalization
% Many models train better if each feature has mean ~0 and std ~1.
% We will do normalization inside the featureInputLayer (see below),
% so we do NOT manually normalize here.

%
%% 3) Convert the one-hot targets to ONE label per observation
% cancerTargets is one-hot encoded: each column corresponds to one sample exactly one entry in the column is 1 (the correct class)
% Example column could be [1;0] or [0;1]. trainNetwork with classificationLayer expects Y as: N x 1 categorical vector (one label per observation)

% Convert one-hot -> class index using max (no toolbox-specific functions needed)
[~, classIdx] = max(cancerTargets, [], 1);   % 1 x 699 (values 1 or 2)
classIdx = classIdx.';                      % 699 x 1

% Give human-friendly names to the classes.
% In most versions of the dataset, row 1 corresponds to Benign and row 2 to Malignant.
% To be safer, we infer which row is more frequent and call it Benign.
counts = sum(cancerTargets, 2);  % 2 x 1
if counts(1) >= counts(2)
    classNames = ["Benign" "Malignant"];  % class 1 = Benign, class 2 = Malignant
else
    classNames = ["Malignant" "Benign"];  % rare case: reversed
end

Y = categorical(classIdx, [1 2], classNames);

% Sanity check: show class distribution
fprintf('Class distribution in full dataset:\n');
disp(countcats(Y));

%
%% 4) Split the data into training and validation sets
% Deep learning workflows typically use: - Training set: used to update the weights - Validation set: used DURING training to monitor generalization A separate test set can be used at the very end (not done here).

trainRatio = 0.80;  % 80% training, 20% validation
[XTrain, YTrain, XVal, YVal] = trainTestSplit(X, Y, trainRatio);

fprintf('After split:\n');
fprintf('  XTrain: %d x %d\n', size(XTrain,1), size(XTrain,2));
fprintf('  XVal  : %d x %d\n', size(XVal,1),   size(XVal,2));
fprintf('  YTrain: %d labels\n', numel(YTrain));
fprintf('  YVal  : %d labels\n\n', numel(YVal));

%
% Note that the features should be spanning the rows, and the samples should span the columns.
% Additionally, the  data type needs to be converted to a categorical.
% Finally, we will only  use one of the target labels for our binary classification t
% Setting up the Network Architecture
% We will train a 5 layer deep feedforward neural network with the following neural network architecture:



%
% Compared to our shallow neural network classifier, we need to specify 4  components in this network architecture, which are described below:
% Input data layer: The function to call for piping in the input dataset differs on the  data type inputted (such as a 2D numerical, tensor, etc.).
% For our  dataset, the sequenceInputLayer() function is needed, requiring the number of features as the parameter.
% Hidden layer: To design a fully connected feedforward neural network, we need to call the fullyConnectedLayer() function, which requires the number of activation nodes as the first parameter. The WeightsInitializer() and BiasInitializer() arguments explain how to initialize the weight and bias parameters. In the code  below, we're specifying using small random numbers from a normal  distribution for symmetry breaking.
% Activation function: Between each hidden layer is the activation function. While we can use  the cross-entropy layer for binary classification, we can also use the softmaxLayer() function, which generalizes to multi-class classification problems.
% Classification layer: The classificationLayer() decodes the probabilities to their predicted labels.
% We'll store these components into the variable layers:
%
%% 5) Define the network architecture
% For tabular feature vectors, use featureInputLayer (NOT sequenceInputLayer). - sequenceInputLayer is intended for sequences/time series where each observation is a sequence of feature vectors. - Here, each observation is just ONE fixed feature vector (9 numbers).

numFeatures = size(XTrain, 2);   % should be 9
numClasses  = numel(categories(Y));  % should be 2

layers = [
    featureInputLayer(numFeatures, "Normalization","zscore", "Name","input")

    fullyConnectedLayer(4, "Name","fc1", ...
        "WeightsInitializer","narrow-normal", ...
        "BiasInitializer","narrow-normal")

    % Activation function:
    % If your MATLAB version does NOT support sigmoidLayer/tanhLayer,
    % replace them with reluLayer (widely available).
    sigmoidLayer("Name","sigmoid1")

   
    fullyConnectedLayer(2, "Name","fc3", ...
        "WeightsInitializer","narrow-normal", ...
        "BiasInitializer","narrow-normal")

    tanhLayer("Name","tanh3")

    % Output layer: number of neurons = number of classes
    fullyConnectedLayer(numClasses, "Name","fc_out", ...
        "WeightsInitializer","narrow-normal", ...
        "BiasInitializer","narrow-normal")

    softmaxLayer("Name","softmax")
    classificationLayer("Name","classOutput")
];

% Visualize the network (nice for classroom demos)
analyzeNetwork(layers);

%
%% 6) Choose training options (hyperparameters)
% - optimizer: adam (good default) - learning rate: controls step size - epochs: how many passes over the training data - validation data: MUST match training data format

options = trainingOptions("sgdm", ...
    "InitialLearnRate", 0.01, ...
    "MaxEpochs", 200, ...
    "MiniBatchSize", 32, ...
    "Shuffle", "every-epoch", ...
    "ValidationData", {XVal, YVal}, ...
    "ValidationFrequency", 30, ...
    "Verbose", true, ...
    "Plots", "training-progress");

%
%% 7) Train the deep neural network

[dnn, trainingInfo] = trainNetwork(XTrain, YTrain, layers, options);

%
%% 8) Make predictions on the validation set

YPred = classify(dnn, XVal);

%
%% 9) Evaluate results with a confusion matrix
% plotconfusion is mostly for the classic NN toolbox and expects one-hot matrices. For categorical labels, use confusionchart.

figure;
confusionchart(YVal, YPred);
title("Deep (5-layer) DNN - Validation Results");

% Simple accuracy (percentage correct)
valAccuracy = mean(YPred == YVal);
fprintf('Validation accuracy: %.2f %%\n', 100*valAccuracy);

%
%% 10) Local helper function (data splitting)

function [XTrain, YTrain, XVal, YVal] = trainTestSplit(X, Y, trainRatio)
% trainTestSplit  Random split into training and validation sets.
%
% INPUTS
%   X          : N x P numeric matrix (N observations, P features)
%   Y          : N x 1 labels (categorical recommended)
%   trainRatio : fraction in [0,1] to allocate to training
%
% OUTPUTS
%   XTrain, YTrain : training set
%   XVal,   YVal   : validation set
%
% NOTE ABOUT cvpartition:
%   HoldOut means "how much to hold out" (validation portion), NOT training size.

    holdoutRatio = 1 - trainRatio;

    % Create a randomized hold-out partition
    cv = cvpartition(size(X,1), "HoldOut", holdoutRatio);

    % cv.test is TRUE for validation observations
    idxVal = cv.test;

    % Training data are the ones NOT in validation
    XTrain = X(~idxVal, :);
    YTrain = Y(~idxVal, :);

    XVal   = X(idxVal, :);
    YVal   = Y(idxVal, :);
end