% GET_EIGEN_BLOCK     Obtain the eigen decomposition of the Hermitian
% matrix, R, for the generating function such that
%
% R=[-T  I]
%   [ I  0]
%
% where T is the positive-definite Toeplitz-like Block matrix
%
% [G,J,rowcol] = GET_EIGEN_BLOCK(M,hyp)
%
% Inputs:   M       : Time-series Inputs of Gaussian process, size nx1
%           hyp     : hyperparameters of the covariance function [a g v]
%
% Outouts:  G       = Generates the Generating function (matrix)
%           J       = J-unitary matrix, J=[I 0;0 -I]
%           rowcol  = block-position 
%
% Note that the covariance function used, is the usual convention, of
%
%       C(Z|a,g,v) = a*exp(-0.5*g*[Zi-Zj]^2) + v
%
% where v is not a*v!
%
%
% (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]=get_eigen_block(M,hyp)

n=size(M,1);
inter=M(2,1)-M(1,1);
b(1,1)=1;
p=2;
for k=2:n
    tmp=(M(k,1)-M(k-1,1))/inter;
    if rdn(tmp,-3)~=1.000
        b(1,p)=k;
        p=p+1;
    end
end
rowcol=b;
clear b;

n=size(M,1);
rk=size(rowcol,2);
A=zeros(rk);
B=zeros(2*n,rk);
B(1:n,1)=-get_row_toep(M,hyp,1);
B(n+1,1)=1;
for k=2:rk
    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)=tmp2-tmp1;
end
for k=1:rk
    for h=1:rk
        p=rowcol(h);
        A(k,h)=B(p,k);
    end
end
M=[A eye(rk);B'*B-A*A zeros(rk)];
[V,D]=eig(M);
x2=-V(1:rk,:);
y2=-V(rk+1:end,:);

tmp0=zeros(size(B));
for k=1:rk
    p=rowcol(k);
    tmp0(p,k)=1;
end
G=B*x2+tmp0*y2;
dd=diag(D);
clear A B tmp0 tmp1 tmp2 tmp3 V D x2 y2;
e=sqrt(sum(G.^2,1));
for k=1:2*rk
    E(:,k)=G(:,k)/e(k);
end
[S,I]=sortrows(-dd);
H=zeros(2*rk);
for k=1:2*rk
    p=I(k);
    H(p,k)=1;
end
for k=1:rk
    a(k,1)=sqrt(-S(k));
    a(rk+k,1)=sqrt(S(rk+k));
end
G=E*H*diag(a);
J=[eye(rk) zeros(rk);zeros(rk) -eye(rk)];

function X=get_row_toep(M,hyp,col)
n=size(M,1);
X=zeros(n,1);
a=hyp(1);
g=hyp(2);
v=hyp(3);
for k=1:n
    if k ~= col
        X(k,1)=a*exp(-0.5*g*(M(k)-M(col)).^2);
    else
        X(k,1)=a*exp(-0.5*g*(M(k)-M(col)).^2)+v;
    end
end

function [x,msg] = rdn(x,n)

%ROUNDN  Rounds input data at specified power of 10
%
%  y = ROUNDN(x) rounds the input data x to the nearest hundredth.
%
%  y = ROUNDN(x,n) rounds the input data x at the specified power
%  of tens position.  For example, n = -2 rounds the input data to
%  the 10E-2 (hundredths) position.
%
%  [y,msg] = ROUNDN(...) returns the text of any error condition
%  encountered in the output variable msg.
%
%  See also ROUND

%  Copyright 1996-2003 The MathWorks, Inc.
%  Written by:  E. Byrns, E. Brown
%   $Revision: 1.9.4.1 $    $Date: 2003/08/01 18:20:03 $

msg = [];   %  Initialize output

if nargin == 0
    error('Incorrect number of arguments')
elseif nargin == 1
    n = -2;
end

%  Test for scalar n

if max(size(n)) ~= 1
   msg = 'Scalar accuracy required';
   if nargout < 2;  error(msg);  end
   return
elseif ~isreal(n)
   warning('Imaginary part of complex N argument ignored')
   n = real(n);
end

%  Compute the exponential factors for rounding at specified
%  power of 10.  Ensure that n is an integer.

factors  = 10 ^ (fix(-n));

%  Set the significant digits for the input data

x = round(x * factors) / factors;