from scipy.integrate import quad as numInt
from scipy.optimize import minimize
import math
import numpy as np
import matplotlib.pyplot as plt
import statLib as sl
# 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 C_PD(p,mu_0,sigma,c_PD):  # Cost of delay, i.e., penalty of default
    return c_PD*B(DD,mu_0/(1 + p/100),sigma)

def C(p,mu_0,sigma,c_O,c_PD): # Objective function
    return C_O(p,c_O) + C_PD(p,mu_0,sigma,c_PD)

# Find optimal value of p
p = 5 # Initial guess
ret = minimize(C,p,args=(mu_0,sigma,c_O,c_PD))
print("Optimal value of p =",ret.x[0])
print("Minimum cost =",ret.fun)


# 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 = False
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()