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

% NOLIN_FG_QUAD3    is an optimisation script for the training of Gaussian
% process prior models. Particularly, the Gaussian Process is 2 Gaussian
% processes with independent explanatory variables. One of them uses a
% Squared Exponential (SE) covariance function, and the other uses a
% Quadratic covariance function.
%
% [logmax,dL] = NOLIN_FG_QUAD3(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,
%                   P1(zi,zj) = a*exp([zi-zj]'*D*[zi-zj])
%                   P2(zi,zj) = wk*Zk^2*Zk^2'
%
% 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);
a=exp(x(1));
v=exp(x(2));
d1=exp(x(3:2+N1)');
d2=exp(x(3+N1:2+N1+N2)');

meanMb=mean(Mb.^2);
Mb=Mb.^2-repmat(meanMb,size(Mb,1),1);

expD=constr_lin(Mb,d2);
expC=constr_cov(Ma,d1);
expA=a*exp(-0.5*sum(expC,3));

Q=diag_add(expA+sum(expD,3),v);
invQ=inv(Q);
invQy=invQ*y;

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

dL=zeros(2+N1+N2,1);
dL(2)=v*(trace(invQ)-invQy'*invQy);
dL(1)=traceaxb(invQ,expA)-invQy'*expA*invQy;
for k=1:N1
    Q=expC(:,:,k).*expA;
    dL(2+k)=-0.5*(traceaxb(invQ,Q)-invQy'*Q*invQy);
end
clear expA expC
for k=1:N2
    dL(2+N1+k)=traceaxb(invQ,expD(:,:,k))-invQy'*expD(:,:,k)*invQy;
end
clear expD


function expC=constr_lin(M,w)
[N0,D0]=size(M);
expC=zeros(N0,N0,D0);
for k=1:D0
    vecM=M(:,k);
    expC(:,:,k)=w(k)*vecM*vecM';
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)
for k=1:size(Q,1)
    Q(k,k)=Q(k,k)+v;
end