###################################################################
#                                                                 #
#  This program plots I_B^{(i)}(t), i = 1, ..., 5 as a function   #
#  of t for a threshold system using a ROBDD representation of    #
#  the system.                                                    #
#                                                                 #
#  Contrary to script_B_3_7 which uses Monte Carlo simulation,    #
#  this script calculates importance analytically. Thus, the      #
#  results are exact I_B^{(i)}(t)-curves.                         #
#                                                                 #
###################################################################

from math import exp
import matplotlib.pyplot as plt
import numpy as np
from c_robdd import ROBDDSystem
from c_threshold import ROBDDThreshold

save_plots = True

# Time parameters
max_t = 120
steps = 120

# 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]

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

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

# Calculate the weight sum
asum = 0
for i in range(num_comps):
    asum += a[i]

# The ROBDD representation of the threshold system
sys = ROBDDSystem(ROBDDThreshold(a, asum, b))

# The reliability function of the system
def reliability(pp):
    result = sys.calculateReliability(pp)
    return result[1]

# 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_03b/birnbaum_03b.pdf")
plt.show()