function [logmax,dL]=nolin_afagvdfdg(x,Ma,Mb,y)

% NOLIN_AFAGVDFDG   is an optimisation script for the training of Gaussian
% process prior models. Particularly, the Gaussian Process is 2 Gaussian
% processes with independent explanatory variables.
%
% [logmax,dL] = NOLIN_AFAGVDFDG(x,Ma,Mb,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) = P1 + P2 + v
%
% where,
%                   P(zi,zj) = a*exp([zi-zj]'*D*[zi-zj])
%
% Inputs:
%       x   = Hyperparameters, [a1 a2 v D1 D2], where D=diag{g1,g2,...,gd}
%       Ma  = Input explanatory variable, of dimension {N0 x D1}
%       Mb  = Input explanatory variable, of dimension {N0 x D2}
%       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. Hessian
% is approximated by finite-differencing.
%
% External function dependencies:-
%   - TRACEAXB  = calculation of the trace of product of two matrices.
%
%
% (C) Multiple Gaussian Processes Toolbox 1.0
% (C) Copyright 2006, Keith Neo <keith.neo@nuim.ie>
%                     http://www.hamilton.ie/keith

N1=size(Ma,2);
N2=size(Mb,2);
a1=exp(x(1));
a2=exp(x(2));
v=exp(x(3));
d1=exp(x(4:3+N1));
d2=exp(x(4+N1:3+N1+N2));

expC=constr_cov(Ma,d1);
expD=constr_cov(Mb,d2);
expA=a1*exp(-0.5*sum(expC,3));
expB=a2*exp(-0.5*sum(expD,3));
Q=diag_add(expA+expB,v);
invQ=inv(Q);
invQy=invQ*y;

dL=zeros(3+N1+N2,1);
dL(3)=v*(trace(invQ)-invQy'*invQy);

for k=1:N1
    P=expC(:,:,k).*expA;
    dL(3+k)=-0.5*(traceaxb(P,invQ)-invQy'*P*invQy);
end
P=invQ*expA;
dL(1)=trace(P)-y'*P*invQy;
clear expC expA

for k=1:N2
    P=expD(:,:,k).*expB;
    dL(3+N1+k)=-0.5*(traceaxb(P,invQ)-invQy'*P*invQy);
end
P=invQ*expB;
dL(2)=trace(P)-y'*P*invQy;
clear expD expB

[P,k]=chol(Q);
if k==0, logdet=2*sum(log(diag(P)));
else logdet=sum(log(svd(Q)));
end
logmax=logdet+y'*invQy;


function Q=diag_add(Q,v)
for k=1:size(Q,1)
    Q(k,k)=Q(k,k)+v;
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]);