% Logistic Regression 2 label, M features
% Uses the Conjugate Gradient Method to Optimize the 
% Objective Function

close all; clear all; clc;


% Load in your data from the MNIST data set
% You will need to have 
% data_train_3vs5.mat needs to have the matrix data_train
% where size(data_train) must be [785, 11552]
% each column corresponds to a vector that is actually a reshaped image
% of size 28x28. The last entry in each column is the class label
% label 1 for a handwritten 3 and 0 for a handwritten 5.
% same holds for data_test_3vs5.mat which must contain a matrix
% data_test that has size(data_test) being [785, 1902]
% It is up to you to create and save these data-sets as described in the
% final exam. Make sure the format is correct before running.
% You are training on over 11 thousand images and testing on nearly 2
% thousand. The training step alone takes several minutes to run.

load data_train_3vs5.mat
load data_test_3vs5.mat

[M,N] = size(data_train); %extract size of data_train

%random shuffle indices
rand_ind = randperm(N);
r_ind = rand_ind(1:25);

%# of Conjugate Gradeint Iterations
NCG = 50;

%store randomly shuffled training data as IM
IM = data_train(1:M-1,:);

%view 25 randomly selected training images
figure(1)
for i=1:25
    im = reshape(IM(:,r_ind(i)),[28,28])';
    subplot(5,5,i); imagesc(im); colormap(gray(256));
    title('training images')
    axis off;
end

fprintf('Preparing to run optimization... \n');
pause(2);
fprintf('Optimizing Now! There will be %0.0d iterations. Be Patient...\n',NCG);

tic

y = data_train(M,:);
X = ones(M,N); 

%store data into X matrix.
X(2:M,:) = data_train(1:M-1,:); 

% initial theta for fminunc
initial_theta = zeros(M,1);

% Set options for fmincg
options = optimset('GradObj', 'on', 'MaxIter', NCG);

% Use unconstrained minimization function: fmincg
% to find approx. theta that minimize the 
% Logistic Regression (LR) Objective Energy
% This function will return theta and the Energy at theta
% The method used is the Conjugate Gradient Method
[theta, energy] = fmincg(@(t)Obj_fun(t,X,y), initial_theta, options);

% Output theta and LR NRG 
%fprintf('LR Energy at theta via Conjugate Gradient: %0.6f \n', energy);
%fprintf('theta: \n');
%fprintf(' %0.6f \n', theta');

p = zeros(N,1); %initialize vect for predicted class label
g = @(z) 1./(1 + exp(-z)); % sigmoid function

for i = 1:N
    p(i) = g(theta'*X(:,i)) > 1/2; % calculate predictions based on decision boundary
end

%output that method finished
fprintf('----------------------------------------------------------------\n');
fprintf('Minimization Method Finished \n')
%fprintf('Values of theta = \n' );
%fprintf(' %0.6f \n', theta);
fprintf('----------------------------------------------------------------\n');

%Compute accuracy on our training set
fprintf('----------------------------------------------------------------\n');
fprintf('Training Accuracy: %.2f %%\n', (sum(p == y')/N)*100 );
fprintf('----------------------------------------------------------------\n');

toc %timer

pause(2)

%%%%%%%%%%%%%%%%%%%%%%%%%% testing phase %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%store training data
X_train = X;
y_train = y;

[M,N] = size(data_test);

y = data_test(M,:);
X = ones(M,N); 

X(2:M,:) = data_test(1:M-1,:); 
y_pos = find(y==1); % find indices where label = 1
y_neg = find(y==0); % ditto for 0

p = zeros(N,1);

for i = 1:N
    p(i) = g(theta'*X(:,i)) > 1/2; % calculate predictions based on decision boundary
end

% Compute accuracy on our testing set
fprintf('----------------------------------------------------------------\n');
fprintf('Testing Accuracy: %.2f %%\n', (sum(p == y')/N)*100 );
fprintf('----------------------------------------------------------------\n');

fprintf('----------------------------------------------------------------\n');
fprintf('Number of Mislabeled: %.0f \n', sum(p ~= y'));
fprintf('Number of Correctly labeled: %.0f \n', sum(p == y'));
fprintf('----------------------------------------------------------------\n');


mis_ind = find(p ~= y'); %find indices of mislabeled data
IMtest = data_test(1:M-1,:); %store data into IMtest

%view first 64 mislabeled data points. 
figure(2)
for i=1:64
    imtest = reshape(IMtest(:,mis_ind(i)),28,28)'; %view 
    subplot(8,8,i); imagesc(imtest); colormap(gray(256));
    title('mislabeled')
    axis off;
end

%find indices for correctly classified images
corr_ind = find(p == y');

%radomly permute the indices so we can view first 64 correctly
%classified images from both classes
r_ind = randperm(length(corr_ind));
corr_ind = corr_ind(r_ind);

%view 64 randomly chosen correctly classified images
figure(3)
for i=1:64
    imtest = reshape(IMtest(:,corr_ind(i)),28,28)';
    subplot(8,8,i); imagesc(imtest); colormap(gray(256));
    title('correct label')
    axis off;
end

% Objective function without gradient returned, for use with Quasi-Newton 
% function LR = Obj_fun(theta,X,y)
%     g = @(z) 1./(1 + exp(-z));
%     M = size(X,2);
%     LR = (1/M)*sum( (-1)*y.*log(g(theta'*X)) - (1-y).*log(1-g(theta'*X)) );
% end


% Objective function and Gradient of it for fmincg
function [LR, gradLR] = Obj_fun(theta,X,y)
    g = @(z) 1./(1 + exp(-z));
    N = size(X,2);
    LR = (1/N)*sum( (-1)*y.*log(g(theta'*X)) - (1-y).*log(1-g(theta'*X)) );
    gradLR = zeros(size(theta));
    len_theta = size(theta,1);
    for j = 1:len_theta
        gradLR(j,:) = (1/N)*( (g(theta'*X)-y)*(X(j,:)') );
    end
end





