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