@ This is a Gauss code for the Computation of a Single PE point estimate of a univariate time series and associated Relative Frequencies @
@ This code accompanies the Aptech blog Permutation entropy which was published on 12/14/2018 @
@ See: https://www.aptech.com/blog/permutation-entropy/ @
@ This Gauss program requires a working copy of GAUSS 18+ @

@ Suggested Citation:
Clower, E. and M. Henry, 2018. "petropy: Gauss module to compute a Single Permutation Entropy point estimate,"
@

NEW;
CLS;
 
/*
** Format:  H = petropy(x, D, tau, method, accu)
**
** Input:
**  x       Column vector, univariate time series data
**
**  D       Scalar, embedding dimension
**
**  tau     Scalar, embedding time delay that determines the time separation between time series values
**
**  method  String treatment of equal values
**          'noise' - add small noise
**          'equal' - allow same rank for equal values
**          'order' - consider order of appearance (first occurence --> lower rank)
**         
**  accu    Scalar, parameter describing accuracy of the values in X
**          by number of decimal places
**
**  Outputs
**  H      Permutation Entropy
**  Hnorm  Normalized Permutation Entropy
**  refreq Relative Frequencies with its plot 
*/

// one-dimensional time series toy example
x = {4, 7, 9, 10, 6, 11, 3};       
x = packr(x);

//embedding dimension
D = 3;

//embedding time delay
tau = 1;

// PE
print;
H = pentropy(x, D, tau, "noise");     
print "Permutation Entropy point estimate =" H[1];     
print;

// RELATIVE FREQUENCIES
rem = {1};                     
RelFreq = delRow(H, rem);
print "Relative Frequencies" RelFreq;
numrows = rows(RelFreq);
xx = seqa(1, 1, numrows);
plotXY(xx, RelFreq);

#include dynargs.dec
proc (1) = pentropy(x, n, tau, ...);
    local n_dynargs, accu, method, M, equal, tmp, shift_map,
    k, j, shift_mat_ind, tmp2, tmp3, shift_mat, ind_mat, ind_vec, sort_ind_mat, ia, ic,
    tmp_sort, permpat_num, refreq, size,tmp0, denom,Hnorm,factor;
    
    n_dynargs = COUNT_DYNARGS;
    
    //Set accuracy of the values in X
    if n_dynargs < 2;
        accu = 4;           // default
    else;
        accu = sysstate(GET_ONE_DYNARG, 2);
    endif;
    
    //Set method
    if n_dynargs < 1;
        method = "order";   // default
    else;
        method = sysstate(GET_ONE_DYNARG, 1);
    endif;
    
    //Get length of x
    M = rows(x);
    
    equal = 0;
    
    if n*log(n)>15;
            errorlogat "Permutation dimension too high";
        end;
    endif;
    
    if (rows(tau)) > 1 and (rows(tau) != n-1);
        errorlogat "Time lag vector has to have n-1 entries";
        end;
    endif;
        
        if ((n-1)*minc(tau) >=M) or maxc(tau) >= M;
        errorlogat  "Too few data points for desired dimension and lags";
    endif;
    
    if lower(method) == "noise";
        print "Method: Add small noise.";
        x = x + rndn(M,1)*10^(-accu-1);
    elseif lower(method) == "equal";
         print "Method: Allow equal ranks.";
        equal = 1;
    elseif lower(method) == "order";
         print "Method: Consider order of occurrence.";
    else;
        errorlogat "Unknown method. Default method 'order' used";
        
    endif;
    print;
    
    if rows(tau) > 1;
        tau = reshape(tau, rows(tau), 1);
        tmp = 0|tau;
        tau = sortc(tmp);
        
        // build n x (M-tau(n))-matrix from shifted values in x
        shift_mat = zeros(n,  M-tau[n]);
        
        for ii(1, n, 1);
            k = tau[ii] + 1;
            j = M - tau[n] + tau[ii];
            shift_mat[ii, .] = x[k:j];
        endfor;
    else;
        
        shift_mat_ind = zeros((M-(n-1)*tau)*n, 1);
        
        tmp0 = seqa(1, tau, n);   // Column vector
        shift_mat_ind[1:n] = tmp0;

        for i(1, (M-(n-1)*tau)-1, 1);   
            shift_mat_ind[(i*n+1):((i+1)*n)] = tmp0 + i;
        endfor;
        
        shift_mat = reshape(x[shift_mat_ind], (M-(n-1)*tau), n)';   
        
    endif;
    
    if equal;
        ind_mat = zeros(size(shift_mat));
        for ii(1, cols(ind_mat), 1);
            tmp = unique(shift_mat[., ii]);
            ind_mat[.,ii] = indnv(shift_mat[., ii], tmp);
        endfor;
    else;
        
        // INDEX ORDERS of smallest to largest
        sort_ind_mat = sortmatind(shift_mat);                       
        
        ind_mat = zeros(rows(sort_ind_mat), cols(sort_ind_mat));    
        
        // Filling the zero matrix with RANK ORDERS
        for ii(1, cols(ind_mat), 1);    
            ind_mat[sort_ind_mat[., ii], ii] = seqa(1, 1, n);       
        endfor;
    endif; 
    
    // assign unique number to each pattern (base-n number system)
    ind_vec = n.^(seqa(0, 1, n))' * (ind_mat-1);    
 
    tmp_sort = sortc(ind_vec',1);   
    tmp2 = unique(tmp_sort);        
    ia = indnv(tmp2, tmp_sort);     
    
    tmp = ia|(cols(ind_vec)+1);     
    
    // Number of admissible Permutation Patterns
    permpat_num = diff(tmp);                           
    
    // RELATIVE FREQUENCIES of PERMUTATIONS p[i]
    refreq = permpat_num/sumc(permpat_num); 
    
    // PE of order D
    H = -sumc(refreq .* log2(refreq));
    
    factor = n!;
    Hnorm = (1/log2(factor)) * H;
    
    retp(Hnorm|refreq);
        
endp;

proc (1) = log2(x);
    retp( log(x)/log(2));
endp;

proc (1) = diff(x);
    local tmp;
    tmp = x[2:rows(x)] - x[1:rows(x)-1];
    retp(tmp);
endp;

proc (1) = sortmatind(x_mat);
    local sort_ind_mat;
    
    sort_ind_mat = zeros(rows(x_mat), cols(x_mat));
    
    for i(1, cols(x_mat), 1);
        sort_ind_mat[.,i] = sortind(x_mat[.,i]);
    endfor;
    retp(sort_ind_mat);
endp;

proc (1) = delRow(x, remove);
   local mask, xout;
   mask = zeros(rows(x), 1);
   mask[remove] = ones(rows(remove), 1);
   xout = delif(x, mask);
   retp(xout);
endp;

proc (1) = getDateRange(X, range);
    local ndigits, idx, X_range;
    ndigits = floor(log(range)) + 1;
    
    range = range .* 10^(14-ndigits);
    
    idx = sumc((x[.,1] .< range[1]) ~ (x[.,1] .< range[2])) + 1;
    
    X_range = X[idx[1]:idx[2],.];
    
    retp(X_range);
endp;