2749.sagews

def runSim(T, N, numYears):
    k=0
    Ic = 100000    # current Inventory
    I1 = 0         # arrive in one week
    I2 = 0         # arrive in two weeks
    Cmiss  = 100
    Chold  = 10
    Corder = 1000000
    C      = 0
    G  = RealDistribution('gaussian', 7000)
    GG = RealDistribution('gaussian', 10000)
    numWeeks = 52*numYears
    while (k < numWeeks):
        I = Ic + I1 + I2
        if (k % 52) < 47:
            D = 100000 + int(G.get_random_element())     # 00..46 Demand
        else:
            D = 150000 + int(GG.get_random_element())    # 47..51 Demand
        #print [k, D, Ic, I1, I2, C]
        Ic = Ic - D
        if Ic >= 0:
            C = C + Ic * Chold
        else:
            C = C - Ic * Cmiss
            Ic = 0
        Ic = Ic + I1
        I1 = I2
        if I < T:
            I2 = N
            C = C + Corder
        else:
            I2 = 0
        k = k + 1
    cost = n(C/numWeeks)
    #print [T, N, cost]
    return cost

def runXSimsDist(T, N, numYears, numRuns):
    C = []
    for i in range(numRuns):
        C.append(runSim(T, N, numYears))
    mean   = sage.stats.basic_stats.mean(C)
    stddev = sage.stats.basic_stats.std(C)
    return (mean, stddev)
    
def runXSims(T, N, numYears, numRuns):
    x = runXSimsDist(T, N, numYears, numRuns)  # get Mean and StdDev
    # http://www.statsoft.com/textbook/distribution-tables/
    # 1 stddev = 0.5 + 0.3413 = 0.8413 = 85%
    # 2 stddev = 0.5 + 0.4772 = 0.9772 = 98%
    # 3 stddev = 0.5 + 0.4987 = 0.9987 = 99.9%
    costTail = x[0] + 2*x[1]   # = mean + 2*stddev
    return costTail      

def f(x, y):
    return runXSims(x, y, 2, 2)  #adjust numYears,numRuns for shorter(2,2)/longer(5,20) run times
    
def plotTail(T, N):
    r = runXSimsDist(T, N, 40, 200)
    mean = r[0]
    stddev = r[1]
    print "T:",T,"N:",N, "G(",mean,",",stddev,") = ", mean + 2*stddev
    G = RealDistribution('gaussian', stddev)
    z = mean
    v = [(x+z, G.cum_distribution_function(x)) for x in [0, 0.05*stddev,..,3*stddev]]
    x1 = z+0
    y1 = G.cum_distribution_function(x1-z)
    x2 = z+2*stddev
    y2 = G.cum_distribution_function(x2-z)
    #g = line(v) + point((x1,y1),color='red') + text(x1+mean, (x1, y1+0.03)) + point((x2,y2),color='red') + text(x2+mean, (x2, y2+0.03))
    g = line(v) + point((x1,y1),color='red') + point((x2,y2),color='red')
    g += line(((z,y2),(x2,y2),(x2,0.5)),color='green')
    return g

def plotZoom():
    p=10
    nc=100
     #c = contour_plot(f, (300000, 440000), (95000, 245000), cmap='jet', plot_points=p, contours=nc)
    c = contour_plot(f, (325000, 425000), (190000, 210000), cmap='jet', plot_points=p, contours=nc)
    c.show(figsize=(10,10))

def plotMatrix():
    p=20
    nc=20
    d1 = density_plot(f, (50000, 2000000), (50000, 1000000), cmap='jet', plot_points=p)
    d2 = density_plot(f, (50000, 800000), (50000, 500000), cmap='jet', plot_points=p)
    d3 = density_plot(f, (325000, 525000), (90000, 390000), cmap='jet', plot_points=p)
    d4 = density_plot(f, (325000, 525000), (90000, 110000), cmap='jet', plot_points=p)
    d5 = density_plot(f, (325000, 425000), (190000, 210000), cmap='jet', plot_points=p)

    c1 = contour_plot(f, (50000, 2000000), (50000, 1000000), cmap='jet', plot_points=p, contours=nc)
    c2 = contour_plot(f, (50000, 800000), (50000, 500000), cmap='jet', plot_points=p, contours=nc)
    c3 = contour_plot(f, (325000, 525000), (90000, 390000), cmap='jet', plot_points=p, contours=nc)
    c4 = contour_plot(f, (325000, 525000), (90000, 110000), cmap='jet', plot_points=p, contours=nc)
    c5 = contour_plot(f, (325000, 425000), (190000, 210000), cmap='jet', plot_points=p, contours=nc)
    g = graphics_array((d1,c1,d2,c2,d3,c3,d4,c4,d5,c5),5,2)
    g.show(figsize=(20,20))

#print(runSim(400000, 160000, 10))
t1 = plotTail(400000, 200000)
t2 = plotTail(395000, 200000)
t3 = plotTail(390000, 200000)
t4 = plotTail(380000, 200000)
show(t1 + t2 + t3 + t4, figsize=10)
plotZoom()
plotMatrix()
plot3d(f, (50000, 800000), (50000, 500000))

T: 400000 N: 200000 G( 1.65225493995192e6 , 17066.0554947630 ) = 1.68638705094145e6 T: 395000 N: 200000 G( 1.65651633629808e6 , 13710.3384432232 ) = 1.68393701318452e6 T: 390000 N: 200000 G( 1.65662767413461e6 , 14167.2420882978 ) = 1.68496215831121e6 T: 380000 N: 200000 G( 1.65725175173077e6 , 14754.3073237698 ) = 1.68676036637831e6