function I=quad_tsallis_fast(X,Y,xcov,ycov,NS)
% function I=quad_tsallis_fast(X,Y,xcov,ycov)
% 
% Quadratic Tsallis mutual information between X and Y.
% Defines mutual information as the sum of the individual
% entropies, minus the joint entropy.
%
% This can be written as I(X;Y) = H(X)+H(Y)-H(X)H(Y) - H(X,Y),
% which simplifies to \int p(x,y)^2 - \int p(x)^2 \int p(y)^2.
% 
% Uses a Nystrom approximation to speed up the computation.
% Specify the number of points to use in the approximation
% with the NS parameter.
%
% Assumes a Parzen density estimate for both variables,
% with spherical, homoscedastic Gaussian kernels.
%
% Currently, this requires that you have the same
% number of samples for X and Y, because it assumes that aligned
% samples are samples from the joint of p(x,y)
%
% D. Wingate 10/25/2006
%

[dimx, nx] = size(X);
[dimy, ny] = size(Y);

if ( nx ~= ny )
  error( 'Number of samples must be equal!' );
  return;
end;

xcov = 2*xcov;
ycov = 2*ycov;

norm_x = 1/sqrt((2*pi*xcov)^dimx);
norm_y = 1/sqrt((2*pi*ycov)^dimy);

nsq = 1/(nx*nx);

% the joint distribution of X and Y.
% we premultiply by the gaussian width here because
% it may be different for X and Y, but our rbf4nn
% function can't cope with that...
XY = [ 1/sqrt(xcov) * X; 1/sqrt(ycov) * Y ];

% get up to NS linearly independent points
%iX = get_lin_ind_pts( X, xcov, 2*xcov, NS );
%iY = get_lin_ind_pts( Y, ycov, 2*ycov, NS );

[iX,bX,bXinv] = get_lin_ind_pts_K( X, xcov, 0.001, NS );
[iY,bY,bYinv] = get_lin_ind_pts_K( Y, ycov, 0.001, NS );
iXY = unique( [iX iY] ); % the union of X and Y

% the subsampled points
sX = X( :, iX );
sY = Y( :, iY );
sXY = XY( :, iXY );

% entropy of X
vX = sum( rbf4nn( sX, X, xcov ), 2 );
SX = norm_x * nsq * ( vX' * bXinv * vX );

% entropy of Y
vY = sum( rbf4nn( sY, Y, ycov ), 2 );
SY = norm_y * nsq * ( vY' * bYinv * vY );

% joint entropy
bXY = rbf4nn( sXY, sXY, 1 );
vXY = sum( rbf4nn( sXY, XY, 1 ), 2 );
SXY = (norm_x*norm_y) * nsq * (vXY' * pinv( bXY ) * vXY );

% Equivalent!
%I = HX + HY + HX*HY - HXY;

I = SXY - SX*SY;

%fprintf( 'HX=%.4f HY=%.4f HX*HY=%.4f HXY=%.4f ->\t %.4f\n', bobX, bobY, ...
%         bobX*bobY, bobXY, I );

return;

%
% ---------------------------------------
%

function [v,Kp,Kpinv] = get_lin_ind_pts_K( pts, sigsq, ald_thresh, maxsize )

  cnt = size(pts,2);
  v = randind( maxsize, cnt );

  K = rbf4nn( pts(:,v), pts(:,v), sigsq );
  dp = kdict_init( K, ald_thresh );

  for i=2:maxsize
    dp = kdict( dp, i );
  end;

  v = v( dp.Dict )';
  Kp = dp.K;
  Kpinv = dp.Kinv;
return;

%
% ---------------------------------------
%

function [v] = get_lin_ind_pts( pts, sigsq, ald_thresh, maxsize )

  dp = dict_init( @rbf4nn, sigsq, ald_thresh, pts( :, 1) );
  v(1) = 1;
  
  cnt = size( pts, 2 );

  for i = 2:cnt
    dp = dict( dp, pts( :, i) );
    if ( dp.addedFlag )
      v( length(v)+1 ) = i;
    end;
    if ( size(dp.Dict,2) >= maxsize )
      break;
    end;
  end;

return;