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. 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
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
j in a defined list
js such that I can reproduce the plot on the right of figure 1 in the main source.
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.