% Logistic Regression 2 label, 2 features
% Uses the Conjugate Gradient Method to Optimize
% the Logistic Regression (LR) Cost Function.
close all; clear all; clc;

% load if your data is saved as a text file
%data = load('your_data.txt')'; % load training data
data = load('your_data.txt')'; % load training data

% You can also store in a .mat file as well. See exam for details on how to
% load your data in this case.

%use the randperm function to randomly permute the indices of your data
%matrix. Create a new data matrix called data_p that is the data matrix
%evaluated at the random indices. (This permutes your data around)
%Break up your data so that 2/3 of your data is for training and 
%1/3 is for testing. i.e. data_train = 2/3 of data_p and 
%data_test = 1/3 of data_p.

%here, you will create the random indices r_ind using 
%the randperm function
%then create data_p, the permuted data entries of data
%i.e. columns are permuted randomly 

%starter code:
[Md,Nd] = size(data);
r_ind = ?? %<-- you fill in
data_p = ?? %<-- you fill in

% if you have 30 data points and already calculated data_p
data_train = data_p(:,1:20); % first 20 data points for training
data_test = data_p(:,21:30); % last 10 data points for testing
%note: you can adjust the amount of data manually for training and testing

% for plotting options
ms = 7; % markersize
xl = 'feature 1'; %xlabel
yl = 'feature 2'; %ylabel
ftitle1 = 'Features: Training';
ftitle2 = 'Features: Testing';
ftitle3 = 'Features: Training and Testing';
l1 = 'label 1'; %class label
l2 = 'label 0'; %ditto

[M,N] = size(data_train);
y = data_train(M,:);
X = ones(M,N); X(2,:) = data_train(1,:); X(3,:) = data_train(2,:); 
x1 = X(2,:); x2 = X(3,:);
y_pos = find(y==1); % find indices where label = 1
y_neg = find(y==0); % ditto for 0

%axis range, 10% past boundary points
x1m = min(x1); x1M = max(x1);
dx1 = 0.10*(abs(x1M-x1m));
x2m = min(x2); x2M = max(x2);
dx2 = 0.10*(abs(x2M-x2m));
x1ma = x1m-dx1; x1Ma = x1M+dx1;
x2ma = x2m-dx2; x2Ma = x2M+dx2;
xrange_db = (x1ma:0.1:x1Ma);

figure(1);
plot(x1(y_pos),x2(y_pos),'r^','MarkerSize',ms,'MarkerFaceColor','r'); hold on; 
plot(x1(y_neg),x2(y_neg),'bo','MarkerSize',ms,'MarkerFaceColor','b'); hold off;
axis([x1ma x1Ma x2ma x2Ma]);
xlabel(xl,'FontWeight','bold'); 
ylabel(yl,'FontWeight','bold');
title(ftitle1)
legend(l1,l2)
pause(2)

g = @(z) 1./(1 + exp(-z)); % sigmoid function

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

x1m = min(X(2,:))-2;
x1M = max(X(2,:))+2;
x2m = -(1/theta(3))*(theta(1) + theta(2)*x1m);
x2M = -(1/theta(3))*(theta(1) + theta(2)*x1M);

% Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 100);


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

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

figure(2);
plot(x1(y_pos),x2(y_pos),'r^','MarkerSize',ms,'MarkerFaceColor','r'); hold on; 
plot(x1(y_neg),x2(y_neg),'bo','MarkerSize',ms,'MarkerFaceColor','b'); hold on; 
db = @(x) (1/theta(3))*(-(theta(2)*x) - theta(1));  
plot(xrange_db, db(xrange_db),'g-','LineWidth',2); hold off; hold off;
axis([x1ma x1Ma x2ma x2Ma]);hold off; hold off;
xlabel(xl,'FontWeight','bold'); 
ylabel(yl,'FontWeight','bold');
title(ftitle1);
legend(l1,l2);

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

%
fprintf('----------------------------------------------------------------\n');
fprintf('Minimization Method Finished \n')
fprintf('Values of theta'' = [%.2f, %.2f, %.2f]\n', theta(1),theta(2),theta(3));
fprintf('----------------------------------------------------------------\n');

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

pause(2)

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

%store training data for visualization
x1t = x1; x2t = x2;
post = y_pos; negt = y_neg;

y = data_test(3,:);
n = length(data_test(1,:));
X = ones(3,n); X(2,:) = data_test(1,:); X(3,:) = data_test(2,:); 
x1 = X(2,:); x2 = X(3,:);
pos = find(y==1); % find indices where label = 1
neg = find(y==0); % ditto for 0

%adjust axis to testing set. 
%uncomment if you desire this
%axis range, 10 percent past boundary points
% x1m = min(x1); x1M = max(x1);
% dx1 = 0.10*(abs(x1M-x1m));
% x2m = min(x2); x2M = max(x2);
% dx2 = 0.10*(abs(x2M-x2m));
% x1ma = x1m-dx1; x1Ma = x1M+dx1;
% x2ma = x2m-dx2; x2Ma = x2M+dx2;

% for plotting options
ms = 7; % markersize
xl = 'Feature 1'; %xlabel
yl = 'Feature 2'; %ylabel
ftitle = 'Features';
l1 = 'label 1'; %class label
l2 = 'label 0'; %ditto

figure(3);
plot(x1(pos),x2(pos),'m^','MarkerSize',ms,'MarkerFaceColor','m'); hold on;
plot(x1(neg),x2(neg),'co','MarkerSize',ms,'MarkerFaceColor','c'); hold off;
xlabel(xl,'FontWeight','bold'); 
ylabel(yl,'FontWeight','bold');
title(ftitle2); axis([x1ma x1Ma x2ma x2Ma]); % axis([xmin xmax ymin ymax])
legend(l1, l2);
pause(2);


figure(4);
plot(x1(pos),x2(pos),'m^','MarkerSize',ms,'MarkerFaceColor','m'); hold on;
plot(x1(neg),x2(neg),'co','MarkerSize',ms,'MarkerFaceColor','c'); hold on;
plot(xrange_db, db(xrange_db), 'g-', 'LineWidth',2);
xlabel(xl,'FontWeight','bold'); 
ylabel(yl,'FontWeight','bold');
title(ftitle2); axis([x1ma x1Ma x2ma x2Ma]); % axis([xmin xmax ymin ymax])
hold off; 
hold off;
legend(l1, l2);

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');

pause(2);

figure(5);
plot(x1(pos),x2(pos),'m^','MarkerSize',ms,'MarkerFaceColor','m'); hold on;
plot(x1(neg),x2(neg),'co','MarkerSize',ms,'MarkerFaceColor','c'); hold on;
plot(x1t(post),x2t(post),'r^','MarkerSize',ms,'MarkerFaceColor','r'); hold on;
plot(x1t(negt),x2t(negt),'bo','MarkerSize',ms,'MarkerFaceColor','b'); hold on;
plot(xrange_db, db(xrange_db), 'g-', 'LineWidth',2);
xlabel(xl,'FontWeight','bold'); 
ylabel(yl,'FontWeight','bold');
title(ftitle3); axis([x1ma x1Ma x2ma x2Ma]); % axis([xmin xmax ymin ymax])
hold off; hold off; hold off; hold off;
legend(l1, l2, 'label 1 (training)', 'label 0 (training)');

% 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 fminunc
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





