function [logmax,dL] = nolin_grad(x,M,y)

% NOLIN_GRAD        is an optimisation script for the training of Gaussian
% process prior models.
%
% [LOGMAX,DL] = NOLIN_GRAD(x,M,y)
%
% This optimisation script uses only Gradient information to perform
% optimisation for Gaussian process prior models, by adapting
% hyperparameters to minimise negative the log-likelihood function,
%
%                   logmax = log(det(Q)) + y'*inv(Q)*y
%
% The prior covariance function used is of the form,
%
%                   C(zi,zj) = a*{exp([zi-zj]'*D*[zi-zj]) + v}
%
% Inputs:
%       x   = Hyperparameters, [a g1 g2 ... gd v], D=diag{g1,g2,...,gd}
%       M   = Input explanatory variable, of dimension {N0 x D0}
%       y   = Outcome (vector) from every input of the explantory variable
%
% Outputs:
%       LOGMAX  = Negative log-likelihood function value
%       DL      = Gradient vector values
%
% This script does not utilise user-supplied Hessian information. To 
% utilise user-supplied Hessian information, choose NOLIN_HESS instead. 
% Note that, Hessian is useful in speeding up optimisation routines.
%
%
% (C) Gaussian Process Toolbox 1.0
% (C) Copyright 2005-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

[N0,L0]=size(M);
a=exp(x(1));
g=exp(x(2:L0+1));
V=exp(x(2+L0));
expC=repmat(permute(M,[1,3,2]),[1,N0,1])-repmat(permute(M',[3,2,1]),[N0,1,1]);
expC=(expC.^2).*repmat(permute(g,[3 2 1]),[N0,N0,1]);
expA=a*exp(-0.5*sum(expC,3));
Q=diag_add(expA,a*V,N0);
invQ=inv(Q);
invQy=invQ*y;
yQy=y'*invQy;
[cholC,pp]=chol(Q);
if pp==0, logdet=2*sum(log(diag(cholC)));
else logdet=sum(log(svd(Q)));
end
logmax=logdet+yQy;
clear cholC;

dL=zeros(L0+2,1);
dL(1,1)=N0-yQy;
dL(L0+2,1)=a*V*(trace(invQ)-invQy'*invQy);
for k=1:L0
    Q=-0.5*expC(:,:,k).*expA;
    dL(1+k,1)=traceAxxB(Q,invQ)-invQy'*Q*invQy;
end

function Q=diag_add(Q,v,N0)
for pp=1:N0
    Q(pp,pp)=Q(pp,pp)+v;
end

function trace_AB=traceAxxB(A,B)
trace_AB=0;
dimAB=size(A,1);
for i=1:dimAB
    trace_AB=trace_AB+A(i,:)*B(:,i);
end;