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

from math import *

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)


# The given point of time
time: float = 25.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 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):
    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
    counter = 0
    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)
                        counter += 1
    print("Number of terms calculated = ", counter)


def printReliabilities():
    rel = 0.0
    for k in range(6):
        phi = 5-k
        rel += reliability[phi]
        print("P(phi(", time, ") >= ", phi, ") = ", rel)
        

calcReliabilities(time)

print("");

printReliabilities()