##########################################################################
#
#    This program plots h(p) for a 2-terminal undirected network system
#    with 7 components, including two bridges (components 4 and 5)
#    The system is described in Exercise 6.6
#
##########################################################################

import matplotlib.pyplot as plt
import numpy as np

def coprod2(p, q):
    return 1-(1-p)*(1-q)


def h(p):
    c1 = coprod2(p, p)
    c2 = coprod2(p*p, p*p)
    c3 = coprod2(p, p*p)
    return p*(c1*c1*c1 + (1-p)*(1-p)*c2) + (1-p)*c3*c3

pp = np.zeros(101)
hh = np.zeros(101)

for j in range(101):
    pp[j] = j / 100
    hh[j] = h(pp[j])
    
plt.plot(pp, hh, label='h(p)')   
plt.xlabel('p')
plt.ylabel('h')
plt.title("Reliability curve")
plt.legend()
plt.show()