function [prt,sdt,prs,sds]=gp4_iden_ssts_matlab(p,M,yy,M1,N,NT,ord,ff)

% GP4_IDEN_SSTS_MATLAB      Performs the calculations of both state-space
% and times-series domains for the predictions and the standard deviations.
% The operations are performed using standard MATLAB functions.
%
%   [prt,sdt,prs,sds] = GP4_IDEN_SSTS_MATLAB(p,M,yy,M1,N,NT,ord,ff)
%
% Inputs:
%       p   = Hyperparameters, [a1 a2 gts g1...gk v]
%       M   = Input state-space explanatory variables
%       yy  = Outcome for every explanatory variables
%       M1  = Input time-series explanatory variables
%       N   = Predicting input state-space explanatory variables
%       NT  = Predicting input time-series explanatory variables
%       ord = Order of reduced state-space size
%       ff  = Fudge Factor, e
%
% Outputs:
%   prt  = Prediction of the time-series domain
%   sdt  = Standard deviation of the time series domain
%   prs  = Prediction of the state-space domain
%   sds  = Standard deviation of the state-space domain
%
%   Note:
%           
%           Q = a1* [ Qf+v(I+eI)      vI    ]
%                   [    vI      Qg+v(I+eI) ]
%
%       such that, Qf = exp(-0.5*Zf^2*d)
%                  Qg = a2*exp(-0.5*Zg^2*dg)
%                   e = machine precision
%
%
% (C) Gaussian Process SSTS Toolbox 1.0
% (C) Copyright 2006-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

D0=size(M,2);
N1=size(M,1);
a1=p(1);
a2=a1*p(2);
g1=p(3);
g2=p(4:D0+3);
v =a1*p(D0+4);
vxp=v*(1+ff);
expA=a1*constr_cov(M1,M1,g1);
expB=a2*constr_cov(M,M,g2);
Qa=diag_add(expA,vxp);
Qb=diag_add(expB,vxp);
T=get_ord(ord,N1);
N2=size(T,1);
Q=[Qa v*T';v*T T*Qb*T'];
clear expA expB Qa Qb;
y=[yy;T*yy];

DexpA=[a1*constr_cov(M1,NT,g1);zeros(N2,size(NT,1))];
DexpB=[zeros(N1,size(N,1));a2*constr_cov(T*M,N,g2)];

invQy=Q\y;
prt=DexpA'*invQy;
prs=DexpB'*invQy;
sdt=sqrt(a1-diag(DexpA'*(Q\DexpA)));
sds=sqrt(a2-diag(DexpB'*(Q\DexpB)));


function expA=constr_cov(M,N,g)
N0=size(M,1);
Nt=size(N,1);
expT=repmat(permute(M,[1,3,2]),[1,Nt,1])-repmat(permute(N',[3,2,1]),[N0,1,1]);
expC=(expT.^2).*repmat(permute(g',[3 2 1]),[N0,Nt,1]);
expA=exp(-0.5*sum(expC,3));

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