#####################################################################
#                                                                   #
#  This program calculates the reliability of a threshold system    #
#  by finding the distribution of:                                  #
#                                                                   #
#     S = (a_{1} X_{1} + ... + a_{n} X_{n}).                        #
#                                                                   #
#                                                                   #
#####################################################################

import numpy as np

# Component reliabilities (p[0] is not in use)
p = [0.0, 0.75, 0.80, 0.82, 0.65, 0.88, 0.91, 0.92, 0.86]

# Component weights (a[0] is not in use)
a = [0, 26, 39, 65, 78, 91, 104, 104, 143]

# Threshold
b = 410

# Number of components
n = len(p) - 1

# Cumulative weights (d[0] is not in use)
d = np.zeros(n+1, dtype=int)

for j in range(1, n + 1):
    d[j] = a[j] + d[j-1]

# After each iteration ps[s] = a_{1} X_{1} + ... + a_{n} X_{n} = s)
ps = np.zeros(d[n] + 1)

# Temporary storage of ps[s]
ph = np.zeros(d[n] + 1)

# Calculate ps[0] = P(X_{1} = 0) and ps[1] = P(X_{1} = 1):
ps[0] = 1 - p[1]
ps[a[1]] = p[1]

# Save these values for the next iteration
ph[0] = ps[0]
ph[a[1]] = ps[a[1]]


for j in range(1, n):
    # Calculate P(a_{1} X_{1} + ... + a_{j} X_{j+1} = s), s = 0, 1, ... , (a[j+1]-1)    
    for s in range(a[j+1]):
        ps[s] = ph[s] * (1 - p[j+1])
        
    # Calculate P(a_{1} X_{1} + ... + a_{j} X_{j+1} = s), s = a[j+1], ..., d[j]:
    for s in range(a[j+1], d[j] + 1):
        ps[s] = ph[s] * (1 - p[j+1]) + ph[s-a[j+1]] * p[j+1]
        
    # Calculate P(a_{1} X_{1} + ... + a_{j} X_{j+1} = s), s = d[j]+1 , ... , d[j+1]:
    for s in range(d[j] + 1, d[j+1] + 1):
        ps[s] = ph[s-a[j+1]] * p[j+1]
    
    # Save these values for the next iteration
    for s in range(d[j+1] + 1):
        ph[s] = ps[s]

# At this stage ps[s] = P(a_{1} X_{1} + ... + a_{j} X_{j+1} = s), s = 0, 1, ... , d[n]

for s in range(d[n] + 1):
    print("P(S = ", s, ") = ", ps[s])

# Calculate the system reliability h = P(S = b) + ... + P(S = d[n])  
h = 0

# for s in range(b, d[n] + 1):
#     h += ps[s]

for s in range(d[n] + 1):
    if s >= b:
        h += ps[s]
    
print("h = ", h)