[y2 idx]=nearest(x,y,type,varargin)

function [y2 idx]=nearest(x,y,type,varargin)
%NEAREST - y2(i)=y(idx(i)) is the nearest neighbor of x(i) in vector y
% Designed for finding pairs quickly in (not equidistant) time series 
% Type parameter:
%  'l': only lower neighbors accepted
%  'u': only upper neighbors accepted
%  'n' or not given: nearest (upper or lower) neighbors selected
%If there is no neighbor for x(i), y2(i)=idx(i)=NaN is returned.
%
%Options:
% 'self': y2(i)==x(i) [zero distance] is accepted as 'neighbor'
% 'unique': if i<>j and y2(i)==y2(j), y2(i)=y2(j)=NaN is returned
%
%Example:
% x=       [ 2   5 7 9 10  12]; 
% y=       [1 3 4 6 8    11];
% [y2 idx2]=nearest(x,y,'l','unique')
% [y3 idx3]=nearest(x,y,'u')
%
% Results:
%  y2=[1 4 6 NaN NaN 11] idx2=[1 3 4 NaN NaN 6]
%  y3=[3 6 8 11  11 NaN] idx3=[2 4 5 6   6 NaN]

%% Transposition
% Internal y is column - y2 has the same dim as input y
if size(y,2)>1 
    y=y';
    transp=true;
else
    transp=false;
end
% x shall be row
if size(x,1)>1 
    x=x';
end

%% Neighbor type selection
if (2==nargin)
    type='n'; %nearest
end

%% Calculation
Y=repmat(y,1,length(x));
X=repmat(x,length(y),1);
D=X-Y; %(i,j) element: x(j) - y(i)>0, if y is a lower neighbor of x
%% Clear not accepted neighbors
switch(type) 
    case 'l'
       D(D<0)=NaN; 
    case 'u'
       D(D>0)=NaN;         
end

if ~ismember('self',varargin)
    D(0==D)=NaN;
end
%% Find the nearest neighbors
D=abs(D);
[dist idx]=min(D,[],1);
y2=y(idx);

if transp
    y2=y2';
end
idx(isnan(dist))=NaN;

%% Check 'unique'
if ismember('unique',varargin)
    tmp=find(diff(idx)==0);
    idx(tmp)=NaN; idx(tmp+1)=NaN;
end

y2(isnan(idx))=NaN;