% ***** REGULARIZED GRADIENT DESCENT QUADRATIC REGRESSION ***** %

% m:  Number of data samples
% n:  Number of factors or dependent variables
% p:  Index value of a single sample within the set of m samples
% q:  Index value of an individual feature within the
%     ith sample's feature (or variable) vector
% x:  Feature vector with n+1 features, i.e. x=[x0 x1 x2 x3] & x0=1
% theta:  Parameter vector with at least n+1 parameters,
%         i.e. theta=[theta0 theta1 theta2 theta3]
% h:  "Hypothesis" or prediction model, i.e. h=theta*transpose(x)
% J:  Cost function dependent on the parameter values in vector theta

close all
clear all


% Data Points
x1=[-2.4  8.2  1.9  12.1  -2.1  9.9  -9  -1   9.2  1];
x2=[ 4.6  2.9  2.7  3.2    4.0  3.4  10   8  -3.4  2];
y=[  0    1     0    1      0    1    0   1   1    0];


% Initialize parameters
theta1(1)=-1.5; % This value changes each iteration
theta2(1)=1.5; % This value changes each iteration
theta3(1)=1.1; % This value changes each iteration

% Initialize "momentum", alpha
alpha=[0.4 0.2 0.2];  % These values do not change during interation

m=length(y);  % Number of data points

g=theta1(1)+theta2(1).*x1+theta3(1).*(x2);
h=1./(1+exp(-g));  % Hypothesized model

%  J(1)=(1/m)*sum((h-y).^2);  % Cost function
J(1)=-(1/m)*sum(y.*log(h)+(1-y).*log(1-h));

e=1;
count=2;  % count is the number of iterations that have occured plus one
lambda=0;  % Regularization constant
frameCount=1;

while (e>0.0001)
    
    theta1(count)=theta1(count-1)-alpha(1).*(2/m).*sum(h-y);
    theta2(count)=theta2(count-1)-alpha(2).*(2/m).*sum((h-y).*x1);
    theta3(count)=theta3(count-1)-alpha(3).*(2/m).*sum((h-y).*x2);
    g=theta1(count)+theta2(count).*x1+theta3(count).*(x2);
    h=1./(1+exp(-g));
    J(count)=-(1/m)*(sum(y.*log(h)+(1-y).*log(1-h))+lambda.*(theta1(count)^2+theta2(count)^2+theta3(count)^2));
    e=abs(J(count)-J(count-1));
    
    figure(1)
    x1_g=[min(x1):((max(x1)-min(x1))/20):max(x1)];
    x2_g=[min(x2):((max(x2)-min(x2))/20):max(x2)];
    [a,b]=meshgrid(x1_g,x2_g);
    g_plot=theta1(count)+theta2(count).*a+theta3(count).*b;
    
    subplot(2,1,1)
    surf(a,b,g_plot)
    title('Decision Boundary Plane','fontsize',11,'fontweight','b')
    axis([-18 18 -18 18 -18 18])
    hold on
    
    for n=1:m
        if y(n)==1
            plot3(x1(n),x2(n),y(n),'ro','MarkerFaceColor','red')
            hold on
        else
            plot3(x1(n),x2(n),y(n),'bo','MarkerFaceColor','blue')
            hold on
        end
        grid on
        hold on
    end

hold off
    
subplot(2,1,2)
semilogx(1:count,J)
title('Cost Function','fontsize',11,'fontweight','b')
xlabel('Epoch Number','fontsize',10)
ylabel('Cost Value','fontsize',10)

axis([1 1000 0 0.9])
grid on

 
% subplot(3,1,3)
% plot3(theta1,theta2,J)
% axis([-7 3 -7 3 -1.2 1.2])
% grid on

    set(gcf,'Renderer','zbuffer')
    M(frameCount)= getframe(gcf);
    frameCount=frameCount+1;

pause(0.1)
hold off    
count=count+1;
    
end

%movie(M)
movie2avi(M,'PlaneRegressionLearner.avi')%,'compression','Cinepak')