%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  File: logisticDemo.m
%  Matlab script file
%  Date: Oct, 08
%
% Demonstrate 2-class logistic regression on simple Gaussian data.

% Dependencies, Toolboxes:
%      iseToolbox/
% Author: ADJ
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Check Path and Constants  %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
%close all;

if ~exist('logisticReg','file')
  addpath('./util');
end

%  Check path.  If these aren't on your path, then you will need
%  to add these toolboxes to your path.  See ~jepson/pub/matlab/startup.m
which histo          % should be in iseToolbox\pyrTools\MEX\

%%%%%%%%  Initialize random number generator %%%%%%%%%%%%%%%%%%%%%%%
% Random number generator seed:
seed = round(sum(1000*clock));
rand('state', seed);
seed0 = seed;

% Set two means. 
%  mn(k, :) are x,y coords of Gaussian mean for each component.
mn = [0 0; 3 3];  

% Set four std devs in principal directions
sigma1 = [2, 2];
sigma2 = [1, 2];

% Set rotation angles of principal directions.
theta1 = 0;
theta2 = pi/3;

% Compute the square root of the covariance matrices for each component.
sqrtCov = zeros(2,2,2);
U = [cos(theta1) -sin(theta1); sin(theta1) cos(theta1)];
sqrtCov(:,:,1) = U*diag(sigma1);

U = [cos(theta2) -sin(theta2); sin(theta2) cos(theta2)];
sqrtCov(:,:,2) = U * diag(sigma2) * U';

% Generate Normally distributed data
nTarg = 1000; nNon = 1000;
testTarg =  sqrtCov(:,:,1) * randn(2, nTarg) + repmat( mn(1,:)', 1, nTarg); 
testNon =  sqrtCov(:,:,2) * randn(2, nNon) + repmat(mn(2,:)', 1, nNon); 

% plot the data set (red non, green targets)
figure(1); clf;
plot(testTarg(1,:), testTarg(2,:), '.g');
hold on;
plot(testNon(1,:), testNon(2,:), '.r');

pause;

%% %%%% Logistic regression %%%%%%%%%%%%%%%%%%%%%%%
% The data matrix has 1's in its first row.
X = [ones(1, nTarg+nNon); [testTarg testNon]];

% The label vector is zeros and ones, (ones for target).
y = [ones(1, nTarg), zeros(1,nNon)];

% Do logistic regression.  
[beta, conv] = logisticReg(X, y);

% Converged?
conv

% Plot result line
figure(1); 
ax = axis;
figCropBox = [ax([1 3]); ax([2 4])];
ep =  cropLineInBox(beta(2:end), beta(1), figCropBox);
figure(1); hold on;
plot(ep(:,1), ep(:,2), 'k');

pause;

% Plot the conditional probability p(Target| X(:,k), beta)
p = logistic(X, beta);
figure(1); clf; hold on;
for k = 1:size(X,2)
  plot(X(2,k), X(3,k), 'ok', 'MarkerSize', 10, 'MarkerFaceColor', p(k)*[1,1,1]);
end
pause;
plot(testTarg(1,:), testTarg(2,:), '.g');
plot(testNon(1,:), testNon(2,:), '.r');
pause;

% True positive and true negative rates.
[sum(p >= 0.5 & y > 0)/nTarg, ...
  sum(p < 0.5 & y == 0)/nNon]

if any(y==1) & any(y == 0)
  [hTarg, hxTarg] = histo(beta'*X(:,y==1), 101);
  hTarg = hTarg/(sum(hTarg) * (hxTarg(2) - hxTarg(1)));
  [hNon, hxNon] = histo(beta'*X(:,y==0), 101);
  hNon = hNon/(sum(hNon) * (hxNon(2) - hxNon(1)));

  figure(2); clf;
  plot(hxTarg, hTarg, 'g');
  hold on;
  plot(hxNon, hNon, 'r');
  ax = axis;
  plot([0 0], ax(3:4), 'k');
  title('Histogram of beta'' X, target/non (g/r)');
  xlabel('beta'' x');
  ylabel('Probability per bin');
end
