###############################################################
#                                                             #
#    Simulate a threshold system with repairable components   #
#                                                             #
###############################################################


from random import seed
from c_event import Weibull, EventQueue, EventComponent, EventSystem, coproduct

MAX_TIME = 2000
SAMPLE_INTERVAL = 10
NUM_SIMULATIONS = 5000

failureAlpha = [1.5, 1.5, 1.5, 1.5, 1.5]      # Shape params for failure dists.
failureBeta = [30.0, 60.0, 20.0, 40.0, 50.0]  # Scale params for failure dists.
repairAlpha = [1.0, 1.0, 1.0, 1.0, 1.0]       # Shape params for repair dists.
repairBeta = [15.0, 30.0, 10.0, 20.0, 25.0]   # Scale params for repair dists.

# Threshold system weights and threshold
w = [2.0, 4.0, 2.3, 3.5, 1.9]
threshold = 5.0

# Number of components
n = len(failureAlpha)

calculateBirnbaum = True 
calculateBarlowProschan = True 


FILE1 = "event02/availability.pdf"
FILE2 = "event02/birnbaum.pdf"
FILE3 = "event02/mean_birnbaum.pdf"
FILE4 = "event02/barlow_proschan.pdf"

SEED_NUM = 135246780
seed(SEED_NUM)

## Create the objects needed ##
eventQueue = EventQueue(MAX_TIME)

eventComps = []
mean_lifecycle = []
availability = []

# Calculate asymptotic component availabilities
for i in range(n):
    failureDist = Weibull(failureAlpha[i], failureBeta[i])
    repairDist = Weibull(repairAlpha[i], repairBeta[i])
    eventComps.append(EventComponent(eventQueue, failureDist, repairDist))
    mean_lifecycle.append(failureDist.mean() + repairDist.mean())
    availability.append(failureDist.mean() / mean_lifecycle[i])

# The structure function is obtained from the definition of a threshold system
def structureFunction(X):
    wsum = 0.0
    for i in range(n):
        wsum += w[i] * X[i]
    if wsum >= threshold:
        return 1
    else:
        return 0

# The reliability function is obtained by pivotal decomposition and s-p-reduction
def reliabilityFunction(p):
    return p[1]*p[3] + (p[1]*(1-p[3]) + (1-p[1])*p[3]) * coproduct(p[0], p[2], p[4]) + \
         (1-p[1])*(1-p[3]) * p[0] * p[2] * p[4]

# The Birnbaum measure is obtained by using the reliabilityFunction()
def birnbaumFunction(p, index):
    tmp = p[index]
    p[index] = 1
    h1 = reliabilityFunction(p)
    p[index] = 0
    h0 = reliabilityFunction(p)
    p[index] = tmp
    return h1 - h0

eventSystem = EventSystem(eventQueue, structureFunction, eventComps, SAMPLE_INTERVAL, calculateBirnbaum, calculateBarlowProschan)


## Run the simulation ##
for m in range(NUM_SIMULATIONS):
    if m % 100 == 0:
        print(".", sep='', end='')
    eventQueue.reinit()
    eventSystem.prepare()
    eventQueue.processAllEvents()

if NUM_SIMULATIONS > 1:
    eventSystem.calculateAverages()
    
eventSystem.plotComboSystemStates(FILE1)

if calculateBirnbaum:
    eventSystem.plotBirnbaumImportance(FILE2)
    eventSystem.plotMeanBirnbaumImportance(FILE3)

if calculateBarlowProschan:
    eventSystem.plotBarlowProschanImportance(FILE4)

print("")
eventSystem.printSystemStatistics()

for i in range(n):
    print("Stationary availability", (i+1), ":", availability[i])
print("------------------------------------")
print("Stationary availability:", reliabilityFunction(availability))
print("------------------------------------")

stationaryBirnbaum = []
stationaryBarlowProschan = []
sum_b_p = 0

for i in range(n):
    stationaryBirnbaum.append(birnbaumFunction(availability, i))
    stationaryBarlowProschan.append(stationaryBirnbaum[i] / mean_lifecycle[i])
    sum_b_p += stationaryBarlowProschan[i]
    
for i in range(n):
    print("Stationary Birnmaum importance", (i+1), ":", stationaryBirnbaum[i])

print("------------------------------------")
    
for i in range(n):
    print("Stationary BarlowProschan", (i+1), ":", stationaryBarlowProschan[i] / sum_b_p)