Daniel@0: function covf = gpcovarf(net, x1, x2)
Daniel@0: %GPCOVARF Calculate the covariance function for a Gaussian Process.
Daniel@0: %
Daniel@0: %	Description
Daniel@0: %
Daniel@0: %	COVF = GPCOVARF(NET, X1, X2) takes  a Gaussian Process data structure
Daniel@0: %	NET together with two matrices X1 and X2 of input vectors,  and
Daniel@0: %	computes the matrix of the covariance function values COVF.
Daniel@0: %
Daniel@0: %	See also
Daniel@0: %	GP, GPCOVAR, GPCOVARP, GPERR, GPGRAD
Daniel@0: %
Daniel@0: 
Daniel@0: %	Copyright (c) Ian T Nabney (1996-2001)
Daniel@0: 
Daniel@0: errstring = consist(net, 'gp', x1);
Daniel@0: if ~isempty(errstring);
Daniel@0:   error(errstring);
Daniel@0: end
Daniel@0: 
Daniel@0: if size(x1, 2) ~= size(x2, 2)
Daniel@0:   error('Number of variables in x1 and x2 must be the same');
Daniel@0: end
Daniel@0: 
Daniel@0: n1 = size(x1, 1);
Daniel@0: n2 = size(x2, 1);
Daniel@0: beta = diag(exp(net.inweights));
Daniel@0: 
Daniel@0: % Compute the weighted squared distances between x1 and x2
Daniel@0: z = (x1.*x1)*beta*ones(net.nin, n2) - 2*x1*beta*x2' ... 
Daniel@0:   + ones(n1, net.nin)*beta*(x2.*x2)';
Daniel@0: 
Daniel@0: switch net.covar_fn
Daniel@0: 
Daniel@0:   case 'sqexp'		% Squared exponential
Daniel@0:     covf = exp(net.fpar(1) - 0.5*z);
Daniel@0: 
Daniel@0:   case 'ratquad'	% Rational quadratic
Daniel@0:     nu = exp(net.fpar(2));
Daniel@0:     covf = exp(net.fpar(1))*((ones(size(z)) + z).^(-nu));
Daniel@0: 
Daniel@0:   otherwise
Daniel@0:     error(['Unknown covariance function ', net.covar_fn]);  
Daniel@0: end