#####################################################################
#                                                                   #
#  This program calculates the reliability of a k-out-of-n system   #
#  by finding the distribution of:                                  #
#                                                                   #
#     S = X_{1} + ... + 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]

# Threshold
b = 6

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


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

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

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

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


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


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

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

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

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