##########################################################################
#                                                                        #
#    This program plots I_B^{(i)}(t) as a function of t for a simple     #
#    2-out-of-3 system with Weibull-distributed component lifetimes.     #
#                                                                        #
##########################################################################

from math import exp
import matplotlib.pyplot as plt
import numpy as np

save_plots = True

# Time parameters
max_t = 120
steps = 120

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

# Number of components in the system
num_comps = len(alpha)

# The reliability function of the system
def reliability(pp):
    return pp[0] * pp[1] + pp[0] * pp[2] + pp[1] * pp[2] - 2 * pp[0] * pp[1] * pp[2]

# The survival function of a Weibull-distribution
def weibull(aa, bb, t):
    if t >= 0:
        return exp(- (t / bb) ** aa)
    else:
        return 1.0

# Component reliability vector
p = np.zeros(num_comps)

# Time values and Birnbaum measures
time = np.linspace(0.0, max_t, steps)
IB = np.zeros((num_comps, steps))

# Calculate the Birnbaum measures
for step in range(steps):
    for i in range(num_comps):
        p[i] = weibull(alpha[i], beta[i], time[step])
    for i in range(num_comps):
        temp = p[i]
        p[i] = 1
        h1 = reliability(p)
        p[i] = 0
        h0 = reliability(p)
        p[i] = temp
        IB[i][step] = h1 - h0
 
# Plot the results
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()
if save_plots:
    plt.savefig("birnbaum_01/birnbaum_01.pdf")
plt.show()