function c = cmi(x, y, z)
% Adapted from mi.m, which is under the following license:
% % % Copyright (c) 2012, Mo Chen
% % % Copyright (c) 2010, Michael Chen
% % % All rights reserved.
% % % 
% % % Redistribution and use in source and binary forms, with or without
% % % modification, are permitted provided that the following conditions are
% % % met:
% % % 
% % %     * Redistributions of source code must retain the above copyright
% % %       notice, this list of conditions and the following disclaimer.
% % %     * Redistributions in binary form must reproduce the above copyright
% % %       notice, this list of conditions and the following disclaimer in
% % %       the documentation and/or other materials provided with the distribution
% % % 
% % % THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% % % AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% % % IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% % % ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% % % LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% % % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% % % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% % % INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% % % CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% % % ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% % % POSSIBILITY OF SUCH DAMAGE.

    assert(numel(x) == numel(y));
    n = numel(x);
    x = reshape(x,1,n);
    y = reshape(y,1,n);
    z = reshape(z,1,n);

    l = min(z);
    z = z - l + 1;
    
    for zIter=min(z):max(z)   
        indices = find(z==zIter);
        xSlice = x(indices);
        ySlice = y(indices);
        
        
        l = min(min(xSlice),min(ySlice));
        xSlice = xSlice-l+1;
        ySlice = ySlice-l+1;
        k = max(max(xSlice),max(ySlice));
       
        
        idx = 1:length(indices);
        Mx = sparse(idx,xSlice,1,length(indices),k,length(indices));
        My = sparse(idx,ySlice,1,length(indices),k,length(indices));
        Pxy = nonzeros(Mx'*My/length(indices)); %joint distribution of x and y
        Hxy = -dot(Pxy,log2(Pxy+eps));
        
        Px = mean(Mx,1);
        Py = mean(My,1);
        
        % entropy of Py and Px
        Hx = -dot(Px,log2(Px+eps));
        Hy = -dot(Py,log2(Py+eps));
        
        % mutual information for this value of z
        q(zIter) = Hx + Hy - Hxy;
        indexCount(zIter) = length(indices);
        
    end
    
    c = dot(q,indexCount/n);
end