Speeding up maximum self-similarity test for heavy tail-exponents

Question

I am trying to reproduce results from a research paper using python. I've checked my method and it works on relatively small sample datasets. However, the code does not run for my actual dataset, which is initially ~11,000 data points and becomes ~70,000 data points after padding zeroes. The algorithm is designed to do a maximum self-similarity test to estimate heavy-tail exponents. The main paper I am trying to use is this , though this secondary source is nearly identical.

Code not shown (but available upon request):

The task is to read the data into lists (list1 -> times in years, list2 -> times in months, ..., list8 -> observed speeds, ... ). Then, the data lists are sliced in the time region of interest. Then, all times are converted into a single list hours, which are then shifted to start from zero. Since there are observations at non-consecutive hours, the maximum speed at each hour is taken, as shown in EX below.

EX:

times = [0, 1, 4, 4, 6, ...] --> [0, 1, 2, 3, 4, 5, 6, ...] # hours
speeds = [280, 370, 290.6, 339.8,  410.1, ...] --> [280, 370, 0, 0, 339.8, 0, 410.1, ...]
# My code does this quickly

Code I would like reviewed:

import random
import numpy as np

def get_data(schema):
    """
    This function returns the speeds and times that correspond
    to the desired data schema.

    schema:
        'random'    --   random Gaussian data sample
        'sample'    --   consecutive integer sample
          'mask'    --   consecutive integer sample
                         with zeroes at some indices
           'CME'    --   solar CME data
    """
    if schema == 'random': 
        mu, sigma = 48, 7
        speeds = [random.gauss(mu, sigma) for index in range(100)]
        # speeds = sorted(speeds)
        times = [idx for idx in range(len(speeds))]
    elif schema == 'sample': 
        speeds = np.linspace(1, 100, 100)
        times = [idx for idx in range(len(speeds))]
    elif schema == 'mask':
        speeds = np.linspace(1, 100, 100)
        times = [idx for idx in range(len(speeds))]
        for idx in range(len(speeds)):
            if idx % 3 == 0 or idx % 2:
                speeds[idx] = 0
    elif schema == 'CME': # code not shown, actual dataset
        speeds = padded_speeds
        times = get_times_for_speeds(speeds)
    else:
        raise ValueError("schema = 'random', 'sample', 'mask', 'CME'")
    return speeds, times

As described in section 2.1 of page 5 of the main paper, there is a formula to calculate X_m (where m = 2**j); I tried to do this using numpy arrays instead of for-loops. The pattern becomes most obvious when choosing the consecutive integer data sample.

speeds, times = get_data('sample') # check that algorithm works on small datasample

def shapeshifter(j, my_array=speeds):
    """
    This function reshapes an array to have
    j columns (j = ncol), and truncates the
    sub-arrays such that the size of the
    array can change.
    """
    my_array = np.array(my_array)
    desired_size_factor = np.prod([n for n in j if n != -1])
    if -1 in j:  # implicit array size
        desired_size = my_array.size // desired_size_factor * desired_size_factor
    else:
        desired_size = desired_size_factor
    return my_array.flat[:desired_size].reshape(j)

def get_D(j, k, style='speed'):
    """
    .- formula on bottom of page 2 of main source
    .- formula 2.5 on page 6 of secondary source
    """
    my_array = shapeshifter((-1, 2**j))
    print("\n shapeshifted data: \n", my_array, "\n") # helps to see the row containing the maximum
    rows = [my_array[index] for index in range(len(my_array))]
    res = []
    if style == 'index':
        for index in range(len(rows)):
            res.append( np.argmax(rows[index]) + index*j)
    elif style == 'speed':
        for index in range(len(rows)):
            res.append( max(rows[index]) )
    else:
        raise ValueError("style = 'index', 'speed'")
    return res[k-1]

To test this step, it helps to see how changing j and k values affects the pattern. Since 2**j is always even (since j is a positive integer), only the evens in the last column of the consecutive integer sample should be the maximums.

j = 1
k = 1
res = get_D(j, k)
print(res)
>> 2

The following may help see the pattern.

j, k = 1, 2 ==> res = 4 # max of [3,4] in [[1,2], [3,4], [5,6], ...]
j, k = 1, 3 ==> res = 6 # max of [5,6] in [[1,2], [3,4], [5,6], ...]
j, k = 2, 1 ==> res = 4 # max of [1,2,3,4] in [[1,2,3,4], [5,6,7,8], ...]
j, k = 3, 1 ==> res = 8 # max of [1,...,8] in [[1,...,8], [9,...,16], ...]
j, k = 2, 2 ==> res = 8 
j, k = 2, 3 ==> res = 12
j, k = 3, 2 ==> res = 16

To get Y(j) described immediately after in the two papers, a new function is defined. The size of each row and the number of rows become relevant, so window parameters are calculated as well.

def window_params(dataset=speeds):
    """
    This function gets the size of each row and number of columns for each j alongside corresponding parameters. 
    """
    ## number of data points (N = 11764)
    numdata = len(dataset)
    ## last term (corresponds to j=13)
    lim = int(np.floor(np.log2(numdata)))
    ## number of rows (corresponds to number of reshapes; j=[1,2,3,...,floor(log_2(N))=13])
    time_sc_index = np.linspace(1, lim, num=lim)
    ## size of window for j^th row of reshaped array (scale=[2,4,8,...,8192])
    window_size = [2**j for j in time_sc_index]
    ## b_j = scale (corrsponds to number of columns and weighted block size)
    block_size = np.floor([numdata/sc for sc in window_size])
    ## return window parameters
    return numdata, time_sc_index, window_size, block_size, lim

numdata, time_sc_index, window_size, block_size, last_j = window_params()


def get_mod_D(j, k, style='speed'):
    """
    This function modifies D(j,k) into a format that is convenient
    as an input to get Y(j). Since the log values are taken, all 
    padded zeroes from the actual dataset must be removed.
    """
    res = np.array(get_D(j, k, style))
    res = np.log2(res[res != 0])
    return res

def get_Y(j, style='speed'):
    """
    This function should be used on datasets that require concatenation.
    """
    nj = int(numdata/2**j)
    ks = np.linspace(1, nj, nj)
    ks = [int(val) for val in ks]
    # res = [get_D(j, k) for k in ks]
    res = np.concatenate([get_mod_D(j, k) for k in ks])
    # res = np.sum([get_mod_D(j, k) for k in ks])
    # print(res)
    return sum(res)/nj

print(get_Y(1))
>> 5.28416276127

Are there any faster ways I could achieve the same result? My goal is to calculate Y(j) for j in a defined list js such that I can reproduce the plot on the right of figure 1 in the main source.

Skippable Aside: My original attempt was to calculate the entire array of D(j,k) arrays using the metalooper(), defined in this post. I couldn't figure out how to compute the sums, but I think it can be adapted for this problem.

Answer 1

I got the code to work for my CME dataset by changing a few things. I still use shapeshifter() but do not use get_D() and get_mod_D(), and get_Y() is altered. New functions are defined to get the array of all D(j, k) for all j and k values. I also pad values with ones instead of zeroes so that values do not need to be checked before taking their logarithmic values.

def looper(ncol, my_array=speeds, style='speed'):
    """
    This function calls the function 'shapeshifter'
    and returns a list of the maximum values of each
    row in 'my_array' for 'ncol' columns. The length
    of each row denotes the size of each window.

    EX:
        ncol = 2 (( via time_sc_index ))
        ==> window_size = 2
        ==> check max( dsample[1], dsample[2] ),
                  max( dsample[3], dsample[4] ),
                  max( dsample[5], dsample[6] ), ...
            for k rows,
                where k = len(my_array)//ncol

    This function will return an error if the reshape
    of the array changes the size of the array.
    """
    my_array = shapeshifter((-1, ncol))
    rows = [my_array[index] for index in range(len(my_array))]
    res = []
    if style == 'index':
        for index in range(len(rows)):
            res.append( np.argmax(rows[index]) + index*ncol)
        # res.append( np.argmax(rows[index], axis=index))
    elif style == 'speed':
        for index in range(len(rows)):
            res.append( max(rows[index]) )
    else:
        raise ValueError("style = 'index', 'speed'")
    return res

def metalooper(window_size, my_array=speeds, style='speed'):
    """
    This function calls 'looper' for each ncol in
    window_size (( block_size ~ b_j )) and returns
    a list of sublists of the maximums in each
    window.
    """
    outer = [looper(win, my_array, style) for win in window_size]
    return outer

max_speeds = metalooper(window_size, speeds, 'speed') ## max speeds

Now to get the Y(j).

def get_Y(j, my_array, style='speed'):
    """
    This function returns the max spectrum for a single j.
    """
    return sum(np.log2(my_array[j-1])) / (numdata/(2**(j)))

def get_all_Y(my_array, style):
    """
    This function returns the full max spectrum for all j's.
    """
    # size = len(my_array)
    # js = np.linspace(1, size, size)
    # js = [int(val) for val in js]
    # return [get_Y(j, my_array, style) for j in js]
    # return [get_Y(j, my_array, style) for j in range(1, len(my_array)+1)]
    ys = [get_Y(j, my_array, style) for j in range(1, len(my_array)+1)]
    js = np.array([int(j+1) for j in range(len(ys))])
    return ys, js

ys, js = get_all_Y(max_speeds, 'speed')

print(ys)
>> [5.2841627612718041, 5.3472605443550068, 5.1868364187260321, 5.3586964954127474, 5.6271880002307695, 3.8399999999999999]

print(js)
>> [1 2 3 4 5 6]