%% Geo-502 Lab 
% Nearest-neighbor and range search
% Phaedon Kyriakidis

%% Matlab functions for neighbor(s) search
% knnsearch   -> search for K nearest neighbors among points of 2 sets
% rangesearch -> search for points (neighbors) within some distance

%% Create some scattered point locations (random coordinates)
rng(1234) % 1234 is the seed
x = rand(100,1)*100; y = rand(100,1)*100;
nEv = numel(x); % number of elements in x

%% Display point locations
plot(x,y,'*')
id = num2str((1:nEv)'); % convert numbers to character strings
text(x,y,id); % add text with id's
title('Event locations')
axis([0 100 0 100]); axis square

%% 1/1: Indices of nearest neighbors and distances between them
[idNN,distNN] = knnsearch([x y],[x y],'K',2);
size(idNN)
disp([idNN(1:3,:) distNN(1:3,:)])
% Array idNN contains the indices (row numbers in array [x y]) 
% of the 2 nearest neighbors for each event, where as array distNN
% contains the corresponding distances.
% 
% If we set K=1, function knnsearch returns the same points as 
% nearest neighbors of themselves and 0 distances; this is not desired. 
% So we need to set K=2, and delete the 1st entries of idNN and distNN 
% to keep the 2nd nearest neighbor only, which is actual the correct 
% nearest neighbor. 
idNN(:,1) = [];
distNN(:,1) = [];

%% 1/2: Display nearest neighbor of some event
ev2plot = 3;
plot(x,y,'*'); % plot event locations
axis equal square
text(x,y,id); % add text with ids
hold on
plot(x(ev2plot),y(ev2plot),'ro'); % event coordinates = row to plot
plot(x(idNN(ev2plot)),y(idNN(ev2plot)),'go'); % nearest neighbor
hold off
line([x(ev2plot) x(idNN(ev2plot))],[y(ev2plot) y(idNN(ev2plot))])
% The length of the line segment is the distance between event #3
% and its nearest event #48
title('Event #3 and its nearest neighbor')

%% 1/3: Histogram of nearest neighbor distances
hist(distNN,20);
figure
cdfplot(distNN);
% From the plot we deduce that the proportion of events within 10 dist.
% units away from their nearest events is about 0.93. We can check this
% as:
sum(distNN <= 10) / nEv

%% 2/1: Range search
searchR = 15;
[idR,distR] = rangesearch([x y],[x y],searchR);
whos idR
idR(1:5)
idR{1}
idR{2}
distR{1}
% Both arrays idR and distR are cell arrays, with each cell containing
% information on the neighbors of each event within a circle of distance
% searchR. Since the number of neighbors within 15 is different from
% one event to another, we need a cell array structure to hold that 
% information.
%
% Same problem as with function knnsearch: the 1st entries of idR 
% are the same indices as those of the events themselves, and the 1st
% entries of distR are the distances between the neighbors and the event
% queried.

%% Delete 1st entries
idR = cellfun(@(ii) ii(2:end),idR,'UniformOutput',false);
idR{1}
distR = cellfun(@(d) d(2:end),distR,'UniformOutput',false);
distR{1}

%% 2/2: Display neighbors within distance from some event
ev2plot = 3;
plot(x,y,'*'); % plot events in point pattern
axis equal square
text(x,y,id); % add text with ids
hold on
plot(x(ev2plot),y(ev2plot),'ro'); % event under consideration
plot(x(idR{ev2plot}),y(idR{ev2plot}),'go'); % neighbors 
drawellipse([x(ev2plot) y(ev2plot)],[searchR searchR],0,500,'g'); % circle
hold off
title('Neighbors of event #3 within circle of radius 15')

%% 2/3: Number of neighbors within distance 
nNeighR = cellfun(@(ii) numel(ii),idR,'UniformOutput',true);
nNeighR(1:5)
scatter(x,y,25,nNeighR,'filled'); colorbar
axis equal square
title('Number of neighbors within distance 15 from each event')
% This resembles a quadrat count procedure, where quadrats are circles
% of radius 15, and colors represent counts within each such circle

%% 2/3: Total # of neighbors within distance
nNeighRtot = sum(nNeighR) 
% sum of number of neighbors within distance 15 from each event = 
% total # of event pairs lying up to 15 distance units away from each other

%% Same result through distance matrix
Dupp = pdist([x y]); 
% This provides the distances in the upper diagonal of the NxN distance
% matrix, since the distance matrix is symmetric.
sum(Dupp<=searchR)
% Here we find which distances are <= 15, and then count them..
% This is the same as computing the number of pairs up to distance 15
% from each other, which is what we already did above.
%
% Alternative procedure, based on the full distance matrix (not just
% its upper diagonal part):
%indDiag = 1:nEv+1:nEv*nEv;
%D(indDiag) = 1000;
%sum(D(:)<=searchR)
%numel(find(D<=searchR))

%% 3/1: Grid specs
GridSpecs = [100 0.5 1; 100 0.5 1]; % Grid specifications
areaCell = prod(GridSpecs(:,3))
xyG = makegrid(GridSpecs); % coordinates of grid nodes
size(xyG)
xGrid = xyG(:,1); yGrid = xyG(:,2);
plot(xGrid,yGrid,'k.'); % plot grid nots
axis equal square
hold on
plot(x,y,'b*') % overlay event locations
hold off
title('Grid nodes and event locations')
% These are the grid nodes at which we will be estimating event intensity

%% 3/2: Neighbors within distance from grid nodes
[idR,distR] = rangesearch([x y],[xGrid yGrid],searchR);
whos idR
size(idR)

%% Display neighbors within distance from some grid node
node2plot = 5050;
idR{node2plot}
plot(x,y,'*'); % event locations
axis equal square
text(x,y,id); % add text with ids
hold on
plot(xGrid(node2plot),yGrid(node2plot),'ro'); % node under consideration
plot(x(idR{node2plot}),y(idR{node2plot}),'go'); % neighbors
% circle of radius 15 centered at the node of interest
drawellipse([xGrid(node2plot) yGrid(node2plot)],[searchR searchR],0,500,'g');
hold off
title('Events within circle of radius 15 from node 5050')

%% 3/3: Number of neighbors within distance from grid nodes
nNeighR = cellfun(@(x) numel(x),idR,'UniformOutput',true);
size(nNeighR)
rastermap(GridSpecs,nNeighR,1,'cont',[]); colorbar
hold on
plot(x,y,'k*'); % event locations
hold off
sum(nNeighR.*areaCell) % Volume under frequency (count) surface
title('Count surface')
% This surface represents the number of events within a circle of 
% radius 15 centered at any grid node. One can see the patterns of 
% the circles used.

%% 3/4: Intensity of point pattern
intensR = nNeighR./(pi*searchR^2);
% The weight given to each event within a circle of radius 15 from any
% grid node = 1/area of circle (equal weighting)
rastermap(GridSpecs,intensR,1,'cont',[]); colorbar
hold on
plot(x,y,'k*'); % event locations
hold off
sum(intensR.*areaCell) % Volume under intensity surface
title('Intensity surface')
% The volume under the intensity surface should be N. Here it is not,
% because some circles (the ones near the borders of the study region)
% extend beyond the region. This discrepancy between the theoretical
% value N and the computed one, is due to edge effects.

%% Display equal weights for events within distance from some grid node
node2plot = 5050;
idR{node2plot}
wR = pi*searchR*searchR;
plot(x,y,'ko'); % event locations
axis equal square
text(x,y,id);
hold on
plot(xGrid(node2plot),yGrid(node2plot),'ro'); % node under consideration
scatter(x(idR{node2plot}),y(idR{node2plot}),25,repmat(wR,nNeighR(node2plot),1),'filled');
colorbar
drawellipse([xGrid(node2plot) yGrid(node2plot)],[searchR searchR],0,500,'g');
hold off
title('Equal weights for events around grid node 5050')

%% Density surface
densR = nNeighR./(nEv*pi*searchR^2);
rastermap(GridSpecs,densR,1,'cont',[]); colorbar
hold on
plot(x,y,'k*'); % event locations
hold off
sum(densR.*areaCell) % Volume under density surface
title('Density surface')
% Volume under the density surface should be 1, but it is not due to 
% edge effects

%% 4/5: Convert distances between neighbors and nodes to weights Quartic
s = searchR;
weightsR = cellfun(@(d) (3/(pi*s^2))*(1-(d/s).^2).^2,distR,'UniformOutput',false);
whos weightsR
weightsR{node2plot}
distR{node2plot}

%% Display weights for events within distance from some grid node
node2plot = 5050;
idR{node2plot}
plot(x,y,'ko'); % event locations
axis equal square
text(x,y,id);
hold on
plot(xGrid(node2plot),yGrid(node2plot),'ro'); % node under consideration
% weights
scatter(x(idR{node2plot}),y(idR{node2plot}),25,weightsR{node2plot},'filled');
colorbar
drawellipse([xGrid(node2plot) yGrid(node2plot)],[searchR searchR],0,500,'g'); % circle
hold off
title('Distance-decay weights around grid node 5050')

%% 4/6: Sum of weights = intensity
intensRkern = cellfun(@(w) sum(w),weightsR,'UniformOutput',true);
rastermap(GridSpecs,intensRkern,1,'cont',[]); colorbar
hold on
plot(x,y,'k*'); % event locations
hold off
sum(intensRkern.*areaCell) % Volume under density surface
title('Kernel-based intensity surface')
% Volume under intensity surface should be N, but it is not due to edge
% effects