% Horacek: [out,outt,outb]=closest(a,v,monotone)
% this function finds index of the closest value of vector 'v' 
% to number 'a' so that  |a-v(i)| is minimal
% if 'a' is also vector, it finds length(a) corresponding values
% if v is a matrix then a has to be vector of length(a)=size(v,1)
% if nargin==3 => don't invert 'a' with 'v'
% if monotone==1 then assuming monoteneity in v => much faster; not done for [out,outt,outb]=...
function [out,outt,outb]=closest(a,v,monotone)
if isempty(a) | isempty(v) out=[]; outt=[]; outb=[]; return; end;
if nargin<3 monotone=0; end;

sz=size(a);
if all(sz>1) 
 for i=1:sz(1) 
  if all(sz==size(v)) vv=v(i,:); else vv=v; end
  if nargout>1 [out(i,:),outt(i,:),outb(i,:)]=closest(a(i,:),vv,monotone);
  else out(i,:)=closest(a(i,:),vv,monotone);
  end
 end
 return
end

%if nargin<3 & length(a(:))>length(v(:)) tmp=v; v=a; a=tmp; end;

%out=[];
s=size(v);
%if     s(1)==1 yy=repmat(v,length(a),1);
%elseif s(2)==1 yy=repmat(v,1,length(a))';
%else yy=v;
%end;
if monotone
 ms=max(s);
 nn=round(sqrt(ms));
 out=closest(a,v(1:nn:end),0); % first approximation
 out=nn*out; %outb=nn*outb; outt=nn*outt; 
 ind8=(-1-2*nn):(2*nn+1);
 for iu=1:length(a)
  if ~isnan(out(iu))
   ind3=out(iu)+ind8; ind3=ind3(ind3>0 & ind3<=ms);
   out(iu)=ind3(closest(a(iu),v(ind3),0)); %outt=ind8(outt1)+outt; outb=ind8(outb1)+outb; 
  end
 end;
 % outt a outb netestovano zda-li v poradku
 return
end;

if 0 & ~all(s~=1) % zkousim to delat maticove, ale je to pomalejsi! protoze musi alokovat hodne pameti
 [tmp,out]=min(abs(repmat(a(:),1,max(s))'-repmat(v(:),1,length(a))));
else
 for j=1:length(a)
  if all(s~=1) y=a(j)-v(j,:); else y=a(j)-v; end;
% if any(s==1)  y=a(j)-v; else y=a(j)-v(j,:); end;
  if nargout>1 temp=find(y<=0); [tmp,ind]=max(y(temp)); outt(j)=min([NaN,temp(ind)]); end;
  if nargout>2 temp=find(y>=0); [tmp,ind]=min(y(temp)); outb(j)=min([NaN,temp(ind)]); end;
  [tmp,ind]=min(abs(y)); out(j)=ind(1);
 end;
end;
indnan=find(isnan(a)); out(indnan)=NaN;
if nargout>1 tmp=find(isnan(outt)); outt(tmp)=out(tmp); outt(indnan)=NaN; end;
if nargout>2 tmp=find(isnan(outb)); outb(tmp)=out(tmp); outb(indnan)=NaN; end;

