###################################################################
#
#    Python script calculating reliability as a function of time
#    of a subsea network for transmitting electronic pulses. 
#    See Section 9.2 of Huseby and Dahl (2021)
#
###################################################################

from math import *

import matplotlib.pyplot as plt
import numpy as np


def mcomb(n, k1, k2, k3) -> int:
    return factorial(n) / (factorial(k1) * factorial(k2) * factorial(k3) * factorial(n-k1-k2-k3))

def binomialDist(n, k, p) -> float:
    return comb(n, k) * p**k * (1-p)**(n-k)

def multinomialDist(n, k1, k2, k3, p1, p2, p3) -> float:
    return mcomb(n, k1, k2, k3) * p1**k1 * p2**k2 * p3**k3 * (1 - p1 - p2 - p3)**(n-k1-k2-k3)


minTime: float = 0.0
maxTime: float = 50.0

# r[0] :             Failure rate pr time unit for the wires.
# r[1], ... , r[8] : Failure rates of the communication links
r = [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01]


# p[0] :             Survival probability for the wires.
# p[1], ... , p[8] : Survival probabilities of the communication links
q = np.zeros(9)

reliability = np.zeros(6)

def initReliabilities():
    for i in range(6):
        reliability[i] = 0.0


def calcLinkProbabilities(t):
    for i in range(9):
        q[i] = exp(-r[i] * t)


def calcConnectionProbabilities():
    q_A: float = q[1] * q[6] * (q[4] + q[3] * q[5] - q[3] * q[4] * q[5]) * (q[7] + q[8] - q[7] * q[8])\
     + q[1] * (1 - q[6]) * (q[4] * q[7] + q[3] * q[5] * q[8] - q[3] * q[4] * q[5] * q[7] * q[8])
    q_B: float = q[2] * q[6] * (q[5] + q[3] * q[4] - q[3] * q[4] * q[5]) * (q[7] + q[8] - q[7] * q[8])\
     + q[2] * (1 - q[6]) * (q[5] * q[8] + q[3] * q[4] * q[7] - q[3] * q[4] * q[5] * q[7] * q[8])   
    q_AB_1_1: float = q[1] * q[2] * (q[4] + q[5] - q[4] * q[5]) * (q[7] + q[8] - q[7] * q[8])
    q_AB_1_0: float = q[1] * q[2] * (q[4] * q[7] + q[5] * q[8] - q[4] * q[5] * q[7] * q[8])
    q_AB_0_1: float = q[1] * q[2] * q[4] * q[5] * (q[7] + q[8] - q[7] * q[8])
    q_AB_0_0: float = q[1] * q[2] * q[4] * q[5] * q[7] * q[8]
    q_AB: float = q_AB_1_1 * q[3] * q[6] + q_AB_1_0 * q[3] * (1 - q[6])\
     + q_AB_0_1 * (1-q[3]) * q[6] + q_AB_0_0 * (1-q[3]) * (1 - q[6])
    return q_AB, q_A, q_B


def calcReliabilities(t):
    initReliabilities()
    calcLinkProbabilities(t)
    result = calcConnectionProbabilities()       # q_AB, q_A, q_B
    p_AB = result[0]                             # p_AB = q_AB
    p_A = result[1] - result[0]                  # p_A = q_A - q_AB
    p_B = result[2] - result[0]                  # p_B = q_B - q_AB
    for y1 in range(7):
        for y2 in range(7):
            for y3 in range(6):
                for y4 in range(6 - y3):
                    for y5 in range(6 - y3 - y4):
                        py1: float = binomialDist(6, y1, q[0])
                        py2: float = binomialDist(6, y2, q[0])
                        py345: float = multinomialDist(5, y3, y4, y5, p_AB, p_A, p_B)
                        w1: int = min(y1, y4)
                        w2: int = min(y2, y5)
                        w3: int = min((y1-w1) + (y2-w2), y3)
                        phi: int = w1 + w2 + w3
                        reliability[phi] += (py1 * py2 * py345)

        


time = np.zeros(101)
rel1 = np.zeros(101)
rel2 = np.zeros(101)
rel3 = np.zeros(101)
rel4 = np.zeros(101)
rel5 = np.zeros(101)

for s in range(101):
    t = minTime + (maxTime - minTime) * s / 100
    time[s] = t
    calcReliabilities(t)
    rel5[s] = reliability[5]
    rel4[s] = rel5[s] + reliability[4]
    rel3[s] = rel4[s] + reliability[3]
    rel2[s] = rel3[s] + reliability[2]
    rel1[s] = rel2[s] + reliability[1]
    
plt.plot(time, rel1, label='P(phi >= 1)')
plt.plot(time, rel2, label='P(phi >= 2)')
plt.plot(time, rel3, label='P(phi >= 3)')
plt.plot(time, rel4, label='P(phi >= 4)')
plt.plot(time, rel5, label='P(phi >= 5)')
plt.xlabel('time')
plt.ylabel('Reliability')
plt.title("Reliabilities as functions of time")
plt.legend()
plt.show()