Test-efficient-entropy-calc.rst

Test efficient calculation of information entropy with masked arrays

# Step 1: load previously calculated results

import pickle
import numpy as np

ua = pickle.load(open("gipps_results.pkl"))

# cast probabilities into numpy array
ua.p_block = np.array(ua.p_block)

def ie_standard(ua):
    # test standard way to calculate information entropy with nested blocks:
    # calculate information entropy and store in self.e_block
    e_block = np.zeros_like(ua.p_block[1])
    for p_block in ua.p_block:
        for i in range(ua.nx):
            for j in range(ua.ny):
                for k in range(ua.nz):
                    if p_block[i, j, k] > 0:
                        e_block[i, j, k] -= p_block[i, j, k] * np.log2(p_block[i, j, k])
    return e_block

%%timeit
e_block_1 = ie_standard(ua)

1 loops, best of 3: 32.8 s per loop

# now test implementation with masked arrays

def ie_masked(ua):
    pm = np.ma.masked_equal(ua.p_block, 0)
    pm2 = - pm * np.ma.log2(pm)
    ie_masked = np.sum(pm2.filled(0), axis = 0)
    return ie_masked

%%timeit
e_block_2 = ie_masked(ua)

1 loops, best of 3: 2.69 s per loop

# Check if result is correct
assert((e_block_1.all() == e_block_2.all()))

import matplotlib.pyplot as plt

%matplotlib inline

plt.imshow(e_block_1[100,:,:].transpose())

<matplotlib.image.AxesImage at 0x10c86cf10>

plt.imshow(e_block_2[100,:,:].transpose())

<matplotlib.image.AxesImage at 0x10c98a490>

np.max(ie_masked)

3.2090368361692407

75 * 75

5625

(1 * 0.05 - 1.0)**2 * 5625

5076.5625

print ua