[or-cvs] r15966: Replaced broken gnuplot code with numpy/pylab code. (torflow/branches/gsoc2008/tools/BTAnalysis)
fallon at seul.org
fallon at seul.org
Wed Jul 16 07:51:44 UTC 2008
Author: fallon
Date: 2008-07-16 03:51:44 -0400 (Wed, 16 Jul 2008)
New Revision: 15966
Modified:
torflow/branches/gsoc2008/tools/BTAnalysis/shufflebt.py
Log:
Replaced broken gnuplot code with numpy/pylab code.
Modified: torflow/branches/gsoc2008/tools/BTAnalysis/shufflebt.py
===================================================================
--- torflow/branches/gsoc2008/tools/BTAnalysis/shufflebt.py 2008-07-16 05:14:18 UTC (rev 15965)
+++ torflow/branches/gsoc2008/tools/BTAnalysis/shufflebt.py 2008-07-16 07:51:44 UTC (rev 15966)
@@ -12,6 +12,9 @@
import math,copy
from scipy.integrate import *
from numpy import trapz
+import numpy
+import pylab
+import matplotlib
class Stats:
def __init__(self,file):
@@ -56,7 +59,20 @@
greatest_idx = v
greatest_val = self.buckets[v]
return greatest_idx
-
+
+
+ def pyhist(self,res,histname):
+ bins = len(self.values) / res
+ print 'bins:',bins
+ x = matplotlib.numerix.arange(1,7000, 0.01)
+ S = pypareto(x,0.918058347945, 1600.0, 32775.0)
+ #pylab.hist(self.values,bins=bins,normed=False, width=1)
+ #(n,bins) = numpy.histogram(self.values,bins=bins,normed=False)
+ #pylab.plot(bins,n )
+ pylab.plot(x,S, 'bo')
+ #pylab.show()
+ pylab.savefig(histname + '.png')
+
# XXX: This doesn't seem to work for small #s of circuits
def makehistogram(self,res,histname):
#res = res /1000.0 # convert ms to s
@@ -143,11 +159,11 @@
#ps = fname+'(x)='+str(N*k*(Xm**k))+'/x**('+str(k+1)+')\n'
return ps
-def pypareto(x, k,Xm,N):
+def pypareto(x, k,Xm):
# gnuplot string for shifted, normalized exponential PDF
# g(x,k,B) = (N * k*(Xm**k)/x**(k+1)))
if x<Xm: return 0
- else: return ((((N*k)*(Xm**k)))/((x)**((k+1))))
+ else: return ((((k)*(Xm**k)))/((x)**((k+1))))
def exp(mean,shift,N,fname):
# gnuplot string for normalized exponential PDF
@@ -244,6 +260,7 @@
p = popen2.Popen4(cmd)
p.wait()
+
if __name__ == "__main__":
sort, truncate,graph,dirname,filenames,k,res = getargs()
@@ -281,21 +298,16 @@
# create histogram from file
s = Stats(newfile)
histfilename = histogram_basefilename(filename,sort,truncate,res,dirname)
-# if not sort and not truncate:
-# histfilename = os.path.join(dirname ,os.path.basename(newfile )+ '.res' + str(res) + '.hist')
-# else:
-# histfilename = newfile + '.res' + str(res) +'.hist'
s.makehistogram(res,histfilename + '.hist')
mean = s.mean()
stddev = s.stddev()
median = s.median()
- #mode = s.mode()/10.0 # relies on s.makehistogram for buckets
mode = s.mode() # relies on s.makehistogram for buckets
parK = s.paretoK(mode)
modeN = s.modeN(mode)
modeMean = s.modeMean(mode)
# verify sanity by integrating scaled distribution:
- modeNint = trapz(map(lambda x: pypareto(x, parK, mode, modeN),
+ modeNint = trapz(map(lambda x: modeN* pypareto(x, parK, mode),
xrange(1,200000)))
print 'Resolution of histogram:',res,'ms'
@@ -306,63 +318,20 @@
# get stats
if graph:
- # create gnuplot file
-# if not sort and not truncate:
-# newfile = os.path.join(dirname, newfile)
-# plotname = newfile + '.plt'
- plotname = histfilename + '.plt'
- ncircuits = str(len(s.values))
-
- xtics = max(s.values) / 10.0
- plotstr = "set terminal png transparent nocrop enhanced size 800,600\nset output '" + histfilename + ".png'\nset style fill solid 1.00 border -1\nset style histogram clustered gap 1 title offset character 0, 0, 0\nset datafile missing '-'\nset title 'Buildtime Distribution Function for "+ ncircuits +" Circuits k=" + str(parK) + "\nset ylabel '# Circuits'\nset xlabel 'time (ms)'\nset xtics " + str(xtics) + " \n"
- plotstr += "set label 'std dev=" + str(stddev) + "' at 25000,100\n"
-
- # FIXME: Hrmm... http://en.wikipedia.org/wiki/Skewness? Seems like a hack
- # Or better: http://en.wikipedia.org/wiki/Gamma_distribution with k=3?
- # Would make sense if this is the sum of 3 paretos for the individual
- # hop distributions.
-
- # Theta estimations
- maxliketheta = s.maxlikelihood(k)
- baytheta = s.bayesian(k)
-
- # N is the value to multipy the probabilities by
- N = len(s.values)
-
- #FIX? Other potential values of N: #circuits that match mode? median? mean?
- #print 'Mean:',mean,'Median:', median,'Mode:', mode
- #i = float("%3.0f" % int(mean * 10)) # crappy way of rounding
- #i = int(mode * 10)
- #N = s.buckets[i] # num. circuits that have buildtimes
- #close to mean/median/mode from histogram
-
- # plotstr += gamma(k,maxliketheta,N, 'maxl')
- # plotstr += gamma(k,baytheta[0],N,'bayplus') # + stddev
- # plotstr += gamma(k,baytheta[1],N,'bayminus') # - stddev
-
- plotstr += pareto(parK,mode,modeN,'pareto')
- plotstr += exp(modeMean*10,mode*10,modeN,'expShifted')
+ # plot histogram
+ # args: values, # bins, normalize y/n, width of bars
+ pylab.hist(s.values,len(s.values) / res, normed=True,width=5)
- #plotstr += pareto(parK,mode*10,modeN,'pareto')
- #plotstr += exp(modeMean*10,mode*10,modeN,'expShifted')
- # plotstr += "plot '" + newfile + ".hist' using 2,\\\n"
- plotstr += "plot '" + histfilename + ".hist' using 1:2 with boxes,\\\n"
-
- plotstr += "pareto(x) title '" + "Shifted Pareto' \\\n"
- #plotstr += "expShifted(x) title '" + "Shifted Exp' \n"
-
-
- f = open(plotname,'w')
- f.write(plotstr)
- f.close()
-
- # plot the file
- p = popen2.Popen4('gnuplot ' + plotname)
- #p = popen2.Popen4('gp4.2 ' + plotname) #peculiarity of fallon's system
-
- p.wait()
- for err in p.fromchild:
- print err
-
-
-
+ #plot Pareto curve
+ X = pylab.arange(mode, max(s.values), 1)
+ Y = map(lambda x: pypareto(x, parK, mode), X)
+ n = len(s.values)
+
+
+ pylab.plot(X,Y,'b-')
+
+ #save figure
+ pylab.savefig(histfilename + '.png')
+ pylab.clf()
+
+
More information about the tor-commits
mailing list