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.


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.

        '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)
        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
        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]) )
        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)
>> 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

>> 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.

  • \$\begingroup\$ I tried padding with 1s instead of 0s since log_2(1) = 0; this way, I do not need to check for nonzero elements in 'get_mod_D(j, k, style)', which has the downside of spending time adding lots of zeroes to each sum. Still too slow.. \$\endgroup\$ – MPath Jun 9 '17 at 8:23

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.

        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]) )
        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
    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')

>> [5.2841627612718041, 5.3472605443550068, 5.1868364187260321, 5.3586964954127474, 5.6271880002307695, 3.8399999999999999]

>> [1 2 3 4 5 6]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.