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

% GP5_FG_LINEART_DERV       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_DERV(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'
%
% Normally, 95% confidence intervals corresponds to 2 times standard
% deviation, sd, for a Gaussian distribution.
%
% 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  = Prediction of the derivative
%       sd  = Standard deviation for the derivative
%
%
% (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);

meanw2=mean(w2);
w2=w2-repmat(meanw2,N0,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

L=ecov(w1,ww1,[a d]);
pr=zeros(n1,D1+D2);
sd=zeros(n1,D1+D2);
for k=1:D1
    Z2=repmat(permute(ww1(:,k)',[3 2 1]),[N0,1,1])-repmat(permute(w1(:,k),[1 3 2]),[1,n1,1]);
    Cy11=-d(k)*Z2.*L;
    pr(:,k)=Cy11'*invQy;
    sd(:,k)=sqrt(a*d(k)-diag(Cy11'*invQ*Cy11));
end
for k=1:D2
    Cy11=repmat(w(k)*w2,[1,n1]);
    pr(:,k+D1)=Cy11'*invQy;
    sd(:,k+D1)=sqrt(w(k)-diag(Cy11'*invQ*Cy11));
end


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