%
% Example on epage XXX
% 
% Problem on epage 517
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');

randn('seed',0); 
rand('seed',0); 

%--
%Estimate for OKWHITE
%--

load okwhite;
X = [okwhx,okwhy];
n = size(X,1);  

% estimate the first order intensity: 
% 

h = 100; 

% Get the convex hull to use in csintkern.
% -- K is the indices to points ON the convex hull
K    = convhull(okwhx, okwhy);
cvh  = [okwhx(K), okwhy(K)]; 
[xl,yl,lamhat] = csintkern(X,cvh,h);
figure; imagesc(xl(:,1),yl(1,:),lamhat); colorbar; 
axis xy; 
map = gray(256);
map = flipud(map); % Flip the colormap so zero is white and max is black.
colormap(map); 
title( 'okwhite first order spatial intensity' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter12/prob_12_5_okwhite_fo', 'epsc' ); 

% estimate the second order properties "Ghat": 
% -- get this function first and plot.
%
% The G function is the distribution function of nearest neighbor distances for each event.
% Find the distances for all points.
dist = pdist(X);
% Convert to a matrix and put large numbers on the diagonal.
D = diag(realmax*ones(1,n)) + squareform(dist);
% Find the smallest distances in each row or col.
mind = min(D);

% create a sampling of nearest neighbor distances: 
dw = (max(mind)-min(mind))/30; 
w  = min(mind):dw:max(mind);
nw = length(w);
ghat = zeros(1,nw);

% Now get the values for ghat.
for i = 1:nw
  ind = find(mind<=w(i));
  ghat(i) = length(ind);
end
ghat = ghat/n;
% Plot the Ghat as a function of w. a fast rise shows evidence of clustering.
gfh=figure; okwg = plot(w,ghat,'g-o'); grid on; 
xlabel('Event-Event Distances - w'); ylabel('Ghat'); 
title( 'Estimated Event-Event Distribution Function' ); 
%%%saveas( gcf, '../../WriteUp/Graphics/Chapter12/prob_12_5_okwhite_so', 'epsc' );

%--
%Estimate for OKBLACK
%--

load okblack;
X = [okblx,okbly];
n = size(X,1);  

% estimate the first order intensity: 
% 

% Get the convex hull to use in csintkern.
% -- K is the indices to points ON the convex hull
K    = convhull(okblx, okbly);
cvh  = [okblx(K), okbly(K)]; 
[xl,yl,lamhat] = csintkern(X,cvh,h);
figure; imagesc(xl(:,1),yl(1,:),lamhat); colorbar; 
axis xy; 
map = gray(256);
map = flipud(map); % Flip the colormap so zero is white and max is black.
colormap(map); 
title( 'okblack first order spatial intensity' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter12/prob_12_5_okblack_fo', 'epsc' ); 

% estimate the second order properties "Ghat": 
% -- get this function first and plot.
%
% The G function is the distribution function of nearest neighbor distances for each event.
% Find the distances for all points.
dist = pdist(X);
% Convert to a matrix and put large numbers on the diagonal.
D = diag(realmax*ones(1,n)) + squareform(dist);
% Find the smallest distances in each row or col.
mind = min(D);

% create a sampling of nearest neighbor distances: 
dw = (max(mind)-min(mind))/30; 
w  = min(mind):dw:max(mind);
nw = length(w);
ghat = zeros(1,nw);

% Now get the values for ghat.
for i = 1:nw
  ind = find(mind<=w(i));
  ghat(i) = length(ind);
end
ghat = ghat/n;
% Plot the Ghat as a function of w. a fast rise shows evidence of clustering.
figure(gfh); hold on;  
okbg=plot(w,ghat,'r-x'); grid on; 
xlabel('Event-Event Distances - w'); ylabel('Ghat'); 
title( 'Estimated Event-Event Distribution Function' ); 
legend( [okwg,okbg], {'ok white thefts', 'ok black thefts'}, 'location', 'southeast' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter12/prob_12_5_so', 'epsc' );

clear functions;