[Author Prev][Author Next][Thread Prev][Thread Next][Author Index][Thread Index]

[or-cvs] r15966: Replaced broken gnuplot code with numpy/pylab code. (torflow/branches/gsoc2008/tools/BTAnalysis)



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()
+
+