%
% An extension of example 9.3
% This illustrates the bivariate histogram. Use the topic 8 and 11 code from Example 7.8.
% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; 
rehash; 
clc; 
clear; addpath( '../../Code/eda_data' ); addpath( '../../Code/eda_toolbox' ); addpath( '../Chapter1' ); 


% Example 9.3
% 
load L1bpm
% Reduce the dimensionality using Isomap.
options.dims = 1:10;    % These are for ISOMAP.
options.display = 0;
[Yiso, Riso, Eiso] = isomap(L1bpm, 'k', 7, options);
% Get the data out.
XX = Yiso.coords{2}';
inds = find(classlab==8 | classlab==11);
x = [XX(inds,:)];
[n,p] = size(x);

% Need bin origins.
bin0 = floor(min(x)); 
% The bin width h, for p = 2:
h = 3.5*std(x)*n^(-0.25); 
% Find the number of bins
nb1 = ceil((max(x(:,1))-bin0(1))/h(1));
nb2 = ceil((max(x(:,2))-bin0(2))/h(2));
% Find the bin edges.
t1 = bin0(1):h(1):(nb1*h(1)+bin0(1));
t2 = bin0(2):h(2):(nb2*h(2)+bin0(2));
[X,Y] = meshgrid(t1,t2);
% Find bin frequencies.
[nr,nc] = size(X);
vu = zeros(nr-1,nc-1);
for i = 1:(nr-1)
   for j = 1:(nc-1)
      xv = [X(i,j) X(i,j+1) X(i+1,j+1) X(i+1,j)];
      yv = [Y(i,j) Y(i,j+1) Y(i+1,j+1) Y(i+1,j)];
      in = inpolygon(x(:,1),x(:,2),xv,yv);
      vu(i,j) = sum(in(:));
   end
end
% Find the proper height of the bins.
Z = vu/(n*h(1)*h(2));
% Plot using bars.
bar3(Z,1,'w')


% Try to get the bin centers.
ct1 = t1 + h(1)/2;  % Columns
ct2 = t2 + h(2)/2;  % Rows
ct1(end) = [];
ct2(end) = [];
for i = 1:length(ct1)
    ct1s{i} = num2str(ct1(i),'%2.1f');
end
for i = 1:length(ct2)
    ct2s{i} = num2str(ct2(i),'%2.1f');
end
set(gca,'XTickLabel',ct1s,'YTickLabel',fliplr(ct2s))

% Plot as a surface plot.
% Get some axes that make sense.
[XX,YY]= meshgrid(linspace(min(x(:,1)),max(x(:,1)),nb1), linspace(min(x(:,2)),max(x(:,2)),nb2));
surf(XX,YY,Z)