function [MI,rMI,bMI]=AdaptiveMI2(count);
%
% Implement a greedy adaptive procedure for estimating mutual
% informtion from a joint distribution count matrix. 
%
% This procedure was empirically found to be unbiased for the
% distributions observed in our data. (Nelken et al 2005).
%
% If using, please cite Nelken et al, J.Comp.Neurosci., 2005. 
% Israel Nelken (C) 2003
% 

h=count;
if(sum(h(:))==0)
  MI=0;
  return
end
nh = sum(h(:)); % June 2005 
if isempty(h)
    MI=0;
else
    ph=h/nh;
    phx=sum(ph,1);
    phy=sum(ph,2);
    approx=phy*phx;
    sph=ph;
    sph(sph==0)=1;
    approx(approx==0)=1;
    rMI = zeros(1,size(h,1)+size(h,2));
    bMI = zeros(1,size(h,1)+size(h,2));    
    last=1;
    rMI(last)=sum(sum(ph.*log2(sph./approx)));
    bMI(last)=bias(h);
    h2u=h;
    while min(size(h2u))>1
        th=h2u;
        [mx,Ix]=min(phx);
        [my,Iy]=min(phy);
        if my<mx
            th=th';
            Ix=Iy;
            phx=phy;
        end
        if Ix==1
            th(:,2)=th(:,2)+th(:,1);
            th=th(:,2:end);
        elseif Ix==length(phx)
            th(:,end-1)=th(:,end-1)+th(:,end);
            th=th(:,1:end-1);
        elseif phx(Ix-1)<phx(Ix+1)
            th(:,Ix-1)=th(:,Ix-1)+th(:,Ix);
            th=[th(:,1:Ix-1) th(:,Ix+1:end)];
        else
            th(:,Ix)=th(:,Ix)+th(:,Ix+1);
            th=[th(:,1:Ix) th(:,Ix+2:end)];
        end
        h2u=th;
        ph=h2u/nh;
        phx=sum(ph,1);
        phy=sum(ph,2);
        approx=phy*phx;
	inds = find(ph>0);
	last = last+1;
        rMI(last)=sum(sum(ph(inds).*log2(ph(inds)./approx(inds))));
        bMI(last)=bias(h2u);
    end
    rMI = rMI(1:last);
    bMI = bMI(1:last);    
    tMI=rMI-bMI;
    MI=max(tMI);
end

% -----------------
function b=bias(h)
  ss=size(h);
  b=(ss(1)-1)*(ss(2)-1)/(2*sum(h(:))*log(2));
return