###################################################################
#                                                                 #
#  This program plots h(p(t)) as a function of t for a threshold  #
#  system using a ROBDD representation of the system.             #
#                                                                 #
#  Contrary to script_B_3_3 which uses Monte Carlo simulation,    #
#  this script calculates reliability analytically. Thus, the     #
#  result is an exact h(p(t))-curve.                              #
#                                                                 #
###################################################################

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 = 175
steps = 175

# 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 Reliability function
time = np.linspace(0.0, max_t, steps)
h = np.zeros(steps)

# Calculate the Reliability function
for step in range(steps):
    for i in range(num_comps):
        p[i] = weibull(alpha[i], beta[i], time[step])
    h[step] = reliability(p)

# Plot the results
fig = plt.figure(figsize = (7, 4))
plt.plot(time, h, label='h(p(t))')
plt.xlabel('Time')
plt.ylabel('Reliability')
plt.title("System reliability as a function of time")
plt.legend()
if save_plots:
    plt.savefig("dynsys_03b/reliability_03b.pdf")
plt.show()