% GET_EIGEN_BLOCK_IDEN3B_MGP     Obtain the eigen decomposition of the
% Hermitian matrix, R, for the generating function such that
%
%      [-T   dQ  y]
%  R = [dQ'  0   0]
%      [ y'  0   0]
%
% where T is the positive-definite Toeplitz-like Block matrix, and dQ is 
% symmetric, not necessary positive-defininite and y is the output data.
% This eigen decomposition function script is used for multiple Gaussian
% Processes, i.e. T = (T1+T2+...)+v, dQ = (dQ1+dQ2+...), such that the
% noise variance, v is identical for all stochastic independent components.
%
% [G,J,rowcol,rowcol_n] = GET_EIGEN_BLOCK_IDEN3B_MGP(M,N,y,hyp)
%
% Inputs:   M       : Time-series Inputs of Gaussian process, size nx1
%           N       : Time-series Inputs variables of prediction, size Ntx1
%           y       : Output data, y of size nx1
%           hyp     : hyperparameters of the covariance function [aa gg vv]
%                     where aa=[a1;a2;..], gg=[g1;g2;..], vv=[v;v;v;..]
%
% Outouts:  G       = Generates the Generating function (matrix)
%           J       = J-unitary matrix, J=[I 0;0 -I]
%           rowcol  = block-position of original explantory variables
%           rowcol_n= block-position of prediction explanatory variables
%
% Note that the covariance function used, is the usual convention, of
%       T(Z|a,g,v) = sum {a * [ exp(-0.5*g*[Zi-Zj]^2) ]} + v
%       dQ(Z|a,g,v)= sum {a * [ exp(-0.5*g*[Zi-tj]^2) ]}
%
%
% (C) Gaussian Process Schur Toolbox 1.0
% (C) Copyright 2005-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

function [G,J,rowcol,rowcol_n]=get_eigen_block_iden3b_mgp(M,N,y,hyp)

[g_row,e_dec]=eigen_tool();

rowcol=g_row(M);
rowcol_n=g_row(N);

n=size(M,1);
Nt=size(N,1);
nn=n+Nt+1;
rk2=size(rowcol,2);
rk2_n=size(rowcol_n,2);
rk=1+rk2+rk2_n;
A=zeros(rk);
B=zeros(nn,rk);
B(1:n,1)=get_row_toep(M,hyp,1);
B(1:n,rk2+1)=get_row_derv_n(M,hyp,1,N);
B(1:n,rk)=-y;
B(n+1:nn-1,1)=get_row_derv_m(M,hyp,1,N);
B(nn,1)=-y(1,1);

for k=2:rk2
    p=rowcol(k);
    tmp1=get_row_toep(M,hyp,p);
    tmp3=get_row_toep(M,hyp,p-1);
    tmp2=[0;tmp3(1:end-1,:)];
    B(1:n,k)=tmp1-tmp2;
    tmp1=get_row_derv_m(M,hyp,p,N);
    tmp3=get_row_derv_m(M,hyp,p-1,N);
    tmp2=[0;tmp3(1:end-1,:)];
    B(n+1:nn-1,k)=tmp1-tmp2;
    B(nn,k)=-y(p);
end
for k=2:rk2_n
    p=rowcol_n(k);
    tmp1=get_row_derv_n(M,hyp,p,N);
    tmp3=get_row_derv_n(M,hyp,p-1,N);
    tmp2=[0;tmp3(1:end-1,:)];
    B(1:n,rk2+k)=tmp1-tmp2;
end

tmp0=zeros(nn,rk);
for h=1:rk2
    p=rowcol(h);
    A(h,:)=B(p,:);
    tmp0(p,h)=1;
end
for h=1:rk2_n
    p=rowcol_n(h);
    A(rk2+h,:)=B(n+p,:);
    tmp0(n+p,rk2+h)=1;
end
A(rk,:)=B(end,:);
tmp0(nn,rk)=1;

[G,J]=e_dec(A,B,tmp0);



function X=get_row_toep(M,hyp,col)
H0=size(hyp,1);
X=zeros(size(M,1),1);
tmp=0;
for k=1:H0
    X=X+hyp(k,1)*exp(-0.5*hyp(k,2)*(M-M(col)).^2);
    tmp=tmp+hyp(k,1);
end
X(col,1)=tmp+hyp(1,3);


function X=get_row_derv_n(M,hyp,col,N)
H0=size(hyp,1);
X=zeros(size(M,1),1);
for k=1:H0
    X=X+hyp(k,1)*exp(-0.5*hyp(k,2)*(M-N(col)).^2);
end


function X=get_row_derv_m(M,hyp,col,N)
H0=size(hyp,1);
X=zeros(size(N,1),1);
for k=1:H0
    X=X+hyp(k,1)*exp(-0.5*hyp(k,2)*(N-M(col)).^2);
end