%--------------------------------------------------------------------------
% TOOL pwlCopula
% J. Ch. Strelen, strelen@informatik.uni-bonn.de
% %On web: 30.9.09
% This program is copyright the original author and the University of
% Bonn and is published here under the GNU General Public.
% License. (See http://www.fsf.org/licenses/licenses.html)
%
% Code is modified to fit the requirements of ResQue eliminating it's
% interactive features.
%--------------------------------------------------------------------------

function out = pwlCopula(numDimen, numRes, seed, accuracy, numResGen, ...
    data, attNames);
% numDimen  No of dimensions/attributes
% numRes    No of resources
% seed      Random seed (1-128)
% accuracy  Accuracy
% numResGen No of resources to generate
% data      Data to use in generation. Arranged as numResxnumDimen matrix
% attNames  List of attribute names.
%Return numResGenxnumDimen matrix

K = accuracy; %e.g., 10, 100, 4000 Accuracy increase with number 
n_by_K = floor (numRes / K);
n = n_by_K * K; %New n based on K & numRes
D = numDimen;
Ds = D;
shift = D;   %random vectors
mult = 1;
min_D_Ds = min(D,Ds);
Z = data(1:n, :)'; %Get only n vectors

%Empirical copula data (sample frequencies), and marginal distributions:
[Z_sorted nU F_emp] = empVF_v3(n,D,Z);

% F_emp wird nicht verwendet! (EN: not used!)
% Z_sorted wird nur für die inverse Trf. gebraucht (und wird gespeichert)
%          (EN: only for the inverse Trf. used (and stored))
Z_sortiert = Z_sorted;  % Zwausg  28.2.  !!!

%copula (marginal) densities:
delta = 1/K;   %Zwausg
j = zeros(1, D);
summand = 1*realpow(K, D)/n;   %Zwausg

%hashing: Build hash table, 1st phase ----------------------------v
ms =  ceil(0.5*n);   %Zwausg
m = ms + n;   %Zwausg
in_1st_part=0;   % # used places used in 1st part of hash table
in_2nd_part=0;   % # used places used in 2nd part of hash table

next_free_entry = ones(1, D - 1)*(ms + 1);  % in overflow part

% Storage for table:
f_value = zeros(D - 1, m);     %f_value( d-1,...) == F(d)_...
j_values = zeros(D - 1, m, D);   % index tuple
coll_ptr = zeros(D - 1, m);     % collision pointer

theta = sqrt(5) * 16;   %Zwausg   %for hash function

j_zw = zeros(1, D);

for i = 1:n
    if mult == 1
        for ds = 1:D
            j(ds) = ceil((nU(ds, i)-0.00001)/n_by_K );
            % actual index tuple;  rounding errors omitted
        end;
    end
    
	% j % Zwausg   15.2.
    for d = 2:D   % for f^(d)
        % calculate hash adress:
        sigma = 0;
        for ds = 1:d-1
            sigma = (sigma + j(ds))/K;
        end;
        h_adr = floor( ms*( sigma*theta - floor(sigma*theta) ) ) + 1;

        % search entry in hash table:
        j_values(d-1, h_adr, 1:d);
        if j_values(d-1, h_adr, 1) == 0  % empty 1st entry
            i_entry = h_adr;
            in_1st_part = in_1st_part+1;   % one more entry
        else  % collision
            j_zw(1:d) =  j_values(d-1, h_adr, 1:d);
            equal_j = j(1:d) == j_zw(1:d);         % compare 2 index tuples
            while  ~ all(equal_j(1:d))
                h_adr_alt = h_adr;
                h_adr = coll_ptr(d-1, h_adr);
                if h_adr ~= 0
                    j_zw(1:d) =  j_values(d-1, h_adr, 1:d);
                    equal_j = j(1:d) == j_zw(1:d);
                else break;
                end;
            end;
            if all(equal_j(1:d)) % entry found
                i_entry = h_adr;
            else % h_adr == 0  % entry j does not yet exist
                i_entry = next_free_entry(d-1);  % in overflow part
                next_free_entry(d-1) = next_free_entry(d-1) + 1;
                coll_ptr(d-1, h_adr_alt) = i_entry; %Collision chain longer
                in_2nd_part=in_2nd_part+1; % one more entry
            end
        end;
                    
        f_value(d-1, i_entry) = f_value(d-1, i_entry) + summand; % density
        j_values(d-1, i_entry, 1:d) = j(1:d);                 % index tuple
    end;
end; %for i=1:n
% End 1st phase ---------------------------------------------------^

% Restructure for Phase2:-----------------------------------------------v
% storage for triple table:
tt_index = zeros(D - 1, n);  % index tuple
tt_f = zeros(D - 1, n);      % density
tt_s = zeros(D - 1, n);      % cumulative density
pair_index = zeros(1, n);    %
pair_f = zeros(1, n);        %

% 2 more fields for the hash table
Beg = zeros(D-1, m);
End = zeros(D-1, m);         % not really used??? omit?

jj = zeros(1, D);
    
for d=2:D   %f2, f3, ...
	ttp = 0;
	for i = 1:ms % for each collision chain
        if j_values(d-1, i, 1) > 0   % entry is used
            h = i; % pointer to chain elements
            % pointer to first chain element of subchain with specifi
            %c j-tuple j
            k = h; 
                        
            k_next = 0;
            for ds = 1:d-1
                j(ds) = j_values(d - 1, k, ds);%first j-tuple in this chain
            end;
            ii = 0;
            while true %for each node in this chain
                % if j == j(d;1:d-1):
                j_zw(1:d-1) = j_values(d-1, h, 1:(d-1));
                for ds = 1:d-1
                    jj(ds) = j(ds);
                end;
                equal_j = jj(1:(d - 1)) ==  j_zw(1:d-1);
                if all(equal_j(1:(d - 1)))
                    % if j == j(d;1:d-1), specific j-tuple
                    ii = ii + 1;
                    pair_index(ii) = j_values(d-1, h, d);
                    pair_f(ii) = f_value(d-1, h);
                    j_values(d-1, h, d) = -1;  % mark entry as extracted
                else % different j-tuple found
                    if k_next == 0 && (j_values(d-1, h, d) > -1)
                        % not yet extracted
                        k_next = h; %next subchain for exraction
                    end;
                end;
                            
                h = coll_ptr(d-1,h);
                if h == 0 % end of chain =
                    % end of sequence with j-tuples == j
                    nmb_pairs = ii;
                    % sort, store in triple table, this_Beg, this_End:
                    this_Beg = ttp+1;
                    [pair_index_sorted place_old] = ...
                        sort(pair_index(1:nmb_pairs));
                    for ii = 1:nmb_pairs
                        ttp = ttp + 1;
                        tt_index(d-1, ttp) = pair_index_sorted(ii);
                        tt_f(d-1, ttp) = pair_f(place_old(ii));
                        if ii==1
                            tt_s(d-1, ttp) = tt_f(d-1, ttp) ;
                        else
                            tt_s(d-1, ttp) = tt_s(d-1, ttp-1) + ...
                                tt_f(d-1, ttp);
                        end;
                    end;
                    this_End = ttp;
                    Beg(d-1,k) = this_Beg;
                    End(d-1,k) = this_End;

                    % reorganise chain:
                    h = k_next;
                    coll_ptr(d-1, k) = k_next; %new chain containes only a
                    %  single entry for each subchain with specific j-tuple
                    k = h;    % first entry of subchain
                    if h ~= 0
                        for ds = 1:d-1
                            j(ds)=j_values(d-1, k, ds); %novel j-tuple in 
                            % this chain = next subchain
                        end;

                        ii = 0;
                        k_next=0;
                    end;
                end; % end of subchain reached
                if h == 0
                    break;  %all chain-members extracted, all subchains
                end;
            end; %while
        end %if
    end %i
end %d
%restructured -------------------------------------------------------^

% Zwausg: f_value j_values coll_ptr tt_index tt_f tt_s pair_index 
% pair_f Beg End
            
% First, generating dependent uniformly distributed random numbers:
trf = zeros(1, min_D_Ds);
do_gen = -1;

while(do_gen ~=0)
    if do_gen < 0
        replication_nbr = 1;
        ZZStrom = 0;
        while (ZZStrom<1) | (ZZStrom>128)
            ZZStrom = seed;
        end
        seed = SetSeed(ZZStrom); % for [LK] random numbers
        rand('state',ZZStrom);   % 1. ZZStrom
        
        n_gen = numResGen; % Sample Size to generate
        do_gen = 1; %replication is set to 1
        max_stat_dev = zeros(1, do_gen);
        maxAutocorrDev = zeros(1, do_gen);
        
        %Check with user on how to handle interpolation of dimensions------
        for d = 1:min_D_Ds
            str = ['Do inverse transform u(d,:) with Linear interpolation or Discrete values for attribute ', ...
                attNames{1}(d, :)];
            % Construct a questdlg with two options
            choice = questdlg(str, 'Inverse transformation', ...
                'Linear interpolation', 'Discrete', 'Linear interpolation');
            % Handle response
            switch choice
                case 'Linear interpolation'
                    trf(d) = 1;
                 case 'Discrete'
                    trf(d) = 2;
            end
        end %d
    end;
                
    u_gen = zeros(D, n_gen);         %copula
    power_K = realpow(K, D-1);
    D_minus_shift = D - shift;
    for i = 1:n_gen
        %------------------------------------------------------------------
        d = 1;
        % i==1 or no AR process but random vectors
        if (d > D_minus_shift) | (i == 1)
            % or d > D-Ds
            [u_gen(1,i) seed] = ZZ(seed);  % u_gen(1,i) = rand(1,1);
            u_gen_1 = u_gen(1,i);    %Zwausg 30.1.
            u1_gen_bar = ceil(u_gen(1,i)*K);   % delta = 1/K
            j(1) = u1_gen_bar;
        else  % time series, AR process
            u_gen(1,i) = u_gen(1+shift ,i-1);
            j(1) = j(1+shift);
        end;
        %------------------------------------------------------------------
        d = 2;
        % Calculate hash table entry adress h ............................v        
        % Given: d and tuple j(1:d-1)
        % Result: Hash table entry address h.  
        %       h>0 -> tuple found,
        %       h==0 -> tuple not in table, hence f == 0 and s == 0

        % Calculate hash adress h ..................................vv
        sigma = 0;
        for dss = 1:d-1
            sigma = (sigma + j(dss))/K;
        end;
        h = floor(ms*( sigma*theta - floor(sigma*theta) ) ) + 1;
        % End calculate hash adress  ..............................^^
                    
        % while j ~= j(d;1:d-1) follow collision chain:
        for dss = 1:d-1
            jj(dss) = j(dss);
        end;
        % j-tuple in the hash table at addr h
        j_zw(1:d-1) = j_values(d-1, h, 1:(d-1));
        equal_j = jj(1:(d-1)) ==  j_zw(1:d-1);
        while h > 0 && ~all(equal_j(1:(d-1)))
            h = coll_ptr(d-1,h);
            j_zw(1:d-1) = j_values(d-1, h, 1:(d-1));
            equal_j = jj(1:(d-1)) ==  j_zw(1:d-1);
        end;
        % End calculate hash table entry adress  .........................^
        
        ttp = Beg(d-1,h);

        % i==1 or no AR process but random vectors
        if (d > D_minus_shift) | (i==1) 
            % or d > D-Ds
            [u_gen_ seed] = ZZ(seed);   % u_gen_ =    rand(1,1);
            u_power_K = u_gen_ * power_K;
            summ = 0;
            j2 = 0;
            while (summ < u_power_K)
                if ttp<1 | ttp>n | d<2 | d>D   % error
                    ttp
                    d
                end;

                j2 = tt_index(d-1,ttp);
                summ_old = summ;
                summ = tt_s(d-1,ttp);
                f_value = tt_f(d-1,ttp);   %  = f2(u1_gen_bar, j2);
                ttp = ttp+1;
            end;
                        
            u2_check = j2-1;
            u2_gen_bar = j2;
            u_gen_rest = u_gen_ - summ_old/power_K;
            % f_value = f2(u1_gen_bar, u2_gen_bar);
            u2_gen_rest = u_gen_rest*realpow(K,D-2) / f_value; 
            u_gen(2,i) = u2_check*delta + u2_gen_rest;
            j(2) = u2_gen_bar;
        else  % AR process
            j2 = j(2 + shift);             % 18.2. wieder reingetan v
            while j2 ~= tt_index(d-1,ttp); % f_value wird nur fuer d=D-Ds
                ttp = ttp + 1;          % gebraucht? Optimierung: Nur dann
            end;                    % berechnen, auch nur dann ttp, s. oben
            %=f2(u1_gen_bar, j2);  18.2. wieder reingetan ^
            f_value = tt_f(d-1, ttp); 
            
            u_gen(d,i) = u_gen(d + shift ,i-1);
            j(d) = j(d + shift);
        end;
                    
        %------------------------------------------------------------------
        for d = 3:D
            % Calculate hash adress:   ................................v
            % Given: d and tuple j(1:d-1)
            % Result: Hash address h.
            %       h>0 -> tuple found,
            %       h==0 -> tuple not in table, hence f == 0 and s == 0
            sigma = 0;
            for dss = 1:d-1
                sigma = (sigma + j(dss))/K;
            end;
            %Zwausg 30.1.
            h = floor( ms*( sigma*theta - floor(sigma*theta) ) ) + 1;
            % j(1:d-1) Zwausg
            for dss=1:d-1
                jj(dss) = j(dss);
            end;
            % j-tuple in the hash table at hashaddr h
            j_zw(1:d-1) =  j_values(d-1, h, 1:(d-1));
            equal_j = jj(1:(d-1)) == j_zw(1:d-1);
            while h > 0 && ~all(equal_j(1:(d-1)))
                h = coll_ptr(d-1,h);        %Zwausg 30.1.
                j_zw(1:d-1) =  j_values(d-1, h, 1:(d-1));
                equal_j = jj(1:(d-1)) ==  j_zw(1:d-1);
            end;
            % End calculate hash adress  ..............................^
            
            ttp = Beg(d-1,h);        %Zwausg 30.1.
            % 18.2. hierher geschoben v
            if (d > D_minus_shift) | (i==1) 
                % i==1 or no AR process but random vectors or d > D-Ds
                [u_gen_ seed] = ZZ(seed);   % u_gen_ =    rand(1,1);
                summ = 0;
                j_d = 0;
                u_gen_f_dm1 = u_gen_* f_value;
                % = f_d-1(u1_gen_bar,...,u_d-1_gen_bar) 
                %18.2. hierher geschoben ^
                
                while (summ < u_gen_f_dm1)
                    j_d = tt_index(d-1,ttp);
                    summ_old = summ;
                    summ = tt_s(d-1,ttp);
                    % = s(d, u1_gen_bar,...,u_d-1_gen_bar, j_d);
                    f_value = tt_f(d-1,ttp);   
                    %  = f(d, u1_gen_bar,...,u_d-1_gen_bar, j_d);
                    ttp = ttp+1;        %Zwausg 30.1.
                end;
                u_d_check = j_d-1;        %Zwausg 30.1.
                u_d_gen_bar = j_d;        %Zwausg 30.1.
                u_gen_rest = delta*(u_gen_f_dm1 - summ_old); %Zwausg 30.1.
                u_d_gen_rest = u_gen_rest / f_value;  
                %Zwausg 30.1. % =  f(d,u1_gen_bar,..., u_d_gen_bar);
                u_gen(d, i) = u_d_check*delta + u_d_gen_rest;
                
                j(d) = u_d_gen_bar;
            else
                % d <= D_minus_shift; AR process neu: 
                %if (d>D_minus_shift) | (i==1)
                jd = j(d + shift);            % 18.2. reingetan v
                while jd ~= tt_index(d-1,ttp);% f_value wird nur fuer d=D-Ds
                    ttp = ttp+1;        % gebraucht? Optimierung: Nur dann
                end;                % berechnen, auch nur dann ttp, s. oben
                f_value = tt_f(d-1,ttp); 
                %=fd(u1_gen_bar, j2);   18.2. reingetan ^
                u_gen(d,i) = u_gen(d + shift, i-1);
                j(d) = j(d + shift);
            end;
        end; % for d
        
        %------------------------------------------------------------------
    end; % for i
    if n_gen < 26
        u_gen;   %Zwausg 15.2.
    end
    
    % Generating, get dependent random numbers using inverse transfom:
    z_gen = zeros(Ds,n_gen);
    
    for d = (D - Ds + 1):D
        ds = mod( d-1,Ds)+1;
        if trf(ds) == 1
            for j = 1:n_gen
                un = u_gen(d,j) * n;
                i = ceil(un)-1;
                if i <= 0
                    if i < 0
                        i = 0;
                    end;
                    % smallest generated value
                    z_gen(ds,j) = Z_sorted(ds,i+1);
                else
                    z_gen(ds,j) = Z_sorted(ds,i) + ...
                        (un-i)*(Z_sorted(ds,i+1)-Z_sorted(ds,i));
                end;
            end;
        else
            for j = 1:n_gen
                un = u_gen(d,j) * n;
                i = ceil(un);
                if i == 0
                    i = 1;
                end;
                z_gen(ds,j) = Z_sorted(ds,i);
            end;
        end;
    end; %d

    size_z_gen = size(z_gen);
    if n_gen < 26
        z_gen;   %Zwausg 15.2.
    end
    
    %Save generated vector
    nSp = Ds * n_gen;
    x = zeros(1, nSp);
    jj = 0;
    for j=1:n_gen
        for d=1:Ds
            jj = jj + 1;
            x(jj) = z_gen(d,j);
        end
    end

    out = [];
    for i = 1:numResGen %Result--------------------------------------------
        startIdx = (i - 1) * numDimen + 1;
        endIdx = startIdx + numDimen - 1;
        tmpx = x(1, startIdx:endIdx);
        out = [out; tmpx];
    end
    
    % Estimate how correct the model is -
    % Compare the generated random numbers and the sample:
%     [maxDev maxAutocorrDev(replication_nbr) Corr_figs] = ...
%     Deviations(n,n_gen,min_D_Ds,Z,z_gen,mult,Corr_figs);
%     max_stat_dev(replication_nbr) = maxDev;

    do_gen = do_gen - 1;
    
    replication_nbr = replication_nbr +1;
end; %Loop for generating: while(do_gen ~=0)


function [maxDev, max_autocorr_dev, Corr_figs] = ...
    Deviations (n,n_gen,D,S,G,mult,Corr_figs);

d_ds = zeros(1,2*D+D*D);

dev_mean = zeros(1,D);
dev_var = zeros(1,D);
std_dev_sample = zeros(1,D);
std_dev_generated = zeros(1,D);
var_measure_sample = zeros(1,D);
var_measure_generated =zeros(1,D);
corr_measure_sample = zeros(D,D);
corr_measure_generated = zeros(D,D);
dev_corr_measure = zeros(D,D);

maxDev=0;
for d=1:D
    mean_sample(d) = mean(S(d,:));
    mean_generated(d) = mean(G(d,:));
    if mean_sample(d) < 0.00001
        nominator=1;
    else nominator = mean_sample(d);
    end;
    dev = abs(mean_sample(d)-mean_generated(d))/ nominator;
    dev_mean(d) = dev;
end;

max_deviation_of_means = max(dev_mean);

if maxDev <  max_deviation_of_means
    maxDev =  max_deviation_of_means;
end;

for d=1:D
    std_dev_sample(d) = sqrt(var(S(d,:)));
    std_dev_generated(d) = sqrt(var(G(d,:)));
    if mean_sample(d)>0.00001 && mean_generated(d)>0.00001
        var_measure_sample(d) = std_dev_sample(d)  / mean_sample(d);
        var_measure_generated(d) =  std_dev_generated(d)  / mean_generated(d);
    else
        var_measure_sample(d) = std_dev_sample(d);
        var_measure_generated(d) =  std_dev_generated(d) ;
    end
    dev = abs(var_measure_sample(d)-var_measure_generated(d));
    dev_var_measure(d) = dev;
end;

max_deviation_of_var_measure = max(dev_var_measure);
disp([ 'max standard deviation measures = ' num2str(max_deviation_of_var_measure)  ]);
if maxDev <  max_deviation_of_var_measure
    maxDev =  max_deviation_of_var_measure;
end;

for d=1:D-1
    for ds=d+1:D
        if std_dev_sample(d)>0.00001 && std_dev_generated(d) > 0.00001
            co = corrcoef(S(d,:),S(ds,:));
            corr_measure_sample(d,ds) = co(1,2);
            co = corrcoef(G(d,:),G(ds,:));
            corr_measure_generated(d,ds) = co(1,2);
        else
            co = cov(S(d,:),S(ds,:));
            corr_measure_sample(d,ds) = co(1,2);
            co = cov(G(d,:),G(ds,:));
            corr_measure_generated(d,ds) = co(1,2);
        end;
        dev_corr_measure(d,ds)  ...
            = abs(   corr_measure_sample(d,ds) - corr_measure_generated(d,ds)  );
    end
end
some_corr_coeff_dev = zeros(16,16);
some_corr_coeff = zeros(16,16);
some_corr_coeff_gen = zeros(16,16);
fib= zeros(17);
d_values = zeros(16);
ds_values = zeros(16,16);

format short

fib(1)=1;
fib(2)=1;
for i=3:16
    fib(i) = fib(i-2)+ fib(i-1);
end;
fib(17)=100000;

i=1;
d=1;
j=1;
while d<D
    j=1;
    ds=d+1;
    d_values(i)=d;
    while ds<=D
        ds_value(i,j)=ds;
        some_corr_coeff_dev(i,j) = dev_corr_measure(d,ds);
        some_corr_coeff(i,j) = corr_measure_sample(d,ds);
        some_corr_coeff_gen(i,j) = corr_measure_generated(d,ds);
        ds=ds+fib(j+1);
        j=j+1;
    end
    d=d+fib(i+1);
    i=i+1;
end
some_corr_coeff = some_corr_coeff (1:i-1,1:j-1)
some_corr_coeff_dev = some_corr_coeff_dev (1:i-1,1:j-1)


max_dev_corr_measure = max(max(dev_corr_measure))

if maxDev <  max_dev_corr_measure
    maxDev =  max_dev_corr_measure;
end;

% For time series:
max_autocorr_dev=0;
if mult>1
    disp('Autocorrelations for time series: ');
    disp('---------------------------------  ');
    
    CorrAnz=200;
    Corr = zeros(D,CorrAnz);
    Corr_gen = zeros(D,CorrAnz);
    Corr_dev = zeros(D,CorrAnz);
    
    for d=1:D
        jj=1;
        while (jj < n/2) && (jj <= CorrAnz) && (jj < n_gen/2)
            lag=jj;
            Corr_zw(:,:,jj) = corrcoef( S(d, 1:n-lag  ),  S(d, 1+lag:n  )  );
            Corr(d,jj)=Corr_zw(1,2,jj);
            Corr_zw(:,:,jj) = corrcoef( G(d, 1:n_gen-lag  ),  G(d, 1+lag:n_gen  )  );
            Corr_gen(d,jj)=Corr_zw(1,2,jj);
            Corr_dev(d,jj) = abs(Corr(d,jj) - Corr_gen(d,jj));
            if jj==1
                max_Corr_dev(d,jj) = Corr_dev(d,jj);
            else
                max_Corr_dev(d,jj) = max(max_Corr_dev(d,jj-1), Corr_dev(d,jj));
            end
            jj=jj+1;
        end
    end
    
    disp('Maximum absolute value of calculated autocorrelations deviations: ');
    max_autocorr_dev = max(max(Corr_dev))
    
    d=1;
    i=1;
    while d<=D
        if Corr_figs==1
            figure ('name', ['Autocorr Dim= '  num2str(d)]);
            plot(Corr(d,:));
        end
        figure ('name', ['Autocorr, max dev. in 1..x, Dim= '  num2str(d)]);
        plot(max_Corr_dev(d,:));
        
        d=d+fib(i);
        i=i+1;
    end
    Corr_figs=0;
end
disp('MAXIMUM STATISTICAL DEVIATION FOR THESE GENERATED VECTORS ');
disp('   (MEANS, STANDARD DEVIATION measures, CORRELATION measures):  ');
maxDev


%****************************************************
function  [U, y]  =ZZ(y);  % alternative random number generator [Law-Kelton]
% faktor=16807;
% modul=1024*1024*2048-1 = 2147483647;
y=mod ((y*16807), 2147483647 );
U= y/2147483647;


%****************************************************
function y=SetSeed(stream);

seeds=[
    1,  ...
    1550655590, 766698560, 1645116277, 1154667137, 1961833381, 460575272,   ...
    1497719440, 901595110, 354421844, 1800697079, 821079954, 1918430133,   ...
    1736896562, 634729389, 23246132, 1008653149, 1689665648, 1628735393,   ...
    550023088, 1267216792, 314116942, 2095003524, 356484144, 2140958617,   ...
    1116852912, 1701937813, 171817360, 360646631, 1652205397, 605448579,   ...
    1223969344, 1821072732, 1237280745, 2125855022, 935058038, 1151620124,   ...
    1970906347, 66562942, 1503754591, 2007872839, 1870986444, 1375265396,   ...
    470646700, 500432100, 347147950, 929595364, 800443847, 5715478,       ...
    2032092669, 996425830, 1884232020, 1821061493, 248900140, 905130946,   ...
    1307421505, 1091346202, 1140141037, 1244153851, 611793617, 808278095,   ...
    2106046913, 1226683843, 90274853, 2147466840, 2140522509, 1146667727,  ...
    1530171233, 319146380, 2077765218, 789950931, 632682054, 1683083109,   ...
    351704670, 171230418, 1986612391, 1411714374, 872179784, 801917373,   ...
    144283230, 1949917822, 113201992, 1965781805, 679060319, 630187802,   ...
    1298843779, 1565131991, 50366922, 144513213, 207666443, 13354949,   ...
    631135695, 958408264, 494931978, 1150735880, 1640573652, 1315013967,   ...
    1254156333, 587564032, 1912367727, 2138169990, 2059764593, 115613893,  ...
    131630606, 1386707532, 2081362360, 1428465736, 1170668648, 931140599,  ...
    197473249, 1381732824, 925311926, 576725369, 198434005, 1296038143,   ...
    653782169, 1503907840, 33491376, 238345926, 1366925216, 1549695660,   ...
    1793656969, 1701980729, 1883864564, 251608257, 642486710, 1115145546,  ...
    1021484058  ] ;
y=seeds(stream);


%************************************************************************
function [Z_sorted, nU] = empVF_time_series2(n,Ds,Z);
% Time series
% D Dimensionen, in jeder eine empirische VF
% Z_sorted(d+D',:)=Z_sorted(d,:), F_emp(d+D',:)=F_emp(d,:)
%Empirical distribution:

step=[1/n:1/n:1]';  % [1/n, 2/n, ... ,1]

for d=1:Ds
    [Z_sorted(d,:) index_alt(d,:)] = sort(Z(d,:));
end

U=zeros(Ds,n);     % beschleunigt die Rechnung erheblich

% same samole values have different values F: i/n, (i+1)/n, ...
for i=1:n
    for d=1:Ds
        nU(d,index_alt(d,i)) =  i;
        % values n*U(d,index_alt), in order to omit rounding errors
    end
end


%******************************************************************
function [Z_sorted, nU, F_emp] = empVF_v3(n,D,Z);
%Empirical distribution:

step=[1/n:1/n:1]';  % [1/n, 2/n, ... ,1]

for d=1:D
    [Z_sorted(d,:) index_alt(d,:)] = sort(Z(d,:));
end
U=zeros(D,n);     % beschleunigt die Rechnung erheblich

%Find  empirical CDF  F_emp
F_emp = zeros(D,n);
for d=1:D
    F_emp(d,:) = step';
end
% Stichprobenwerte mehrfach have different values F i/n, (i+1)/n, ...
for i=1:n
    for d=1:D
        nU(d,index_alt(d,i)) =  i;
        % values n*U(d,index_alt) in order to omit rounding errors
    end
end