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

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

• 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..
– user127168
Jun 9, 2017 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.

EX:
ncol = 2 (( via time_sc_index ))
==> window_size = 2
==> check max( dsample, dsample ),
max( dsample, dsample ),
max( dsample, dsample ), ...
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]