clear
clc
tic

layer1size = 784;%this does not change unless there are different input images
layer2size = 16;%pretty arbitrary
layer3size = 16;%also fairly arbitrary
layer4size = 10;%this does not change, there will always be 10 outputs for 0-9
loopAmount = 5;%how many times the program should go through all the data. Each loops is ~40 minutes
trainSize = 60000;%number of data points to use in training
batchSize = 100;%number of data points to use in each batch

%open necessary files
fileID = fopen('train-images.idx3-ubyte');
trainImages = fread(fileID);
fileID = fopen('train-labels.idx1-ubyte');
trainLabels = fread(fileID);
fileID = fopen('t10k-images.idx3-ubyte');
testImages = fread(fileID);
fileID = fopen('t10k-labels.idx1-ubyte');
testLabels = fread(fileID);

trainImages = trainImages(17:size(trainImages,1),:);
trainLabels = trainLabels(9:size(trainLabels,1),:);

%randomizes the training images
trainImages2 = trainImages;
randomPositions = randperm(trainSize);
for index = 1:trainSize
    trainImages(784*(index)-783:784*(index)) = ...
        trainImages2(784*randomPositions(index)-783:784*randomPositions(index));
end
%reassociates the labels with the images
trainLabels2 = trainLabels;
for index = 1:trainSize
    trainLabels(index) = trainLabels2(randomPositions(index));
end

%{
%viewable format
testOutput = zeros((trainSize-1)*28,28);
for k = 0:(trainSize-1)
    for n = 1:28
        for m = 1:28
            testOutput(n+28*k,m) = trainImages(-27+28*n+m+k*28^2);
        end
    end
end
%}
randSize = 10;

step1WeightMatrix = -randSize/2+randSize*rand(784,layer2size);
step1bias = -randSize/2+randSize*rand(1,layer2size);

step2WeightMatrix = -randSize/2+randSize*rand(layer2size,layer3size);
step2bias = -randSize/2+randSize*rand(1,layer3size);

step3WeightMatrix = -randSize/2+randSize*rand(layer3size,10);
step3bias = -randSize/2+randSize*rand(1,10);

step3WeightAdjust = zeros(size(step3WeightMatrix,1)*batchSize,size(step3WeightMatrix,2));
step3WeightAdjustAverage = zeros(size(step3WeightMatrix,1),size(step3WeightMatrix,2));
step3BiasAdjust = zeros(size(step3bias,1)*batchSize,size(step3bias,2));
step3BiasAdjustAverage = zeros(size(step3bias,1),size(step3bias,2));
step3layer3backProp = zeros(1,layer3size);

step2WeightAdjust = zeros(size(step2WeightMatrix,1)*batchSize,size(step2WeightMatrix,2));
step2WeightAdjustAverage = zeros(size(step2WeightMatrix,1),size(step2WeightMatrix,2));
step2BiasAdjust = zeros(size(step2bias,1)*batchSize,size(step2bias,2));
step2BiasAdjustAverage = zeros(size(step2bias,1),size(step2bias,2));
step2layer2backProp = zeros(1,layer2size);

step1WeightAdjust = zeros(size(step1WeightMatrix,1)*batchSize,size(step1WeightMatrix,2));
step1WeightAdjustAverage = zeros(size(step1WeightMatrix,1),size(step1WeightMatrix,2));
step1BiasAdjust = zeros(size(step1bias,1)*batchSize,size(step1bias,2));
step1BiasAdjustAverage = zeros(size(step1bias,1),size(step1bias,2));

cost = zeros(1,10);
costTotal = 0;
resultsArray = zeros(loopAmount*trainSize,2);

for n = 1:loopAmount
    for trainingIndex=0:(trainSize/batchSize)-1
        for batchIndex = 0:batchSize-1

            currentData(1,:) = trainImages(1+(batchIndex + trainingIndex*batchSize)*784:...
                ((batchIndex + trainingIndex*batchSize)+1)*784);%makes currentData the next up data set
            currentData(1,:) = currentData(1,:)/256;%scales data down
            currentLayer2 = currentData*step1WeightMatrix - step1bias; %applies weights and biases
            currentLayer2 = 1./(1+exp(-currentLayer2)); %makes values 0-1 (sigmoid/logistic function)

            currentLayer3 = currentLayer2*step2WeightMatrix - step2bias; %applies weights and biases
            currentLayer3 = 1./(1+exp(-currentLayer3)); %makes values 0-1

            currentLayer4 = currentLayer3*step3WeightMatrix - step3bias; %applies weights and biases
            currentLayer4 = 1./(1+exp(-currentLayer4)); %makes values 0-1

            costMatrix = zeros(1,10); %find the cost of the image
            costMatrix(1,trainLabels((batchIndex + trainingIndex*batchSize)+1)+1) = 1;
            cost(1,:) = (currentLayer4(1,:)-costMatrix(1,:)).^2;
            costTotal = sum(cost);

            %step 3 training
            %weight adjustment
            for columns = batchIndex*layer3size+1:batchIndex*layer3size+layer3size %1-15
                for rows = 1:layer4size%1-10
                    step3WeightAdjust(columns,rows) = currentLayer3(1,columns-batchIndex*layer3size)...
                        *sigmoidDerivative(currentLayer3*step3WeightMatrix(:,rows)+step3bias(rows))...
                        *2*(currentLayer4(1,rows)-costMatrix(1,rows));
                end
            end
            %bias adjustment
            for rows = 1:layer4size%1-10
                step3BiasAdjust(batchIndex+1,rows) = sigmoidDerivative(currentLayer3*step3WeightMatrix(:,rows)+step3bias(rows))...
                    *2*(currentLayer4(1,rows)-costMatrix(1,rows));
            end
            %value adjustment
            for columns = 1:layer3size
                for sums = 1:layer4size
                    step3layer3backProp(1,columns) = step3layer3backProp(1,columns)+ ...
                        step3WeightMatrix(columns,sums)*sigmoidDerivative(currentLayer3*step3WeightMatrix(:,sums)+step3bias(sums))...
                        *2*(currentLayer4(1,sums)-costMatrix(1,sums));
                end
            end
            %step3layer3backProp = -step3layer3backProp;
            %step 2 training
            %weight adjustment
            for columns = batchIndex*layer2size+1:batchIndex*layer2size+layer2size %1-20
                for rows = 1:layer3size %1-15
                    step2WeightAdjust(columns,rows) = currentLayer2(1,columns-batchIndex*layer2size)...
                        *sigmoidDerivative(currentLayer2*step2WeightMatrix(:,rows)+step2bias(rows))...
                        *step3layer3backProp(1,rows);
                end
            end
            %bias adjustment
            for rows = 1:layer3size%1-15
                step2BiasAdjust(batchIndex+1,rows) = sigmoidDerivative(currentLayer2...
                    *step2WeightMatrix(:,rows)+step2bias(rows))*step3layer3backProp(1,rows);
            end
            %value adjustment
            for columns = 1:layer2size%1-20
                for sums = 1:layer3size%1-15
                    step2layer2backProp(1,columns) = step2layer2backProp(1,columns)+ ...
                        step2WeightMatrix(columns,sums)*sigmoidDerivative(currentLayer2*step2WeightMatrix(:,sums)+step2bias(sums))...
                        *step3layer3backProp(1,sums);
                end
            end
            %step2layer2backProp = -step2layer2backProp;
            %step 1 training
            %weight adjustment
            for columns = batchIndex*layer1size+1:batchIndex*layer1size+layer1size %1-784
                for rows = 1:layer2size %1-20
                    step1WeightAdjust(columns,rows) = currentData(1,columns-batchIndex*layer1size)...
                        *sigmoidDerivative(currentData*step1WeightMatrix(:,rows)+step1bias(rows))...
                        *step2layer2backProp(1,rows);
                end
            end
            %bias ajustment
            for rows = 1:layer2size%1-20
                step1BiasAdjust(batchIndex+1,rows) = sigmoidDerivative(currentData...
                    *step1WeightMatrix(:,rows)+step1bias(rows))*step2layer2backProp(1,rows);
            end
            %update results array
            a= max(costMatrix);
            b= max(currentLayer4);
            if find(costMatrix == a) == find(currentLayer4 == b)
                resultsArray((n-1)*trainSize+trainingIndex*batchSize+batchIndex+1,1) = 1;
            end
            resultsArray((n-1)*trainSize+trainingIndex*batchSize+batchIndex+1,2) = sum(cost);
        end

        %average weightadjust, biasadjust step 1
        for index = 1:size(step1WeightAdjustAverage,1)
            for matrixIterator = 0:batchSize-1
                step1WeightAdjustAverage(index,:) = step1WeightAdjustAverage(index,:)+...
                    step1WeightAdjust(index+matrixIterator*size(step1WeightAdjustAverage,1),:);
            end
        end
        step1WeightAdjustAverage = step1WeightAdjustAverage./batchSize;
        for index = 1:size(step1BiasAdjustAverage,1)
            for matrixIterator = 0:batchSize-1
                step1BiasAdjustAverage(index,:) = step1BiasAdjustAverage(index,:)+...
                    step1BiasAdjust(index+matrixIterator*size(step1BiasAdjustAverage,1),:);
            end
        end
        step1BiasAdjustAverage = step1BiasAdjustAverage./batchSize;
        %average weightadjust, biasadjust step 2
        for index = 1:size(step2WeightAdjustAverage,1)
            for matrixIterator = 0:batchSize-1
                step2WeightAdjustAverage(index,:) = step2WeightAdjustAverage(index,:)+...
                    step2WeightAdjust(index+matrixIterator*size(step2WeightAdjustAverage,1),:);
            end
        end
        step2WeightAdjustAverage = step2WeightAdjustAverage./batchSize;
        for index = 1:size(step2BiasAdjustAverage,1)
            for matrixIterator = 0:batchSize-1
                step2BiasAdjustAverage(index,:) = step2BiasAdjustAverage(index,:)+...
                    step2BiasAdjust(index+matrixIterator*size(step2BiasAdjustAverage,1),:);
            end
        end
        step2BiasAdjustAverage = step2BiasAdjustAverage./batchSize;
        %average weightadjust, biasadjust step 3
        for index = 1:size(step3WeightAdjustAverage,1)
            for matrixIterator = 0:batchSize-1
                step3WeightAdjustAverage(index,:) = step3WeightAdjustAverage(index,:)+...
                    step3WeightAdjust(index+matrixIterator*size(step3WeightAdjustAverage,1),:);
            end
        end
        step3WeightAdjustAverage = step3WeightAdjustAverage./batchSize;
        for index = 1:size(step3BiasAdjustAverage,1)
            for matrixIterator = 0:batchSize-1
                step3BiasAdjustAverage(index,:) = step3BiasAdjustAverage(index,:)+...
                    step3BiasAdjust(index+matrixIterator*size(step3BiasAdjustAverage,1),:);
            end
        end
        step3BiasAdjustAverage = step3BiasAdjustAverage./batchSize;
        %update biases, weights
        step1bias = step1bias - step1BiasAdjustAverage;
        step1WeightMatrix = step1WeightMatrix - step1WeightAdjustAverage;

        step2bias = step2bias - step2BiasAdjustAverage;
        step2WeightMatrix = step2WeightMatrix - step2WeightAdjustAverage;

        step3bias = step3bias - step3BiasAdjustAverage;
        step3WeightMatrix = step3WeightMatrix - step3WeightAdjustAverage;

    end
end
toc

function y = sigmoidDerivative(x)
z = exp(-x)./(1+exp(-x)).^2;
if isnan(z)
    y = 0;%arbitrary fix if 0/0 (nan)
elseif z < 10^(-7) && z > -10^(-7)
    y = 0;
else
    y = z;
end
end