%%======================================================================
%% STEP 0: Initialise constants and parameters
%
%  Here we define and initialise some constants which allow your code
%  to be used more generally on any arbitrary input. 
%  We also initialise some parameters used for tuning the model.

inputSize = 28 * 28+1; % Size of input vector (MNIST images are 28x28)
numClasses = 2;     % Number of classes (MNIST images fall into 10 classes)

% lambda = 1e-4; % Weight decay parameter

%%======================================================================
%% STEP 1: Load data
%
%  In this section, we load the input and output data.
%  For softmax regression on MNIST pixels, 
%  the input data is the images, and 
%  the output data is the labels.
%

% Change the filenames if you've saved the files under different names
% On some platforms, the files might be saved as 
% train-images.idx3-ubyte / train-labels.idx1-ubyte

images = loadMNISTImages('mnist/train-images-idx3-ubyte');
labels = loadMNISTLabels('mnist/train-labels-idx1-ubyte');


index=(labels==0|labels==1);

images=images(:,index);
labels=labels(index);

inputData = [images;ones(1,size(images,2))];


% Randomly initialise theta


%%======================================================================
%% STEP 2: Implement softmaxCost
%
%  Implement softmaxCost in softmaxCost.m. 

% [cost, grad] = logisticCost(theta, inputSize,inputData, labels);
                                     


%%======================================================================
%% STEP 4: Learning parameters
%
%  Once you have verified that your gradients are correct, 
%  you can start training your softmax regression code using softmaxTrain
%  (which uses minFunc).

options.maxIter = 100;
options.alpha = 0.1;
options.method = 'Grad';
theta = logisticTrain( inputData, labels,options);
                          
% Although we only use 100 iterations here to train a classifier for the 
% MNIST data set, in practice, training for more iterations is usually
% beneficial.

%%======================================================================
%% STEP 5: Testing
%
%  You should now test your model against the test images.
%  To do this, you will first need to write softmaxPredict
%  (in softmaxPredict.m), which should return predictions
%  given a softmax model and the input data.

images = loadMNISTImages('mnist/t10k-images-idx3-ubyte');
labels = loadMNISTLabels('mnist/t10k-labels-idx1-ubyte');

index=(labels==0|labels==1);
images=images(:,index);
labels=labels(index);

inputData = [images;ones(1,size(images,2))];

% You will have to implement softmaxPredict in softmaxPredict.m
[pred] = logisticPredict(theta, inputData);

acc = mean(labels(:) == pred(:));
fprintf('Accuracy: %0.3f%%\n', acc * 100);

% Accuracy is the proportion of correctly classified images
% After 100 iterations, the results for our implementation were:
%
% Accuracy: 92.200%
%
% If your values are too low (accuracy less than 0.91), you should check 
% your code for errors, and make sure you are training on the 
% entire data set of 60000 28x28 training images 
% (unless you modified the loading code, this should be the case)

function [modelTheta] = logisticTrain(inputData, labels,option)

if ~exist('options', 'var')
    options = struct;
end

if ~isfield(options, 'maxIter')
    options.maxIter = 400;
end

if ~isfield(options, 'method')
    options.method = 'Newton';
end

if ~isfield(options, 'alpha')
    options.method = 0.01;
end
theta = 0.005 * randn(size(inputData,1),1);
iter=1;

maxIter=option.maxIter;
alpha=option.alpha;
method=option.method;
fprintf('iter\tStep Length\n');
lastSteps=0;
while iter<=maxIter
    
    h=sigmoid(theta'*inputData);
%     cost=sum(labels'.*log(h)+(1-labels').*log(1-h),2)/size(inputData,2);
    grad=inputData*(labels'-h)';
    
    if strcmp(method,'Grad')>0
        steps=alpha.*grad;
%     else
%         H = inputData* diag(h) * diag(1-h) * inputData';
%         steps=-alpha.*H\grad;
    end
    theta=theta+steps;
    stepLength=sum(steps.^2)/size(steps,1);
    fprintf('%d\t%f\n',iter,stepLength);
    if abs(stepLength)<1e-9
        break;
    end
    iter=iter+1;
end
modelTheta=theta;

    function z=sigmoid(x)
        z=1./(1+exp(-1.*x));
    end
end

function [pred] = logisticPredict(theta, data)

% softmaxModel - model trained using softmaxTrain
% data - the N x M input matrix, where each column data(:, i) corresponds to
%        a single test set
%
% Your code should produce the prediction matrix 
% pred, where pred(i) is argmax_c P(y(c) | x(i)).

%% ---------- YOUR CODE HERE --------------------------------------
%  Instructions: Compute pred using theta assuming that the labels start 
pred=theta'*data>0.5;
% ---------------------------------------------------------------------
end