Python loops (numpy)

Version's name: Python loops (numpy) ; a version of the Python loops program.
Repository: [home] and version downloads: [.zip] [.tar.gz] [.tar.bz2] [.tar]
Patterns and behaviours: Inefficient Python loops ·
Implemented best practices: Usage of Numba and Numpy to improve Python's serial efficiency ·

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.

The following experiments have been registered: