###################################################################
#                                                                 #
#  This program plots I_B^{(i)}(t) as a function of t for a       #
#  threshold 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

# Time parameters
max_t = 120
steps = 120

# Number of simulations
num_sims = 100000

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

# Threshold system weights and threshold
w = [2.0, 4.0, 2.3, 3.5, 1.9]
threshold = 5.0

# Number of components and minimal paths
num_comps = 5

# Birnbaum measures
IB = np.zeros((num_comps, steps))

# A component is identified as critical when t is in [L, U)
L = np.zeros(num_comps)
U = np.zeros(num_comps)

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

# The structure function of the threshold system
def phi(xx):
    wsum = 0.0
    for i in range(num_comps):
        wsum += w[i] * xx[i]
    if wsum >= threshold:
        return 1
    else:
        return 0

# 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 simulated lifetimes of the components
t = np.zeros(num_comps)

# Time steps
time = np.linspace(0.0, max_t, steps)

# Run simulations
for k in range(num_sims):
    for i in range(num_comps):
        t[i] = weibullvariate(beta[i], alpha[i])
    for i in range(num_comps):
        temp = t[i]
        t[i] = np.inf
        U[i] = sys_lifetime(t)
        t[i] = 0
        L[i] = sys_lifetime(t)
        t[i] = temp
    for step in range(steps):
        for i in range(num_comps):
            if L[i] <= time[step] and time[step] < U[i]:
                IB[i][step] += 1

# Calculate average importance values
for step in range(steps):
    for i in range(num_comps):
        IB[i][step] = IB[i][step] / num_sims

fig = plt.figure(figsize = (7, 4))
for i in range(num_comps):
    plt.plot(time, IB[i], label='IB(' + str(i+1) + ", t)")
plt.xlabel('Time')
plt.ylabel('Importance')
plt.title("Importance as a function of time")
plt.legend()
plt.show()