function [pr,sd]=gp4_iden_ssts_ts(x,M,yz,M1,N,ord,ff)

% GP4_IDEN_SSTS_TS      Performs the calculations of the prediction with 
% standard deviations of the time-series domain
% using state-space information, with incorporation of fast algorithms. The
% generalised Schur algorithm is used here.
%
%   [pr,sd] = GP4_IDEN_SSTS_TS(p,M,yy,M1,N,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
%       ord = Order of reduced state-space size
%       ff  = Fudge Factor, e
%
% Outputs:
%   pr  = Prediction of the time-series domain
%   sd  = Standard deviation of the time series 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

[N0,D0]=size(M);
a1 = x(1);
a2 = a1*x(2);
g1 = x(3);
g2 = x(4:3+D0);
v = a1*x(4+D0);
Vxp=v*(1+ff);

G=get_eigen_block(M1,[a1 g1 Vxp]);
[invQ,logdetQf,trQf,invFz]=schur_invssts_c(G,ord,yz);
clear logdetQf trQf G
ss_M=getord(ord,M);
Tz=getord(ord,yz);
T0=size(ss_M,1);
expA2=expacov(ss_M,ss_M,[a2 g2]);
QG=diag_add(expA2,Vxp,T0)-(v*v)*invQ;
MatGF=inv(QG);
Y=MatGF*(Tz-v*getord(ord,invFz));
ct=toep2_genct(M1,[a1 g1 Vxp]);
X=invord(ord,Y,N0);
X=toep2_invcy(ct,X);
X=invFz-v*X;
% ct(1)=a1;
% pr=toe_times_y(ct,X);

pr=zeros(size(N));
for k=1:length(N)
    tmp=exp((-0.5*g1)*(M1'-N(k)).^2);
    pr(k)=a1*(tmp*X);
end
G=get_eigen_block_ssts_tssd(M1,N,[a1 g1 Vxp/a1]);
[invP,NQf]=schur_invssts_cz(G,ord,N0);
sd=sqrt(a1-invP-(v*v)*diag(NQf'*MatGF*NQf));

% tmp1=Vxp-Vxp*Vxp*diag(invQ)-v*v*diag(MatGF);
% tmp2=v*v*(2*Vxp*diag(invQ*MatGF)-Vxp*Vxp*diag(invQ*MatGF*invQ));
% sd=sqrt(tmp1+tmp2);
% if ord>0
%     tmpsd=zeros(T0,ord+1);
%     tmpsd(:,1)=sd;
%     for k=1:ord
%         invQ=schur_invssts_cx(G,ord,k);
%         ss_M=getordx(ord,M,k);
%         expA2=expacov(ss_M,ss_M,[a2 g2]);
%         QG=diag_add(expA2,Vxp,T0)-(v*v)*invQ;
%         MatGF=inv(QG);
%         tmp1=Vxp-Vxp*Vxp*diag(invQ)-v*v*diag(MatGF);
%         tmp2=v*v*(2*Vxp*diag(invQ*MatGF)-Vxp*Vxp*diag(invQ*MatGF*invQ));
%         tmpsd(:,k+1)=sqrt(tmp1+tmp2);
%     end
%     sd=reshape(tmpsd',(ord+1)*T0,1);
% end


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

function z=getordx(x,y,p)
if x==0
    z=y;
else
    x1=x+1;
    D0=ceil(size(y,1)/x1);
    z=zeros(D0,size(y,2));
    for k=1:D0
        pp=x1*k-x+p;
        z(k,:)=y(pp,:);
    end
end

function z=getord(x,y)
if x==0
    z=y;
else
    x1=x+1;
    D0=ceil(size(y,1)/x1);
    z=zeros(D0,size(y,2));
    for k=1:D0
        pp=x1*k-x;
        z(k,:)=y(pp,:);
    end
end


function z=invord(x,y,N1)
Y0=size(y);
if x==0 && Y0(1)==N1
    z=y;
else
    yx1=Y0(1)*(x+1);
    z=zeros(N1,Y0(2));
    for k=1:Y0(2)
        tmpZ=[y(:,k)';zeros(x,Y0(1))];
        tmpZ=reshape(tmpZ,yx1,1);
        z(:,k)=tmpZ(1:N1);
    end
end