function [minlambda,maxlambda,mine,maxe, H, R2,integralx,integraly]=plot_sdle(input_data, max_orig_scale,doplot,lambdaskip)

% Postprocessing script to get Scale-dependent Lyapunov exponent (SDLE) curve
% Copyrighted by Jing Hu and Jianbo Gao, ECE department, University of Florida

% input_data: the data you want to apply SDLE to. Actually this is the
% output file from the computation of the Time-dependent exponent (TDE)
% curves. The C and Fortran codes for computing TDE can be downloaded from 
% http://www.gao.ece.ufl.edu/GCTH_Wileybook/programs/lambda_k_curves/
% The output file of TDE curves has the following format:
% the first row stores the number of pairs of vectors within each shell
% the rest rows store the information of \Lambda(k) curves
% each column corresponds to one shell used for computing TDE curves.
% Another input parameter is max_orig_scale, which is the max scale
% chosen for computing the TDE curve, the actual shell size is
% 2^(max_orig_scale).

% close all the open figure windows
%close all;

if (nargin<3)
    doplot=1;
end
if (nargin<4)
    lambdaskip=1;
end

% get the number of pairs of vectors within each shell
num_of_pair = input_data(1, :);

row_num = size(input_data, 1) - 1;
shell_num = size(input_data, 2);

% to determine the scale of each shell
orig_scale = max_orig_scale : max_orig_scale + shell_num - 1;
data_total = input_data(2:end, :);

% You may adjust the start point for plotting the sdle curves here. The
% purpose is to exclude points V_i and V_j (involved in the computation 
% of the time-dependent exponent curves) that are very close in time, 
% since those points may be on the same orbit. A rule of thumb suggested by
% us is to choose start_pt greater than the embedding window size used for 
% computing the time-dependent exponent curves. For more details, please 
% refer to our book: "Multiscale Analysis of Complex Time Series:
% Integration of Chaos and Random Fractal Theory, and Beyond" 
% http://www.wiley.com/WileyCDA/WileyTitle/productCd-0471654701,descCd-authorInfo.html 
% By default it is set to 2 (corresponding to the embedding parameters m = 2, L = 1).
start_pt = 2; 

% You may set the end point for plotting the sdle curves. The default value
% is the total length of the data.
end_pt = size(input_data, 1) - 2;

% exclude shells with too few pairs of vectors. You may change the threhold
% value depending on your data
pair_num_thred = 1;

% note even though here I labeled it as max shell, it actually corresponds
% to the smallest scale resolvable
max_shell = max(find(num_of_pair > pair_num_thred));
minlambda=+inf;
maxlambda=-inf;
mine=+inf;
maxe=-inf;
masxxx=[];
masyyy=[];
if (doplot)
    figure;
    h2=axes;
    figure;
    h3=axes;
end

for i = 1 : max_shell
    data = data_total(:, i);
    xxx = log10(exp(log(2^(-orig_scale(i)/2)) + data(1:end - 1)));
    yyy = diff(data);
    minlambda=min([minlambda;yyy(lambdaskip:length(yyy))]);
    maxlambda=max([maxlambda;yyy(lambdaskip:length(yyy))]);
    mine=min([mine;xxx(lambdaskip:length(xxx))]);
    maxe=max([maxe;xxx(lambdaskip:length(xxx))]);
    masxxx=[masxxx;xxx];
    masyyy=[masyyy;yyy];
    integral2=xxx(1)+cumsum(yyy);

    if (doplot)
        axes(h2);
        plot(xxx(start_pt:end_pt), yyy(start_pt:end_pt), '.-');
        hold on;
        axes(h3);
        plot(start_pt:end_pt, integral2(start_pt:end_pt), '.-');
        hold on;

    end
    
end
integralx=start_pt:end_pt;
integraly=integral2(start_pt:end_pt);
coefs=polyfit(masxxx,masyyy,1);
a1=coefs(1);
a0=coefs(2);
yyy_teor=a0+a1*masxxx;

u2=masyyy-yyy_teor;
R2=(var(masyyy)-var(u2))/var(masyyy);
if (doplot)
    axes(h2);
    plot(masxxx,yyy_teor,'r');
    title(['SLDE, H=' num2str(-a1) ', R^2=' num2str(R2)]);
    set(gca, 'FontSize', 14);
    ylabel('\lambda(\epsilon)', 'FontSize', 14);
    xlabel('log_{10} \epsilon', 'FontSize', 14);
    axes(h3);
    set(gca, 'FontSize', 14);
    ylabel('ln(\epsilon_t)', 'FontSize', 14);
    xlabel('t', 'FontSize', 14);
    while (strcmp(questdlg('Провести линию наклона не по всем точкам?','','Да', 'Нет','Да'),'Да')==1)
        axes(h2);
        [mouseX,mouseY]=ginput(2);
        xmin=min(mouseX(1),mouseX(2));
        xmax=max(mouseX(1),mouseX(2));
        nesessary_points_coords=find(masxxx>xmin);
        nesessary_pointsx=masxxx(nesessary_points_coords);
        nesessary_pointsy=masyyy(nesessary_points_coords);
        nesessary_points_coords2=find(nesessary_pointsx<xmax);
        nesessary_pointsx=nesessary_pointsx(nesessary_points_coords2);
        nesessary_pointsy=nesessary_pointsy(nesessary_points_coords2);
        coefs=polyfit(nesessary_pointsx,nesessary_pointsy,1);
        a1=coefs(1);
        a0=coefs(2);
        yyy_teor=a0+a1*masxxx;
        yyy_teor2=a0+a1*nesessary_pointsx;
        u2=nesessary_pointsy-yyy_teor2;
        R2=(var(nesessary_pointsy)-var(u2))/var(nesessary_pointsy);
        axes(h2);
        plot(masxxx,yyy_teor,'r');
        title(['SLDE, H=' num2str(-a1) ', R^2=' num2str(R2)]);
        
    end
end
H=-a1;