###################################################################
#                                                                 #
#  This program plots h(p(t)) as a function of t for a network    #
#  system using Monte Carlo simulations of the component          #
#  lifetimes and calculating the resulting system lifetime.       #
#                                                                 #
###################################################################

from math import *
from random import *
import matplotlib.pyplot as plt
import numpy as np

# Number of simulations
num_sims = 100000

# Weibull-parameters for the components
alpha = [2.0, 2.5, 3.0, 1.0, 1.5, 2.0, 2.5]
beta = [50.0, 60.0, 70.0, 40.0, 50.0, 55.0, 65.0]

# Number of components
num_comps = 7

# Component state vector
x = np.zeros(num_comps)

# The coproduct function
def coprod(x1, x2):
    return x1 + x2 - x1 * x2

# The structure function given that component 3 is functioning
def phi_3_1(xx):
    return coprod(xx[2], xx[4]) * coprod(xx[0], xx[1]) * coprod(xx[5], xx[6]) \
        + (1 - coprod(xx[2], xx[4])) * coprod(xx[0] * xx[5], xx[1] * xx[6])

# The structure function given that component 3 is failed
def phi_3_0(xx):
    return coprod(xx[0] * xx[2], xx[1]) * coprod(xx[4] * xx[5], xx[6])

# The structure function of the network system (using pivotal decomposition)
def phi(xx):
    return xx[3] * phi_3_1(xx) + (1 - xx[3]) * phi_3_0(xx)

# Calculate the system lifetime given the component lifetimes
def sys_lifetime(tt):
    # Initialize component states
    for i in range(num_comps):
        x[i] = 1
    # Initialize system state and lifetime
    sys_state = phi(x)
    sys_life = 0.0
    # Determine the ordering of the failure times
    order = tt.argsort()
    # Initialize failure counter
    c = 0
    # Switch off one component at a time in the order of the failure times
    while (sys_state == 1) and (c < num_comps):
        x[order[c]] = 0
        sys_life = tt[order[c]]
        sys_state = phi(x)
        c += 1
    if sys_state == 0:
        return sys_life
    else: # This happens only for trivial systems where phi(0,...,0) = 1.
        return np.inf 

# The upper percentile levels
q = np.zeros(num_sims)

# The simulated system lifetimes
s = np.zeros(num_sims)

# The simulated lifetimes of the components
t = np.zeros(num_comps)

for k in range(num_sims):
    q[k] = 1.0 - k / num_sims
    for i in range(num_comps):
        t[i] = weibullvariate(beta[i], alpha[i])
    s[k] = sys_lifetime(t)

s.sort()

fig = plt.figure(figsize = (7, 4))
plt.plot(s, q, label='h(p(t))')
plt.xlabel('Time')
plt.ylabel('Reliability')
plt.title("System reliability as a function of time")
plt.legend()
plt.show()