function kl = kl_gaussian( mu_1, Sigma_1, mu_2, Sigma_2 )

% function kl_gaussian( mu_1, Sigma_1, mu_2, Sigma_2 )
%
% KL divergence of two single Gaussian processes with full covariance matrix.
%
% "Kullback-Leibler Divergences of Normal, Gamma, Dirichlet and Wishart Densities."
% W.D. Penny
% Technical report, Wellcome Department of Cognitive Neurology, 2001
% 2001
%
% KL( N( mu_1, Sigma_1 ) || N( mu_2, Sigma_2 ) )
  
  if nargin < 4
    error( [ mfilename ' needs 4 input arguments.' ] );
  end
  
  % Check sizes
  
  a_size = size( mu_1 );
  
  if ( numel( a_size ) ~= 2 ) | ( a_size( 2 ) ~= 1 )
    error( [ mfilename ': input argument "mu_1" must be N x 1.' ] );
  end
  
  N = a_size( 1 );
  
  a_size = size( Sigma_1 );
  
  if ( numel( a_size ) ~= 2 ) | ( any( a_size( 1:2 ) ~= [ N N ] ) )
    error( [ mfilename ': input argument "Sigma_1" must be N x N, where "mu_1" is N x 1.' ] );
  end
  
  a_size = size( mu_2 );
  
  if ( numel( a_size ) ~= 2 ) | ( any( a_size( 1:2 ) ~= [ N 1 ] ) )
    error( [ mfilename ': input argument "mu_2" must be N x 1, where "mu_1" is also N x 1.' ] );
  end
  
  a_size = size( Sigma_2 );
  
  if ( numel( a_size ) ~= 2 ) | any( a_size( 1:2 ) ~= [ N N ] )
    error( [ mfilename ': input argument "Sigma_2" must be N x N, where "mu_1" is N x 1.' ] );
  end
  
  % Do it!
  
  delta_mu = mu_2 - mu_1;
  
  inv_Sigma_2 = inv( Sigma_2 );
  Sigma       = Sigma_1 * inv_Sigma_2;
  
  kl = 0.5 * ( delta_mu.' * inv_Sigma_2 * delta_mu + trace( Sigma ) - N - log( det( Sigma ) ) );