function detectHandWriting_01

% Import MNIST data
[X,Y] = importTrainingData(60000);
[Xtest,Ytest] = importTestingData(10000);

% make sure everything looks okay
Xsample = reshape(Xtest(:,1),20,20);
figure('Position',[38   239   729   629])
imagesc(Xsample)
colormap gray

% Training Notes:
% We train the network in stages.
% During stage k, we perform N cycles of gradient descent based on test 
% data from index i1 to index i2, where [i1 i2 N] = testSchedule(k,:)

nodeNums = [400 15 16 10]; % number of nodes in each layer
alpha = 0.002;             % learning rate
iend = 500:500:size(X,2);
istart = 1 + [0 iend(1:end-1)];
testSchedule = [istart' iend' 100*ones(size(istart'))];
testSchedule([1 2],3) = 500;
% results in a final performance score of 90.72%

nodeNums = [400 15 16 10];
alpha = 0.002; % learning rate
iend = 1000:1000:size(X,2);
istart = 1 + [0 iend(1:end-1)];
testSchedule = [istart' iend' 150*ones(size(istart'))];
testSchedule([1 2],3) = 500;
% results in a final performance score of 91.64% (temp high of 92.5%)

nodeNums = [400 15 16 10];
alpha = 0.002; % learning rate
iend = 1500:1500:size(X,2);
istart = 1 + [0 iend(1:end-1)];
testSchedule = [istart' iend' 200*ones(size(istart'))];
testSchedule([1 2],3) = 500;
% numberOfParameters: 6441
% results in a final performance score of 92.18%

nodeNums = [400 20 25 10];
alpha = 0.002; % learning rate
iend = 1500:1500:size(X,2);
istart = 1 + [0 iend(1:end-1)];
testSchedule = [istart' iend' 200*ones(size(istart'))];
testSchedule([1 2],3) = 500;
% numberOfParameters: 8805
% results in a final performance score of 93.3%

nodeNums = [400 30 25 10];
alpha = 0.002; % learning rate
iend = 1500:1500:size(X,2);
istart = 1 + [0 iend(1:end-1)];
testSchedule = [istart' iend' 200*ones(size(istart'))];
testSchedule([1 2],3) = 500;
% numberOfParameters: 13065
% results in a final performance score of 94.5%

[A,b] = trainModel(X,Y,Xtest,Ytest,nodeNums,alpha,testSchedule);

% save('PKneuralnet20240519.mat','A','b')

function [A,b] = trainModel(X,Y,Xtest,Ytest,nodeNums,alpha,testSchedule)
%
%input:
% nodeNums is 1xn for a network with n layers, with nodeNums(i) the number
%          of neurons in layer i. With ni and nf the initial and final
%          elements of nodeNums,
% X is an nixM array of input training data
% Y is an nfxM array of output training data
% Xtest is an nixM array of input testing data
% Ytest is an nfxM array of output testing data
%
%ouput:
% A 1x(n-1) cell array with A{i} a double array of weights
% b 1x(n-1) cell array with B{i} a double array of biases
%
% The different elements in A (and in b) generally have different sizes

% initialize A and b
rng(1)
n = numel(nodeNums)-1;
A = cell(1,n);
b = cell(1,n);
numberOfParameters = 0;
for k = 1:n
    A{k} = rand(nodeNums(k+1),nodeNums(k))-0.5;
    b{k} = rand(nodeNums(k+1),1)-0.5;
    numberOfParameters = numberOfParameters + numel(A{k}) + numel(b{k});
end
disp(['numberOfParameters: ' num2str(numberOfParameters)])

ntotal = sum(testSchedule(:,3));
F = NaN(1,ntotal);
Fd = NaN(10,ntotal);
G = NaN(1,ntotal);
count = 1;

for j = 1:size(testSchedule,1)

    istart =  testSchedule(j,1);
    iend =    testSchedule(j,2);
    I = istart:iend;
    for k = 1:testSchedule(j,3)

        [ftotal,fdigits] = performanceCheck(A,b,Xtest,Ytest);
        F(count) = ftotal;
        Fd(:,count) = fdigits;
        G(count) = costFunction(A,b,X(:,I),Y(:,I));
        count = count + 1;
        
        % evaluate the gradient
        [gradA,gradb] = evaluateGradient(A,b,X(:,I),Y(:,I));

        % adjust A and b in the direction opposite the gradient
        for i = 1:numel(A)
            A{i} = A{i} - alpha*gradA{i};
            b{i} = b{i} - alpha*gradb{i};
        end
    end
    disp([count-1 F(count-1)])

end

figure('Position',[137 4 1007 911],'Color',[1 1 1])
line('xdata',1:numel(F),'ydata',F,'color',[1 0 0])
grid on
title(F(count-1))

figure('Position',[137 4 1007 911],'Color',[1 1 1])
for k = 1:10
    subplot(10,1,k)
    Fdk = Fd(k,:);
    line('xdata',1:numel(Fdk),'ydata',Fdk,'color',[1 0 0])
    title(Fdk(end))
end
grid on

function [ftotal,fdigits] = performanceCheck(A,b,X,Y)
% return fraction of test samples that are predicted correctly by the model
% do this for each 

totalRight = zeros(10,1);
totalAttempted = zeros(10,1);
for i = 1:size(X,2)

    u = X(:,i);
    for k = 1:numel(A)
        u = 1./(1+exp(-(A{k}*u+b{k})));
    end
    [~,imax] = max(u);

    j = find(Y(:,i)); % assigned index
    totalAttempted(j) = totalAttempted(j) + 1;
    if j == imax
        totalRight(j) = totalRight(j) + 1;
    end
end

ftotal = sum(totalRight)/sum(totalAttempted);
fdigits = totalRight./totalAttempted;

function G = costFunction(A,b,X,Y)
% return the cost function associated with the model (A,b) and the training
% data (X,Y)

G = 0;
for i = 1:size(X,2)

    u = X(:,i);
    for k = 1:numel(A)
        u = 1./(1+exp(-(A{k}*u+b{k})));
    end
    e = u - Y(:,i);
    G = G + e'*e;

end

function [gradA,gradb] = evaluateGradient(A,b,X,Y)
% evaluage the gradient based on the given training data (X,Y)

% initialize gradA and gradb
gradA = cell(size(A));
gradb = cell(size(b));
for i = 1:numel(A)
    gradA{i} = zeros(size(A{i}));
    gradb{i} = zeros(size(b{i}));
end

for i = 1:size(X,2) % for each set of training data

    % determine u and w at each layer
    ui = {X(:,i)};
    wi = {[]};
    for k = 1:numel(A)
        ui{k+1} = 1./(1+exp(-(A{k}*ui{k}+b{k})));
        wi{k+1} = ui{k+1}.*(1-ui{k+1});
    end

    % initialize g
    g = 2*(ui{end}-Y(:,i)).*wi{end};

    % upate gradient arrays
    for k = numel(A):-1:1
        gradA{k} = gradA{k} + g*ui{k}';
        gradb{k} = gradb{k} + g;
        if k==1
            break;
        end
        g = (g'*A{k})'.*wi{k};
    end

end

function [X,Y] = importTestingData(M)
% M<=10000 is the number of data sets to import.

if M<2
    M=2;
elseif M>10000
    M=10000;
end

% testing data (10k images)
imageFile = 'MNIST/t10k-images.idx3-ubyte';
labelFile = 'MNIST/t10k-labels.idx1-ubyte';
istart = 1; % index (starting from 1) of initial image
iend = M;   % index (starting from 1) of final image
[Xraw,Yraw] = readMNIST(imageFile, labelFile, iend, istart - 1);
% Xraw is 20x20xM
% Yraw is Mx1

% build X
% Xraw is 20x20xM, with Xraw(:,:,i) the ith image
X = reshape(Xraw(:),20*20,M); % X is 400xM, each image a column vector
% We can recover Xraw(:,:,i) with reshape(X(:,i),20,20) 

% build Y
% Yraw is an Mx1 array of integers (0 through 9)
% Y is a 10xM array of doubles, with values 0 or 1
% Y(:,i) is zeros except for element Yraw(i)+1, which is one
Y = zeros(10,M);
I = Yraw + 1 + 10*(0:(M-1))';
Y(I) = 1;

function [X,Y] = importTrainingData(M)
% M<=60000 is the number of data sets to import.

if M<2
    M=2;
elseif M>60000
    M=60000;
end

% import training data (60k images)
imageFile = 'MNIST/train-images.idx3-ubyte';
labelFile = 'MNIST/train-labels.idx1-ubyte';
istart = 1; % index (starting from 1) of initial image
iend = M;   % index (starting from 1) of final image
[Xraw,Yraw] = readMNIST(imageFile, labelFile, iend, istart - 1);
% Xraw is 20x20xM
% Yraw is Mx1

% build X
% Xraw is 20x20xM, with Xraw(:,:,i) the ith image
X = reshape(Xraw(:),20*20,M); % X is 400xM, each image a column vector
% We can recover Xraw(:,:,i) with reshape(X(:,i),20,20) 

% build Y
% Yraw is an Mx1 array of integers (0 through 9)
% Y is a 10xM array of doubles, with values 0 or 1
% Y(:,i) is zeros except for element Yraw(i)+1, which is one
Y = zeros(10,M);
I = Yraw + 1 + 10*(0:(M-1))';
Y(I) = 1;

function [imgs,labels] = readMNIST(imgFile, labelFile, readDigits, offset)
% readMNIST by Siddharth Hegde
%
% Description:
% Read digits and labels from raw MNIST data files
% File format as specified on http://yann.lecun.com/exdb/mnist/
% Note: The 4 pixel padding around the digits will be remove
%       Pixel values will be normalised to the [0...1] range
%
% Usage:
% [imgs labels] = readMNIST(imgFile, labelFile, readDigits, offset)
%
% Parameters:
% imgFile = name of the image file
% labelFile = name of the label file
% readDigits = number of digits to be read
% offset = skips the first offset number of digits before reading starts
%
% Returns:
% imgs = 20 x 20 x readDigits sized matrix of digits
% labels = readDigits x 1 matrix containing labels for each digit
%

% Read digits
fid = fopen(imgFile, 'r', 'b');
header = fread(fid, 1, 'int32');
if header ~= 2051
    error('Invalid image file header');
end
count = fread(fid, 1, 'int32');
if count < readDigits+offset
    error('Trying to read too many digits');
end

h = fread(fid, 1, 'int32');
w = fread(fid, 1, 'int32');

if offset > 0
    fseek(fid, w*h*offset, 'cof');
end

imgs = zeros([h w readDigits]);

for i=1:readDigits
    for y=1:h
        imgs(y,:,i) = fread(fid, w, 'uint8');
    end
end

fclose(fid);

% Read digit labels
fid = fopen(labelFile, 'r', 'b');
header = fread(fid, 1, 'int32');
if header ~= 2049
    error('Invalid label file header');
end
count = fread(fid, 1, 'int32');
if count < readDigits+offset
    error('Trying to read too many digits');
end

if offset > 0
    fseek(fid, offset, 'cof');
end

labels = fread(fid, readDigits, 'uint8');
fclose(fid);

% Calc avg digit and count
imgs = trimDigits(imgs, 4);
imgs = normalizePixValue(imgs);

function digits = trimDigits(digitsIn, border)
dSize = size(digitsIn);
digits = zeros([dSize(1)-(border*2) dSize(2)-(border*2) dSize(3)]);
for i=1:dSize(3)
    digits(:,:,i) = digitsIn(border+1:dSize(1)-border, border+1:dSize(2)-border, i);
end

function digits = normalizePixValue(digits)
digits = double(digits);
for i=1:size(digits, 3)
    digits(:,:,i) = digits(:,:,i)./255.0;
end
