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

% NOLIN_GRAD_OPTAG      is an optimisation script for the training of
% Gaussian process prior models.
%
% [LOGMAX,DL] = NOLIN_GRAD_OPTAG(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)) + N0*log(y'*inv(Q)*y)
%
% The prior covariance function used is of the form,
%
%                   C(zi,zj) = exp([zi-zj]'*D*[zi-zj]) + v
%
% Inputs:
%       x   = Hyperparameters, [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_OPTAG
% 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);
g=exp(x(1:L0));
V =exp(x(L0+1));

expC=constr_cov(M,g);
expA=exp(-0.5*sum(expC,3));
Q=diag_add(expA,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+N0*log(yQy);
clear cholC

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

function expC=constr_cov(M,g)
N0=size(M,1);
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]);

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;