function output_result = lambda(x, m, L, k, w, max_shell, num_of_shell);

% input parameters:
%x: input data
%m: embedding dimension 
%L: delay time
%k: number of iterations
%w: minimum interval, related to the decorrelation time. Please refer to
%pages 245, 249-250 of the book "Multiscale analysis of complex time series --- 
%Integration of Chaos and Random Fractal Theory, and Beyond" for more
%details
%max_shell: maximum shell 2^(-max_shell)       
%num_of_shell: number of shells

% output results are stored in the variale named "output_result", it 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.

% load data from the file
%load input_file.txt;

% normalize the input data to the unit interval [0 1]
x_min = min(x);
x_max = max(x);
x = (x - x_min)/(x_max - x_min); % стандартизация в интервал 0..1

% initialization of shells
r = zeros(1, num_of_shell);
for i = 1 : num_of_shell
    r(i) = power(2, -(max_shell + i - 1)); % подсчет массива r(i) =2^(-max_shell+i-1)
end

% initialization
result = zeros(k, num_of_shell);
counter = zeros(k, num_of_shell);
pair_counter = zeros(1, num_of_shell);

num_of_vector = size(x, 1) - (m - 1)*L;

% main loop
for i = 1 : num_of_vector - w - k
    for j = i + w : num_of_vector - k
        % distance ||x(i) - x(j)||^2
        dist = sum((x(i : L : i + L*(m - 1)) - x(j : L : j + L*(m - 1))).^2);
        
        if dist == 0
            continue;
        end
        
        % decide ||x(i)-x(j)||^2 belong to which shell
        il = -1;
        for n = 1 : num_of_shell
            if (dist > r(n))
                break;
            else
                il = n;
            end
        end
        
        if il > 0
            pair_counter(il) = pair_counter(il) + 1;
            for n = 1 : k
                dist_n = sum((x(i + n : L : i + n + L*(m - 1)) - x(j + n : L : j + n + L*(m - 1))).^2);
                if dist_n ~= 0
                    result(n, il) = result(n, il) + log(dist_n/dist)/2;
                    counter(n, il) = counter(n, il) + 1;
                end
            end
        end
    end
end

% average the results
for i = 1 : k
    for j = 1 : num_of_shell
        if counter(i, j) ~= 0
            result(i, j) = result(i, j)/counter(i, j);
        end
    end
end

output_result = zeros(k+1, num_of_shell);
output_result = [pair_counter; result];