function [Hxy, Hyx] = Hi(X,Y,m,tau,nperm)

% Hi Calculates a time resolved measure of directional couplings as
% introduced in 
% Andrzejak RG, Ledberg A, Deco G (2006): Detecting event-related time-dependent directional couplings. New Journal of Physics 8, 6
%
% 
% 
% usage: [...]=Hi(X,Y,m,tau);
% Inputs:
%          X is a #trials x # time-points data matrix 
%          Y is a #trials x # time-points data matrix 
%          m is the number of points in the embedding window
%          TAU is the embedding delay
%          nperm is the number of permutations to use
% Outputs: 
%         [Hxy,Hyx] = Hi(X,Y,m,tau);
%         returns the measures from x to y and y to x as two vectors
%         
%
% 
%
% 14 October 2005
% Ralph G. Andrzejak, Anders Ledberg, 
% Univesritat Pompeu Fabra

% Check that enough inputs are given

if (nargin < 2)
  'error: too few arguments' 
  'usage: '
  help Hi
  return
end

% set default values for unspecified arguments
if (nargin <= 2)
  m=2;
  tau=1;
  nperm=100;
end
if (nargin <= 3)
  tau=1;
  nperm=100;
end
if (nargin <= 4)
  nperm=100;
end


% check that data matrics are of same size

% note that you should use the != operator instead of ~= below if you want to run
% this code under GNU Octave

if(size(X)~=size(Y))
  'error: input data are not of same size'
  help permute_Hi
  return  
end


ntrials=size(X,1);
ntimepts=size(X,2);

% the number of time points in the output is equal 
% to ntimepts-(m-1)*tau
Hxy=zeros(ntimepts-(m-1)*tau,nperm + 1);
Hyx=zeros(ntimepts-(m-1)*tau,nperm + 1);

% initialize some temporary variables

index=1:tau:(m-1)*tau+1;
distX=zeros(ntrials,ntrials);
distY=zeros(ntrials,ntrials);


% initalize permutations
permindx=zeros(ntrials,100);
for i=1:100
  permindx(:,i)=randperm(ntrials)';
end
% here's the main loop over time 

for t=1:ntimepts-(m-1)*tau
  % first we extract the data for the current window
  currX=X(:,index);
  currY=Y(:,index);
  % then we generate all the distance for X and Y
  for i=1:ntrials
    currTrial=currX(i,:);
    distX(i,:)=sum( ((currX-repmat(currTrial,ntrials,1)).^2)' );
    currTrial=currY(i,:);
    distY(i,:)=sum( ((currY-repmat(currTrial,ntrials,1)).^2)' );
  
  end
  % store raw distances for later
  rawdistX=distX;
  rawdistY=distY;

  % first we calculate the avearge distance for all trials
  Rx=mean(distX);
  Ry=mean(distY);

  % to find the Rxy and Ryx we first add a large value to the diagonal (which
  % are 0 otherwise). 
  distX=eye(ntrials)*1e10+distX;
  distY=eye(ntrials)*1e10+distY;
  
  [minX,indxX]=min(distX);
  [minY,indxY]=min(distY);
  Rxy=diag(distX(:,indxY));
  Ryx=diag(distY(:,indxX));
  Hxy(t,1)=mean(log(Rx./Rxy'));
  Hyx(t,1)=mean(log(Ry./Ryx'));

  % here's the loop over permutations 
  % notice that we don't need to recalculate the distances
  for p=1:nperm
    % permute Y
    distY=rawdistY(permindx(:,p),permindx(:,p));
    distX=rawdistX;

    Rx=mean(distX);
    Ry=mean(distY);

    distX=eye(ntrials)*1e10+distX;
    distY=eye(ntrials)*1e10+distY;
  
    [minX,indxX]=min(distX);
    [minY,indxY]=min(distY);
    Rxy=diag(distX(:,indxY));
    Ryx=diag(distY(:,indxX));
    Hxy(t,1+p)=mean(log(Rx./Rxy'));
    Hyx(t,1+p)=mean(log(Ry./Ryx'));
  end
  % finally move the index one step
  index = index +1;
end


% note use 'endfunction' if you are running this code under GNU Octave


