clear all;
close all;
clc;


% get the data file
% assigning to the variable

load lorenz_x;
x = lorenz_x;

% defining some constants
% data length: N
% signular value decomposition size, which is taken as much greater than the dimension (number of Lyapunov exponents) expected: svd_size
% number of embedded data(dynamical) vectors points: M
% number of evolve steps: a 
%
N = length(x); 
svd_size = 10;
M = N - svd_size;
a = 5;

% the threshold for the svd decomposition: svd_thr
% the reinitialization after: jump
%
svd_thr = 100;
jump = 1;

% basic embedding done
% 
for i = 1:M
    X1(i, 1:svd_size) = x(i: (i + svd_size - 1));
end;

% singular value decomposition step
% to concentrate on the correction dimension, 
% we will support those many number of columns which have 
% singular values higher than the noise floor level
% Matrix of singular vectors (eigenvectors): S
% Matrix of singular vectors (convariance matrix): C
% Matrix of singular values (diagonal elements): sigma
%
[S, sigma, C] = svd(X1);

% debug command
%
%plot(diag(sigma));

% since the singular values cotnained in sigma are 
% the root mean square projections onto the basis vectors, 
% the number of non-zero values in signma should give the 
% number of degress of freedom in embedded vector variables. 
% Noise on the data, however, will generate the spurious degress 
% of freedom, evident as extra non-zero values in sigma. 
% Hence, svd_thr is defined as some higher value than zero, 
% to account for this spurious behavior
%
D_embed = sum(diag(sigma) > svd_thr);

% Now the embedded data vector is embedded to the "correct" dimensions
% 
X = X1 * C(:, 1:D_embed);

% other variables required for the computation are described here
% during QR decomposition of the Tangent maps, initial value of Q 
% is taken as identity matrix
% we will get only D_embed number of Lyapunov exponents in the 
% lypunov spectrum
% 
Q = eye(D_embed);
lyp_est = zeros(1, D_embed);

% the count over which the averaging is done needs to be stored
% 
n = 0;

% starting with the core code of computation of Lyapunov exponent
% 
for i = 1: jump: (M - (a + 1))

% to have non-intersecting attractors, 
% the dimensions for correct computation should be at least given by the 
% condition below
% In practice, we might not follow the condition 'strictly'.
    num_nearest_neighbours = 2 * D_embed + 1;

% for a single step calculation
% 
    for j = 1: (M - (a + 1))

       % to include all points which are collected within the neighbourhood 
       % without considering the distance from the point to itself :)
       % self-distance carries no meaning, it will always be the shortest :)
       % 
       if (j ~= i)
            d(j) = sqrt(sum((X(i, :) - X(j, :)).^2));
        end;
    end;

    % getting the distance and the index sorted
    % 
    [Y, I] = sort(d);

    % pick of the nearest neighbours only
    %
    I = I(1:num_nearest_neighbours);
    
    % note the maximum distance considered 
    % right now we are not using it, but, this is part of the tweaking 
    % that can be done to make the algorithm perform "correct"
    %
    epsilon = Y(num_nearest_neighbours);
    
    % to calculate the trajectory matrix
    % 
    for k = 1:num_nearest_neighbours
        B_i(k, :)=X(I(k), :) - X(i, :);
        B_i_plus_a(k, :) = X(I(k)+a, :) - X(i + a, :);
    end;
    T=(pinv(B_i) * B_i_plus_a)';
    [Q,R] = qr(T*Q);
    if (sum((diag(R) > 0)) == D_embed)
        lyp_est = lyp_est + log(diag(R))';
        n = n + 1;
    end;
end;

% expected results are:
% lambda1 = 1.50,   lambda2 = ZERO,   lambda3 = -22.46
% but for 1000 points we are getting the values as
% lambda1 = 0.1807, lambda2 = 0.0015, lambda3 = -1.0934
% n = total number of iterations (which satisfied the step 121) = 553;
% 
lyp_est = lyp_est / n

% in such a case, we have following things to do
% 1. increase the data length
% 2. remove the transient at the start
% 3. 