function filename = logisticRegression(varagin)
%
% trains the logistic regression classifier
% varagin = number of iterations
%
% logisticRegression takes a variable integer input that determines how 
% "accurate" the classifier will be. Larger numbers means the function runs
% longer. Practically, the default value will get you pretty close to the 
% best one. Your answer will change a lot with smaller numbers, but once
% it gets larger than, say, 20, you  get pretty much the same result every 
% time.
% 
% "t" (short for theta) is the output of logisticRegression. It's a row 
% vector that contains the parameters for the classifier. (It's also the 
% row vector that contains the equation for the separating line / hyperplane.)
%
% The equation for the separating line is 
% t(1) + t(2)*dist + t(3) *freq + t(4)*orientation + t(5)*dEab (or spatial dEab) = 0.
%
% See also CLASSIFYLR, PREPAREDATA
%

%% Check type
if (nargin > 1)
    error('Too many input arguments.');
elseif (nargin > 0)
    if(isnumeric(varagin)) num_iter = round(varagin);
    else error('You must specify an integer number of iterations!'); end
else num_iter = 20; end

%% Prepare the data
[xe xs xh y] = prepareData('LR');

norm_xe = [];
norm_xs = [];
norm_xh = [];

for M = 1:size(xe,1),
    norm_xe = [norm_xe; norm(xe(M,:))];
    norm_xs = [norm_xs; norm(xs(M,:))];
    norm_xh = [norm_xh; norm(xh(M,:))];
end

xe = xe./repmat(norm_xe,1,size(xe,2));
xs = xs./repmat(norm_xs,1,size(xs,2));
xh = xh./repmat(norm_xh,1,size(xh,2));

%% Use logistic regression to find dividing hyperplane for binary (0,1)
%% classification
xe = [ones(size(xe,1),1) xe];
xs = [ones(size(xs,1),1) xs];
xh = [ones(size(xh,1),1) xh];

data = {xe xs xh};
type = {'CIELAB' 'sCIELAB' 'hasler'};

for J=1:3,
    x = data{J};
    t = zeros(1,size(x,2)); % initial search point
    h_t = 1./(1 + exp(-(t*x')'));
    grad_l = (y - h_t)'*x;
    % Hessian
    temp = repmat((h_t.*(1-h_t)),1, size(x,2));
    H = -(x'*(temp.*x));


    N = num_iter; % number of iterations
    for I = 1:N
        h_t = 1./(1 + exp(-(t*x')'));
        grad_l = (y - h_t)'*x;
        temp = repmat((h_t.*(1-h_t)),1, size(x,2));
        H = -(x'*(temp.*x));
        t = t - (inv(H)*grad_l')';
    end

    t % solution. t(1) is the intercept term.
    eval(['LRmodels.' type{J} '_theta = t;']);
end

save LRmodels.mat LRmodels
filename = 'LRmodels.mat';
%% Real Code Ends Here
%% Make Pretty Plots
% % modify the assignment to temp_x to display different vectors--either x
% % (CIELAB error) or xs (sCIELAB error)
% 
for K = 1:3,

    figure()
    x=data{K};
    eval(['t = LRmodels.' type{K} '_theta;']);
    index = find(y==0);
    temp_x = x(index,:);

    % plot the data; change the index into temp_x to see comparisons among
    % different dimensions; temp_x(:,1) is constant so it's useless to plot
    % against it (just ignore it)
    plot(temp_x(:,5), temp_x(:,2), 'bo');
    hold on

    % index = find(y==1);
    % temp_x = x(index,:);
    % plot(temp_x(:,5), temp_x(:,2),'g^');
    % hold on

    index = find(y==1);
    temp_x = x(index,:);
    plot(temp_x(:,5), temp_x(:,2),'rx');
    hold on

    % index = find(y==3);
    % temp_x = x(index,:);
    % plot(temp_x(:,5), temp_x(:,2),'rx');

    grid on

    v = axis;
    x1 = v(3):.01:v(4);
    %[x1 x2]= meshgrid(v(3):v(4), v(5):v(6));
    %x1 = [v(5):.01:v(6)];
    %x2 = [v(3):.01:v(4)];
    % Since we're looking for h_t = 0.5, that is equivalent to (t*x')' = 0. So
    % we want t(1) + t(2)*x1 + t(3)*x2 + t(4)*x3 + t(5)*x4 = 0. Ignoring the x2 
    % and x3 dimensions, and using the x4 vector above, we get the x1 vector below:
    x4 = -(t(1)+t(2)*x1)./t(5);

    plot(x4,x1)
    hold off
    title(type{K});
    xlabel('Metric');
    ylabel('Distance (cm)');
end