function M=edgemap(edges)
% edgemap - Create arbitrary discretisation map from bin edges
%
% edgemap :: [[N+1]] -> dmap(N).
%
% See also: dmap

N=length(edges)-1;
M=dmap(N,@(t)map(edges,t),@(t)revmap(edges,t));

function I=map(edges,X)
	I=X; % make sure I is same shape as X
	min=edges(1);
	max=edges(end);
	for k=1:numel(X)
		x=X(k);
		if     x<min,  I(k)=-inf; 
		elseif x>=max, I(k)=inf;
		else I(k)=bsearch_it(edges,x); end
	end

function X=revmap(edges,I)
	X=[edges(I);edges(I+1)];
	

% bsearch :: real, [[N]] -> 1..N.

% recursive version.
function i=bsearch(y,x)
	l=length(y);
	if l==2, i=1;
	else
		half=ceil(l/2);
		if x<y(half), i=bsearch(x,y(1:half));
		else          i=bsearch(x,y(half:end))+half-1;
		end
	end
	
% iterative version
function i=bsearch_it(y,x)
	off=0;
	l=length(y);
	while l>2
		half=ceil(l/2);
		if x<y(half)
			y=y(1:half);
		else 
			y=y(half:end);
			off=off+half-1; 
		end
		l=length(y);
	end
	i=off+1;
