##########################################################################
#
#    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
#
#    This is an extended version of bridge2.py which also includes
#    calculations of the bounds given in Corollary 6.2.6 (l1(p) and u1(p))
#    and Corollary 6.2.8 (l2(p) and u2(p))
#
#    Finally, we calculate the *best* bounds by taking the largest lower
#    bound and smallest upper bound for each p.
#
##########################################################################

import matplotlib.pyplot as plt
import numpy as np

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

def coprod3(p, q, r):
    return 1-(1-p)*(1-q)*(1-r)

def coprod_k(p, k):
    return 1-(1-p)**k


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

def l1(p):
    return p*p

def u1(p):
    return coprod2(p, p)

def l2(p):
    c2 = coprod2(p, p)
    c3 = coprod_k(p, 3)
    c4 = coprod_k(p, 4)
    return c2 * c2 * c3 * c3 * c4 * c4

def u2(p):
    c2 = p*p
    c3 = coprod_k(p**3, 3)
    c4 = coprod_k(p**4, 3)
    return coprod3(c2, c3, c4)


pp = np.zeros(101)
hh = np.zeros(101)
ll1 = np.zeros(101)
uu1 = np.zeros(101)
ll2 = np.zeros(101)
uu2 = np.zeros(101)
ll = np.zeros(101)
uu = np.zeros(101)


for j in range(101):
    pp[j] = j / 100
    hh[j] = h(pp[j])
    ll1[j] = l1(pp[j])
    uu1[j] = u1(pp[j])
    ll2[j] = l2(pp[j])
    uu2[j] = u2(pp[j])
    ll[j] = max(ll1[j], ll2[j])
    uu[j] = min(uu1[j], uu2[j])
    
    
plt.plot(pp, hh, label='h(p)')   
plt.xlabel('p')
plt.ylabel('h')
plt.title("Reliability curve")
plt.legend()
plt.show()
    
plt.plot(pp, hh, label='h(p)')   
plt.plot(pp, ll1, label='l1(p)')   
plt.plot(pp, uu1, label='u1(p)')   
plt.plot(pp, ll2, label='l2(p)')   
plt.plot(pp, uu2, label='u2(p)')   
plt.xlabel('p')
plt.ylabel('h')
plt.title("Reliability curve with bounds")
plt.legend()
plt.show()

    
plt.plot(pp, hh, label='h(p)')   
plt.plot(pp, ll, label='l*(p)')   
plt.plot(pp, uu, label='u*(p)')   
plt.xlabel('p')
plt.ylabel('h')
plt.title("Reliability curve with best bounds")
plt.legend()
plt.show()