% EIGEN_TOOL    A compilation of Toolbox functions for the use of
% Generalised Schur Algorithm application in Gaussian processes.
%
% [g_row,e_dec,g_toe,g_dqg,g_to3,g_drn,g_drm,g_hes] = EIGEN_TOOL()
%
% The various functions are:
%
%   [rr]    = g_row (M)
%   [G,J]   = e_dec (A,B,T)
%   [X]     = g_toe (M,hyp,col)
%   [X]     = g_dqg (M,hyp,col)
%   [X]     = g_to3 (M,hy3,col)
%   [X]     = g_drn (M,hyp,col,N)
%   [X]     = g_drm (M,hyp,col,N)
%   [X]     = g_hes (M,hyp,col)
%   
% Inputs:   M       : Time-series Inputs or explanatory variables, size nx1
%           A       : Symmetric Matrix A, to form the Eigenschur matrix
%           B       : Matrix B, to form the Eigenschur matrix
%           T       : I-unitary Matrix T, to form the Eigenschur matrix
%           hyp     : Hyperparameters of the covariance function [g v]
%           hy3     : Hyperparameters of the covariance function [a g v]
%           col     : Column position of a vector
%           N       : Posterior Time-series Inputs, size ntx1
%
% Outouts:  G       = Generates the Generating function (matrix)
%           J       = J-unitary matrix, J=[I 0;0 -I]
%           rr      = Block-positions of a quasi-Toeplitz matrix
%           X       = Vector of specific covariance function
%
% Note that the covariance function used, is the usual convention, of
%       C1(Z|a,g,v) = exp(-0.5*g*[Zi-Zj]^2) + v
%       C2(Z|a,g,v) = a * { exp(-0.5*g*[Zi-Zj]^2) + v },   for <g_to3> only
%
% It follows that the Eigenschur matrix, Mx is
%
%   Mx = [  A       I]
%        [B'*B-AA   0]
%
%
% (C) Gaussian Process Schur Toolbox 1.0
% (C) Copyright 2005-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

function [g_row,e_dec,g_toe,g_dqg,g_to3,g_drn,g_drm,g_hes]=eigen_tool()

g_row = @get_rowcolumn;
e_dec = @eig_decomp;
g_toe = @get_row_toep;
g_dqg = @get_row_derv;
g_to3 = @get_row_toep3;
g_drn = @get_row_derv_n;
g_drm = @get_row_derv_m;
g_hes = @get_row_derv2;



function b=get_rowcolumn(M)
n=length(M);
b=1;
p=2;
inter=M(2,1)-M(1,1);
TT=(M(3:n,1)-M(2:n-1,1))/inter;
for k=1:n-2
    if rdn(TT(k),-3)~=1.000
        b(1,p)=k+2;
        p=p+1;
    end
end


function [EG,J]=eig_decomp(A,B,T0)
rk=length(A);
n1=size(B,1);
M=[A eye(rk);B'*B-A*A zeros(rk)];
[V,D]=eig(M);
dd=real(diag(D));   % Removes complex numeric components
V=real(V)+imag(V);
x2=V(1:rk,:);
y2=V(rk+1:2*rk,:);
G=B*x2+T0*y2;
gg=abs(G'*G);
dz=diag(gg);
p=0;
p=sum(dz<1e-307);
q=sum(sum(gg,2)-dz);
if p>0 && q>1e-7
    disp('WARNING.......');
    disp('Orthogonality unsure, possible off-diagonal blocks.');
    disp('Eigen Stabilisation on Eigenschur matrix not performed.');
else
    for k=1:2*rk    % Ensure orthogonality of G
        vt2=zeros(n1,1);
        for h=1:(k-1)
            vt2=vt2+((G(:,h)'*G(:,k))/(G(:,h)'*G(:,h)))*G(:,h);
        end
        G(:,k)=G(:,k)-vt2;
    end
end
e=sqrt(sum(G.^2,1));

[S,I]=sortrows(-dd);    % Reordering of eigen-values
a=sqrt(abs(S));
EG=zeros(n1,2*rk);
for k=1:2*rk
    p=I(k);
    EG(:,k)=G(:,p)*(e(p)\a(k));
end
J=[eye(rk) zeros(rk);zeros(rk) -eye(rk)];


function X=get_row_toep(M,hyp,col)
X=exp(-0.5*hyp(1)*(M-M(col)).^2);
X(col,1)=1+hyp(2);


function X=get_row_derv(M,hyp,col)
tmp=-0.5*hyp(1)*(M-M(col)).^2;
X=exp(tmp).*tmp;


function X=get_row_toep3(M,hyp,col)
X=hyp(1)*exp(-0.5*hyp(2)*(M-M(col)).^2);
X(col,1)=hyp(1)*(1+hyp(3));


function X=get_row_derv_n(M,hyp,col,N)
X=hyp(1)*exp(-0.5*hyp(2)*(M-N(col)).^2);


function X=get_row_derv_m(M,hyp,col,N)
X=hyp(1)*exp(-0.5*hyp(2)*(N-M(col)).^2);


function X=get_row_derv2(M,hyp,col)
tmp=-0.5*hyp(1)*(M-M(col)).^2;
X=tmp.*exp(tmp).*(1+tmp);


function [x,msg] = rdn(x,n)
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;