% Logistic Regression via Gradient Descent
close all; clear all; clc;

%data = load('exam_scores_train2.txt')'; % load regression data
data = load('exam_scores_training.txt')'; % load regression data


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

y = data(3,:);
M = length(data(1,:));
X = ones(3,M); X(2,:) = data(1,:); X(3,:) = data(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

g = @(z) 1./(1 + exp(-z)); % sigmoid function
theta = zeros(3,1); % intialize theta vector
theta(1) = -50; theta(2) = 0.1; theta(3)= 5; % initial guess for theta's
%theta(1) = 0; theta(2) = 0; theta(3)= 0; % initial guess for theta's

%h_theta = @(theta, x) g(theta'*x); % define h_theta
h_theta = @(Theta, xx) g(Theta'*xx); % define h_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);

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]);
xlabel(xl,'FontWeight','bold'); 
ylabel(yl,'FontWeight','bold');
title(ftitle1)
legend(l1,l2)
pause


Alpha = 2e-4; % learning rate, time-step for gradient descent
% scale alpha to number of training examples

%return 
N_max = 50;
%eta = 0

for j = 1:N_max
    for i = 1:M    
    
    %A = y - h_theta(theta,x)';
    %rh_sum = 0;
    %for i=1:m
    %    rh_sum = ( y(i) - h_theta(theta,x(:,i)) )*x(:,i);
    %end
    %theta_prev = theta;
    %theta = theta + alpha*(rh_sum); %batch
    %theta = theta + alpha*( (y(i) - h_theta(theta,x(:,i)))*x(:,i) );
    
    theta(1) = theta(1) + Alpha*( (y(i) - h_theta(theta, X(:,i)))*X(1,i));
    theta(2) = theta(2) + Alpha*( (y(i) - h_theta(theta, X(:,i)))*X(2,i));
    theta(3) = theta(3) + Alpha*( (y(i) - h_theta(theta, X(:,i)))*X(3,i));
    
    %i 
    %theta
    %norm(theta_prev-theta,'fro');
    end
    if rem(j,1) == 0
        theta 
        sprintf('iteration # %d',j)
        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; 
        %x2m = -(1/theta(3))*(theta(1) + theta(2)*x1m);
        %x2M = -(1/theta(3))*(theta(1) + theta(2)*x1M);
        db = @(x) (1/theta(3))*(-(theta(2)*x) - theta(1));  
        plot(xrange_db, db(xrange_db),'g-','LineWidth',2); hold off; hold off;
        %plot([x1m x1M], [x2m x2M],'g-','LineWidth',2); hold off; hold off;
        axis([x1ma x1Ma x2ma x2Ma]);
        xlabel(xl,'FontWeight','bold'); 
        ylabel(yl,'FontWeight','bold');
        title(ftitle1);
        legend(l1,l2);
        pause(0.25)
       
    end
end

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; 
%x2m = -(1/theta(3))*(theta(1) + theta(2)*x1m);
%x2M = -(1/theta(3))*(theta(1) + theta(2)*x1M);
%plot([x1m x1M], [x2m x2M],'g-','LineWidth',2);
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(M,1); %initialize vect for predicted class label
g = @(z) 1./(1 + exp(-z)); % sigmoid function

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

%
fprintf('----------------------------------------------------------------\n');
fprintf('Stochastic Gradient Descent 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')/M)*100 );
fprintf('----------------------------------------------------------------\n');

pause

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

data = load('exam_scores_testing.txt')'; % load regression testing data
%data = load('exam_scores.txt')'; % load regression testing data


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

y = data(3,:);
m = length(data(1,:));
X = ones(3,m); X(2,:) = data(1,:); X(3,:) = data(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 = 'Exam 1'; %xlabel
yl = 'Exam 2'; %ylabel
ftitle = 'Exam Scores';
l1 = 'passed'; %class label
l2 = 'not passed'; %ditto

figure(6);
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;


figure(7);
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(M,1);

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

for i = 1:M
    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')/M)*100 );
fprintf('----------------------------------------------------------------\n');

pause;

figure(8);
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, 'passed (training)', 'not passed (training)');




