###################################################################
#                                                                 #
#  This program estimates the Natvig importance measure for       #
#  components in a system with given minimal path sets using      # 
#  Monte Carlo simulations.                                       #
#                                                                 #
###################################################################

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]
beta = [50.0, 60.0, 70.0, 40.0, 50.0]

# Number of components and minimal paths
num_comps = 5
num_paths = 4


# Minimal path sets
paths = [{1,4}, {1,3,5}, {2,3,4}, {2,5}]

# Calculate the system lifetime given the component lifetimes
def sys_lifetime(tt):
    sys_life = 0.0
    for j in range(num_paths):
        path_life = np.inf
        for i in paths[j]:
            if tt[i-1] < path_life:
                path_life = tt[i-1]
        if path_life > sys_life:
            sys_life = path_life
    return sys_life

# The improvement in system lifetimes
# due to minimal repairs
z = np.zeros(num_comps)

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

# The simulated lifetimes of the components
# sampled from the conditional distribution
# given that they exceed the corresponding t
v = np.zeros(num_comps)

for k in range(num_sims):
    for i in range(num_comps):
        t[i] = weibullvariate(beta[i], alpha[i])
        v[i] = 0.0
        while v[i] < t[i]:
            v[i] = weibullvariate(beta[i], alpha[i])
    s = sys_lifetime(t)
    for i in range(num_comps):
        temp = t[i]
        t[i] = v[i]
        z[i] += (sys_lifetime(t) - s)
        t[i] = temp

# Calculate the unstandardized Natvig measure
z_sum = 0.0
z_imp = np.zeros(num_comps)
for i in range(num_comps):
    z_imp[i] = z[i] / num_sims
    z_sum += z_imp[i]

# Standardize the measure so that it adds up to 1
nat_imp = np.zeros(num_comps)
for i in range(num_comps):
    nat_imp[i] = z_imp[i] / z_sum

comp = np.linspace(1, num_comps, num_comps)

fig = plt.figure(figsize = (7, 4))
plt.bar(comp, nat_imp, color ='gray', width = 0.4)
plt.xlabel("Components")
plt.ylabel("Importance")
plt.title("Natvig importance")
plt.show()