clear all;
close all;
clc;

load x;
x=x(8000:12000);

N=length(x);
svd_size=15;
a=5;

svd_thr=100;
jump=5;
lyp_exp=[];

%for svd_size=15:5:90
for a=1:2:19
    M=N-(svd_size+1);
    X1=[]; d=[];
    for i=1:M
        X1(i,1:svd_size)=x(i:i+svd_size-1);
    end;
    disp('svd starting');
    [S,sigma,C]=svd(X1);
    %plot(diag(sigma));
    D_embed=3;
    %D_embed=sum(diag(sigma)>svd_thr);
    X=X1*C(:,1:D_embed);
    Q=eye(D_embed);
    %lyp_est=zeros(1,D_embed);
    lyp_est=[];
    svd_size
    disp('svd over');
    n=0;troubles=0;
    for i=1:jump:M-(a+1)
        num_nearest_neighbours=20;%2*D_embed+1;
        for j=1:M-(a+1)
            if (j~=i)
                d(j)=sqrt(sum((X(i,:)-X(j,:)).^2));
            end;
        end;
        d(i)=max(d)+100;
        [Y,I]=sort(d);
        I=I(1:num_nearest_neighbours);
        epsilon=Y(num_nearest_neighbours);
        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,:);
            d1(k)=sqrt(sum((B_i(k,:)).^2));
            d2(k)=sqrt(sum((B_i_plus_a(k,:)).^2));
            if (2*epsilon<d2(k))
                %            display('Here trouble');
                troubles=troubles+1;
                %            stophere;
            end;
        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))'/(0.01*a)];
            n=n+1
            if n>5000
                break;
            end;
        end;
    end;
    lyp_exp=[lyp_exp; mean(lyp_est)];
    clear X1,X;
end;
plot([1:2:19],lyp_exp)
%lyp_est=lyp_est/n