###############################################################
#                                                             #
#    Simulate a bridge 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.

# Number of components
n = len(failureAlpha)

calculateBirnbaum = True 
calculateBarlowProschan = True 


FILE1 = "event01/availability.pdf"
FILE2 = "event01/birnbaum.pdf"
FILE3 = "event01/mean_birnbaum.pdf"
FILE4 = "event01/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 by pivotal decomposition and s-p-reduction
def structureFunction(X):
    return X[2] * (coproduct(X[0], X[1])) * (coproduct(X[3], X[4])) \
        + (1 - X[2]) * coproduct(X[0] * X[3], X[1] * X[4])
    
# The reliability function is obtained by pivotal decomposition and s-p-reduction
def reliabilityFunction(p):
    return p[2] * (coproduct(p[0], p[1])) * (coproduct(p[3], p[4])) \
        + (1 - p[2]) * coproduct(p[0] * p[3], p[1] * 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)