from scipy.integrate import quad as numInt
import math
import numpy as np
import matplotlib.pyplot as plt
import statLib as sl  # statLib contiains statistical functions 
# Raw data
data = [25, 35, 28, 41, 22, 31, 28, 25]
# Model parameters:
DD = 30 # Due data = 30 days from now on
c_O = 1000
c_PD = 10000
mu_0 = sl.average(data)
sigma = sl.sd(data)

print("Estimate of mu_0 =",mu_0)
print("Estimate of sigma =",sigma)

def C_O(p,c_O): # Cost of overtime
    return p*c_O

def B(r, mu , sigma): # Backorder function, used in later problems
    return sigma * sl.phi((r - mu) / sigma) + (mu - r) * (1 - sl.Phi((r - mu) / sigma))

def B_integrand(x,DD,mu,sigma):  # Integrand in the backorder function, numerical integration
    return (x - DD)*sl.pdfNormal(x, mu, sigma)

def B_numInt(DD,mu,sigma): # Numerical integration to solve the backorder function
    high = mu + sigma*6
    I, err = numInt(B_integrand, DD, high, args=(DD,mu,sigma))
    return I
    
def C_PD(p,mu_0,sigma,c_PD):  # Cost of delay, i.e., penalty of default
    return c_PD*B_numInt(DD,mu_0/(1 + p/100),sigma)

print("E(# of days delayed), num.int =",B_numInt(DD,mu_0,sigma))
print("E(# of days delayed), backorder funcition =",B(DD,mu_0,sigma))

# Create plot
xpoints = np.zeros(21)
ypoints_tot = np.zeros(21)
ypoints_O = np.zeros(21)
ypoints_PD = np.zeros(21)

p = 2.5
for i in range(21):     
    xpoints[i] = p
    ypoints_O[i] = C_O(p,c_O)
    ypoints_PD[i] =  C_PD(p,mu_0,sigma,c_PD)
    ypoints_tot[i] = C_O(p,c_O) + C_PD(p,mu_0,sigma,c_PD)
    p += .25
    
onlyCost = True
if onlyCost == False:
    plt.plot(xpoints,ypoints_O,"-.b",label ="Overtime")
    plt.plot(xpoints,ypoints_PD,"-..r",label ="Penalty")
plt.plot(xpoints,ypoints_tot,"-g",label ="Total")
plt.legend(loc="lower right")
plt.show()