function mi = mikl( x, y, type, k )
%MIKL: Kozachenko-Leonenko estimates of the continuous mutual information
% [1] or its offset version [2]
%
% USAGE: mi = mikl( rand( 100, 3 ), rand( 100, 3 ) )
%
% INPUT: x and y: matrices with sample in each rows dimensions as columns
%        type   : type of nearest: raw Kozachenko-Leonenko estimator, 'kl',
%                 or offset Kozachenko-Leonenko
%                 estimator,'klo'(default)
%        k      : nearest-neightbour criterion: k = 1 (default) for closest
%                 neightbour, k = 2 for second closet neighbour, ... See
%                 ref [3]
%
% OUTPUT: mi: estimated mutual information
%
% The efficient mex C++ file was created by Martin Sanz (MS: anade
% direccion email)
%
% ADD HERE THE COMMONS CREATIVE LICENCES (GPL3?)
%
% REFERENCES:
% [1] L. F. Kozachenko and N. N. Leonenko. Sample estimate of the entropy
%     of a random vector. Problems of Information Transmission,
%     23(2):95?101, 1987. [Translated from Problemy Peredachi Informatsii,
%     23(2): 9?16, 1987.]
%
% [2] I. Marin-Franch and D. H. Foster. Estimating information from image
%     colors: An application to digital cameras and natural scenes. IEEE
%     Transactions on Pattern Analysis and Machine Intelligence, 35(1):78?
%     91, 2013.
%
% [3] M. N. Goria, N. N. Leonenko, V. V. Mergel, and P. L. Novi Inverardi.
%     A new class of random vec- tor entropy estimators and its
%     applications in testing statistical hypotheses. Journal of
%     Nonparametic Statistics, 17(3):277?297, 2005.
%
% Ivan Marin-Franch 01/06/2014 (EU date format)

%%%% NOTICE THAT WE CAN CALCULATE, TECHNICALLY, THE INFORMATION BETWEEN A
%%%% 3-DIMENSIONAL RANDOM VARIABLE AND A 6-DIMENSIONAL ONE

% check input
if nargin < 3
    error( 'input must be x, y, and the type of estimator ''kl'' for the original estimator or ''klo'' for the offset version' );
end

if nargin < 4
    k = 1;
end

assert( k > 0, 'k-nearest neighbour must be larger than 0' )

% check input: correct type of estimator?
if ~strcmp( type, 'kl' ) && ~strcmp( type, 'klo' )
    error( 'incorrect type of estimator' )
end

% check input: equal sizes?
if size( x, 1 ) ~= size( y, 1 )
    error( 'variables x and y should have the same size (number of rows)' )
end

% get the joint distribution
xy = [ x y ];

% if type is 'klo' (classic estimator) then calculate differential
% entropies with x, y and xy directly. If type is 'kloffset', then
% calculate mutual information as if they were Gaussian and prepare the
% data to calculate the offset differential entropies
miG = 0;
if strcmp( type, 'klo' )
    cx  = cov( x );
    cy  = cov( y );
    cxy = cov( xy );
    x   = x  / sqrtm( cx  );
    y   = y  / sqrtm( cy  );
    xy  = xy / sqrtm( cxy );
    miG = 1 / 2 * log2( det( cx ) * det( cy ) / det( cxy ) );
end

% Compute continuous mutual information between x and y, or between
% transformed x and transformed y
mi = miG                             ...
     + log2( exp( entkl_internal( x, k  ) ) ) ...
     + log2( exp( entkl_internal( y, k  ) ) ) ...
     - log2( exp( entkl_internal( xy, k ) ) );