In [1]:
import numpy as np
import matplotlib.pyplot as plt
import analysisUtils.analysisUtils as au
import pylab
import scipy.stats as st
import GPy

In [2]:
# read test file
prefix = 'sampling/'
testFile = prefix + 'test.csv'

data = np.genfromtxt(testFile,dtype=float, delimiter=';', skip_header=0,names=True)
print data['column_B']


[ 17.  20.  23.]

In [3]:
cps = [50, 200, 500]
sps = [5, 10, 15, 25, 35]

costsColumn = 'Top_level_costs_regio_central'
runtimeColumn = 'Total_runtime_regio_central'


for cp in cps:
    # load central file 
    
    dataCentral = np.genfromtxt(prefix+'cp_'+str(cp)+'_.csv', dtype=float, delimiter=';', skip_header=0,names=True)
    costsCentral = au.readCol(dataCentral,  'Top_level_costs_central') 
    meanCostsCentral = np.mean(costsCentral)
    
    runtimeCentral = au.readCol(dataCentral, 'Total_runtime_central')
    meanTimeCentral = np.mean(runtimeCentral)
        
    print 'For ' + str(cp) + ' plants: Opt. Costs ', meanCostsCentral, ' Runtime Central ', meanTimeCentral
    
    allSpCosts = []
    for sp in sps:
        fileName = prefix + 'sps_'+str(sp)+ '_cp_'+str(cp)+'_.csv'
        
        data = np.genfromtxt(fileName,dtype=float, delimiter=';', skip_header=0,names=True)
        costs = au.readCol(data, costsColumn)
        meanCosts = np.mean(costs)
        overShoot = meanCosts / meanCostsCentral
        additionalCosts = meanCosts - meanCostsCentral
        
        runtime = au.readCol(data, runtimeColumn) 
        meanTime = np.mean(runtime)
        
        allSpCosts += [additionalCosts / cp]
        print '  '+str(sp) + ' Sps: ' + str(meanCosts) + ' (+ ' + str(overShoot) + ' ; + ' + str(additionalCosts) + ' € ) / ' + str(meanTime)
        
    y_pos = np.arange(len(sps))
    plt.bar(y_pos, allSpCosts, align='center', alpha=0.4)
    plt.xticks(y_pos, sps)
    plt.title(str(cp) + ' plants ')
    plt.figure()


For 50 plants: Opt. Costs  2180.70860152  Runtime Central  6.9299751808
  5 Sps: 2187.73445677 (+ 1.00322182214 ; + 7.02585524524 € ) / 27.2947242019
  10 Sps: 2184.30859122 (+ 1.00165083483 ; + 3.59998970299 € ) / 41.4176581546
  15 Sps: 2183.45375808 (+ 1.00125883695 ; + 2.74515655707 € ) / 52.8120872461
  25 Sps: 2182.05306845 (+ 1.00061652755 ; + 1.34446693238 € ) / 78.1143778371
  35 Sps: 2181.89756481 (+ 1.00054521878 ; + 1.18896328445 € ) / 106.864988682
For 200 plants: Opt. Costs  8587.10311621  Runtime Central  56.3830966762
  5 Sps: 8717.41081394 (+ 1.0151748146 ; + 130.307697733 € ) / 107.054980488
  10 Sps: 8671.0749875 (+ 1.00977883579 ; + 83.971871293 € ) / 172.88823309
  15 Sps: 8713.46094086 (+ 1.01471483723 ; + 126.35782465 € ) / 225.948511393
  25 Sps: 8618.24372105 (+ 1.00362643891 ; + 31.1406048427 € ) / 329.837883109
  35 Sps: 8616.31100852 (+ 1.00340136737 ; + 29.2078923145 € ) / 460.119622181
For 500 plants: Opt. Costs  21157.7182061  Runtime Central  208.723564949
  5 Sps: 21393.5325794 (+ 1.01114554845 ; + 235.814373296 € ) / 244.003133418
  10 Sps: 21428.2394027 (+ 1.01278593438 ; + 270.521196574 € ) / 357.594805563
  15 Sps: 21339.9226693 (+ 1.00861172559 ; + 182.204463221 € ) / 470.31848035
  25 Sps: 21180.4995945 (+ 1.00107674127 ; + 22.7813884284 € ) / 694.298222776
  35 Sps: 21176.5500575 (+ 1.00089007006 ; + 18.831851436 € ) / 934.632445593

In [12]:
cps = [50, 200, 500]
sps = [5, 10, 15, 25, 35]

costsColumn = 'Top_level_costs_regio_central'
runtimeColumn = 'Total_runtime_regio_central'

allOffSetCosts = { 5 : [],
                   10: [],
                   15 : [],
                   25 : [],
                   35 : [] }# for sp = 5 a list (50, 200, 500), for sp = 10 a list (50, 200, 500) ...

allOffSetCostsStd = { 5 : [],
                   10: [],
                   15: [],
                   25 : [],
                   35 : [] }# for sp = 5 a list (50, 200, 500), for sp = 10 a list (50, 200, 500) ...

allRuntimes = { 5 : [],
                   10: [],
                   15 : [],
                   25 : [],
                   35 : [] }# for sp = 5 a list (50, 200, 500), for sp = 10 a list (50, 200, 500) ...
for cp in cps:
    # load central file 
    
    dataCentral = np.genfromtxt(prefix+'cp_'+str(cp)+'_.csv', dtype=float, delimiter=';', skip_header=0,names=True)
    costsCentral = au.readCol(dataCentral,  'Top_level_costs_central') 
    meanCostsCentral = np.mean(costsCentral)
    
    runtimeCentral = au.readCol(dataCentral, 'Total_runtime_central')
    meanTimeCentral = np.mean(runtimeCentral)
        
    print 'For ' + str(cp) + ' plants: Opt. Costs ', meanCostsCentral, ' Runtime Central ', meanTimeCentral
    
    allSpCosts = []
    for sp in sps:
        fileName = prefix + 'sps_'+str(sp)+ '_cp_'+str(cp)+'_.csv'
        
        data = np.genfromtxt(fileName,dtype=float, delimiter=';', skip_header=0,names=True)
        costs = au.readCol(data, costsColumn)
        
        allExtraCosts = costs - costsCentral
        meanExtraCosts = np.mean(allExtraCosts)
        overShoot = meanExtraCosts / meanCostsCentral
        stdExtraCosts = np.sqrt(np.std(allExtraCosts))
        stdExtraCosts = 0.25 * np.std(allExtraCosts)
      
        allOffSetCosts[sp] += [meanExtraCosts / cp]
        allOffSetCostsStd[sp] += [stdExtraCosts / cp]
                                  
        runtime = au.readCol(data, runtimeColumn) 
        meanTime = np.mean(runtime)

        print '  '+str(sp) + ' Sps: ' + str(meanExtraCosts) + ' (+ ' + str(overShoot) + ' ; + ' + str(meanExtraCosts) + ' € ) / ' + str(meanTime)
        
    #y_pos = np.arange(len(sps))
    #plt.bar(y_pos, allSpCosts, align='center', alpha=0.4)
    #plt.xticks(y_pos, sps)
    #plt.title(str(cp) + ' plants ')
    #plt.figure()


For 50 plants: Opt. Costs  2180.70860152  Runtime Central  6.9299751808
  5 Sps: 7.02585524524 (+ 0.00322182213632 ; + 7.02585524524 € ) / 27.2947242019
  10 Sps: 3.599989703 (+ 0.0016508348252 ; + 3.599989703 € ) / 41.4176581546
  15 Sps: 2.74515655707 (+ 0.00125883694647 ; + 2.74515655707 € ) / 52.8120872461
  25 Sps: 1.34446693238 (+ 0.000616527550469 ; + 1.34446693238 € ) / 78.1143778371
  35 Sps: 1.18896328445 (+ 0.000545218780549 ; + 1.18896328445 € ) / 106.864988682
For 200 plants: Opt. Costs  8587.10311621  Runtime Central  56.3830966762
  5 Sps: 130.307697733 (+ 0.0151748145992 ; + 130.307697733 € ) / 107.054980488
  10 Sps: 83.971871293 (+ 0.00977883579091 ; + 83.971871293 € ) / 172.88823309
  15 Sps: 126.35782465 (+ 0.0147148372321 ; + 126.35782465 € ) / 225.948511393
  25 Sps: 31.1406048427 (+ 0.00362643890742 ; + 31.1406048427 € ) / 329.837883109
  35 Sps: 29.2078923145 (+ 0.00340136736676 ; + 29.2078923145 € ) / 460.119622181
For 500 plants: Opt. Costs  21157.7182061  Runtime Central  208.723564949
  5 Sps: 235.814373296 (+ 0.0111455484471 ; + 235.814373296 € ) / 244.003133418
  10 Sps: 270.521196574 (+ 0.0127859343781 ; + 270.521196574 € ) / 357.594805563
  15 Sps: 182.204463221 (+ 0.0086117255862 ; + 182.204463221 € ) / 470.31848035
  25 Sps: 22.7813884284 (+ 0.00107674127269 ; + 22.7813884284 € ) / 694.298222776
  35 Sps: 18.831851436 (+ 0.000890070056352 ; + 18.831851436 € ) / 934.632445593

In [13]:
N = len(cps)
ind = np.arange(N)
width = 0.15

# find gray values

colors =  { 5 : '0.95', 10: '0.8', 15 : '0.65', 25 : '0.45', 35 : '0.25' }
hatches = { 5 : '//', 10: '\\\\', 15 : 'x', 25 : 'o', 35 : '-' }
fig, ax = plt.subplots()
i = 0
for sp in sps:
    rects = ax.bar(ind+i*width, allOffSetCosts[sp], width, color=colors[sp], yerr = allOffSetCostsStd[sp], ecolor='black', hatch = hatches[sp], label=str(sp) + ' sps')
    i+=1
    
ax.set_ylabel('Extra costs per plant (EUR)') 
ax.set_xticks(ind+width)
ax.set_xticklabels(cps)
ax.set_xlabel('Problem sizes (\# plants)')
plt.legend(loc=2)
savefig('SamplingCosts.pdf')
plt.show()



In [6]:
from pylab import arange,pi,sin,cos,sqrt
fig_width_pt = 400  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inch
golden_mean = (sqrt(5)-1.0)/2.0         # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt  # width in inches
fig_height = fig_width*golden_mean      # height in inches
fig_size =  [fig_width,fig_height]
params = {'backend': 'ps',
          'axes.labelsize': 12,
          'text.fontsize': 12,
          'legend.fontsize': 12,
          'xtick.labelsize': 10,
          'ytick.labelsize': 10,
          'text.usetex': True,
          'figure.figsize': fig_size}
pylab.rcParams.update(params)

In [7]:
cps = [50, 200, 500]
sps = [5, 10, 15, 25, 35]

costsColumn = 'Top_level_costs_regio_central'
runtimeColumn = 'Total_runtime_regio_central'


topString = "\\begin{tabular*}{\\textwidth}{@{\extracolsep{\\fill} }l"
cols = 5
singleCol = "d{3}"
colString = singleCol * cols
topString = topString + colString + "}"
print topString
print "\\toprule"

colHeaders = ['\# sampling points', 'costs per step', 'rel. overhead', 'abs. overhead', 'runtime']
headers = ""
first = True 
for h in colHeaders:
    if(first):
        first = False
    else:
        headers = headers + " & "
    headers = headers + "\multicolumn{1}{c}{"+h+"}"

print headers + "\\\\"

for cp in cps:
    # load central file 
    
    dataCentral = np.genfromtxt(prefix+'cp_'+str(cp)+'_.csv', dtype=float, delimiter=';', skip_header=0,names=True)
    costsCentral = au.readCol(dataCentral,  'Top_level_costs_central') 
    meanCostsCentral = np.mean(costsCentral)
    stdCostsCentral = np.std(costsCentral)
    
    runtimeCentral = au.readCol(dataCentral, 'Total_runtime_central')
    meanTimeCentral = np.mean(runtimeCentral)
    stdTimeCentral = np.std(runtimeCentral)
    print '\midrule'    
    print '\multicolumn{'+str(cols)+'}{l}{'+str(cp) + ' power plants : } \\\\' # (Mean Opt. Costs ', np.around(meanCostsCentral, 2),' ):  \\\\' # Opt. Costs ', meanCostsCentral, ' Runtime Central ', meanTimeCentral
    print '\midrule'    
    allSpCosts = []
    for sp in sps:
        fileName = prefix + 'sps_'+str(sp)+ '_cp_'+str(cp)+'_.csv'
        
        data = np.genfromtxt(fileName,dtype=float, delimiter=';', skip_header=0,names=True)
        costs = au.readCol(data, costsColumn)
        
        allExtraCosts = costs - costsCentral
        meanExtraCosts = np.mean(allExtraCosts)
        
        allOverShoots = ((costs / costsCentral ) - 1) * 100
        stdOvershoot = np.std(allOverShoots)
        overShoot = meanExtraCosts / meanCostsCentral
        stdExtraCosts = np.sqrt(np.std(allExtraCosts))
            
        meanCosts = np.mean(costs)
        stdCosts = np.std(costs)
        overShoot = ((meanCosts / meanCostsCentral) - 1) * 100
        additionalCosts = meanCosts - meanCostsCentral
        
        runtime = au.readCol(data, runtimeColumn) 
        meanTime = np.mean(runtime)
        stdTime = np.std(runtime)
        
        allSpCosts += [additionalCosts / cp]
        print '  '+str(sp) + ' sps & ' + str(np.around(meanCosts, 2)) + '~$\euro$ ~('+str(np.around(stdCosts, 2))+')& + ' + str(np.around(overShoot, 2)) + ' \% ~('+str(np.around(stdOvershoot, 2))+')& + ' + str(np.around(additionalCosts, 2)) + ' ~$\euro$ ~('+str(np.around(stdExtraCosts, 2))+') & ' + str(np.around(meanTime, 2)) + " ~\\text{sec} ~("+str(np.around(stdTime, 2)) +") \\\\"
    print '  \\emph{Optimum} & ' + au.boldify(str(np.around(meanCostsCentral, 2))) + "~$\euro$ ~("+str(np.around(stdCostsCentral, 2))+ ") \\\\"
print '\\bottomrule'
print '\end{tabular*}'


\begin{tabular*}{\textwidth}{@{\extracolsep{\fill} }ld{3}d{3}d{3}d{3}d{3}}
\toprule
\multicolumn{1}{c}{\# sampling points} & \multicolumn{1}{c}{costs per step} & \multicolumn{1}{c}{rel. overhead} & \multicolumn{1}{c}{abs. overhead} & \multicolumn{1}{c}{runtime}\\
\midrule
\multicolumn{5}{l}{50 power plants : } \\
\midrule
  5 sps & 2187.73~$\euro$ ~(453.57)& + 0.32 \% ~(0.29)& + 7.03 ~$\euro$ ~(2.37) & 27.29 ~\text{sec} ~(2.82) \\
  10 sps & 2184.31~$\euro$ ~(453.87)& + 0.17 \% ~(0.16)& + 3.6 ~$\euro$ ~(1.79) & 41.42 ~\text{sec} ~(4.16) \\
  15 sps & 2183.45~$\euro$ ~(453.72)& + 0.13 \% ~(0.15)& + 2.75 ~$\euro$ ~(1.77) & 52.81 ~\text{sec} ~(2.52) \\
  25 sps & 2182.05~$\euro$ ~(454.01)& + 0.06 \% ~(0.06)& + 1.34 ~$\euro$ ~(1.1) & 78.11 ~\text{sec} ~(2.92) \\
  35 sps & 2181.9~$\euro$ ~(453.86)& + 0.05 \% ~(0.07)& + 1.19 ~$\euro$ ~(1.15) & 106.86 ~\text{sec} ~(3.67) \\
  \emph{Optimum} & \textbf{2180}.\textbf{71}~$\euro$ ~(453.91) \\
\midrule
\multicolumn{5}{l}{200 power plants : } \\
\midrule
  5 sps & 8717.41~$\euro$ ~(1014.05)& + 1.52 \% ~(1.53)& + 130.31 ~$\euro$ ~(12.28) & 107.05 ~\text{sec} ~(3.61) \\
  10 sps & 8671.07~$\euro$ ~(975.07)& + 0.98 \% ~(1.36)& + 83.97 ~$\euro$ ~(11.2) & 172.89 ~\text{sec} ~(33.44) \\
  15 sps & 8713.46~$\euro$ ~(995.13)& + 1.47 \% ~(2.22)& + 126.36 ~$\euro$ ~(14.02) & 225.95 ~\text{sec} ~(43.85) \\
  25 sps & 8618.24~$\euro$ ~(938.57)& + 0.36 \% ~(0.83)& + 31.14 ~$\euro$ ~(8.27) & 329.84 ~\text{sec} ~(20.76) \\
  35 sps & 8616.31~$\euro$ ~(952.83)& + 0.34 \% ~(0.7)& + 29.21 ~$\euro$ ~(8.14) & 460.12 ~\text{sec} ~(51.72) \\
  \emph{Optimum} & \textbf{8587}.\textbf{1}~$\euro$ ~(938.36) \\
\midrule
\multicolumn{5}{l}{500 power plants : } \\
\midrule
  5 sps & 21393.53~$\euro$ ~(1885.9)& + 1.11 \% ~(1.08)& + 235.81 ~$\euro$ ~(14.97) & 244.0 ~\text{sec} ~(11.28) \\
  10 sps & 21428.24~$\euro$ ~(1868.18)& + 1.28 \% ~(1.16)& + 270.52 ~$\euro$ ~(15.51) & 357.59 ~\text{sec} ~(10.71) \\
  15 sps & 21339.92~$\euro$ ~(1870.11)& + 0.86 \% ~(0.99)& + 182.2 ~$\euro$ ~(14.37) & 470.32 ~\text{sec} ~(18.24) \\
  25 sps & 21180.5~$\euro$ ~(1853.91)& + 0.11 \% ~(0.04)& + 22.78 ~$\euro$ ~(2.86) & 694.3 ~\text{sec} ~(43.68) \\
  35 sps & 21176.55~$\euro$ ~(1853.2)& + 0.09 \% ~(0.09)& + 18.83 ~$\euro$ ~(4.28) & 934.63 ~\text{sec} ~(26.98) \\
  \emph{Optimum} & \textbf{21157}.\textbf{72}~$\euro$ ~(1853.51) \\
\bottomrule
\end{tabular*}

In [8]:
cps = [50, 200, 500]
sps = [5, 10, 15, 25, 35]

costsColumn = 'Top_level_costs_regio_central'
runtimeColumn = 'Total_runtime_regio_central'


for cp in cps:
    # load central file 
    
    allSpCosts = []
    
    
    allCosts = { 5 : [],
                 15 : [],
                 25 : [],
                 35 : [] }
    print '--------------------------'
    print 'Considering ', str(cp), ' plants'
    for sp in sps:
        fileName = prefix + 'sps_'+str(sp)+ '_cp_'+str(cp)+'_.csv'
        
        data = np.genfromtxt(fileName,dtype=float, delimiter=';', skip_header=0,names=True)
        costs = au.readCol(data, costsColumn)
        
        allCosts[sp] = costs
    
    for sp in sps:
        print 'Comparing ',sp,' to optimal solution: '
        dataCentral = np.genfromtxt(prefix+'cp_'+str(cp)+'_.csv', dtype=float, delimiter=';', skip_header=0,names=True)
        costsCentral = au.readCol(dataCentral,  'Top_level_costs_central') 
        costsSp = np.array(allCosts[sp])
        [t, prob] = st.ttest_rel(np.array(allCosts[sp]), np.array(costsCentral))
        print 't: ', t, ' prob: ', prob

        if prob < 0.01:
            print '    SIGNIFICANT'
        else:
            print '    insignficant'
        
        for otherSp in sps:
            if otherSp > sp:
                costsOtherSp = np.array(allCosts[otherSp])
                np.array(allCosts[sp])
                print 'testing : ', sp, ' and ', otherSp
                #print '  var ', sp, ' : ', np.var(costsSp)
                #print '  var ', otherSp, ' : ', np.var(np.array(allCosts[otherSp]))
                #print '  bartlett, ', st.bartlett(costsOtherSp, costsSp)                               
                
                [t, prob] = st.ttest_rel(np.array(allCosts[sp]), np.array(allCosts[otherSp]))
                print '  t: ', t, ' prob: ', prob

                if prob < 0.01:
                    print '    SIGNIFICANT'
                else:
                    print '    insignficant'
                       
    print '------------------------'


--------------------------
Considering  50  plants
Comparing  5  to optimal solution: 
t:  19.786969805  prob:  5.27648947667e-53
    SIGNIFICANT
testing :  5  and  10
  t:  9.76691653978  prob:  2.76298312738e-19
    SIGNIFICANT
testing :  5  and  15
  t:  11.3538893479  prob:  2.37423459545e-24
    SIGNIFICANT
testing :  5  and  25
  t:  16.8457852209  prob:  5.13141065328e-43
    SIGNIFICANT
testing :  5  and  35
  t:  17.207765701  prob:  2.94643959458e-44
    SIGNIFICANT
Comparing  10  to optimal solution: 
t:  17.6846161866  prob:  6.89184483104e-46
    SIGNIFICANT
testing :  10  and  15
  t:  3.72047251117  prob:  0.000245795213229
    SIGNIFICANT
testing :  10  and  25
  t:  12.4086095886  prob:  7.57614762219e-28
    SIGNIFICANT
testing :  10  and  35
  t:  12.4532018339  prob:  5.36978094445e-28
    SIGNIFICANT
Comparing  15  to optimal solution: 
t:  13.8799897435  prob:  7.85999976866e-33
    SIGNIFICANT
testing :  15  and  25
  t:  8.17593321325  prob:  1.50153151191e-14
    SIGNIFICANT
testing :  15  and  35
  t:  8.56217999465  prob:  1.16370660937e-15
    SIGNIFICANT
Comparing  25  to optimal solution: 
t:  17.3870231506  prob:  7.17192086156e-45
    SIGNIFICANT
testing :  25  and  35
  t:  2.14324667387  prob:  0.033062623846
    insignficant
Comparing  35  to optimal solution: 
t:  14.0704732941  prob:  1.75491409521e-33
    SIGNIFICANT
------------------------
--------------------------
Considering  200  plants
Comparing  5  to optimal solution: 
t:  13.6339866948  prob:  5.43012418784e-32
    SIGNIFICANT
testing :  5  and  10
  t:  5.40155122682  prob:  1.5425184275e-07
    SIGNIFICANT
testing :  5  and  15
  t:  0.324318430943  prob:  0.745969242726
    insignficant
testing :  5  and  25
  t:  9.13352386796  prob:  2.37221624022e-17
    SIGNIFICANT
testing :  5  and  35
  t:  10.2336419659  prob:  9.57740912093e-21
    SIGNIFICANT
Comparing  10  to optimal solution: 
t:  10.563364853  prob:  8.58499628291e-22
    SIGNIFICANT
testing :  10  and  15
  t:  -3.72870592735  prob:  0.000238315946795
    SIGNIFICANT
testing :  10  and  25
  t:  6.04799367263  prob:  5.34102578915e-09
    SIGNIFICANT
testing :  10  and  35
  t:  7.06907282722  prob:  1.56487827312e-11
    SIGNIFICANT
Comparing  15  to optimal solution: 
t:  10.1441141571  prob:  1.83426280649e-20
    SIGNIFICANT
testing :  15  and  25
  t:  7.24659672895  prob:  5.34882983706e-12
    SIGNIFICANT
testing :  15  and  35
  t:  7.80962640523  prob:  1.59821185259e-13
    SIGNIFICANT
Comparing  25  to optimal solution: 
t:  7.1890919001  prob:  7.58697215766e-12
    SIGNIFICANT
testing :  25  and  35
  t:  0.305705658439  prob:  0.760084185501
    insignficant
Comparing  35  to optimal solution: 
t:  6.9529372591  prob:  3.12995422954e-11
    SIGNIFICANT
------------------------
--------------------------
Considering  500  plants
Comparing  5  to optimal solution: 
t:  16.5966162656  prob:  3.67777431003e-42
    SIGNIFICANT
testing :  5  and  10
  t:  -1.74896287426  prob:  0.0815296153437
    insignficant
testing :  5  and  15
  t:  3.02438184958  prob:  0.00275212569962
    SIGNIFICANT
testing :  5  and  25
  t:  14.9713525358  prob:  1.4262084272e-36
    SIGNIFICANT
testing :  5  and  35
  t:  15.1577158658  prob:  3.26234488841e-37
    SIGNIFICANT
Comparing  10  to optimal solution: 
t:  17.7452141101  prob:  4.27985290947e-46
    SIGNIFICANT
testing :  10  and  15
  t:  5.21864662852  prob:  3.79884793523e-07
    SIGNIFICANT
testing :  10  and  25
  t:  16.2098120383  prob:  7.84882055062e-41
    SIGNIFICANT
testing :  10  and  35
  t:  16.3903056181  prob:  1.88094269576e-41
    SIGNIFICANT
Comparing  15  to optimal solution: 
t:  13.9289947623  prob:  5.34557340332e-33
    SIGNIFICANT
testing :  15  and  25
  t:  12.2418695543  prob:  2.73748023598e-27
    SIGNIFICANT
testing :  15  and  35
  t:  12.4191132723  prob:  6.98628672999e-28
    SIGNIFICANT
Comparing  25  to optimal solution: 
t:  44.1022312874  prob:  1.18041024984e-119
    SIGNIFICANT
testing :  25  and  35
  t:  3.3580717626  prob:  0.000907824600767
    SIGNIFICANT
Comparing  35  to optimal solution: 
t:  16.2049322796  prob:  8.1579624869e-41
    SIGNIFICANT
------------------------

In [5]:
cps = [500]
sps = [5, 10, 15, 25, 35]

costsColumn = 'Top_level_costs_regio_central'
runtimeColumn = 'Total_runtime_regio_central'

allOffSetCosts = { 5 : [],
                   10: [],
                   15 : [],
                   25 : [],
                   35 : [] }# for sp = 5 a list (50, 200, 500), for sp = 10 a list (50, 200, 500) ...

allOffSetCostsStd = { 5 : [],
                   10: [],
                   15: [],
                   25 : [],
                   35 : [] }# for sp = 5 a list (50, 200, 500), for sp = 10 a list (50, 200, 500) ...

allRuntimes = { 5 : [],
                   10: [],
                   15 : [],
                   25 : [],
                   35 : [] }# for sp = 5 a list (50, 200, 500), for sp = 10 a list (50, 200, 500) ...
for cp in cps:
    # load central file 
    
    dataCentral = np.genfromtxt(prefix+'cp_'+str(cp)+'_.csv', dtype=float, delimiter=';', skip_header=0,names=True)
    costsCentral = au.readCol(dataCentral,  'Top_level_costs_central') 
    meanCostsCentral = np.mean(costsCentral)
    
    runtimeCentral = au.readCol(dataCentral, 'Total_runtime_central')
    meanTimeCentral = np.mean(runtimeCentral)
        
    print 'For ' + str(cp) + ' plants: Opt. Costs ', meanCostsCentral, ' Runtime Central ', meanTimeCentral
    
    allSpCosts = []
    allExcessCosts = []
    
    for sp in sps:
        fileName = prefix + 'sps_'+str(sp)+ '_cp_'+str(cp)+'_.csv'
        
        data = np.genfromtxt(fileName,dtype=float, delimiter=';', skip_header=0,names=True)
        costs = au.readCol(data, costsColumn)
        
        allExtraCosts = costs - costsCentral
        meanExtraCosts = np.mean(allExtraCosts)
        overShoot = meanExtraCosts / meanCostsCentral
        stdExtraCosts = np.sqrt(np.std(allExtraCosts))
        stdExtraCosts = 0.25 * np.std(allExtraCosts)
      
        allOffSetCosts[sp] += [meanExtraCosts / cp]
        allOffSetCostsStd[sp] += [stdExtraCosts / cp]
                                  
        allExcessCosts += [meanExtraCosts]
        
        runtime = au.readCol(data, runtimeColumn) 
        meanTime = np.mean(runtime)

        print '  '+str(sp) + ' Sps: ' + str(meanExtraCosts) + ' (+ ' + str(overShoot) + ' ; + ' + str(meanExtraCosts) + ' € ) / ' + str(meanTime)
    kernel = GPy.kern.RBF(input_dim=1, variance=1., lengthscale=1.)
    m = GPy.models.GPRegression(np.atleast_2d(sps).T, np.atleast_2d(allExcessCosts).T, kernel)
    m.optimize_restarts(num_restarts=10)
    m.plot()
    xlabel('\# sampling points')
    ylabel('Extra costs per plant (EUR)')
    #title('Extra costs for \#sampling points')
    savefig('sampling-extra-costs'+str(cp)+'.pdf')
    figure()


For 500 plants: Opt. Costs  21157.7182061  Runtime Central  208.723564949
  5 Sps: 235.814373296 (+ 0.0111455484471 ; + 235.814373296 € ) / 244.003133418
  10 Sps: 270.521196574 (+ 0.0127859343781 ; + 270.521196574 € ) / 357.594805563
  15 Sps: 182.204463221 (+ 0.0086117255862 ; + 182.204463221 € ) / 470.31848035
  25 Sps: 22.7813884284 (+ 0.00107674127269 ; + 22.7813884284 € ) / 694.298222776
  35 Sps: 18.831851436 (+ 0.000890070056352 ; + 18.831851436 € ) / 934.632445593
Optimization restart 1/10, f = 29.2793853059
Optimization restart 2/10, f = 29.3367334222
Optimization restart 3/10, f = 33.0749381241
Optimization restart 4/10, f = 33.0860770604
Optimization restart 5/10, f = 29.3367334218
Optimization restart 6/10, f = 31.9926001141
Optimization restart 7/10, f = 29.2793770319
Optimization restart 8/10, f = 31.9971015787
Optimization restart 9/10, f = 33.0831693011
Optimization restart 10/10, f = 29.2793776557

In [10]:
from pylab import arange,pi,sin,cos,sqrt
fig_width_pt = 400  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inch
golden_mean = (sqrt(5)-1.0)/2.0         # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt  # width in inches
fig_height = fig_width*golden_mean      # height in inches
fig_size =  [fig_width,fig_height]
params = {'backend': 'ps',
          'axes.labelsize': 10,
          'text.fontsize': 10,
          'legend.fontsize': 10,
          'xtick.labelsize': 8,
          'ytick.labelsize': 8,
          'text.usetex': True,
          'figure.figsize': fig_size}
pylab.rcParams.update(params)

In [10]:
cps = [500]
sps = [5, 10, 15, 25, 35]

costsColumn = 'Top_level_costs_regio_central'
runtimeColumn = 'Total_runtime_regio_central'

allOffSetCosts = { 5 : [],
                   10: [],
                   15 : [],
                   25 : [],
                   35 : [] }# for sp = 5 a list (50, 200, 500), for sp = 10 a list (50, 200, 500) ...

allOffSetCostsStd = { 5 : [],
                   10: [],
                   15: [],
                   25 : [],
                   35 : [] }# for sp = 5 a list (50, 200, 500), for sp = 10 a list (50, 200, 500) ...

allRuntimes = { 5 : [],
                   10: [],
                   15 : [],
                   25 : [],
                   35 : [] }# for sp = 5 a list (50, 200, 500), for sp = 10 a list (50, 200, 500) ...
for cp in cps:
    # load central file 
    
    dataCentral = np.genfromtxt(prefix+'cp_'+str(cp)+'_.csv', dtype=float, delimiter=';', skip_header=0,names=True)
    costsCentral = au.readCol(dataCentral,  'Top_level_costs_central') 
    meanCostsCentral = np.mean(costsCentral)
    
    runtimeCentral = au.readCol(dataCentral, 'Total_runtime_central')
    meanTimeCentral = np.mean(runtimeCentral)
        
    print 'For ' + str(cp) + ' plants: Opt. Costs ', meanCostsCentral, ' Runtime Central ', meanTimeCentral
    
    allSpCosts = []
    allExcessCosts = []
    
    allInputTimes = []
    for sp in sps:
        fileName = prefix + 'sps_'+str(sp)+ '_cp_'+str(cp)+'_.csv'
        
        data = np.genfromtxt(fileName,dtype=float, delimiter=';', skip_header=0,names=True)
        costs = au.readCol(data, costsColumn)
        
        allExtraCosts = costs - costsCentral
        meanExtraCosts = np.mean(allExtraCosts)
        overShoot = meanExtraCosts / meanCostsCentral
        stdExtraCosts = np.sqrt(np.std(allExtraCosts))
        stdExtraCosts = 0.25 * np.std(allExtraCosts)
      
        allOffSetCosts[sp] += [meanExtraCosts / cp]
        allOffSetCostsStd[sp] += [stdExtraCosts / cp]
                                  
        allExcessCosts += [meanExtraCosts]
        
        
        runtime = au.readCol(data, runtimeColumn) 
        meanTime = np.mean(runtime)
        
        allInputTimes += [meanTime]

        print '  '+str(sp) + ' Sps: ' + str(meanExtraCosts) + ' (+ ' + str(overShoot) + ' ; + ' + str(meanExtraCosts) + ' € ) / ' + str(meanTime)
    kernel = GPy.kern.Brownian(input_dim=1)
    m = GPy.models.GPRegression(np.atleast_2d(allInputTimes).T, np.atleast_2d(allExcessCosts).T, kernel)
    m.optimize_restarts(num_restarts=10)
    m.plot()
    xlabel('\# sampling points')
    ylabel('Extra costs per plant (EUR)')
    #title('Extra costs for \#sampling points')
    savefig('sampling-extra-time'+str(cp)+'.pdf')
    figure()


For 500 plants: Opt. Costs  21157.7182061  Runtime Central  208.723564949
  5 Sps: 235.814373296 (+ 0.0111455484471 ; + 235.814373296 € ) / 244.003133418
  10 Sps: 270.521196574 (+ 0.0127859343781 ; + 270.521196574 € ) / 357.594805563
  15 Sps: 182.204463221 (+ 0.0086117255862 ; + 182.204463221 € ) / 470.31848035
  25 Sps: 22.7813884284 (+ 0.00107674127269 ; + 22.7813884284 € ) / 694.298222776
  35 Sps: 18.831851436 (+ 0.000890070056352 ; + 18.831851436 € ) / 934.632445593
Optimization restart 1/10, f = 31.1032556492
Optimization restart 2/10, f = 31.1032566267
Optimization restart 3/10, f = 31.1032506195
Optimization restart 4/10, f = 31.1032563493
Optimization restart 5/10, f = 31.1032555507
Optimization restart 6/10, f = 31.1032519085
Optimization restart 7/10, f = 31.1032524311
Optimization restart 8/10, f = 31.1032560239
Optimization restart 9/10, f = 31.1032496756
Optimization restart 10/10, f = 31.1032572718

In [ ]: