function Q = expacov(M,N,p)

% EXPACOV       Computes the covariance matrix of a squared exponential
% covariance function, C(zi,tj).
%
% [Q] = EXPACOV(M,N,p)
%
% This script is optimised to calculate the covariance matrix, given the
% covariance function of the form,
%
%           C(zi,zj) = a*exp([zi-tj]'*D*[zi-tj])
%
% such that a and D are hyperparameters of the covariance function, and zi
% is the i-th entry of vector z of M, and tj is the j-th entry of vector t
% of N.
%
% It returns Q in matrix form.
%
% Inputs:
%       M   = Input explanatory variable of dimension {n1 x d}
%       N   = Input explanatory variable of dimension {n2 x d}
%       p   = Hyperparameters, [a g1 g2 ... gd], D=diag{g1,g2,...,gd}
%
% Output:
%       Q   = Resulting covariance matrix
%
%
% (C) Gaussian Process Toolbox 1.0
% (C) Copyright 2005-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

[N0,D0]=size(M);
Nt=size(N,1);
a=p(1);
g=p(2:D0+1);
Z2=(repmat(permute(M,[1,3,2]),[1,Nt,1])-repmat(permute(N',[3,2 1]),[N0,1,1])).^2;
Q=a*exp(-0.5*sum(Z2.*repmat(permute(g',[3,2,1]),[N0,Nt,1]),3));