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

% GP5_FG_DERV           Obtains the derivative of the posterior given the
% joint probability distribution of the Gaussian process.
%
% [pr,sd] = GP5_FG_DERV(w1,w2,ka,kd,kv,y,ww1,ww2)
%
% This function computes the prediction and standard deviation of the
% derivative observation of the Gaussian process, given the prior
% information. 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])
%
% 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}
%       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  = 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);

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;
clear Lf Lg

L=ecov(w1,ww1,[ka(1) kd(1:D1)]);
pr=zeros(n1,D1+D2);
sdev=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=-kd(k)*Z2.*L;
    pr(:,k)=Cy11'*invQy;
    sdev(:,k)=sqrt(kd(k)*ka(1)-diag(Cy11'*invQ*Cy11));
end
L=ecov(w2,ww2,[ka(2) kd(D1+1:D1+D2)]);
for k=1:D2
    Z2=repmat(permute(ww2(:,k)',[3,2 1]),[N0,1,1])-repmat(permute(w2(:,k),[1,3,2]),[1,n1,1]);
    Cy11=-kd(k+D1)*Z2.*L;
    pr(:,k+D1)=Cy11'*invQy;
    sdev(:,k+D1)=sqrt(kd(k+D1)*ka(2)-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=diag_add(Q,v)
for i=1:size(Q,1)
    Q(i,i)=Q(i,i)+v;
end