The numpy version of the Python loops kernel uses Numpy’s functions and vectorization feature. This way, we minimize the usage of Python’s interpreter, relying instead on the low-level implementaion of Numpy library.
The following code snippet shows the change introduced to the code:
def compute(prob, vents, load_orig, w, float_th, i_size, exposure_time):
i_sel = vents[:, :]
i_sel = np.transpose(i_sel)
load_iarea = np.transpose(load_orig[:exposure_time], axes=(0,2,1)) # Transposes last 2 dimensions
load_iarea = load_iarea.reshape(load_iarea.shape[0], load_iarea.shape[1] * load_iarea.shape[2]) # Flattens last 2 dimensions
load_iarea = load_iarea[:, i_sel] * 1000.0
load_iarea = load_iarea[:, :, :, np.newaxis] # Adapts matrix's dimensions for broadcasting
persist_th = np.sum(np.greater(load_iarea, float_th), axis=0)
prob[i_size, :, :, :] += (persist_th / exposure_time) * w
return prob
def compute_step_1(sizes, weight_list, num_scen_found, load_orig, prob, vents, float_th, exposure_time):
for size in sizes:
if size == 'E':
continue
i_size = sizes.index(size)
num_scen = num_scen_found[i_size-1]
weights = weight_list[i_size-1]
for i_scen in range(num_scen):
w = weights[i_scen]
if w == 0:
# weight of this scenario is 0, pass
continue
prob = compute(prob, vents, load_orig[i_scen], w, float_th, i_size, exposure_time)
return prob
def compute_step_2(alpha, beta, prob, h):
alpha[:, :, 1:, 0] += np.transpose(h * prob[1:, :, :, 0], (1, 2, 0)) # Initializes thresholds 0
beta[:, :, 1:, 0] += np.transpose(h * (1 - prob[1:, :, :, 0]), (1, 2, 0)) # Initializes thresholds 0
tai = prob[1:, :, :, 1:]
tai1 = prob[1:, :, :, :-1]
condition = tai1 > 0 # Evaluates condition
i_ctrue = np.where(condition == True) # Retrieves indices from prob_model where the condition is True
i_cfalse = np.where(condition == False) # Retrieves indices from prob_model where the condiion is False
f = tai[i_ctrue[0], i_ctrue[1], i_ctrue[2], i_ctrue[3]] / tai1[i_ctrue[0], i_ctrue[1], i_ctrue[2], i_ctrue[3]] # Computes fact only for those elements where the condition is true
# Computes alpha and beta for those elements where the condition was true
alpha[i_ctrue[1], i_ctrue[2], i_ctrue[0] + 1, i_ctrue[3] + 1] += h * f
beta[i_ctrue[1], i_ctrue[2], i_ctrue[0] + 1, i_ctrue[3] + 1] += h * (1 - f)
# Computes alpha and beta for those elements where the condition was false
alpha[i_cfalse[1], i_cfalse[2], i_cfalse[0] + 1, i_cfalse[3] + 1] += 0
beta[i_cfalse[1], i_cfalse[2], i_cfalse[0] + 1, i_cfalse[3] + 1] += h
As you can notice, we have replaced the 3 innermost levels of for-loops in compute_step_1
by Numpy operations. In compute_step_2
we have removed all loop levels, and reimplemented the if-else
clause with a mask.