function [I,SX,x,W,F] = ngausshermi(f,n,m,P,param)
% NGAUSSHERMI - ND scaled Gauss-Hermite quadrature (cubature) rule
%
% Syntax:
%   [I,x,W,F] = ngausshermi(f,p,m,P,param)
%
% In:
%   f - Function f(x,param) as inline, name or reference
%   n - Polynomial order
%   m - Mean of the d-dimensional Gaussian distribution
%   P - Covariance of the Gaussian distribution
%   param - Optional parameters for the function
%
% Out:
%   I - The integral value
%   x - Evaluation points
%   W - Weights
%   F - Function values
%
% Description:
%   Approximates a Gaussian integral using the Gauss-Hermite method
%   in multiple dimensions:
%     int f(x) N(x | m,P) dx

% History:
%   2009 - Initial version (ssarkka)
%   2010 - Partially rewritten (asolin)

% Copyright (c) 2010 Simo Sarkka, Arno Solin
%
% This software is distributed under the GNU General Public 
% Licence (version 2 or later); please refer to the file 
% Licence.txt, included with the software, for details.

%% THEORY OF OPERATION:
%
% Consider the multidimensional integral:
%
%  int f(x1,...,xn) exp(-x1^2-...-xn^2) dx1...dxn
%
% If we form quadrature for dimension 1, then
% for dimension 2 and so on, we get
%
% = sum w_i1 int f(x1^i1,...,xn) exp(-x2^2-...-xn^2) dx2...dxn
% = sum w_i1 w_i1 int f(x1^i1,x2^i2,...,xn) exp(-x3^2-...-xn^2) dx3...dxn
% = ...
% = sum w_i1 ... w_in f(x1^i1,...,xn^in) 
%
% Now let X = [x1,...,xn] and consider
%
%  1/(2pi)^{n/2}/sqrt(|P|) int f(X) exp(-1/2 (X-M)' iP (X-M)) dX
%
% Let P = L L' and change variable to
%
%   X = sqrt(2) L Y + M  <=>  Y = 1/sqrt(2) iL (X-M)
%
% then dX = sqrt(2)^n |L] dY = sqrt(2)^n sqrt(|P]) dY i.e. we get
%
%  1/(2pi)^{n/2}/sqrt(|P|) int f(X) exp(-1/2 (X-M)' iP (X-M)) dX
%   = 1/(pi)^{n/2} int f(sqrt(2) L Y + M) exp(-Y'Y) dY
%   = int g(Y) exp(-Y'Y) dY
%
% which is of the previous form if we define
%
%   g(Y) = 1/(pi)^{n/2} f(sqrt(2) L Y + M)
% 

%% The Gauss-Hermite cubature rule

  % The hermite polynomial of order n
  p = hermitepolynomial(n);
  
  % Evaluation points
  x = roots(p);
  
  % 1D coefficients
  Wc = pow2(n-1)*factorial(n)*sqrt(pi)/n^2;
  p2 = hermitepolynomial(n-1);
  W  = zeros(1,n);
  for i=1:n
     W(i)=Wc*polyval(p2,x(i)).^-2;
  end
    
  d = size(m,1);
  if d == 1
    x = x';
  end
  
  % Generate all n^d collections of indexes by
  % transforming numbers 0...n^d-1) into n-base system
  % and by adding 1 to each digit
  num = 0:(n^d-1);
  ind = zeros(d,n^d);
  for i=1:d
      ind(i,:) = rem(num,n)+1;
      num = floor(num / n);
  end
  
  % Form the sigma points and weights
  L = chol(P)';
  SX = sqrt(2)*L*x(ind)+repmat(m,1,size(ind,2));
  W = prod(W(ind),1); % ND weights
  
  % Evaluate the function at the sigma points
  if ischar(f) || strcmp(class(f),'function_handle')
      if nargin < 5
         F = feval(f,SX);
      else
         F = feval(f,SX,param);
      end
  elseif isnumeric(f)
         F = f*SX;
  else
      if nargin < 5
         F = f(SX);
      else
         F = f(SX,param);
      end
  end

  % Evaluate the integral
  I = 1/sqrt(pi)^d*sum(F.*repmat(W,size(F,1),1),2);