import numpy as np

h = np.array([0.1, 0.2, 0.3, 0.0, 0.01, 0.01, 0.1, 0.28])

w0 = [0.3, 0.0, 0.01, 0.01, 0.1, 0.28, 0.1, 0.2]
w1 = [0.01, 0.01, 0.1, 0.28, 0.1, 0.2, 0.3, 0.0]
w2 = [0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125]

w = np.array([w0, w1, w2])

eps = 0.2

h0 = h
hw0 = h0 * w[0]
sum_hw0 = np.sum(hw0)

h1 = (h0 + eps*hw0/sum_hw0) /(1+eps)
hw1 = h1 * w[1]
sum_hw1 = np.sum(hw1)

h2 = (h1 + eps*hw1/sum_hw1) /(1+eps)
hw2 = h2 * w[2]
sum_hw2 = np.sum(hw2)

h3 = (h2 + eps*hw2/sum_hw2) /(1+eps)

hp0 = h0
hpw0 = hp0 * w[0]
sum_hpw0 = np.sum(hpw0)
norm_hp0 = np.sum(hp0)

hp1 = (sum_hpw0 * hp0 + eps * norm_hp0 * hpw0)
hpw1 = hp1 * w[1]
sum_hpw1 = np.sum(hpw1)
norm_hp1 = np.sum(hp1)

hp2 = (sum_hpw1*hp1 + eps * norm_hp1 * hpw1)
hpw2 = hp2 * w[2]
sum_hpw2 = np.sum(hpw2)
norm_hp2 = np.sum(hp2)

hp3 = (sum_hpw2*hp2 + eps * norm_hp2 * hpw2)


# Show that hp are just multiple of h
print(hp1/h1)
print(hp2/h2)
print(hp3/h3)

# This should be the output of HW model
print(hp0, hp1, hp2, hp3)