function [pr,sd]=gp5_fg_lineart(w1,w2,a,d,w,v,y,ww1,ww2)

% GP5_FG_LINEART        provides a prediction procedure for obtaining the 
% posterior of the Gaussian process prior model, with 2 indendent Gaussian 
% processes of independent explanatory variables, i.e. y~f(u,v)+g(w)
%
% [pr,sd] = GP5_FG_LINEART(w1,w2,a,d,w,v,y,ww1,ww2)
%
% 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*Zk'
%
% Inputs:
%       w1  = Explanatory variable, of dimension {N0 x D1}
%       w2  = Explanatory variable, of dimension {N0 x D2}
%       a   = Hyperparameter of 'a' = [a1 a2]
%       d   = Hyperparameter of 'd' = [D1 D2], where D=diag{g1,g2,...,gd}
%       w   = Hyperparameter of 'w' = {w1,w2,...,wd}
%       v   = Hyperparameter of 'v' (noise variance)
%       y   = Outcome (vector) from every input of the explantory variable
%       ww1 = Predicted explanatory variable, of dimension {N0 x D1}
%       ww2 = Predicted explanatory variable, of dimension {N0 x D2}
%
% Outputs:
%       pr  = Predictions
%       sd  = Standard deviation of the predictions
%
%
% (C) Multiple Gaussian Processes Toolbox 1.0
% (C) Copyright 2006, Keith Neo <keith.neo@nuim.ie>
%                     http://www.hamilton.ie/keith

n1=size(ww1,1);
if (size(w1,1)~=size(w2,1)) || (n1~=size(ww2,1))
    disp('Incorrect dimension');
    return
end

meanw2=mean(w2);
w2=w2-repmat(meanw2,size(w2,1),1);
ww2=ww2-repmat(meanw2,n1,1);

Lf=ecov(w1,w1,[a d]);
Lg=elin(w2,w2,w);
invQ=inv(diag_add(Lf+Lg,v));
invQy=invQ*y;
clear Lf Lg

Lf=ecov(w1,ww1,[a d]);
Lg=elin(w2,ww2,w);

pr=[Lf'*invQy Lg'*invQy];

stdf=a-diag(Lf'*invQ*Lf);
stdg=sum(ww2.^2.*repmat(w,n1,1),2)-diag(Lg'*invQ*Lg);
sd=sqrt([stdf stdg]);


function Q=ecov(M,N,p)
[N0,D0]=size(M);
Nt=size(N,1);
g=p(2:D0+1);
Z2=(repmat(permute(M,[1,3,2]),[1,Nt,1])-repmat(permute(N',[3,2 1]),[N0,1,1])).^2;
Q=p(1)*exp(-0.5*sum(Z2.*repmat(permute(g',[3,2,1]),[N0,Nt,1]),3));

function Q=elin(M,N,p)
[N0,D0]=size(M);
Nt=size(N,1);
Q=zeros(N0,Nt);
for k=1:D0
    Q=p(k)*M(:,k)*N(:,k)'+Q;
end

function Q=diag_add(Q,v)
for i=1:size(Q,1)
    Q(i,i)=Q(i,i)+v;
end