<jupyter_start><jupyter_text>Workflow debug<jupyter_code>import time
import numba as nb
import numpy as np
year_hist = np.arange(2000, 2019).astype(int) # 2000 --> 2018
# importing GEFS historical fcst
GEFS_history = {}

for year in year_hist:
    # random number to replace real data for debugging
    # size = (365, 64, 64)
    apcp_ = np.random.normal(0.0, scale=1.0, size=(365, 64, 64))
    pwat_ = np.random.normal(0.0, scale=1.0, size=(365, 64, 64))
    
    GEFS_history['apcp_{}'.format(year)] = apcp_
    GEFS_history['pwat_{}'.format(year)] = pwat_
# importing new fcst
apcp_new = np.random.normal(0.0, scale=1.0, size=(365, 32, 32))
pwat_new = np.random.normal(0.0, scale=1.0, size=(365, 32, 32))<jupyter_output><empty_output><jupyter_text>Objective Given new forecasts with 32-by-32 grids, identify the **top 50 similar** 32-by-32 patterns within the historical forecast with 64-by-64 grids and 20 years through a **sliding window**.* One historical day can have only one window.**pseudo code*** loop over new forecast days* loop over old forecast days* sliding window on each day, found the "single best" window* Comparing the "single best" with the existing "top 50"* If sigle best is good, then insert into the "new top 50"<jupyter_code>N_fcst_day = 365 # new forecast has 365 days
EN = 50 # top 50 similar
size = 32 # size of new forecast
detect_size = 64 # size of historical forecast

# number of windows on both x and y directions
N_window = detect_size-size
@nb.jit(nopython=True)
def shift_one(arr, num=1, fill_value=999):
    result = np.empty_like(arr)
    result[:num] = fill_value
    result[num:] = arr[:-num]
    return result
<EMAIL>(nopython=True)
def window_slider(apcp_cadidate, pwat_cadidate, apcp_base, pwat_base, 
                  hist_day, N_window, size, 
                  record_single_day, indx_single_day, indy_single_day):
    '''
    '''
    apcp_base_sub = apcp_base[hist_day, ...]
    pwat_base_sub = pwat_base[hist_day, ...]
    # windows
    for i in range(0, N_window):
        for j in range(0, N_window):
            apcp_window = apcp_base_sub[i:i+size, j:j+size]
            pwat_window = pwat_base_sub[i:i+size, j:j+size]
            # the rule of similar (cannot be changed)
            record_temp = 0.7*np.sum(np.abs(apcp_cadidate-apcp_window))+0.3*np.sum(np.abs(pwat_cadidate-pwat_window))
            # identify single best
            if record_temp < record_single_day:
                # good match
                record_single_day = record_temp
                indx_single_day = i
                indy_single_day = j
                
    return record_single_day, indx_single_day, indy_single_day<jupyter_output><empty_output><jupyter_text>The part that is slow<jupyter_code># global allocations
record_full = np.ones([N_fcst_day, EN], dtype=np.float16)*9999.9
record_full_day = np.ones([N_fcst_day, EN], dtype=np.int)*999
record_full_year = np.ones([N_fcst_day, EN], dtype=np.int)*999
record_full_indx = np.zeros([N_fcst_day, EN], dtype=np.int)
record_full_indy = np.zeros([N_fcst_day, EN], dtype=np.int)

# local allocation
record_day = np.empty([EN], dtype=np.int)
record_year = np.empty([EN], dtype=np.int)
record_indx = np.empty([EN], dtype=np.int)
record_indy = np.empty([EN], dtype=np.int)
# -------------------- #

indx_single_day = 999
indy_single_day = 999

# loop over days of new forecasts
for fcst_day in range(N_fcst_day):
    start_time = time.time()
    apcp_cadidate = apcp_new[fcst_day, ...]
    pwat_cadidate = pwat_new[fcst_day, ...]
    # initializing the top 50 records
    record = np.ones([EN], dtype=np.float32)*9999
    # -------------------- #
    # loop over historical forecast by years
    for year in year_hist:
        # -------------------- #
        # current year information
        apcp_base = GEFS_history['apcp_{}'.format(year)]
        pwat_base = GEFS_history['pwat_{}'.format(year)]
        # -------------------- #
        days_history = np.arange(365) # <----- ignoring leap year for testing
        
        # loop over historical forecast days
        for hist_day in days_history:
            record_single_day = 9999.9
            # ---------- Sliding window (JIT) ---------- #
            # update single day record
            # update indx
            # update indy
            record_single_day, indx_single_day, indy_single_day = window_slider(
                apcp_cadidate, 
                pwat_cadidate, 
                apcp_base, 
                pwat_base, 
                hist_day, 
                N_window, 
                size, 
                record_single_day, 
                indx_single_day, 
                indy_single_day)
            # ------------------------------------------ #
            # compare single best with global best        
            if record_single_day < record[-1]:
                # if the current day can be added to the "top 50"
                analog_ind = np.searchsorted(record, record_single_day)
                # shift one to the left and free space
                record[analog_ind:] = shift_one(record[analog_ind:])
                record_day[analog_ind:] = shift_one(record_day[analog_ind:])
                record_year[analog_ind:] = shift_one(record_year[analog_ind:])
                record_indx[analog_ind:] = shift_one(record_indx[analog_ind:])
                record_indy[analog_ind:] = shift_one(record_indy[analog_ind:])
                # insert the current day
                record[analog_ind] = record_single_day
                record_day[analog_ind] = hist_day
                record_year[analog_ind] = year
                record_indx[analog_ind] = indx_single_day
                record_indy[analog_ind] = indy_single_day
                
    record_full[fcst_day, :] = record
    record_full_day[fcst_day, :] = record_day
    record_full_year[fcst_day, :] = record_year
    record_full_indx[fcst_day, :] = record_indx
    record_full_indy[fcst_day, :] = record_indy

    print("--- %s seconds ---" % (time.time() - start_time))<jupyter_output>--- 36.706010580062866 seconds ---
--- 36.18890118598938 seconds ---
--- 36.18970227241516 seconds ---
--- 36.08813261985779 seconds ---
{'fraction_non_alphanumeric': 0.12710765239948119, 'fraction_numerical': 0.03469520103761349, 'mean_word_length': 2.446368715083799, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '56096', 'license': 'Unlicense', 'n_tokens_mistral': 2247, 'n_tokens_neox': 2007, 'n_words': 459, 'path': 'Python/Workflow_debug.ipynb', 'repo_name': 'yingkaisha/notebook_code_vault', 'size': '6169'}