function [pr,sd]=gp5_fg(w1,w2,ka,kd,kv,y,ww1,ww2)

% GP5_FG        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(w1,w2,ka,kd,kv,y,ww1,ww2)
%
% 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:
%       w1  = Explanatory variable, of dimension {N0 x D1}
%       w2  = Explanatory variable, of dimension {N0 x D2}
%       ka  = Hyperparameter of 'a' = [a1 a2]
%       kd  = Hyperparameter of 'd' = [D1 D2], where D=diag{g1,g2,...,gd}
%       kv  = 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,D1]=size(ww1);
n0=size(w1,1);
if (n0~=size(w2,1)) || (n1~=size(ww2,1))
    disp('Incorrect dimension');
    return
end
D2=size(ww2,2);
n2=n0+n0;

Lf=ecov(w1,w1,[ka(1) kd(1:D1)]);
Lg=ecov(w2,w2,[ka(2) kd(D1+1:D1+D2)]);
invQ=inv(diag_add(Lf+Lg,kv));
invQy=invQ*y;

fg=sum((Lf-Lg)'*invQy)/n2;

Jf=ecov(w1,ww1,[ka(1) kd(1:D1)]);
Jg=ecov(w2,ww2,[ka(2) kd(D1+1:D1+D2)]);
pr=[Jf'*invQy-repmat(fg,n1,1) Jg'*invQy+repmat(fg,n1,1)];
% fg=f-g;
% pr=[f g]+repmat(sum([-fg fg],1)/n2,n1,1);

tmpf=invQ*Lf;
tmpg=invQ*Lg;
Cff22=sum(sum(diag_add(Lg,kv)*tmpf))/(n2^2);
Cgg22=sum(sum(diag_add(Lf,kv)*tmpg))/(n2^2);
Cgf22=-sum(sum(Lg'*tmpf))/(n2^2);
Cfg22=-sum(sum(Lf'*tmpg))/(n2^2);

vf=Jf'*invQ;
vg=Jg'*invQ;
Cff21=sum(diag_add(Lg,kv)*vf',1)/n2;
Cgg21=sum(diag_add(Lf,kv)*vg',1)/n2;
Cgf21=-sum(tmpg'*Jf,1)/n2;
Cfg21=-sum(tmpf'*Jg,1)/n2;

% Cff12=sum(vf*diag_add(Lg,kv),2)/n2;
% Cgg12=sum(vg*diag_add(Lf,kv),2)/n2;
% Cgf12=-sum(Jg'*tmpf,2)/n2;
% Cfg12=-sum(Jf'*tmpg,2)/n2;
Cff12=Cff21';
Cgg12=Cgg21';
Cgf12=Cfg21';
Cfg12=Cgf21';

Cff11=repmat(ka(1),n1,1)-diag(vf*Jf);
Cgg11=repmat(ka(2),n1,1)-diag(vg*Jg);
% Cgf11=-vg*Jf
% Cfg11=-vf*Jg;

v1=Cff11-Cff21'-Cff12+repmat(Cff22,n1,1);
v2=repmat(Cgg22,n1,1);
v3=Cgf21'-repmat(Cgf22,n1,1);
v4=Cfg12-repmat(Cfg22,n1,1);
sd_f=v1+v2+v3+v4;

v1=repmat(Cff22,n1,1);
v2=Cgg11-Cgg21'-Cgg12+repmat(Cgg22,n1,1);
v3=Cgf12-repmat(Cgf22,n1,1);
v4=Cfg21'-repmat(Cfg22,n1,1);
sd_g=v1+v2+v3+v4;

sd=sqrt([sd_f sd_g]);


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