Index: trunk/UTIL/PYTHON/powerlaw/ddhist.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/ddhist.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/ddhist.py	(revision 1053)
@@ -0,0 +1,140 @@
+#! /usr/bin/env python
+import matplotlib.pyplot as mpl
+import ppplot
+import numpy as np
+import scipy.optimize as sciopt
+
+## SAVED
+namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = False ; typefit=1
+namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = True ; typefit=2
+namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = False ; typefit=1
+namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = True ; typefit=1
+namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt" ; drop = False ; typefit=1 #bof. tous les 10 est mieux.
+
+##########################################################
+#namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm1_1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm2_1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm1_1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt"
+namefile = "/home/aymeric/Big_Data/LES_dd/press_ustm_exomars.nc1.txt"
+namefile = "/home/aymeric/Big_Data/LES_dd/sav/press_ustm_exomars.ncm2_1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_1.txt"
+drop = False ; typefit=1
+#drop = True ; typefit=1
+#drop = True ; typefit=2
+##########################################################
+## define bins (we expect fit do not depend too much on this -- to be checked)
+##nbins = 7 ; www = 3.
+#nbins = 10 ; www = 2.5
+nbins = 15 ; www = 2. 
+##nbins = 20 ; www = 1.75
+##nbins = 30 ; www = 1.5
+##nbins = 50 ; www = 1.2
+##########################################################
+
+# -- functions for fitting
+# -- http://wiki.scipy.org/Cookbook/FittingData
+# power law
+def fitfunc(x,a,b): return a*(x**(-b))
+# power law + constant
+def fitfunc3(x,a,b,c): return a*(x**(-b)) + c * (x**(-0.5))
+# exponential law
+def fitfunc2(x,a,b): return a*np.exp(-b*x)
+
+# load data
+data = np.loadtxt(namefile,delimiter=";")
+t = data[:,0] ; s = data[:,1] ; d = data[:,2]
+i = data[:,3] ; j = data[:,4] # int?
+
+# a way to guess resolution
+dx = np.min(s)/2. ; print "resolution: ", dx
+
+# choose a variable
+if drop: var = d
+else: var = s
+
+# restrictions
+restrict = (s > 0) # initialization (True everywhere)
+#restrict = restrict*(np.abs(i-j) <= 6.*dx) # condition sur i,j (width,height)
+#restrict = restrict*(d > 0.9) # limit on drop for size (casse la power law? seulement si keep smaller devils)
+#restrict = restrict*(d > 0.5)
+#restrict = restrict*(s >= 3.*dx) # remove lowest sizes (detection limit) 
+var = var[restrict]
+
+## define bins
+zebins = [np.min(var)]
+for i in range(0,nbins):  zebins.append(zebins[i]*(www**0.5))
+zebins = np.array(zebins)
+middle = 0.5*(zebins[1:] + zebins[:-1])
+binwidth = zebins[1:] - zebins[:-1]
+
+# plot histogram
+yeah = mpl.hist(var,log=True,bins=zebins,normed=True,color='white') ; mpl.xscale('log')
+#yeah = mpl.hist(var,bins=zebins)#,normed=True)
+
+# print info
+print "min %5.2e // max %5.2e" % (np.min(var),np.max(var))
+for iii in range(len(zebins)-1):
+    print "%5.2e in [%5.2e %5.2e]" % (yeah[0][iii],zebins[iii],zebins[iii+1])
+
+# fitting function to superimpose
+if typefit == 1: xx = sciopt.curve_fit(fitfunc, middle, yeah[0])
+elif typefit == 2: xx = sciopt.curve_fit(fitfunc2, middle, yeah[0])
+elif typefit == 3: xx = sciopt.curve_fit(fitfunc3, middle, yeah[0])
+print "exponent",xx[0][1]
+print "variance %",100.*xx[1][1][1]/xx[0][1]
+
+# plot obtained fit along with actual points
+if typefit == 1: func = fitfunc(middle,xx[0][0],xx[0][1])
+elif typefit == 2: func = fitfunc2(middle,xx[0][0],xx[0][1])
+elif typefit == 3: func = fitfunc3(middle,xx[0][0],xx[0][1],xx[0][2]) 
+func[func<0]=np.nan
+ind = yeah[0]>0
+mpl.plot(middle[ind],yeah[0][ind],'k.')
+mpl.plot(middle[ind],func[ind],'r-')
+mpl.plot(middle[ind],func[ind],'r.')
+
+# print fit vs. actual populations
+total = var.shape[0]
+fit = func*binwidth*total
+pop = yeah[0]*binwidth*total
+for iii in range(len(fit)):
+    if fit[iii] > 0.99 or pop[iii] != 0:
+        print "fit %5.0f actual %5.0f" % (fit[iii],pop[iii])
+
+# plot settings
+ax = mpl.gca() ; ax.set_xbound(lower=np.min(zebins),upper=np.max(zebins))
+if drop: mpl.xlabel("Pressure drop (Pa)") ; ax.set_xbound(lower=1.e-1,upper=10.)
+else: mpl.xlabel("Vortex size (m)") ; ax.set_xbound(lower=10.,upper=5000.)
+mpl.ylabel('Population density $n / N w_{bin}$')
+mpl.title("Statistics on "+str(total)+" detected vortices")
+
+# show plot and end
+mpl.show()
+exit()
+
+################################################################################
+################################################################################
+################################################################################
+
+## UNCOMMENT
+##import plfit,plpva,plplot
+
+#[alpha, xmin, L] = plfit.plfit(var,'xmin',0.3)
+[alpha, xmin, L] = plfit.plfit(var)
+print alpha,xmin
+
+#a = plpva.plpva(var,0.75,'xmin',0.75)
+#print a
+
+
+h = plplot.plplot(var,xmin,alpha)
+
+ppplot.save(mode="png")
+
+#mpl.loglog(h[0], h[1], 'k--',linewidth=2)
+
+#mpl.show()
+
+
+
Index: trunk/UTIL/PYTHON/powerlaw/ddstat.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/ddstat.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/ddstat.py	(revision 1053)
@@ -0,0 +1,54 @@
+#! /usr/bin/env python
+import ppplot
+import numpy as np
+
+##########################################################
+namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm1_2.txt"
+namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm2_2.txt"
+##########################################################
+namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm1_2.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_2.txt"
+##########################################################
+#namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_2.txt"
+##########################################################
+
+# load data
+data = np.loadtxt(namefile,delimiter=";")
+t = data[:,0] ; n = data[:,1] ; s = data[:,2] ; d = data[:,3]
+
+# remove size and drop point when no vortex detected
+d[np.where(n==0)] = np.nan
+s[np.where(n==0)] = np.nan
+
+## PLOTS
+
+number = ppplot.plot1d()
+number.f = n
+number.x = t
+number.linestyle = ''
+number.marker = '.'
+number.color = 'b'
+number.xlabel = "Local time (hour)"
+number.ylabel = "Detected vortices"
+number.makeshow()
+
+drop = ppplot.plot1d()
+drop.f = d
+drop.x = t
+drop.linestyle = ''
+drop.marker = '.'
+drop.color = 'r'
+drop.fmt = "%.1f"
+drop.xlabel = "Local time (hour)"
+drop.ylabel = "Maximum drop of detected vortices (Pa)"
+drop.makeshow()
+
+size = ppplot.plot1d()
+size.f = s
+size.x = t
+size.linestyle = ''
+size.marker = '.'
+size.color = 'g'
+size.xlabel = "Local time (hour)"
+size.ylabel = "Maximum size of detected vortices (m)"
+size.makeshow()
Index: trunk/UTIL/PYTHON/powerlaw/dropsize.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/dropsize.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/dropsize.py	(revision 1053)
@@ -0,0 +1,32 @@
+#! /usr/bin/env python
+import matplotlib.pyplot as mpl
+import ppplot
+import numpy as np
+
+##########################################################
+namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm1_1.txt"
+namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm2_1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm1_1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm2_1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/press_ustm_exomars.nc1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_1.txt"
+##########################################################
+
+# read data np.loadtxt?
+## import pylab as plb data = plb.loadtxt('data.dat') x = data[:,0] y= data[:,1]
+f = open(namefile,'r')
+t = [] ; s = [] ; d = [] ; i = [] ; j = []
+for line in f:
+    ct,cs,cd,ci,cj = line.strip().split(';')
+    t.append(float(ct))
+    s.append(float(cs)) ; d.append(float(cd))
+    i.append(int(ci)) ; j.append(int(cj))
+f.close()
+# ensure numpy array
+t = np.array(t) ; s = np.array(s) ; d = np.array(d) ; i = np.array(i) ; j = np.array(j)
+
+drop = ppplot.plot1d() ; drop.field = d ; drop.absc = s
+drop.lstyle = '' ; drop.marker = '.' ; drop.color = 'r' ; drop.fmt = "%.1f"
+drop.xlabel = "Vortex size (m)" ; drop.ylabel = "Pressure drop (Pa)"
+drop.make() ; ppplot.save()
+
Index: trunk/UTIL/PYTHON/powerlaw/find_devils.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/find_devils.py	(revision 946)
+++ 	(revision )
@@ -1,329 +1,0 @@
-#! /usr/bin/env python
-
-def detsize( xx, res=1, thres=3, loga=False ):
-    import numpy as np
-    import math
-    size = []
-    sizecalc = 1
-    diff = np.asarray( np.roll(xx,-1) - xx )
-    for i in diff:
-        if abs(i) > 1:
-            if sizecalc >= thres: 
-                if loga: addthis = math.log(sizecalc*res)
-                else:    addthis = sizecalc*res
-                size.append(addthis)
-            sizecalc = 1
-        else:
-            sizecalc += 1
-    return size
-
-def getsize(filename):
-
-    import numpy as np
-    from scipy.ndimage.measurements import minimum_position
-    from netCDF4 import Dataset
-    import matplotlib.pyplot as plt
-    import myplot as myp
-   
-    ### LOAD NETCDF DATA
-    nc = Dataset(filename) 
-    psfc = nc.variables["PSFC"]
-    
-    ### LOOP on TIME
-    ### NB: a same event could be counted several times...
-    shape = np.array(psfc).shape
-    allsizesx = []
-    allsizesy = []
-    depression = []
-    stride = 1 #5
-    stride = 20
-    #stride = 50
-    #stride = 100
-    start = 0
-    start = stride
-    for i in range(start,shape[0],stride):
-
-        psfc2d = np.array ( psfc [ i, : , : ] )
-    
-        ############### CRITERION
-        ave = np.mean(psfc2d,dtype=np.float64)  ## dtype otherwise inaccuracy
-
-        #limdp = -0.2 ## on en loupe pas mal
-        #where = np.where(psfc2d - ave < limdp)  ## comme le papier Phoenix
-    
-        std = np.std(psfc2d,dtype=np.float64)   ## dtype otherwise inaccuracy
-        fac = 4.   ## how many sigmas. not too low, otherwise vortices are not caught. 4 good choice.
-                   ## 2.5 clearly too low, 3.5 not too bad, 4 probably good
-        fac = 3.5
-        fac = 3.2
-        ##fac = 2.5
-        #fac = 3. ## final choice
-        #fac = 2.5
-        lim = ave - fac*std
-        where = np.where(psfc2d < lim)
-        ############### END CRITERION
-
-        depression = np.append(depression,np.ravel(psfc2d[where])-ave)
-   
-        lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine
-        lab[where] = 1.  ## do not treat points close to 'mean' (background) pressure
-
-        draw = False
-        #draw = True
-        if draw:
-        ##################################################################################
-            vmin = -0.3
-            vmax =  0.0
-            ndiv = 3
-            palette = plt.get_cmap(name="YlGnBu") 
-            what_I_plot = psfc2d-ave
-            zevmin, zevmax = myp.calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
-            what_I_plot = myp.bounds(what_I_plot,zevmin,zevmax)
-            zelevels = np.linspace(zevmin,zevmax)
-            fig = plt.figure(figsize=(16,8)) 
-            subv,subh = myp.definesubplot(2,fig)
-            plt.subplot(subv,subh,1) 
-            plt.contourf(what_I_plot,zelevels,cmap=palette) 
-            plt.colorbar(fraction=0.05,pad=0.03,format="%.1f",\
-                           ticks=np.linspace(zevmin,zevmax,ndiv+1),\
-                           extend='both',spacing='proportional')
-            plt.subplot(subv,subh,2)
-            palette = plt.get_cmap(name="hot") #"binary")
-            #lab[np.where(lab > 0.99)] = np.NaN
-            #plt.pcolor(lab,cmap=palette)
-            plt.contourf(lab,2,cmap=palette)
-            plt.show() 
-        ##################################################################################
-  
-        xx = []
-        yy = []
-        while 1 in lab:
-            p = minimum_position(psfc2d,labels=lab)
-            lab[p] = 0 ## once a minimum has been found in a grid point, do not search here again.
-            if p[0] not in xx: xx.append(p[0]) ## if x coordinate not yet in the list add it
-            if p[1] not in yy: yy.append(p[1]) ## if y coordinate not yet in the list add it
-        xx.sort()
-        yy.sort()
-        ### now xx and yy are sorted arrays containing grid points with pressure minimum
-       
-        ######## DETERMINE SIZE OF STRUCTURES
-        ######## this is rather brute-force...
-        sizex = detsize( xx, res = 10, loga=False, thres=2 )
-        sizey = detsize( yy, res = 10, loga=False, thres=2 )
-        #sizex = detsize( xx, res = 10, loga=False, thres=3 )
-        #sizey = detsize( yy, res = 10, loga=False, thres=3 )
-        sizex = detsize( xx, res = 15, loga=False, thres=2 )
-        sizey = detsize( yy, res = 15, loga=False, thres=2 )
-        ###
-        print sizex, sizey
-        #if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex  ### un peu limite dans certains cas
-        if (len(sizex) > len(sizey))     : sizey = sizex        ### plus fidele mais petit souci lorsque PBC
-        elif (len(sizex) == len(sizey))  : 
-          if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex
-          else                                  : sizex = sizey
-        else                            : sizex = sizey
-        allsizesx = np.append(allsizesx,sizex)
-        allsizesy = np.append(allsizesy,sizey)
-        print i, ' on ', shape[0], ' caught ', len(sizex), ' vortices ', sizex
-        ########
-  
-    allsizesx.sort()
-    allsizesy.sort()
-    
-    return allsizesx, allsizesy, depression
-
-#########################################################################
-#########################################################################
-
-import matplotlib.pyplot as plt
-import pickle
-import numpy as np
-import matplotlib.mlab as mlab
-import mymath as mym
-import myplot as myp
-
-import plfit
-import plplot
-import randht
-import plpva
-
-save = True
-#save = False
-pression = False
-#pression = True
-
-filename = "/home/aymeric/Big_Data/psfc_f18.nc"
-
-if save:
-    allsizesx, allsizesy, depression = getsize(filename)
-    ### sauvegarde texte pour inspection
-    mym.writeascii(allsizesx,'allsizex.txt')
-    mym.writeascii(allsizesy,'allsizey.txt')
-    mym.writeascii(depression,'alldepression.txt')
-    ### sauvegarde binaire pour utilisation python
-    myfile = open('allsizex.bin', 'wb')
-    pickle.dump(allsizesx, myfile)
-    myfile.close()
-    myfile = open('allsizey.bin', 'wb')
-    pickle.dump(allsizesy, myfile)
-    myfile.close()
-    myfile = open('alldepression.bin', 'wb')
-    pickle.dump(depression, myfile)
-    myfile.close()
-
-### load files
-myfile = open('allsizex.bin', 'r')
-allsizesx = pickle.load(myfile)
-myfile = open('allsizey.bin', 'r')
-allsizesy = pickle.load(myfile)
-myfile = open('alldepression.bin', 'r')
-depression = pickle.load(myfile)
-depression = np.array(abs(depression))#*1000.
-
-### sizes
-#plothist = np.append(allsizesx,allsizesy)
-plothist = allsizesx
-if pression: plothist = depression
-plothist.sort()
-print 'mean ', np.mean(plothist,dtype=np.float64)  
-print 'std ', np.std(plothist,dtype=np.float64)
-print 'max ', np.max(plothist)
-print 'min ', np.min(plothist)
-print 'len ', len(plothist)
-
-
-### MAKE BINS
-nbins = 100
-zebins = [2.0]
-#nbins = 8
-#zebins = [19.0]
-#nbins = 15
-#zebins = [11.] ##12 non mais donne un peu la meme chose
-#zebins = [20.] ##20 non car trop pres du premier
-nbins = 100
-zebins = [2./np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
-nbins = 200
-
-if pression: zebins = [0.3]
-
-for i in range(0,nbins):  zebins.append(zebins[i]*np.sqrt(2))
-zebins = np.array(zebins)
-#### select reasonable bins for DD
-if not pression:
-    zebins = zebins [ zebins > 15. ] 
-    #zebins = zebins [ zebins > 20. ]
-    zebins = zebins [ zebins > 25. ]
-    zebins = zebins [ zebins < 1000. ]
-else:
-    zebins = zebins [ zebins < 10. ]
-print 'corrected bins ',zebins
-
-#### HISTOGRAM
-plt.figure(1)
-plt.hist(    plothist,\
-             log=True,\
-             bins=zebins,\
-#             cumulative=-1,\
-             normed=True,\
-        )
-plt.xscale('log')
-if pression:    plt.xlabel('Pressure (Pa)')
-else:           plt.xlabel('Diameter (m)')
-plt.ylabel('Population (normalized)')
-if pression: prefix="p"
-else:        prefix=""
-myp.makeplotres(prefix+"histogram",res=200,disp=False)
-plt.close(1)
-
-### COMPARED HISTOGRAMS
-### --- FIT WITH POWER LAW
-if pression: [alpha, xmin, L] = plfit.plfit(plothist,'xmin',0.3)
-else:        [alpha, xmin, L] = plfit.plfit(plothist,'limit',20.)
-print alpha,xmin
-
-#a = plpva.plpva(plothist,0.75,'xmin',0.75)
-#print a
-
-#### DEUXIEME ROUTINE
-####IL FAUT UTILISER LE DISCRET POUR LA TAILLE !!!
-#if pression:   myplfit = plfit.plfit(plothist,verbose=True,xmin=0.75)
-#else:          myplfit = plfit.plfit(plothist,verbose=True,xmin=20.)
-#myplfit.plotppf()
-#plt.show()
-#exit()
-
-
-#plt.figure(1)
-#h = plplot.plplot(plothist,xmin,alpha)
-#myp.makeplotres(prefix+"fit",res=200,disp=False)
-plt.figure(2)
-### --- POWER LAW (factor does not really matter)
-power = (xmin/2.2)*np.array(randht.randht(10000,'powerlaw',alpha))
-#power = (xmin/2.2)*np.array(randht.randht(10000,'cutoff',alpha,10.)) ##marche pas si trop grand
-print 'mean ', np.mean(power,dtype=np.float64)
-### --- EXPONENTIAL LAW
-expo = randht.randht(10000,'exponential',1./(np.mean(power,dtype=np.float64)*1.00))
-print 'mean ', np.mean(expo,dtype=np.float64)
-### --- PLOT
-plt.hist(    [plothist,power,expo],\
-             label=['LES vortices','Power law '+'{:.1f}'.format(alpha),'Exponential law'],\
-             log=True,\
-             bins=zebins,\
-#            cumulative=-1,\
-             normed=True,\
-        )
-plt.legend()
-plt.xscale('log')
-if pression:    plt.xlabel('Pressure (Pa)')
-else:           plt.xlabel('Diameter (m)')
-plt.ylabel('Population (normalized)')
-myp.makeplotres(prefix+"comparison",res=200,disp=False)
-plt.close(2)
-
-########################
-########################
-zebins = [30.,42.,60.,84.,120.,170.,240.,340.]
-plothist = []
-plothist = np.append(plothist,30 *np.ones(306))
-plothist = np.append(plothist,42 *np.ones(58) )
-plothist = np.append(plothist,60 *np.ones(66) )
-plothist = np.append(plothist,84 *np.ones(41) )
-plothist = np.append(plothist,120*np.ones(19) )
-plothist = np.append(plothist,170*np.ones(9)  )
-plothist = np.append(plothist,240*np.ones(2)  )
-plothist = np.append(plothist,340*np.ones(1)  )
-
-#zebins = [50.,71.,100.,141.,200.,282.,400.]
-#plothist = []
-#plothist = np.append(plothist,50. *np.ones(36))
-#plothist = np.append(plothist,71. *np.ones(18) )
-#plothist = np.append(plothist,100. *np.ones(12) )
-#plothist = np.append(plothist,141. *np.ones(6) )
-#plothist = np.append(plothist,200.*np.ones(4) )
-#plothist = np.append(plothist,282.*np.ones(1)  )
-#plothist = np.append(plothist,400.*np.ones(2)  )
-
-exit()
-
-plt.figure(3)
-[alpha, xmin, L] = plfit.plfit(plothist,'xmin',30)#50.)
-print alpha,xmin
-#a = plpva.plpva(plothist,30,'xmin',30)
-h = plplot.plplot(plothist,xmin,alpha)
-plt.loglog(h[0], h[1], 'k--',linewidth=2)
-plt.hist(    plothist,\
-             log=True,\
-             bins=zebins,\
-#             cumulative=-1,\
-             normed=True,\
-        )
-plt.xscale('log')
-plt.xlabel('Pressure (micro Pa)')
-plt.ylabel('Population (normalized)')
-myp.makeplotres("data",res=200,disp=False)
-
-#plt.figure(4)
-#[alpha, xmin, L] = plfit.plfit(plothist,'xmin',50.) #,'xmin',30.)
-#print alpha,xmin
-#h = plplot.plplot(plothist,xmin,alpha)
-#myp.makeplotres("datafit",res=200,disp=False)
Index: trunk/UTIL/PYTHON/powerlaw/finddd.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/finddd.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/finddd.py	(revision 1053)
@@ -0,0 +1,340 @@
+#! /usr/bin/env python
+from ppclass import pp
+from ppplot import plot2d,save,writeascii
+import numpy as np
+from scipy.ndimage.measurements import minimum_position
+import matplotlib.pyplot as mpl
+import ppcompute
+
+###############################################################################
+# FIND DUST DEVILS
+# filefile --> file
+# dx --> horizontal resolution
+# dt_out --> frequency of outputs (s) 
+# lt_start --> starting local time
+# halolim --> limit for halos
+# method --> method used to detect vortices
+# plotplot --> plot or not
+# save --> save results or not
+###############################################################################
+def finddd(filefile,\
+           timelist,\
+           dx=50.,\
+           dt_out=50.,\
+           lt_start=8.,\
+           halolim=30.,\
+           method=1,\
+           plotplot=False,\
+           save=True):
+
+    ###############################################################################
+    ########################## FOR METHOD 1 FOR METHOD 1 ##########################
+    ## FACLIST : how many sigmas below mean we start to consider this could be a vortex
+    ## ... < 3.0 dubious low-intensity minima are caught
+    ## ... 3 ~ 3.1 is a bit too low but helps? because size of pressure drop could be an underestimate of actual dust devil
+    ## ... 3.2 ~ 3.3 is probably right, this is the one we choose for exploration [3.25]
+    ## ... NB: values 3.2 to 3.6 yields the same number of devils but measured sizes could vary a bit
+    ## ... > 3.6 makes strong minima disappear... especially >3.8
+    ## ... EVENTUALLY 3.5-4 is better for problematic cases 
+    ###############################################################################
+    faclist = [3.75]
+
+    ## (see below) NEIGHBOR_FAC is the multiple of std used to evaluate size
+    ## --> 1: limit not discriminative enough. plus does not separate neighbouring vortices.
+    ##        ... but interesting: gives an exponential law (because vortices are artificially merged?)
+    ## --> 2.7: very good for method 1. corresponds usually to ~0.3
+    ## --> 2: so-so. do not know what to think.
+    neighbor_fac = 2.7
+
+    ###############################################################################
+    if save:
+        myfile1 = open(filefile+'m'+str(method)+'_'+'1.txt', 'w')
+        myfile2 = open(filefile+'m'+str(method)+'_'+'2.txt', 'w')
+    ###############################################################################
+    
+    ## mean and std calculations
+    ## -- std is used in both methods for limits
+    ## -- mean is only used in method 1
+    print "calculate mean and std, please wait."
+    ## -- get time series of 2D surface pressure
+    psfc = pp(file=filefile,var="PSFC").getf()
+    ## -- calculate mean and standard deviation
+    ## -- ... calculating std at all time is not right!
+    ## -- ... for mean value though, similar results with both methods
+    mean = np.mean(psfc,dtype=np.float64)
+    std = np.std(psfc,dtype=np.float64)
+    ## some information about inferred limits
+    print "**************************************************************"
+    print "MEAN",mean
+    print "STD",std
+    print "LIMIT FOR PRESSURE MINIMUM:",-np.array(faclist)*std
+    print "LIMIT FOR EVALUATING SIZE OF A GIVEN LOCAL MINIMUM",-neighbor_fac*std
+    print "**************************************************************"
+    
+    ## LOOP ON TIME
+    for time in timelist:
+    
+     ## get 2D surface pressure at a given time
+     ## (this is actually so quick we don't use psfc above)
+     psfc2d = pp(file=filefile,var="PSFC",t=time).getf()
+     #ustm = pp(file=filefile,var="USTM",t=time).getf()
+
+     ## MAIN ANALYSIS. LOOP ON FAC.
+     for fac in faclist:
+     #fac = 3.75
+     #for method in [2,1]:
+  
+      ## initialize arrays
+      tabij = [] ; tabsize = [] ; tabdrop = []
+      tabijcenter = [] ; tabijvortex = [] ; tabijnotconv = [] ; tabdim = []
+    
+      ################ FIND RELEVANT POINTS TO BE ANALYZED
+      ## lab is 1 for points to be treated by minimum_position routine
+      ## we set elements at 1 where pressure is under mean-fac*std
+      ## otherwise we set to 0 because this means we are close enough to mean pressure (background)
+      lab = np.zeros(psfc2d.shape)
+      if method == 1:
+          # method 1: standard deviation
+          lab[np.where(psfc2d < mean-fac*std)] = 1
+      else:
+          # method 2: polynomial fit
+          # ... tried smooth but too difficult, not accurate and too expensive
+          deg = 5 #plutot bien (deg 10 ~pareil) mais loupe les gros (~conv cell)
+          #deg = 2 #pas mal mais false positive (pareil 3-4 mm si un peu mieux)
+          #        #(OK now with fixing the false positive bug)
+          nx = psfc2d.shape[1] ; ny = psfc2d.shape[0]
+          xxx = np.array(range(nx)) ; yyy = np.array(range(ny))
+          anopsfc2d = psfc2d*0. ; polypsfc2d = psfc2d*0.
+          for iii in range(0,nx,1):
+             poly = np.poly1d(np.polyfit(yyy,psfc2d[iii,:],deg))
+             polypsfc2d[iii,:] = poly(yyy)
+          for jjj in range(0,ny,1):
+             poly = np.poly1d(np.polyfit(xxx,psfc2d[:,jjj],deg))
+             polypsfc2d[:,jjj] = 0.5*polypsfc2d[:,jjj] + 0.5*poly(xxx)
+          ## smooth a little to avoid 'crosses' (plus, this removes PBC problems)
+          polypsfc2d = ppcompute.smooth2diter(polypsfc2d,n=deg)
+          # compute anomaly and find points to be explored
+          anopsfc2d = psfc2d - polypsfc2d
+          limlim = fac*std ## same as method 1
+          lab[np.where(anopsfc2d < -limlim)] = 1
+    
+      ## while there are still points to be analyzed...
+      while 1 in lab:
+        ## ... get the point with the minimum field values
+        if method == 1:
+            ij = minimum_position(psfc2d,labels=lab)
+        else:
+            ij = minimum_position(anopsfc2d,labels=lab)
+        ## ... store the indexes of the point in tabij
+        tabij.append(ij)
+        ## ... remove the point from labels to be further explored by minimum_position
+        lab[ij] = 0
+    
+      ################ GET SIZES BASED ON THOSE FOUND POINTS
+      ## reslab is the same as lab
+      ## except for scanning purpose we keep the information
+      ## about how went the detection
+      ## --> and we set a lower fac
+      ## --> above a high fac is good to catch only strong vortices
+      ## --> but here a casual fac=3 is better to get accurate sizes
+      ## --> or even lower as shown by plotting reslab 
+      reslab = np.zeros(psfc2d.shape)
+      if method == 1:
+          reslab[np.where(psfc2d < mean-neighbor_fac*std)] = 1
+      else:
+          reslab[np.where(anopsfc2d < -neighbor_fac*std)] = 1
+     
+      ## initialize halomax and while loop
+      ## HALOMAX : maximum halo defined around a minima to evaluate the size
+      ## ... halomax must be large enough to encompass a vortex
+      ## ... but not too large otherwise neighboring vortex are caught
+      halomax = 3 ; notconv = 9999 ; yorgl = 9999
+    
+      ## WHILE LOOP on HALOMAX exploration
+      while ( notconv > 0 and halomax < halolim and yorgl != 0 ):
+       # now browse through all points caught in previous loop
+       for ij in tabij:
+        ## ... OK. take each indexes found before with minimum_position
+        i,j = ij[0],ij[1]
+        ## EITHER
+        ## ... if reslab is already 0, we do not have to do anything
+        ## ... because this means point is already part of detected vortex
+        ## OR
+        ## ... if the ij couple is already in a previously detected vortex
+        ## ... we don't actually need to consider it again
+        ## ... this is necessary otherwise (sometimes a lot) of false positives
+        if reslab[i,j] <= 0 or ij in tabijvortex:
+          pass
+        else:
+          ## ... then define a growing halo around this point
+          ## ... we make the halo grow until convergence of 
+          ## ... the number of under-limit points (i.e. with lab=1) within halo
+          ## ... which means the whole vortex is encompassed
+          ## ... we start with a halo of 1 point around minimum point
+          ## ... and we end when we reached halomax which likely means dubious case
+          halo = 1 ; nmesh = -9999. ; prevmesh = 9999. ; notconverged = False
+          while nmesh != prevmesh and halo <= halomax:
+              ## ... save the number of vortex points calculated at previous iteration
+              prevmesh = nmesh
+              ## ... define a halo around the minimum point
+              minx,maxx,miny,maxy = i-halo,i+halo+1,j-halo,j+halo+1
+              ## ... treat the boundary case (TBD: periodic boundary conditions)
+              if minx < 0: minx = 0
+              if miny < 0: miny = 0
+              if maxx > psfc2d.shape[0]: maxx = psfc2d.shape[0]
+              if maxy > psfc2d.shape[1]: maxy = psfc2d.shape[1]
+              ## ... define the patch, made of the halo of points
+              patch = reslab[minx:maxx,miny:maxy]       
+              ## ... count how many 1 are inside this patch
+              ## ... these are points for which value is >0
+              ## ... because close to mean is 0 and already caught is <0
+              ## ... and not converged are >1 so still to be scanned
+              nmesh = len(np.where(patch >= 1)[0])
+              ## ... in case widening halo adds 1-2 points only
+              ## ... we consider vortex is encompassed
+              ## ... this helps to separate neighboring vortices
+              ## ... and only yields marginal inaccuracies
+              if halo > 3 and abs(nmesh-prevmesh)<=2: prevmesh = nmesh
+              ## ... if with halo=1 we caught a 1-point vortex, end here because spurious
+              if halo == 1 and nmesh == 1: prevmesh = nmesh
+              ## ... with the last halo=halomax test, we can detect not-converged cases
+              if halo == halomax and nmesh != prevmesh: notconverged = True
+              ## ... increment halo size before the next iteration
+              halo = halo + 1
+          ## ... now we got the grid points encompassed by the vortex in nmesh
+          ## ... then to continue scanning we store results in reslab
+          if notconverged:
+            ## multiply reslab by 2. 
+            ## --> if it is >1, point could be part of another vortex found later
+            ## --> if it is <0, point is already within a caught vortex, so reslab should remain <0
+            ## --> if it is =0, this should remain 0 because close to average
+            reslab[i,j] = reslab[i,j]*2
+            tabijnotconv.append(ij)
+          else:
+            ## OK. this is most likely an actual vortex. we get the drop.
+            ## we multiply by mesh area, then square to get approx. size of vortex
+            if method == 1:
+                drop = -psfc2d[i,j]+mean
+            else:
+                drop = -anopsfc2d[i,j]
+            size = int(np.sqrt(nmesh*dx*dx))
+            ## if vortex is too small (i.e. too close to mesh grid resolution) we don't store information
+            ## ... we just remove the patch from labels to be further explored
+            facdx = 2.
+            if size < facdx*dx:
+                reslab[minx:maxx,miny:maxy] = 0
+            ## otherwise it is a VORTEX! we store info in arrays
+            else:
+                ## we put the vortex points into tabijvortex so that 
+                ## a (i,j) couple into a vortex is not considered in further explorations
+                ## -- also we evaluate the x-size and y-size of the vortex (max distance to center)
+                ## -- this can be useful to detect not-so-round fake vortices (convective gusts?)
+                maxw = -9999 ; maxh = -9999 ; maxu = -9999
+                for iii in range(minx,maxx):
+                 for jjj in range(miny,maxy):
+                  if reslab[iii,jjj] != 0:
+                     tabijvortex.append((iii,jjj))
+                     width = np.abs(iii-ij[0])
+                     if width > maxw: maxw = width
+                     height = np.abs(jjj-ij[1])
+                     if height > maxh: maxh = height
+                     #ustar = ustm[iii,jjj]
+                     #if ustar > maxu: maxu = ustar
+                #print size,drop,ustar
+                ## store info in dedicated arrays
+                tabdim.append((maxw*dx,maxh*dx))
+                tabsize.append(size)
+                tabdrop.append(drop)
+                tabijcenter.append(ij)
+                #print "... VORTEX!!!! size %.0f drop %.1f coord %.0f %.0f" % (size,drop,i,j)
+                ## we remove the patch from labels to be further explored
+                ## we multiply reslab by -1 to plot detection maps
+                reslab[minx:maxx,miny:maxy] = patch*-1
+     
+       ## count how many points are not converged and left to be analyzed
+       notconv = len(np.where(reslab > 1)[0])
+       yorgl = len(np.where(reslab == 1)[0])
+    
+       ## increment halomax
+       ## to speed-up the increment is slightly increasing with considered halomax
+       halomax = halomax + halomax / 2
+    
+      ## just for simpler plots.
+      reslab[reslab > 2] = 4
+      reslab[reslab < -2] = -4
+
+      ## give some info to the user
+      if len(tabsize) > 0:
+        nvortex = len(tabsize)
+        maxsize = np.max(tabsize)
+        maxdrop = np.max(tabdrop)
+      else:
+        nvortex = 0
+        maxsize = 0
+        maxdrop = 0.
+      notconv = len(np.where(reslab > 1)[0])
+      print "t=%3.0f / n=%2.0f / s_max=%4.0f / d_max=%4.1f / halo_out=%3.0f / notconvp=%3.1f" \
+            % (time,nvortex,maxsize,maxdrop,halomax,100.*notconv/float(reslab.size))        
+    
+      ## save results in a text file
+      if save:
+          # convert t in local time
+          ttt = lt_start + time*dt_out/3700.      
+          # write files
+          myfile2.write( "%5.2f ; %5.0f ; %5.0f ; %7.2f\n" % (ttt,nvortex,maxsize,maxdrop) )
+          for iii in range(len(tabsize)):
+              myfile1.write( "%5.2f ; %5.0f ; %7.2f ; %5.0f ; %5.0f\n" \
+              % (ttt,tabsize[iii],tabdrop[iii],tabdim[iii][0],tabdim[iii][1]) )
+    
+      #### PLOT PLOT PLOT PLOT
+      if nvortex>0 and plotplot:
+       mpl.figure(figsize=(12,8))
+       myplot = plot2d()
+       myplot.absc = np.array(range(psfc2d.shape[1]))*dx/1000.
+       myplot.ordi = np.array(range(psfc2d.shape[0]))*dx/1000.
+       myplot.title = str(nvortex)+" vortices found (indicated diameter / pressure drop)"
+       myplot.xlabel = "x distance (km)"
+       myplot.ylabel = "y distance (km)"
+       if method > 0:
+       #if method == 1:
+           #myplot.field = ustm 
+           myplot.field = psfc2d
+           #myplot.field = polypsfc2d
+           myplot.vmin = -2.*std 
+           myplot.vmax = +2.*std
+           myplot.vmin = mean - 6.*std
+           myplot.vmax = mean + 6.*std
+       else:
+           myplot.field = anopsfc2d
+           myplot.vmin = -1.5
+           myplot.vmax = 0.5
+       myplot.fmt = "%.1f"
+       myplot.div = 20
+       myplot.colorb = "spectral"
+       myplot.make()
+      
+       ### ANNOTATIONS
+       for iii in range(len(tabsize)):
+        ij = tabijcenter[iii]
+        coord1 = ij[1]*dx/1000.
+        coord2 = ij[0]*dx/1000.
+        txt = "%.0f/%.1f" % (tabsize[iii],tabdrop[iii])
+        mpl.annotate(txt,xy=(coord1,coord2),
+             xytext=(-10,-30),textcoords='offset points',ha='center',va='bottom',\
+             bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.3),\
+             arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.3',color='red'),\
+             size='small')
+    
+       ###show detection
+       #lev = [-4,-2,-1,0,1,2,4]
+       #lev = [-1,0]
+       #mpl.contourf(myplot.absc,myplot.ordi,reslab,alpha=0.9,cmap=mpl.cm.get_cmap("binary_r"),levels=lev)
+    
+       ### SHOW OR SAVE IN FILE
+       mpl.show()
+       #save(mode="png",filename="detectm"+"_"+str(time)+"_"+str(method),folder="detect/",includedate=False)
+    
+    ## close data files
+    if save:
+        myfile1.close()
+        myfile2.close()
Index: trunk/UTIL/PYTHON/powerlaw/performdd.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/performdd.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/performdd.py	(revision 1053)
@@ -0,0 +1,59 @@
+#! /usr/bin/env python
+from finddd import finddd
+
+####################
+####################
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",range(0,591,1))
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",range(0,591,1),method=2)
+####################
+####################
+
+
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",range(0,591,10),method=2)
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",range(0,591,10),method=2,plotplot=True)
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",[300.,350.,400.],method=2,plotplot=True)
+#### checks
+#plotplot = True ; save = False
+#timelist = [270.,302.,305.,352.,358.,412.,423.] #joli
+#timelist = [300.,310.,320.]
+#timelist = [210.,250.,280.,400.]
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",range(350,450,5),method=2,plotplot=True,save=False)
+
+####################
+####################
+finddd("/home/aymeric/Big_Data/LES_dd/press_ustm_exomars.nc",range(0,54,1),dx=12.,halolim=100,method=2)
+#timelist = [28.,48.,49.]
+#plotplot = True
+#save = False
+####################
+####################
+
+####################
+####################
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.nc",\
+#           range(0,336,1),dx=10.,lt_start=9.,halolim=150,method=2)
+###halolim 50 idem
+####################
+####################
+
+
+
+###############################################################################
+#dx = 100. # horizontal resolution
+#lt_out = 10. # frequency of outputs (s) --- wrong?????
+#lt_start = 6. # starting local time 
+#filefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.nc"
+#halolim = 20 # limit for halos
+#timelist = range(0,4440,1) # times in file
+#timelist = range(0,4440,100) # times in file
+#method = 1 # method used to detect vortices
+#plotplot = False # plot or not
+#save = True # save results or not
+###############################################################################
+
+
+##namefile = "/home/aymeric/Big_Data/LES_dd/press_ustm_exomars.nc2.txt"
+##dt_out = 50.
+##lt_start = 8.
+
+
Index: trunk/UTIL/PYTHON/powerlaw/plfit.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/plfit.py	(revision 946)
+++ 	(revision )
@@ -1,429 +1,0 @@
-from math import *
-
-# function [alpha, xmin, L]=plfit(x, varargin)
-# PLFIT fits a power-law distributional model to data.
-#    Source: http://www.santafe.edu/~aaronc/powerlaws/
-# 
-#    PLFIT(x) estimates x_min and alpha according to the goodness-of-fit
-#    based method described in Clauset, Shalizi, Newman (2007). x is a 
-#    vector of observations of some quantity to which we wish to fit the 
-#    power-law distribution p(x) ~  x^-alpha for x >= xmin.
-#    PLFIT automatically detects whether x is composed of real or integer
-#    values, and applies the appropriate method. For discrete data, if
-#    min(x) > 1000, PLFIT uses the continuous approximation, which is 
-#    a reliable in this regime.
-#   
-#    The fitting procedure works as follows:
-#    1) For each possible choice of x_min, we estimate alpha via the 
-#       method of maximum likelihood, and calculate the Kolmogorov-Smirnov
-#       goodness-of-fit statistic D.
-#    2) We then select as our estimate of x_min, the value that gives the
-#       minimum value D over all values of x_min.
-#
-#    Note that this procedure gives no estimate of the uncertainty of the 
-#    fitted parameters, nor of the validity of the fit.
-#
-#    Example:
-#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
-#       [alpha, xmin, L] = plfit(x)
-#   or  a = plfit(x)
-#
-#    The output 'alpha' is the maximum likelihood estimate of the scaling
-#    exponent, 'xmin' is the estimate of the lower bound of the power-law
-#    behavior, and L is the log-likelihood of the data x>=xmin under the
-#    fitted power law.
-#    
-#    For more information, try 'type plfit'
-#
-#    See also PLVAR, PLPVA
-
-# Version 1.0.10 (2010 January)
-# Copyright (C) 2008-2011 Aaron Clauset (Santa Fe Institute)
-
-# Ported to Python by Joel Ornstein (2011 July)
-# (joel_ornstein@hmc.edu)
-
-# Distributed under GPL 2.0
-# http://www.gnu.org/copyleft/gpl.html
-# PLFIT comes with ABSOLUTELY NO WARRANTY
-#
-# 
-# The 'zeta' helper function is modified from the open-source library 'mpmath'
-#   mpmath: a Python library for arbitrary-precision floating-point arithmetic
-#   http://code.google.com/p/mpmath/
-#   version 0.17 (February 2011) by Fredrik Johansson and others
-# 
-
-# Notes:
-# 
-# 1. In order to implement the integer-based methods in Matlab, the numeric
-#    maximization of the log-likelihood function was used. This requires
-#    that we specify the range of scaling parameters considered. We set
-#    this range to be 1.50 to 3.50 at 0.01 intervals by default. 
-#    This range can be set by the user like so,
-#    
-#       a = plfit(x,'range',[1.50,3.50,0.01])
-#    
-# 2. PLFIT can be told to limit the range of values considered as estimates
-#    for xmin in three ways. First, it can be instructed to sample these
-#    possible values like so,
-#    
-#       a = plfit(x,'sample',100)
-#    
-#    which uses 100 uniformly distributed values on the sorted list of
-#    unique values in the data set. Second, it can simply omit all
-#    candidates above a hard limit, like so
-#    
-#       a = plfit(x,'limit',3.4)
-#    
-#    Finally, it can be forced to use a fixed value, like so
-#    
-#       a = plfit(x,'xmin',3.4)
-#    
-#    In the case of discrete data, it rounds the limit to the nearest
-#    integer.
-# 
-# 3. When the input sample size is small (e.g., < 100), the continuous 
-#    estimator is slightly biased (toward larger values of alpha). To
-#    explicitly use an experimental finite-size correction, call PLFIT like
-#    so
-#    
-#       a = plfit(x,'finite')
-#    
-#    which does a small-size correction to alpha.
-#
-# 4. For continuous data, PLFIT can return erroneously large estimates of 
-#    alpha when xmin is so large that the number of obs x >= xmin is very 
-#    small. To prevent this, we can truncate the search over xmin values 
-#    before the finite-size bias becomes significant by calling PLFIT as
-#    
-#       a = plfit(x,'nosmall')
-#    
-#    which skips values xmin with finite size bias > 0.1.
-
-def plfit(x, *varargin):
-    vec     = []
-    sample  = []
-    xminx   = []
-    limit   = []
-    finite  = False
-    nosmall = False
-    nowarn  = False
-
-    # parse command-line parameters trap for bad input
-    i=0 
-    while i<len(varargin): 
-        argok = 1 
-        if type(varargin[i])==str: 
-            if varargin[i]=='range':
-                Range = varargin[i+1]
-                if Range[1]>Range[0]:
-                    argok=0
-                    vec=[]
-                try:
-                    vec=map(lambda X:X*float(Range[2])+Range[0],\
-                            range(int((Range[1]-Range[0])/Range[2])))
-                    
-                    
-                except:
-                    argok=0
-                    vec=[]
-                    
-
-                if Range[0]>=Range[1]:
-                    argok=0
-                    vec=[]
-                    i-=1
-
-                i+=1
-                    
-                            
-            elif varargin[i]== 'sample':
-                sample  = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'limit':
-                limit   = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'xmin':
-                xminx   = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'finite':       finite  = True    
-            elif varargin[i]==  'nowarn':       nowarn  = True    
-            elif varargin[i]==  'nosmall':      nosmall = True    
-            else: argok=0 
-        
-      
-        if not argok:
-            print '(PLFIT) Ignoring invalid argument #',i+1 
-      
-        i = i+1 
-
-    if vec!=[] and (type(vec)!=list or min(vec)<=1):
-        print '(PLFIT) Error: ''range'' argument must contain a vector or minimum <= 1. using default.\n'                        
-              
-        vec = []
-    
-    if sample!=[] and sample<2:
-        print'(PLFIT) Error: ''sample'' argument must be a positive integer > 1. using default.\n'
-        sample = []
-    
-    if limit!=[] and limit<min(x):
-        print'(PLFIT) Error: ''limit'' argument must be a positive value >= 1. using default.\n'
-        limit = []
-    
-    if xminx!=[] and xminx>=max(x):
-        print'(PLFIT) Error: ''xmin'' argument must be a positive value < max(x). using default behavior.\n'
-        xminx = []
-    
-
-
-    # select method (discrete or continuous) for fitting
-    if     reduce(lambda X,Y:X==True and floor(Y)==float(Y),x,True): f_dattype = 'INTS'
-    elif reduce(lambda X,Y:X==True and (type(Y)==int or type(Y)==float or type(Y)==long),x,True):    f_dattype = 'REAL'
-    #else:                 f_dattype = 'UNKN'
-    ##AS
-    else:                 f_dattype = 'REAL'  
-
-    if f_dattype=='INTS' and min(x) > 1000 and len(x)>100:
-        f_dattype = 'REAL'
-   
-    print f_dattype
-
-    # estimate xmin and alpha, accordingly
-        
-    if f_dattype== 'REAL':
-        xmins = unique(x)
-        xmins.sort()
-        xmins = xmins[0:-1]
-        if xminx!=[]:
-            
-            xmins = [min(filter(lambda X: X>=xminx,xmins))]
-            
-        
-        if limit!=[]:
-            xmins=filter(lambda X: X<=limit,xmins)
-            if xmins==[]: xmins = [min(x)]
-            
-        if sample!=[]:
-            step = float(len(xmins))/(sample-1)
-            index_curr=0
-            new_xmins=[]
-            for i in range (0,sample):
-                if round(index_curr)==len(xmins): index_curr-=1
-                new_xmins.append(xmins[int(round(index_curr))])
-                index_curr+=step
-            xmins = unique(new_xmins)
-            xmins.sort()
-            
-            
-        
-        dat   = []
-        z     = sorted(x)
-        
-        for xm in range(0,len(xmins)):
-            xmin = xmins[xm]
-            z    = filter(lambda X:X>=xmin,z)
-            
-            n    = len(z)
-            # estimate alpha using direct MLE
-
-            a    = float(n) / sum(map(lambda X: log(float(X)/xmin),z))
-            if nosmall:
-                if (a-1)/sqrt(n) > 0.1 and dat!=[]:
-                    xm = len(xmins)+1
-                    break
-                
-            
-            # compute KS statistic
-            #cx   = map(lambda X:float(X)/n,range(0,n))
-            cf   = map(lambda X:1-pow((float(xmin)/X),a),z)
-            dat.append( max( map(lambda X: abs(cf[X]-float(X)/n),range(0,n))))
-        D     = min(dat)
-        xmin  = xmins[dat.index(D)]
-        z     = filter(lambda X:X>=xmin,x)
-        z.sort()
-        n     = len(z) 
-        alpha = 1 + n / sum(map(lambda X: log(float(X)/xmin),z))
-        if finite: alpha = alpha*float(n-1)/n+1./n  # finite-size correction
-        if n < 50 and not finite and not nowarn:
-            print '(PLFIT) Warning: finite-size bias may be present.\n'
-        
-        L = n*log((alpha-1)/xmin) - alpha*sum(map(lambda X: log(float(X)/xmin),z))
-    elif f_dattype== 'INTS':
-        
-        x=map(int,x)
-        if vec==[]:
-            for X in range(150,351):
-                vec.append(X/100.)    # covers range of most practical 
-                                    # scaling parameters
-        zvec = map(zeta, vec)
-        
-        xmins = unique(x)
-        xmins.sort()
-        xmins = xmins[0:-1]
-        if xminx!=[]:
-            xmins = [min(filter(lambda X: X>=xminx,xmins))]
-        
-        if limit!=[]:
-            limit = round(limit)
-            xmins=filter(lambda X: X<=limit,xmins)
-            if xmins==[]: xmins = [min(x)]
-        
-        if sample!=[]:
-            step = float(len(xmins))/(sample-1)
-            index_curr=0
-            new_xmins=[]
-            for i in range (0,sample):
-                if round(index_curr)==len(xmins): index_curr-=1
-                new_xmins.append(xmins[int(round(index_curr))])
-                index_curr+=step
-            xmins = unique(new_xmins)
-            xmins.sort()
-        
-        if xmins==[]:
-            print '(PLFIT) Error: x must contain at least two unique values.\n'
-            alpha = 'Not a Number'
-            xmin = x[0]
-            D = 'Not a Number'
-            return [alpha,xmin,D]
-        
-        xmax   = max(x)
-        
-        z      = x
-        z.sort()
-        datA=[]
-        datB=[]
-
-        for xm in range(0,len(xmins)):
-            xmin = xmins[xm]
-            z    = filter(lambda X:X>=xmin,z)
-            n    = len(z)
-            # estimate alpha via direct maximization of likelihood function
-
-            # force iterative calculation 
-            L       = []
-            slogz   = sum(map(log,z))
-            xminvec = map(float,range(1,xmin))
-            for k in range(0,len(vec)):
-                L.append(-vec[k]*float(slogz) - float(n)*log(float(zvec[k]) - sum(map(lambda X:pow(float(X),-vec[k]),xminvec))))
-            
-            
-            I = L.index(max(L))
-            # compute KS statistic
-            fit = reduce(lambda X,Y: X+[Y+X[-1]],\
-                         (map(lambda X: pow(X,-vec[I])/(float(zvec[I])-sum(map(lambda X: pow(X,-vec[I]),map(float,range(1,xmin))))),range(xmin,xmax+1))),[0])[1:]
-            cdi=[]
-            for XM in range(xmin,xmax+1):
-                cdi.append(len(filter(lambda X: floor(X)<=XM,z))/float(n))
-            
-            datA.append(max( map(lambda X: abs(fit[X] - cdi[X]),range(0,xmax-xmin+1))))
-            datB.append(vec[I])
-        # select the index for the minimum value of D
-        I = datA.index(min(datA))
-        xmin  = xmins[I]
-        z     = filter(lambda X:X>=xmin,x)
-        n     = len(z)
-        alpha = datB[I]
-        if finite: alpha = alpha*(n-1.)/n+1./n  # finite-size correction
-        if n < 50 and not finite and not nowarn:
-            print '(PLFIT) Warning: finite-size bias may be present.\n'
-        
-        L     = -alpha*sum(map(log,z)) - n*log(zvec[vec.index(max(filter(lambda X:X<=alpha,vec)))] - \
-                                              sum(map(lambda X: pow(X,-alpha),range(1,xmin))))
-    else:
-        print '(PLFIT) Error: x must contain only reals or only integers.\n'
-        alpha = []
-        xmin  = []
-        L     = []
-
-    return [alpha,xmin,L]
-
-
-# helper functions (unique and zeta)
-
-
-def unique(seq): 
-    # not order preserving 
-    set = {} 
-    map(set.__setitem__, seq, []) 
-    return set.keys()
-
-def _polyval(coeffs, x):
-    p = coeffs[0]
-    for c in coeffs[1:]:
-        p = c + x*p
-    return p
-
-_zeta_int = [\
--0.5,
-0.0,
-1.6449340668482264365,1.2020569031595942854,1.0823232337111381915,
-1.0369277551433699263,1.0173430619844491397,1.0083492773819228268,
-1.0040773561979443394,1.0020083928260822144,1.0009945751278180853,
-1.0004941886041194646,1.0002460865533080483,1.0001227133475784891,
-1.0000612481350587048,1.0000305882363070205,1.0000152822594086519,
-1.0000076371976378998,1.0000038172932649998,1.0000019082127165539,
-1.0000009539620338728,1.0000004769329867878,1.0000002384505027277,
-1.0000001192199259653,1.0000000596081890513,1.0000000298035035147,
-1.0000000149015548284]
-
-_zeta_P = [-3.50000000087575873, -0.701274355654678147,
--0.0672313458590012612, -0.00398731457954257841,
--0.000160948723019303141, -4.67633010038383371e-6,
--1.02078104417700585e-7, -1.68030037095896287e-9,
--1.85231868742346722e-11][::-1]
-
-_zeta_Q = [1.00000000000000000, -0.936552848762465319,
--0.0588835413263763741, -0.00441498861482948666,
--0.000143416758067432622, -5.10691659585090782e-6,
--9.58813053268913799e-8, -1.72963791443181972e-9,
--1.83527919681474132e-11][::-1]
-
-_zeta_1 = [3.03768838606128127e-10, -1.21924525236601262e-8,
-2.01201845887608893e-7, -1.53917240683468381e-6,
--5.09890411005967954e-7, 0.000122464707271619326,
--0.000905721539353130232, -0.00239315326074843037,
-0.084239750013159168, 0.418938517907442414, 0.500000001921884009]
-
-_zeta_0 = [-3.46092485016748794e-10, -6.42610089468292485e-9,
-1.76409071536679773e-7, -1.47141263991560698e-6, -6.38880222546167613e-7,
-0.000122641099800668209, -0.000905894913516772796, -0.00239303348507992713,
-0.0842396947501199816, 0.418938533204660256, 0.500000000000000052]
-
-def zeta(s):
-    """
-    Riemann zeta function, real argument
-    """
-    if not isinstance(s, (float, int)):
-        try:
-            s = float(s)
-        except (ValueError, TypeError):
-            try:
-                s = complex(s)
-                if not s.imag:
-                    return complex(zeta(s.real))
-            except (ValueError, TypeError):
-                pass
-            raise NotImplementedError
-    if s == 1:
-        raise ValueError("zeta(1) pole")
-    if s >= 27:
-        return 1.0 + 2.0**(-s) + 3.0**(-s)
-    n = int(s)
-    if n == s:
-        if n >= 0:
-            return _zeta_int[n]
-        if not (n % 2):
-            return 0.0
-    if s <= 0.0:
-        return 0
-    if s <= 2.0:
-        if s <= 1.0:
-            return _polyval(_zeta_0,s)/(s-1)
-        return _polyval(_zeta_1,s)/(s-1)
-    z = _polyval(_zeta_P,s) / _polyval(_zeta_Q,s)
-    return 1.0 + 2.0**(-s) + 3.0**(-s) + 4.0**(-s)*z
-
-
-
-    
Index: trunk/UTIL/PYTHON/powerlaw/plot_granulo.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/plot_granulo.py	(revision 946)
+++ 	(revision )
@@ -1,96 +1,0 @@
-#! /usr/bin/env python
-
-import numpy as np
-from scipy import ndimage
-import matplotlib.pyplot as plt
-
-from myplot import getfield
-from netCDF4 import Dataset
-from scipy import ndimage
-
-def disk_structure(n):
-    struct = np.zeros((2 * n + 1, 2 * n + 1))
-    x, y = np.indices((2 * n + 1, 2 * n + 1))
-    mask = (x - n)**2 + (y - n)**2 <= n**2
-    struct[mask] = 1
-    return struct.astype(np.bool)
-
-
-def granulometry(data, sizes=None):
-    s = max(data.shape)
-    if sizes == None:
-        sizes = range(1, s/2, 2)
-    granulo = [ndimage.binary_opening(data, \
-            structure=disk_structure(n)).sum() for n in sizes]
-    return granulo
-
-
-#np.random.seed(1)
-#n = 10
-#l = 256
-#im = np.zeros((l, l))
-#points = l*np.random.random((2, n**2))
-#im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
-#im = ndimage.gaussian_filter(im, sigma=l/(4.*n))
-#mask = im > im.mean()
-
-namefile = "/home/aymeric/Big_Data/psfc_f18.nc"
-nc  = Dataset(namefile)
-psfc_full=getfield(nc,'PSFC')
-
-#mask = psfc_full[250,100:200,100:200]
-mask = psfc_full[250,:,:]
-mask = ndimage.laplace(mask)
-#mask = mask > 0.3
-
-print mask.std()
-print mask.mean()
-
-lim = mask.std() * 12.
-
-ind1 = np.where(mask > lim)
-ind2 = np.where(mask < lim)
-mask[ind1] = 1
-mask[ind2] = 0
-
-print mask
-
-
-labeled_array, num_features = ndimage.label(mask)
-
-print num_features
-
-plt.pcolor(labeled_array)
-plt.show()
-exit()
-
-
-
-
-mask = mask < mask.mean() - 0.3
-
-#mask = mask > 0.3
-
-
-#plt.pcolor(mask)
-#plt.show()
-
-granulo = granulometry(mask, sizes=[2,5,20])
-
-print granulo
-
-plt.figure(figsize=(6, 2.2))
-
-plt.subplot(121)
-plt.imshow(mask, cmap=plt.cm.gray)
-opened = ndimage.binary_opening(mask, structure=disk_structure(10))
-opened_more = ndimage.binary_opening(mask, structure=disk_structure(14))
-plt.contour(opened, [0.5], colors='b', linewidths=2)
-plt.contour(opened_more, [0.5], colors='r', linewidths=2)
-plt.axis('off')
-plt.subplot(122)
-#plt.plot(np.arange(2, 19, 4), granulo, 'ok', ms=8)
-
-
-plt.subplots_adjust(wspace=0.02, hspace=0.15, top=0.95, bottom=0.15, left=0, right=0.95)
-plt.show()
Index: trunk/UTIL/PYTHON/powerlaw/plplot.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/plplot.py	(revision 946)
+++ 	(revision )
@@ -1,201 +1,0 @@
-import matplotlib.pyplot as plt
-from math import *
-
-# function h=plplot(x, xmin, alpha)
-# PLPLOT visualizes a power-law distributional model with empirical data.
-#    Source: http://www.santafe.edu/~aaronc/powerlaws/
-#
-#    PLPLOT  REQUIRES  the use of the free library: matplotlib
-#
-#    PLPLOT(x, xmin, alpha) plots (on log axes) the data contained in x 
-#    and a power-law distribution of the form p(x) ~ x^-alpha for 
-#    x >= xmin. For additional customization, PLPLOT returns a pair of 
-#    handles, one to the empirical and one to the fitted data series. By 
-#    default, the empirical data is plotted as 'bo' and the fitted form is
-#    plotted as 'k--'. PLPLOT automatically detects whether x is composed 
-#    of real or integer values, and applies the appropriate plotting 
-#    method. For discrete data, if min(x) > 50, PLFIT uses the continuous 
-#    approximation, which is a reliable in this regime.
-#
-#    Example:
-#       xmin  = 47;
-#       alpha = 2.71;
-#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
-#       h = plplot(x,xmin,alpha);
-#
-#    For more information, try 'type plplot'
-#
-#    See also PLFIT, PLVAR, PLPVA
-
-# Version 1.0   (2008 February)
-# Ported to python by Joel Ornstein(2011 July)
-# (joel_ornstein@hmc.edu)
-
-# Copyright (C) 2008-2011 Aaron Clauset (Santa Fe Institute)
-# Distributed under GPL 2.0
-# http://www.gnu.org/copyleft/gpl.html
-# PLFIT comes with ABSOLUTELY NO WARRANTY
-# 
-# The 'zeta' helper function is modified from the open-source library 'mpmath'
-#   mpmath: a Python library for arbitrary-precision floating-point arithmetic
-#   http://code.google.com/p/mpmath/
-#   version 0.17 (February 2011) by Fredrik Johansson and others
-# 
-# No Notes
-#
-
-def plplot(x,xmin,alpha):
-        # select method (discrete or continuous) for fitting
-    if     reduce(lambda X,Y:X==True and floor(Y)==float(Y),x,True): f_dattype = 'INTS'
-    elif reduce(lambda X,Y:X==True and (type(Y)==int or type(Y)==float or type(Y)==long),x,True):    f_dattype = 'REAL'
-    else:                 f_dattype = 'UNKN'
-    
-    if f_dattype=='INTS' and min(x) > 1000 and len(x)>100:
-        f_dattype = 'REAL'
-    plt.close()
-    plt.ion()
-    h=[[],[]]
-    # estimate xmin and alpha, accordingly
-    if f_dattype== 'REAL':
-        n = len(x)
-        c1 = sorted(x)
-        c2 = map(lambda X:X/float(n),range(n,0,-1))
-        q = sorted(filter(lambda X:X>=xmin,x))
-        cf = map(lambda X:pow(float(X)/xmin,1.-alpha),q)
-        cf = map(lambda X:X*float(c2[c1.index(q[0])]),cf)
-
-        h[0]=plt.loglog(c1, c2, 'bo',markersize=8,markerfacecolor=[1,1,1],markeredgecolor=[0,0,1])
-        h[1]=plt.loglog(q, cf, 'k--',linewidth=2)
-        
-        xr1 = pow(10,floor(log(min(x),10)))
-        xr2 = pow(10,ceil(log(min(x),10)))
-        yr1 = pow(10,floor(log(1./n,10)))
-        yr2 = 1
-        
-
-        plt.axhspan(ymin=yr1,ymax=yr2,xmin=xr1,xmax=xr2)
-        plt.ylabel('Pr(X >= x)',fontsize=16);
-        plt.xlabel('x',fontsize=16)
-        plt.draw()
-        
-    elif f_dattype== 'INTS':
-        n = len(x)
-        q = sorted(unique(x))
-        c=[]
-        for Q in q:
-            c.append(len(filter(lambda X: floor(X)==Q,x))/float(n))
-        c1 = q+[q[-1]+1]
-        c2 = map(lambda Z: 1.-Z,reduce(lambda X,Y: X+[Y+X[-1]],c,[0]))
-        c2 = filter(lambda X:float(X)>=pow(10,-10.),c2)
-        c1 = c1[0:len(c2)]
-        cf = map(lambda X:pow(X,-alpha)/(float(zeta(alpha)) - sum(map(lambda Y:pow(Y,-alpha),range(1,xmin)))),range(xmin,q[-1]+1))
-        cf1 = range(xmin,q[-1]+2)
-        cf2 = map(lambda Z: 1.-Z,reduce(lambda X,Y: X+[Y+X[-1]],cf,[0]))
-        cf2 = map(lambda X: X*float(c2[c1.index(xmin)]),cf2)
-
-        h[0]=plt.loglog(c1, c2, 'bo',markersize=8,markerfacecolor=[1,1,1],markeredgecolor=[0,0,1])
-        h[1]=plt.loglog(cf1, cf2, 'k--',linewidth=2)
-        
-        xr1 = pow(10,floor(log(min(x),10)))
-        xr2 = pow(10,ceil(log(min(x),10)))
-        yr1 = pow(10,floor(log(1./n,10)))
-        yr2 = 1
-
-
-        plt.axhspan(ymin=yr1,ymax=yr2,xmin=xr1,xmax=xr2)
-        plt.ylabel('Pr(X >= x)',fontsize=16);
-        plt.xlabel('x',fontsize=16)
-        plt.draw()
-                 
-          
-
-    return h
-
-# helper functions (unique and zeta)
-
-
-def unique(seq): 
-    # not order preserving 
-    set = {} 
-    map(set.__setitem__, seq, []) 
-    return set.keys()
-
-def _polyval(coeffs, x):
-    p = coeffs[0]
-    for c in coeffs[1:]:
-        p = c + x*p
-    return p
-
-_zeta_int = [\
--0.5,
-0.0,
-1.6449340668482264365,1.2020569031595942854,1.0823232337111381915,
-1.0369277551433699263,1.0173430619844491397,1.0083492773819228268,
-1.0040773561979443394,1.0020083928260822144,1.0009945751278180853,
-1.0004941886041194646,1.0002460865533080483,1.0001227133475784891,
-1.0000612481350587048,1.0000305882363070205,1.0000152822594086519,
-1.0000076371976378998,1.0000038172932649998,1.0000019082127165539,
-1.0000009539620338728,1.0000004769329867878,1.0000002384505027277,
-1.0000001192199259653,1.0000000596081890513,1.0000000298035035147,
-1.0000000149015548284]
-
-_zeta_P = [-3.50000000087575873, -0.701274355654678147,
--0.0672313458590012612, -0.00398731457954257841,
--0.000160948723019303141, -4.67633010038383371e-6,
--1.02078104417700585e-7, -1.68030037095896287e-9,
--1.85231868742346722e-11][::-1]
-
-_zeta_Q = [1.00000000000000000, -0.936552848762465319,
--0.0588835413263763741, -0.00441498861482948666,
--0.000143416758067432622, -5.10691659585090782e-6,
--9.58813053268913799e-8, -1.72963791443181972e-9,
--1.83527919681474132e-11][::-1]
-
-_zeta_1 = [3.03768838606128127e-10, -1.21924525236601262e-8,
-2.01201845887608893e-7, -1.53917240683468381e-6,
--5.09890411005967954e-7, 0.000122464707271619326,
--0.000905721539353130232, -0.00239315326074843037,
-0.084239750013159168, 0.418938517907442414, 0.500000001921884009]
-
-_zeta_0 = [-3.46092485016748794e-10, -6.42610089468292485e-9,
-1.76409071536679773e-7, -1.47141263991560698e-6, -6.38880222546167613e-7,
-0.000122641099800668209, -0.000905894913516772796, -0.00239303348507992713,
-0.0842396947501199816, 0.418938533204660256, 0.500000000000000052]
-
-def zeta(s):
-    """
-    Riemann zeta function, real argument
-    """
-    if not isinstance(s, (float, int)):
-        try:
-            s = float(s)
-        except (ValueError, TypeError):
-            try:
-                s = complex(s)
-                if not s.imag:
-                    return complex(zeta(s.real))
-            except (ValueError, TypeError):
-                pass
-            raise NotImplementedError
-    if s == 1:
-        raise ValueError("zeta(1) pole")
-    if s >= 27:
-        return 1.0 + 2.0**(-s) + 3.0**(-s)
-    n = int(s)
-    if n == s:
-        if n >= 0:
-            return _zeta_int[n]
-        if not (n % 2):
-            return 0.0
-    if s <= 0.0:
-        return 0
-    if s <= 2.0:
-        if s <= 1.0:
-            return _polyval(_zeta_0,s)/(s-1)
-        return _polyval(_zeta_1,s)/(s-1)
-    z = _polyval(_zeta_P,s) / _polyval(_zeta_Q,s)
-    return 1.0 + 2.0**(-s) + 3.0**(-s) + 4.0**(-s)*z
-
-
-
-    
Index: trunk/UTIL/PYTHON/powerlaw/plpva.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/plpva.py	(revision 946)
+++ 	(revision )
@@ -1,476 +1,0 @@
-from math import *
-from random import *
-
-# function [alpha, xmin, n]=plpva(x, xmin, varargin)
-# PLPVA calculates the p-value for the given power-law fit to some data.
-#    Source: http://www.santafe.edu/~aaronc/powerlaws/
-# 
-#    PLPVA(x, xmin) takes data x and given lower cutoff for the power-law
-#    behavior xmin and computes the corresponding p-value for the
-#    Kolmogorov-Smirnov test, according to the method described in 
-#    Clauset, Shalizi, Newman (2007).
-#    PLPVA automatically detects whether x is composed of real or integer
-#    values, and applies the appropriate method. For discrete data, if
-#    min(x) > 1000, PLPVA uses the continuous approximation, which is 
-#    a reliable in this regime.
-#   
-#    The fitting procedure works as follows:
-#    1) For each possible choice of x_min, we estimate alpha via the 
-#       method of maximum likelihood, and calculate the Kolmogorov-Smirnov
-#       goodness-of-fit statistic D.
-#    2) We then select as our estimate of x_min, the value that gives the
-#       minimum value D over all values of x_min.
-#
-#    Note that this procedure gives no estimate of the uncertainty of the 
-#    fitted parameters, nor of the validity of the fit.
-#
-#    Example:
-#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
-#       [p, gof] = plpva(x, 1);
-#
-#    For more information, try 'type plpva'
-#
-#    See also PLFIT, PLVAR
-
-# Version 1.0.7 (2009 October)
-# Copyright (C) 2008-2011 Aaron Clauset (Santa Fe Institute)
-
-# Ported to Python by Joel Ornstein (2011 July)
-#(joel_ornstein@hmc.edu)
-
-# Distributed under GPL 2.0
-# http://www.gnu.org/copyleft/gpl.html
-# PLPVA comes with ABSOLUTELY NO WARRANTY
-#
-# 
-# The 'zeta' helper function is modified from the open-source library 'mpmath'
-#   mpmath: a Python library for arbitrary-precision floating-point arithmetic
-#   http://code.google.com/p/mpmath/
-#   version 0.17 (February 2011) by Fredrik Johansson and others
-# 
-
-# Notes:
-# 
-# 1. In order to implement the integer-based methods in Matlab, the numeric
-#    maximization of the log-likelihood function was used. This requires
-#    that we specify the range of scaling parameters considered. We set
-#    this range to be 1.50 to 3.50 at 0.01 intervals by default. 
-#    This range can be set by the user like so,
-#
-#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
-#       a = plpva(x,1,'range',[1.50,3.50,0.01])
-#    
-# 2. PLPVA can be told to limit the range of values considered as estimates
-#    for xmin in two ways. First, it can be instructed to sample these
-#    possible values like so,
-#    
-#       a = plpva(x,1,'sample',100);
-#    
-#    which uses 100 uniformly distributed values on the sorted list of
-#    unique values in the data set. Second, it can simply omit all
-#    candidates above a hard limit, like so
-#    
-#       a = plpva(x,1,'limit',3.4);
-#    
-#    Finally, it can be forced to use a fixed value, like so
-#    
-#       a = plpva(x,1,'xmin',1);
-#    
-#    In the case of discrete data, it rounds the limit to the nearest
-#    integer.
-# 
-# 3. The default number of semiparametric repetitions of the fitting
-# procedure is 1000. This number can be changed like so
-#    
-#       p = plpva(x, 1,'reps',10000);
-# 
-# 4. To silence the textual output to the screen, do this
-#    
-#       p = plpva(x, 1,'reps',10000,'silent');
-# 
-
-def plpva(x,xmin, *varargin):
-    vec     = []
-    sample  = []
-    xminx   = []
-    limit   = []
-    Bt      = []
-    quiet   = False
-    
-
-    # parse command-line parameters trap for bad input
-    i=0 
-    while i<len(varargin): 
-        argok = 1 
-        if type(varargin[i])==str: 
-            if varargin[i]=='range':
-                Range = varargin[i+1]
-                if Range[1]>Range[0]:
-                    argok=0
-                    vec=[]
-                try:
-                    vec=map(lambda X:X*float(Range[2])+Range[0],\
-                            range(int((Range[1]-Range[0])/Range[2])))
-                    
-                    
-                except:
-                    argok=0
-                    vec=[]
-                    
-
-                if Range[0]>=Range[1]:
-                    argok=0
-                    vec=[]
-                    i-=1
-
-                i+=1
-                    
-                            
-            elif varargin[i]== 'sample':
-                sample  = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'limit':
-                limit   = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'xmin':
-                xminx   = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'reps':
-                Bt   = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'silent':       quiet  = True    
- 
-            else: argok=0 
-        
-      
-        if not argok:
-            print '(PLPVA) Ignoring invalid argument #',i+1 
-      
-        i = i+1 
-
-    if vec!=[] and (type(vec)!=list or min(vec)<=1):
-        print '(PLPVA) Error: ''range'' argument must contain a vector or minimum <= 1. using default.\n'                        
-              
-        vec = []
-    
-    if sample!=[] and sample<2:
-        print'(PLPVA) Error: ''sample'' argument must be a positive integer > 1. using default.\n'
-        sample = []
-    
-    if limit!=[] and limit<min(x):
-        print'(PLPVA) Error: ''limit'' argument must be a positive value >= 1. using default.\n'
-        limit = []
-    
-    if xminx!=[] and xminx>=max(x):
-        print'(PLPVA) Error: ''xmin'' argument must be a positive value < max(x). using default behavior.\n'
-        xminx = []
-    if Bt!=[] and Bt<2:
-        print '(PLPVA) Error: ''reps'' argument must be a positive value > 1; using default.\n'
-        Bt = [];
-
-
-    # select method (discrete or continuous) for fitting
-    if     reduce(lambda X,Y:X==True and floor(Y)==float(Y),x,True): f_dattype = 'INTS'
-    elif reduce(lambda X,Y:X==True and (type(Y)==int or type(Y)==float or type(Y)==long),x,True):    f_dattype = 'REAL'
-    else:                 f_dattype = 'UNKN'
-    
-    if f_dattype=='INTS' and min(x) > 1000 and len(x)>100:
-        f_dattype = 'REAL'
-    
-    N=len(x)
-    
-    if Bt==[]: Bt=1000
-    nof = []
-    if not quiet:
-        print 'Power-law Distribution, p-value calculation\n'
-        print '   Copyright 2007-2010 Aaron Clauset\n'
-        print '   Warning: This can be a slow calculation; please be patient.\n'
-        print '   n    = ',len(x),'\n   xmin = ',xmin,'\n   reps = ',Bt
-    # estimate xmin and alpha, accordingly
-    if f_dattype== 'REAL':
-        # compute D for the empirical distribution
-        z = filter(lambda X:X>=xmin,x)
-        z.sort()
-        nz = len(z)
-        y = filter(lambda X:X<xmin,x)
-        ny = len(y)
-        alpha = 1 + float(nz) / sum(map(lambda X: log(float(X)/xmin),z))
-        cz    = map(lambda X:X/float(nz),range(0,nz))
-        cf    = map(lambda X:1-pow((xmin/float(X)),(alpha-1)),z)
-        gof   = max( map(lambda X:abs(cz[X] - cf[X]),range(0,nz)))
-        pz    = nz/float(N)
-        
-        # compute distribution of gofs from semi-parametric bootstrap
-        # of entire data set with fit
-        for B in range(0,Bt):
-            # semi-parametric bootstrap of data
-            n1=0
-            for i in range(0,N):
-                if random()>pz: n1+=1
-
-            q1=[]
-            for i in range(0,n1):
-                q1.append(y[int(floor(ny*random()))])
-            n2 = N-n1;
-            q2=[]
-            for i in range (0,n2):
-                q2.append(random())
-            q2 = map(lambda X:xmin*pow(1.-X,(-1./(alpha-1.))),q2)
-
-            q  = q1+q2
-            q.sort()
-
-            
-
-            # estimate xmin and alpha via GoF-method
-            qmins = unique(q);
-            qmins.sort()
-            qmins = qmins[0:-1];
-            
-            if xminx!=[]:
-                
-                qmins = [min(filter(lambda X: X>=xminx,qmins))]
-                
-            
-            if limit!=[]:
-                qmins=filter(lambda X: X<=limit,qmins)
-                if qmins==[]: qmins=[min(q)]
-                
-            if sample!=[]:
-                step = float(len(qmins))/(sample-1)
-                index_curr=0
-                new_qmins=[]
-                for i in range (0,sample):
-                    if round(index_curr)==len(qmins): index_curr-=1
-                    new_qmins.append(qmins[int(round(index_curr))])
-                    index_curr+=step
-                qmins = unique(new_qmins)
-                qmins.sort()
-            
-            
-        
-            dat   = []
-            
-            for qm in range(0,len(qmins)):
-                qmin = qmins[qm]
-                zq   = filter(lambda X:X>=qmin,q)
-                nq   = len(zq)
-               
-                a    = float(nq)/sum(map(lambda X: log(float(X)/qmin),zq))
-                cf   = map(lambda X:1-pow((float(qmin)/X),a),zq)
-                
-                dat.append( max( map(lambda X: abs(cf[X]-float(X)/nq),range(0,nq))))
-            if not quiet:
-                print '['+str(B+1)+']\tp = ',round(len(filter(lambda X:X>=gof,nof))/float(B+1),3)
-            nof.append(min(dat))
-        p=len(filter(lambda X:X>=gof,nof))/float(Bt+1) 
-
-    elif f_dattype== 'INTS':
-        x=map(int,x)
-        if vec==[]:
-            for X in range(150,351):
-                vec.append(X/100.)    # covers range of most practical 
-                                    # scaling parameters
-        zvec = map(zeta, vec)
-        z = filter (lambda X:X>=xmin,x)
-        nz = len(z)
-        y = filter (lambda X:X<xmin,x)
-        ny = len(y)
-        xmax=max(z)
-        L       = []
-        for k in range(0,len(vec)):
-            L.append(-vec[k]*float(sum(map(log,z))) - float(nz)*log(float(zvec[k]) - sum(map(lambda X:pow(float(X),-vec[k]),range(1,xmin)))))
-        I = L.index(max(L))
-        alpha=vec[I]
-        fit = reduce(lambda X,Y: X+[Y+X[-1]],\
-                     (map(lambda X: pow(X,-alpha)/(float(zvec[I])-sum(map(lambda X: pow(X,-alpha),map(float,range(1,xmin))))),range(xmin,xmax+1))),[0])[1:]
-        cdi=[]
-        for XM in range(xmin,xmax+1):
-            cdi.append(len(filter(lambda X: floor(X)<=XM,z))/float(nz))
-        
-        gof= max( map(lambda X: abs(fit[X] - cdi[X]),range(0,xmax-xmin+1)))
-        pz=nz/float(N)
-        mmax = 20*xmax
-        pdf=[]
-        for i in range(0,xmin):
-            pdf.append(0)
-        pdf += map(lambda X: pow(X,-alpha)/(float(zvec[I])-sum(map(lambda X: pow(X,-alpha),map(float,range(1,xmin))))),range(xmin,mmax+1))
-        cdf = reduce(lambda X,Y: X+[Y+X[-1]],pdf,[0])[1:]+[1]
-
-        # compute distribution of gofs from semi-parametric bootstrap
-        # of entire data set with fit
-        for B in range(0,Bt):
-            # semi-parametric bootstrap of data
-            n1=0
-            for i in range(0,N):
-                if random()>pz: n1+=1
-
-            q1=[]
-            for i in range(0,n1):
-                q1.append( y[int(floor(ny*random()))])
-            n2 = N-n1;
-            r2=[]
-            for i in range (0,n2):
-                r2.append(random())
-            r2.sort()
-            c=0
-            k=0
-            q2=[]
-            for i in range(xmin,mmax+2):
-                while c<n2 and r2[c]<=cdf[i]:c+=1
-                for k in range(k,c):q2.append(i)
-                k=c
-                if k>=n2: break
-
-            q  = q1+q2
-            qmins = unique(q)
-            qmins.sort()
-            qmins = qmins[0:-1]
-
-            if xminx!=[]:
-                
-                qmins = [min(filter(lambda X: X>=xminx,qmins))]
-                
-            
-            if limit!=[]:
-                qmins=filter(lambda X: X<=limit,qmins)
-                if qmins==[]: qmins=[min(q)]
-                
-            if sample!=[]:
-                step = float(len(qmins))/(sample-1)
-                index_curr=0
-                new_qmins=[]
-                for i in range (0,sample):
-                    if round(index_curr)==len(qmins): index_curr-=1
-                    new_qmins.append(qmins[int(round(index_curr))])
-                    index_curr+=step
-                qmins = unique(new_qmins)
-                qmins.sort()
-        
-            qmax   = max(q)
-            dat=[]
-            zq=q
-
-
-            for qm in range(0,len(qmins)):
-                qmin = qmins[qm]
-                zq    = filter(lambda X:X>=qmin,zq)
-                nq    = len(zq)
-                if nq>1:
-                    L=[]
-                    slogzq = sum(map(log,zq))
-                    qminvec = range (1,qmin)
-                    for k in range (1,len(vec)):
-                        L.append(-vec[k]*float(slogzq) - float(nq)*log(float(zvec[k]) - sum(map(lambda X:pow(float(X),-vec[k]),qminvec))))
-            
-                    I = L.index(max(L))
-                    fit = reduce(lambda X,Y: X+[Y+X[-1]],\
-                                 (map(lambda X: pow(X,-vec[I])/(float(zvec[I])-sum(map(lambda X: pow(X,-vec[I]),map(float,range(1,qmin))))),range(qmin,qmax+1))),[0])[1:]
-                    cdi=[]
-                    for QM in range(qmin,qmax+1):
-                        cdi.append(len(filter(lambda X: floor(X)<=QM,zq))/float(nq))
-                    
-                    dat.append(max( map(lambda X: abs(fit[X] - cdi[X]),range(0,qmax-qmin+1))))
-                else: dat[qm]=float('-inf')
-            if not quiet:
-                print '['+str(B+1)+']\tp = ',round(len(filter(lambda X:X>=gof,nof))/float(B+1),3)
-            nof.append(min(dat))
-        
-        p=len(filter(lambda X:X>=gof,nof))/float(Bt)
-    else:
-        print '(PLPVA) Error: x must contain only reals or only integers.\n'
-        p = []
-        gof  = []
-        
-
-    return [p,gof]
-
-
-# helper functions (unique and zeta)
-
-
-def unique(seq): 
-    # not order preserving 
-    set = {} 
-    map(set.__setitem__, seq, []) 
-    return set.keys()
-
-def _polyval(coeffs, x):
-    p = coeffs[0]
-    for c in coeffs[1:]:
-        p = c + x*p
-    return p
-
-_zeta_int = [\
--0.5,
-0.0,
-1.6449340668482264365,1.2020569031595942854,1.0823232337111381915,
-1.0369277551433699263,1.0173430619844491397,1.0083492773819228268,
-1.0040773561979443394,1.0020083928260822144,1.0009945751278180853,
-1.0004941886041194646,1.0002460865533080483,1.0001227133475784891,
-1.0000612481350587048,1.0000305882363070205,1.0000152822594086519,
-1.0000076371976378998,1.0000038172932649998,1.0000019082127165539,
-1.0000009539620338728,1.0000004769329867878,1.0000002384505027277,
-1.0000001192199259653,1.0000000596081890513,1.0000000298035035147,
-1.0000000149015548284]
-
-_zeta_P = [-3.50000000087575873, -0.701274355654678147,
--0.0672313458590012612, -0.00398731457954257841,
--0.000160948723019303141, -4.67633010038383371e-6,
--1.02078104417700585e-7, -1.68030037095896287e-9,
--1.85231868742346722e-11][::-1]
-
-_zeta_Q = [1.00000000000000000, -0.936552848762465319,
--0.0588835413263763741, -0.00441498861482948666,
--0.000143416758067432622, -5.10691659585090782e-6,
--9.58813053268913799e-8, -1.72963791443181972e-9,
--1.83527919681474132e-11][::-1]
-
-_zeta_1 = [3.03768838606128127e-10, -1.21924525236601262e-8,
-2.01201845887608893e-7, -1.53917240683468381e-6,
--5.09890411005967954e-7, 0.000122464707271619326,
--0.000905721539353130232, -0.00239315326074843037,
-0.084239750013159168, 0.418938517907442414, 0.500000001921884009]
-
-_zeta_0 = [-3.46092485016748794e-10, -6.42610089468292485e-9,
-1.76409071536679773e-7, -1.47141263991560698e-6, -6.38880222546167613e-7,
-0.000122641099800668209, -0.000905894913516772796, -0.00239303348507992713,
-0.0842396947501199816, 0.418938533204660256, 0.500000000000000052]
-
-def zeta(s):
-    """
-    Riemann zeta function, real argument
-    """
-    if not isinstance(s, (float, int)):
-        try:
-            s = float(s)
-        except (ValueError, TypeError):
-            try:
-                s = complex(s)
-                if not s.imag:
-                    return complex(zeta(s.real))
-            except (ValueError, TypeError):
-                pass
-            raise NotImplementedError
-    if s == 1:
-        raise ValueError("zeta(1) pole")
-    if s >= 27:
-        return 1.0 + 2.0**(-s) + 3.0**(-s)
-    n = int(s)
-    if n == s:
-        if n >= 0:
-            return _zeta_int[n]
-        if not (n % 2):
-            return 0.0
-    if s <= 0.0:
-        return 0
-    if s <= 2.0:
-        if s <= 1.0:
-            return _polyval(_zeta_0,s)/(s-1)
-        return _polyval(_zeta_1,s)/(s-1)
-    z = _polyval(_zeta_P,s) / _polyval(_zeta_Q,s)
-    return 1.0 + 2.0**(-s) + 3.0**(-s) + 4.0**(-s)*z
-
-
-
-    
Index: trunk/UTIL/PYTHON/powerlaw/plvar.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/plvar.py	(revision 946)
+++ 	(revision )
@@ -1,444 +1,0 @@
-from math import *
-from random import *
-
-# function [alpha, xmin, n]=plvar(x, varargin)
-# PLVAR estimates the uncertainty in the estimated power-law parameters.
-#    Source: http://www.santafe.edu/~aaronc/powerlaws/
-# 
-#    PLVAR(x) takes a vector of observations x and returns estimated
-#    uncertainties in the estimated power-law parameters, based on the
-#    nonparametric approach described in Clauset, Shalizi, Newman (2007). 
-#    PLVAR automatically detects whether x is composed of real or integer
-#    values, and applies the appropriate method. For discrete data, if
-#    min(x) > 1000, PLVAR uses the continuous approximation, which is 
-#    a reliable in this regime.
-#   
-#    The fitting procedure works as follows:
-#    1) For each possible choice of x_min, we estimate alpha via the 
-#       method of maximum likelihood, and calculate the Kolmogorov-Smirnov
-#       goodness-of-fit statistic D.
-#    2) We then select as our estimate of x_min, the value that gives the
-#       minimum value D over all values of x_min.
-#
-#    Note that this procedure gives no estimate of the validity of the fit.
-#
-#    Example:
-#
-#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
-#       [alpha, xmin, ntail] = plvar(x);
-#
-#    For more information, try 'type plvar'
-#
-#    See also PLFIT, PLPVA
-
-
-# Version 1.0.8 (2010 April)
-# Copyright (C) 2008-2011 Aaron Clauset (Santa Fe Institute)
-
-# Ported to Python by Joel Ornstein (2011 July)
-#(joel_ornstein@hmc.edu)
-
-# Distributed under GPL 2.0
-# http://www.gnu.org/copyleft/gpl.html
-# PLVAR comes with ABSOLUTELY NO WARRANTY
-#
-# 
-# The 'zeta' helper function is modified from the open-source library 'mpmath'
-#   mpmath: a Python library for arbitrary-precision floating-point arithmetic
-#   http://code.google.com/p/mpmath/
-#   version 0.17 (February 2011) by Fredrik Johansson and others
-# 
-
-# Notes:
-# 
-# 1. In order to implement the integer-based methods in Matlab, the numeric
-#    maximization of the log-likelihood function was used. This requires
-#    that we specify the range of scaling parameters considered. We set
-#    this range to be 1.50 to 3.50 at 0.01 intervals by default. 
-#    This range can be set by the user like so,
-#    
-#       a = plvar(x,'range',[1.50,3.50,0.01])
-#    
-#    
-# 2. PLVAR can be told to limit the range of values considered as estimates
-#    for xmin in three ways. First, it can be instructed to sample these
-#    possible values like so,
-#    
-#       a = plvar(x,'sample',100);
-#    
-#    which uses 100 uniformly distributed values on the sorted list of
-#    unique values in the data set. Second, it can simply omit all
-#    candidates above a hard limit, like so
-#    
-#       a = plvar(x,'limit',3.4);
-#    
-#    Finally, it can be forced to use a fixed value, like so
-#    
-#       a = plvar(x,'xmin',3.4);
-#    
-#    In the case of discrete data, it rounds the limit to the nearest
-#    integer.
-# 
-# 3. The default number of nonparametric repetitions of the fitting
-# procedure is 1000. This number can be changed like so
-#    
-#       a = plvar(x,'reps',10000);
-#    
-# 4. To silence the textual output to the screen, do this
-#    
-#       p = plvar(x,'silent');
-# 
-#
-
-def plvar(x,*varargin):
-    vec     = []
-    sample  = []
-    xminx   = []
-    limit   = []
-    Bt      = []
-    quiet   = False
-    
-
-    # parse command-line parameters trap for bad input
-    i=0 
-    while i<len(varargin): 
-        argok = 1 
-        if type(varargin[i])==str: 
-            if varargin[i]=='range':
-                Range = varargin[i+1]
-                if Range[1]>Range[0]:
-                    argok=0
-                    vec=[]
-                try:
-                    vec=map(lambda X:X*float(Range[2])+Range[0],\
-                            range(int((Range[1]-Range[0])/Range[2])))
-                    
-                    
-                except:
-                    argok=0
-                    vec=[]
-                    
-
-                if Range[0]>=Range[1]:
-                    argok=0
-                    vec=[]
-                    i-=1
-
-                i+=1
-                    
-                            
-            elif varargin[i]== 'sample':
-                sample  = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'limit':
-                limit   = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'xmin':
-                xminx   = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'reps':
-                Bt   = varargin[i+1]
-                i = i + 1
-            elif varargin[i]==  'silent':       quiet  = True    
- 
-            else: argok=0 
-        
-      
-        if not argok:
-            print '(PLVAR) Ignoring invalid argument #',i+1 
-      
-        i = i+1 
-
-    if vec!=[] and (type(vec)!=list or min(vec)<=1):
-        print '(PLVAR) Error: ''range'' argument must contain a vector or minimum <= 1. using default.\n'                        
-              
-        vec = []
-    
-    if sample!=[] and sample<2:
-        print'(PLVAR) Error: ''sample'' argument must be a positive integer > 1. using default.\n'
-        sample = []
-    
-    if limit!=[] and limit<min(x):
-        print'(PLVAR) Error: ''limit'' argument must be a positive value >= 1. using default.\n'
-        limit = []
-    
-    if xminx!=[] and xminx>=max(x):
-        print'(PLVAR) Error: ''xmin'' argument must be a positive value < max(x). using default behavior.\n'
-        xminx = []
-    if Bt!=[] and Bt<2:
-        print '(PLVAR) Error: ''reps'' argument must be a positive value > 1; using default.\n'
-        Bt = [];
-
-
-    # select method (discrete or continuous) for fitting
-    if     reduce(lambda X,Y:X==True and floor(Y)==float(Y),x,True): f_dattype = 'INTS'
-    elif reduce(lambda X,Y:X==True and (type(Y)==int or type(Y)==float or type(Y)==long),x,True):    f_dattype = 'REAL'
-    else:                 f_dattype = 'UNKN'
-    
-    if f_dattype=='INTS' and min(x) > 1000 and len(x)>100:
-        f_dattype = 'REAL'
-    
-    N=len(x)
-    
-    if Bt==[]: Bt=1000
-    bofA = []
-    bofB = []
-    bofC = []
-    if not quiet:
-        print 'Power-law Distribution, parameter error calculation\n'
-        print '   Copyright 2007-2009 Aaron Clauset\n'
-        print '   Warning: This can be a slow calculation; please be patient.\n'
-        print '   n    = ',len(x),'\n   reps = ',Bt
-    # estimate xmin and alpha, accordingly
-    if f_dattype== 'REAL':
-
-        for B in range(0,Bt):
-            #  bootstrap resample
-            y = []
-            for i in range(0,N):
-                y.append(x[int(floor(N*random()))])
-            ymins = unique(y)
-            ymins.sort()
-            ymins=ymins[0:-1]
-            
-            if xminx!=[]:
-                
-                ymins = [min(filter(lambda X: X>=xminx,ymins))]
-                
-            
-            if limit!=[]:
-                qmins=filter(lambda X: X<=limit,qmins)
-                if qmins==[]: qmins=[min(y)]
-                
-            if sample!=[]:
-                step = float(len(ymins))/(sample-1)
-                index_curr=0
-                new_ymins=[]
-                for i in range (0,sample):
-                    if round(index_curr)==len(ymins): index_curr-=1
-                    new_ymins.append(ymins[int(round(index_curr))])
-                    index_curr+=step
-                ymins = unique(new_ymins)
-                ymins.sort()
-
-            z = sorted(y)
-
-            dat   = []
-            
-            for xm in range(0,len(ymins)):
-                xmin = ymins[xm]
-                z   = filter(lambda X:X>=xmin,z)
-                n   = len(z)
-                # estimate alpha using direct MLE
-                a    = float(n)/sum(map(lambda X: log(float(X)/xmin),z))
-                # compute KS statistic
-                cf   = map(lambda X:1-pow((float(xmin)/X),a),z)
-                
-                dat.append( max( map(lambda X: abs(cf[X]-float(X)/n),range(0,n))))
-            ymin = ymins[dat.index(min(dat))]
-            z = filter(lambda X: X>=ymin,y)
-            n = len(z)
-            alpha    = 1+float(n)/sum(map(lambda X: log(float(X)/ymin),z))
-            bofA.append(n)
-            bofB.append(ymin)
-            bofC.append(alpha)
-            # store distribution of estimated parameter values 
-            if not quiet:
-                print '['+str(B+1)+']\tntail = ',round(mean(bofA),3),' (',round(std(bofA),3),')','\txmin = ',\
-                      round(mean(bofB),3),' (',round(std(bofB),3),')','\talpha = ',round(mean(bofC),3),' (',round(std(bofC),3),')'
-
-        n = std(bofA)
-        xmin = std(bofB)
-        alpha = std(bofC)
-
-    elif f_dattype== 'INTS':
-        x=map(int,x)
-        if vec==[]:
-            for X in range(150,351):
-                vec.append(X/100.)    # covers range of most practical 
-                                    # scaling parameters
-        zvec = map(zeta, vec)
-
-        for B in range(0,Bt):
-            #  bootstrap resample
-            y = []
-            for i in range(0,N):
-                y.append(x[int(floor(N*random()))])
-            ymins = unique(y)
-            ymins.sort()
-            ymins=ymins[0:-1]
-            
-            if xminx!=[]:
-                
-                ymins = [min(filter(lambda X: X>=xminx,ymins))]
-                
-            
-            if limit!=[]:
-                qmins=filter(lambda X: X<=limit,qmins)
-                if qmins==[]: qmins=[min(y)]
-                
-            if sample!=[]:
-                step = float(len(ymins))/(sample-1)
-                index_curr=0
-                new_ymins=[]
-                for i in range (0,sample):
-                    if round(index_curr)==len(ymins): index_curr-=1
-                    new_ymins.append(ymins[int(round(index_curr))])
-                    index_curr+=step
-                ymins = unique(new_ymins)
-                ymins.sort()
-                
-            ymax = max(y)
-            z = sorted(y)
-
-            datA = []
-            datB = []
-            for xm in range(0,len(ymins)):
-                xmin = ymins[xm]
-                z = filter(lambda X:X>=xmin,z)
-                n = len(z)
-                L = []
-                slogz = sum(map(log,z))
-                xminvec = range (1,xmin)
-                for k in range (1,len(vec)):
-                    L.append(-vec[k]*float(slogz) - float(n)*log(float(zvec[k]) - sum(map(lambda X:pow(float(X),-vec[k]),xminvec))))
-        
-                I = L.index(max(L))
-
-                # compute KS statistic
-                fit = reduce(lambda X,Y: X+[Y+X[-1]],\
-                             (map(lambda X: pow(X,-vec[I])/(float(zvec[I])-sum(map(lambda X: pow(X,-vec[I]),map(float,range(1,xmin))))),range(xmin,ymax+1))),[0])[1:]
-                cdi=[]
-                for XM in range(xmin,ymax+1):
-                    cdi.append(len(filter(lambda X: floor(X)<=XM,z))/float(n))
-                
-                datA.append(max( map(lambda X: abs(fit[X] - cdi[X]),range(0,ymax-xmin+1))))
-                datB.append(vec[I])
-                
-            I = datA.index(min(datA))
-            ymin = ymins[I]
-            n = len(filter(lambda X:X>=ymin,y))
-            alpha = datB[I]
-
-            bofA.append(n)
-            bofB.append(ymin)
-            bofC.append(alpha)
-            # store distribution of estimated parameter values 
-            if not quiet:
-                print '['+str(B+1)+']\tntail = ',round(mean(bofA),3),' (',round(std(bofA),3),')','\txmin = ',\
-                      round(mean(bofB),3),' (',round(std(bofB),3),')','\talpha = ',round(mean(bofC),3),' (',round(std(bofC),3),')'
-        n = std(bofA)
-        xmin = std(bofB)
-        alpha = std(bofC)
-        
-    else:
-        print '(PLVAR) Error: x must contain only reals or only integers.\n'
-        n = []
-        xmin = []
-        alpha = []
-        
-
-    return [alpha,xmin,n]
-
-
-# helper functions (unique and zeta)
-def mean(L):
-    try:
-        return float(sum(L))/len(L)
-    except:
-        
-        return 0
-
-def std(L):
-    try:
-        u = mean(L)
-        return sqrt((1./(len(L)-1))*sum(map(lambda X: pow(X-u,2),L)))
-    except:
-        return 0
-        
-
-def unique(seq): 
-    # not order preserving 
-    set = {} 
-    map(set.__setitem__, seq, []) 
-    return set.keys()
-
-def _polyval(coeffs, x):
-    p = coeffs[0]
-    for c in coeffs[1:]:
-        p = c + x*p
-    return p
-
-_zeta_int = [\
--0.5,
-0.0,
-1.6449340668482264365,1.2020569031595942854,1.0823232337111381915,
-1.0369277551433699263,1.0173430619844491397,1.0083492773819228268,
-1.0040773561979443394,1.0020083928260822144,1.0009945751278180853,
-1.0004941886041194646,1.0002460865533080483,1.0001227133475784891,
-1.0000612481350587048,1.0000305882363070205,1.0000152822594086519,
-1.0000076371976378998,1.0000038172932649998,1.0000019082127165539,
-1.0000009539620338728,1.0000004769329867878,1.0000002384505027277,
-1.0000001192199259653,1.0000000596081890513,1.0000000298035035147,
-1.0000000149015548284]
-
-_zeta_P = [-3.50000000087575873, -0.701274355654678147,
--0.0672313458590012612, -0.00398731457954257841,
--0.000160948723019303141, -4.67633010038383371e-6,
--1.02078104417700585e-7, -1.68030037095896287e-9,
--1.85231868742346722e-11][::-1]
-
-_zeta_Q = [1.00000000000000000, -0.936552848762465319,
--0.0588835413263763741, -0.00441498861482948666,
--0.000143416758067432622, -5.10691659585090782e-6,
--9.58813053268913799e-8, -1.72963791443181972e-9,
--1.83527919681474132e-11][::-1]
-
-_zeta_1 = [3.03768838606128127e-10, -1.21924525236601262e-8,
-2.01201845887608893e-7, -1.53917240683468381e-6,
--5.09890411005967954e-7, 0.000122464707271619326,
--0.000905721539353130232, -0.00239315326074843037,
-0.084239750013159168, 0.418938517907442414, 0.500000001921884009]
-
-_zeta_0 = [-3.46092485016748794e-10, -6.42610089468292485e-9,
-1.76409071536679773e-7, -1.47141263991560698e-6, -6.38880222546167613e-7,
-0.000122641099800668209, -0.000905894913516772796, -0.00239303348507992713,
-0.0842396947501199816, 0.418938533204660256, 0.500000000000000052]
-
-def zeta(s):
-    """
-    Riemann zeta function, real argument
-    """
-    if not isinstance(s, (float, int)):
-        try:
-            s = float(s)
-        except (ValueError, TypeError):
-            try:
-                s = complex(s)
-                if not s.imag:
-                    return complex(zeta(s.real))
-            except (ValueError, TypeError):
-                pass
-            raise NotImplementedError
-    if s == 1:
-        raise ValueError("zeta(1) pole")
-    if s >= 27:
-        return 1.0 + 2.0**(-s) + 3.0**(-s)
-    n = int(s)
-    if n == s:
-        if n >= 0:
-            return _zeta_int[n]
-        if not (n % 2):
-            return 0.0
-    if s <= 0.0:
-        return 0
-    if s <= 2.0:
-        if s <= 1.0:
-            return _polyval(_zeta_0,s)/(s-1)
-        return _polyval(_zeta_1,s)/(s-1)
-    z = _polyval(_zeta_P,s) / _polyval(_zeta_Q,s)
-    return 1.0 + 2.0**(-s) + 3.0**(-s) + 4.0**(-s)*z
-
-
-
-    
Index: trunk/UTIL/PYTHON/powerlaw/previous_tests/find_devils.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/previous_tests/find_devils.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/previous_tests/find_devils.py	(revision 1053)
@@ -0,0 +1,300 @@
+#! /usr/bin/env python
+
+def detsize( xx, res=1, thres=3, loga=False ):
+    import numpy as np
+    import math
+    size = []
+    sizecalc = 1
+    diff = np.asarray( np.roll(xx,-1) - xx )
+    for i in diff:
+        if abs(i) > 1:
+            if sizecalc >= thres: 
+                if loga: addthis = math.log(sizecalc*res)
+                else:    addthis = sizecalc*res
+                size.append(addthis)
+            sizecalc = 1
+        else:
+            sizecalc += 1
+    return size
+
+def getsize(filename):
+
+    import numpy as np
+    from scipy.ndimage.measurements import minimum_position
+    from scipy import ndimage
+    from netCDF4 import Dataset
+    import matplotlib.pyplot as plt
+    import myplot as myp
+   
+    ### LOAD NETCDF DATA
+    nc = Dataset(filename) 
+    psfc = nc.variables["PSFC"]
+    print "yeah"
+    
+    ### LOOP on TIME
+    ### NB: a same event could be counted several times...
+    shape = np.array(psfc).shape
+    allsizesx = []
+    allsizesy = []
+    depression = []
+    stride = 1 #5
+    stride = 20
+    #stride = 50
+    stride = 100
+    start = 0
+    start = stride
+    for i in range(start,shape[0],stride):
+
+        psfc2d = np.array ( psfc [ i, : , : ] )
+    
+        ############### CRITERION
+        ave = np.mean(psfc2d,dtype=np.float64)  ## dtype otherwise inaccuracy
+
+        #limdp = -0.2 ## on en loupe pas mal
+        #where = np.where(psfc2d - ave < limdp)  ## comme le papier Phoenix
+    
+        std = np.std(psfc2d,dtype=np.float64)   ## dtype otherwise inaccuracy
+        fac = 4.   ## how many sigmas. not too low, otherwise vortices are not caught. 4 good choice.
+                   ## 2.5 clearly too low, 3.5 not too bad, 4 probably good
+        fac = 3.5
+        fac = 3.2
+        ##fac = 2.5
+        #fac = 3. ## final choice
+        #fac = 2.5
+        lim = ave - fac*std
+        where = np.where(psfc2d < lim)
+        ############### END CRITERION
+
+        depression = np.append(depression,np.ravel(psfc2d[where])-ave)
+   
+        ## lab is 0 or 1
+        lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine
+        lab[where] = 1.  ## do not treat points close to 'mean' (background) pressure
+  
+        xx = []
+        yy = []
+        while 1 in lab:
+            p = minimum_position(psfc2d,labels=lab)
+            lab[p] = 0 ## once a minimum has been found in a grid point, do not search here again.
+            if p[0] not in xx: xx.append(p[0]) ## if x coordinate not yet in the list add it
+            if p[1] not in yy: yy.append(p[1]) ## if y coordinate not yet in the list add it
+        xx.sort()
+        yy.sort()
+        ### now xx and yy are sorted arrays containing grid points with pressure minimum
+       
+        ######## DETERMINE SIZE OF STRUCTURES
+        ######## this is rather brute-force...
+        sizex = detsize( xx, res = 10, loga=False, thres=2 )
+        sizey = detsize( yy, res = 10, loga=False, thres=2 )
+        #sizex = detsize( xx, res = 10, loga=False, thres=3 )
+        #sizey = detsize( yy, res = 10, loga=False, thres=3 )
+        sizex = detsize( xx, res = 15, loga=False, thres=2 )
+        sizey = detsize( yy, res = 15, loga=False, thres=2 )
+        ###
+        print sizex, sizey
+        #if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex  ### un peu limite dans certains cas
+        if (len(sizex) > len(sizey))     : sizey = sizex        ### plus fidele mais petit souci lorsque PBC
+        elif (len(sizex) == len(sizey))  : 
+          if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex
+          else                                  : sizex = sizey
+        else                            : sizex = sizey
+        allsizesx = np.append(allsizesx,sizex)
+        allsizesy = np.append(allsizesy,sizey)
+        print i, ' on ', shape[0], ' caught ', len(sizex), ' vortices ', sizex
+        ########
+  
+    allsizesx.sort()
+    allsizesy.sort()
+    
+    return allsizesx, allsizesy, depression
+
+#########################################################################
+#########################################################################
+
+import matplotlib.pyplot as plt
+import pickle
+import numpy as np
+import matplotlib.mlab as mlab
+import mymath as mym
+import myplot as myp
+
+import plfit
+import plplot
+import randht
+import plpva
+
+save = True
+#save = False
+pression = False
+pression = True
+
+filename = "/home/aymeric/Big_Data/psfc_f18.nc"
+
+if save:
+    ### getsize
+    allsizesx, allsizesy, depression = getsize(filename)
+    ### sauvegarde texte pour inspection
+    mym.writeascii(allsizesx,'allsizex.txt')
+    mym.writeascii(allsizesy,'allsizey.txt')
+    mym.writeascii(depression,'alldepression.txt')
+    ### sauvegarde binaire pour utilisation python
+    myfile = open('allsizex.bin', 'wb') ; pickle.dump(allsizesx, myfile) ; myfile.close()
+    myfile = open('allsizey.bin', 'wb') ; pickle.dump(allsizesy, myfile) ; myfile.close()
+    myfile = open('alldepression.bin', 'wb') ; pickle.dump(depression, myfile) ; myfile.close()
+
+### load files
+myfile = open('allsizex.bin', 'r')
+allsizesx = pickle.load(myfile)
+myfile = open('allsizey.bin', 'r')
+allsizesy = pickle.load(myfile)
+myfile = open('alldepression.bin', 'r')
+depression = pickle.load(myfile)
+depression = np.array(abs(depression))#*1000.
+
+### sizes
+#plothist = np.append(allsizesx,allsizesy)
+plothist = allsizesx
+if pression: plothist = depression
+plothist.sort()
+print 'mean ', np.mean(plothist,dtype=np.float64)  
+print 'std ', np.std(plothist,dtype=np.float64)
+print 'max ', np.max(plothist)
+print 'min ', np.min(plothist)
+print 'len ', len(plothist)
+
+
+### MAKE BINS
+nbins = 100
+zebins = [2.0]
+#nbins = 8
+#zebins = [19.0]
+#nbins = 15
+#zebins = [11.] ##12 non mais donne un peu la meme chose
+#zebins = [20.] ##20 non car trop pres du premier
+nbins = 100
+zebins = [2./np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
+nbins = 200
+
+if pression: zebins = [0.3]
+
+for i in range(0,nbins):  zebins.append(zebins[i]*np.sqrt(2))
+zebins = np.array(zebins)
+#### select reasonable bins for DD
+if not pression:
+    zebins = zebins [ zebins > 15. ] 
+    #zebins = zebins [ zebins > 20. ]
+    zebins = zebins [ zebins > 25. ]
+    zebins = zebins [ zebins < 1000. ]
+else:
+    zebins = zebins [ zebins < 10. ]
+print 'corrected bins ',zebins
+
+#### HISTOGRAM
+plt.figure(1)
+plt.hist(    plothist,\
+             log=True,\
+             bins=zebins,\
+#             cumulative=-1,\
+             normed=True,\
+        )
+plt.xscale('log')
+if pression:    plt.xlabel('Pressure (Pa)')
+else:           plt.xlabel('Diameter (m)')
+plt.ylabel('Population (normalized)')
+if pression: prefix="p"
+else:        prefix=""
+myp.makeplotres(prefix+"histogram",res=200,disp=False)
+plt.close(1)
+
+### COMPARED HISTOGRAMS
+### --- FIT WITH POWER LAW
+if pression: [alpha, xmin, L] = plfit.plfit(plothist,'xmin',0.3)
+else:        [alpha, xmin, L] = plfit.plfit(plothist,'limit',20.)
+print alpha,xmin
+
+#a = plpva.plpva(plothist,0.75,'xmin',0.75)
+#print a
+
+#### DEUXIEME ROUTINE
+####IL FAUT UTILISER LE DISCRET POUR LA TAILLE !!!
+#if pression:   myplfit = plfit.plfit(plothist,verbose=True,xmin=0.75)
+#else:          myplfit = plfit.plfit(plothist,verbose=True,xmin=20.)
+#myplfit.plotppf()
+#plt.show()
+#exit()
+
+
+#plt.figure(1)
+#h = plplot.plplot(plothist,xmin,alpha)
+#myp.makeplotres(prefix+"fit",res=200,disp=False)
+plt.figure(2)
+### --- POWER LAW (factor does not really matter)
+power = (xmin/2.2)*np.array(randht.randht(10000,'powerlaw',alpha))
+#power = (xmin/2.2)*np.array(randht.randht(10000,'cutoff',alpha,10.)) ##marche pas si trop grand
+print 'mean ', np.mean(power,dtype=np.float64)
+### --- EXPONENTIAL LAW
+expo = randht.randht(10000,'exponential',1./(np.mean(power,dtype=np.float64)*1.00))
+print 'mean ', np.mean(expo,dtype=np.float64)
+### --- PLOT
+plt.hist(    [plothist,power,expo],\
+             label=['LES vortices','Power law '+'{:.1f}'.format(alpha),'Exponential law'],\
+             log=True,\
+             bins=zebins,\
+#            cumulative=-1,\
+             normed=True,\
+        )
+plt.legend()
+plt.xscale('log')
+if pression:    plt.xlabel('Pressure (Pa)')
+else:           plt.xlabel('Diameter (m)')
+plt.ylabel('Population (normalized)')
+myp.makeplotres(prefix+"comparison",res=200,disp=False)
+plt.close(2)
+
+########################
+########################
+zebins = [30.,42.,60.,84.,120.,170.,240.,340.]
+plothist = []
+plothist = np.append(plothist,30 *np.ones(306))
+plothist = np.append(plothist,42 *np.ones(58) )
+plothist = np.append(plothist,60 *np.ones(66) )
+plothist = np.append(plothist,84 *np.ones(41) )
+plothist = np.append(plothist,120*np.ones(19) )
+plothist = np.append(plothist,170*np.ones(9)  )
+plothist = np.append(plothist,240*np.ones(2)  )
+plothist = np.append(plothist,340*np.ones(1)  )
+
+#zebins = [50.,71.,100.,141.,200.,282.,400.]
+#plothist = []
+#plothist = np.append(plothist,50. *np.ones(36))
+#plothist = np.append(plothist,71. *np.ones(18) )
+#plothist = np.append(plothist,100. *np.ones(12) )
+#plothist = np.append(plothist,141. *np.ones(6) )
+#plothist = np.append(plothist,200.*np.ones(4) )
+#plothist = np.append(plothist,282.*np.ones(1)  )
+#plothist = np.append(plothist,400.*np.ones(2)  )
+
+exit()
+
+plt.figure(3)
+[alpha, xmin, L] = plfit.plfit(plothist,'xmin',30)#50.)
+print alpha,xmin
+#a = plpva.plpva(plothist,30,'xmin',30)
+h = plplot.plplot(plothist,xmin,alpha)
+plt.loglog(h[0], h[1], 'k--',linewidth=2)
+plt.hist(    plothist,\
+             log=True,\
+             bins=zebins,\
+#             cumulative=-1,\
+             normed=True,\
+        )
+plt.xscale('log')
+plt.xlabel('Pressure (micro Pa)')
+plt.ylabel('Population (normalized)')
+myp.makeplotres("data",res=200,disp=False)
+
+#plt.figure(4)
+#[alpha, xmin, L] = plfit.plfit(plothist,'xmin',50.) #,'xmin',30.)
+#print alpha,xmin
+#h = plplot.plplot(plothist,xmin,alpha)
+#myp.makeplotres("datafit",res=200,disp=False)
Index: trunk/UTIL/PYTHON/powerlaw/previous_tests/plfit.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/previous_tests/plfit.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/previous_tests/plfit.py	(revision 1053)
@@ -0,0 +1,429 @@
+from math import *
+
+# function [alpha, xmin, L]=plfit(x, varargin)
+# PLFIT fits a power-law distributional model to data.
+#    Source: http://www.santafe.edu/~aaronc/powerlaws/
+# 
+#    PLFIT(x) estimates x_min and alpha according to the goodness-of-fit
+#    based method described in Clauset, Shalizi, Newman (2007). x is a 
+#    vector of observations of some quantity to which we wish to fit the 
+#    power-law distribution p(x) ~  x^-alpha for x >= xmin.
+#    PLFIT automatically detects whether x is composed of real or integer
+#    values, and applies the appropriate method. For discrete data, if
+#    min(x) > 1000, PLFIT uses the continuous approximation, which is 
+#    a reliable in this regime.
+#   
+#    The fitting procedure works as follows:
+#    1) For each possible choice of x_min, we estimate alpha via the 
+#       method of maximum likelihood, and calculate the Kolmogorov-Smirnov
+#       goodness-of-fit statistic D.
+#    2) We then select as our estimate of x_min, the value that gives the
+#       minimum value D over all values of x_min.
+#
+#    Note that this procedure gives no estimate of the uncertainty of the 
+#    fitted parameters, nor of the validity of the fit.
+#
+#    Example:
+#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
+#       [alpha, xmin, L] = plfit(x)
+#   or  a = plfit(x)
+#
+#    The output 'alpha' is the maximum likelihood estimate of the scaling
+#    exponent, 'xmin' is the estimate of the lower bound of the power-law
+#    behavior, and L is the log-likelihood of the data x>=xmin under the
+#    fitted power law.
+#    
+#    For more information, try 'type plfit'
+#
+#    See also PLVAR, PLPVA
+
+# Version 1.0.10 (2010 January)
+# Copyright (C) 2008-2011 Aaron Clauset (Santa Fe Institute)
+
+# Ported to Python by Joel Ornstein (2011 July)
+# (joel_ornstein@hmc.edu)
+
+# Distributed under GPL 2.0
+# http://www.gnu.org/copyleft/gpl.html
+# PLFIT comes with ABSOLUTELY NO WARRANTY
+#
+# 
+# The 'zeta' helper function is modified from the open-source library 'mpmath'
+#   mpmath: a Python library for arbitrary-precision floating-point arithmetic
+#   http://code.google.com/p/mpmath/
+#   version 0.17 (February 2011) by Fredrik Johansson and others
+# 
+
+# Notes:
+# 
+# 1. In order to implement the integer-based methods in Matlab, the numeric
+#    maximization of the log-likelihood function was used. This requires
+#    that we specify the range of scaling parameters considered. We set
+#    this range to be 1.50 to 3.50 at 0.01 intervals by default. 
+#    This range can be set by the user like so,
+#    
+#       a = plfit(x,'range',[1.50,3.50,0.01])
+#    
+# 2. PLFIT can be told to limit the range of values considered as estimates
+#    for xmin in three ways. First, it can be instructed to sample these
+#    possible values like so,
+#    
+#       a = plfit(x,'sample',100)
+#    
+#    which uses 100 uniformly distributed values on the sorted list of
+#    unique values in the data set. Second, it can simply omit all
+#    candidates above a hard limit, like so
+#    
+#       a = plfit(x,'limit',3.4)
+#    
+#    Finally, it can be forced to use a fixed value, like so
+#    
+#       a = plfit(x,'xmin',3.4)
+#    
+#    In the case of discrete data, it rounds the limit to the nearest
+#    integer.
+# 
+# 3. When the input sample size is small (e.g., < 100), the continuous 
+#    estimator is slightly biased (toward larger values of alpha). To
+#    explicitly use an experimental finite-size correction, call PLFIT like
+#    so
+#    
+#       a = plfit(x,'finite')
+#    
+#    which does a small-size correction to alpha.
+#
+# 4. For continuous data, PLFIT can return erroneously large estimates of 
+#    alpha when xmin is so large that the number of obs x >= xmin is very 
+#    small. To prevent this, we can truncate the search over xmin values 
+#    before the finite-size bias becomes significant by calling PLFIT as
+#    
+#       a = plfit(x,'nosmall')
+#    
+#    which skips values xmin with finite size bias > 0.1.
+
+def plfit(x, *varargin):
+    vec     = []
+    sample  = []
+    xminx   = []
+    limit   = []
+    finite  = False
+    nosmall = False
+    nowarn  = False
+
+    # parse command-line parameters trap for bad input
+    i=0 
+    while i<len(varargin): 
+        argok = 1 
+        if type(varargin[i])==str: 
+            if varargin[i]=='range':
+                Range = varargin[i+1]
+                if Range[1]>Range[0]:
+                    argok=0
+                    vec=[]
+                try:
+                    vec=map(lambda X:X*float(Range[2])+Range[0],\
+                            range(int((Range[1]-Range[0])/Range[2])))
+                    
+                    
+                except:
+                    argok=0
+                    vec=[]
+                    
+
+                if Range[0]>=Range[1]:
+                    argok=0
+                    vec=[]
+                    i-=1
+
+                i+=1
+                    
+                            
+            elif varargin[i]== 'sample':
+                sample  = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'limit':
+                limit   = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'xmin':
+                xminx   = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'finite':       finite  = True    
+            elif varargin[i]==  'nowarn':       nowarn  = True    
+            elif varargin[i]==  'nosmall':      nosmall = True    
+            else: argok=0 
+        
+      
+        if not argok:
+            print '(PLFIT) Ignoring invalid argument #',i+1 
+      
+        i = i+1 
+
+    if vec!=[] and (type(vec)!=list or min(vec)<=1):
+        print '(PLFIT) Error: ''range'' argument must contain a vector or minimum <= 1. using default.\n'                        
+              
+        vec = []
+    
+    if sample!=[] and sample<2:
+        print'(PLFIT) Error: ''sample'' argument must be a positive integer > 1. using default.\n'
+        sample = []
+    
+    if limit!=[] and limit<min(x):
+        print'(PLFIT) Error: ''limit'' argument must be a positive value >= 1. using default.\n'
+        limit = []
+    
+    if xminx!=[] and xminx>=max(x):
+        print'(PLFIT) Error: ''xmin'' argument must be a positive value < max(x). using default behavior.\n'
+        xminx = []
+    
+
+
+    # select method (discrete or continuous) for fitting
+    if     reduce(lambda X,Y:X==True and floor(Y)==float(Y),x,True): f_dattype = 'INTS'
+    elif reduce(lambda X,Y:X==True and (type(Y)==int or type(Y)==float or type(Y)==long),x,True):    f_dattype = 'REAL'
+    #else:                 f_dattype = 'UNKN'
+    ##AS
+    else:                 f_dattype = 'REAL'  
+
+    if f_dattype=='INTS' and min(x) > 1000 and len(x)>100:
+        f_dattype = 'REAL'
+   
+    print f_dattype
+
+    # estimate xmin and alpha, accordingly
+        
+    if f_dattype== 'REAL':
+        xmins = unique(x)
+        xmins.sort()
+        xmins = xmins[0:-1]
+        if xminx!=[]:
+            
+            xmins = [min(filter(lambda X: X>=xminx,xmins))]
+            
+        
+        if limit!=[]:
+            xmins=filter(lambda X: X<=limit,xmins)
+            if xmins==[]: xmins = [min(x)]
+            
+        if sample!=[]:
+            step = float(len(xmins))/(sample-1)
+            index_curr=0
+            new_xmins=[]
+            for i in range (0,sample):
+                if round(index_curr)==len(xmins): index_curr-=1
+                new_xmins.append(xmins[int(round(index_curr))])
+                index_curr+=step
+            xmins = unique(new_xmins)
+            xmins.sort()
+            
+            
+        
+        dat   = []
+        z     = sorted(x)
+        
+        for xm in range(0,len(xmins)):
+            xmin = xmins[xm]
+            z    = filter(lambda X:X>=xmin,z)
+            
+            n    = len(z)
+            # estimate alpha using direct MLE
+
+            a    = float(n) / sum(map(lambda X: log(float(X)/xmin),z))
+            if nosmall:
+                if (a-1)/sqrt(n) > 0.1 and dat!=[]:
+                    xm = len(xmins)+1
+                    break
+                
+            
+            # compute KS statistic
+            #cx   = map(lambda X:float(X)/n,range(0,n))
+            cf   = map(lambda X:1-pow((float(xmin)/X),a),z)
+            dat.append( max( map(lambda X: abs(cf[X]-float(X)/n),range(0,n))))
+        D     = min(dat)
+        xmin  = xmins[dat.index(D)]
+        z     = filter(lambda X:X>=xmin,x)
+        z.sort()
+        n     = len(z) 
+        alpha = 1 + n / sum(map(lambda X: log(float(X)/xmin),z))
+        if finite: alpha = alpha*float(n-1)/n+1./n  # finite-size correction
+        if n < 50 and not finite and not nowarn:
+            print '(PLFIT) Warning: finite-size bias may be present.\n'
+        
+        L = n*log((alpha-1)/xmin) - alpha*sum(map(lambda X: log(float(X)/xmin),z))
+    elif f_dattype== 'INTS':
+        
+        x=map(int,x)
+        if vec==[]:
+            for X in range(150,351):
+                vec.append(X/100.)    # covers range of most practical 
+                                    # scaling parameters
+        zvec = map(zeta, vec)
+        
+        xmins = unique(x)
+        xmins.sort()
+        xmins = xmins[0:-1]
+        if xminx!=[]:
+            xmins = [min(filter(lambda X: X>=xminx,xmins))]
+        
+        if limit!=[]:
+            limit = round(limit)
+            xmins=filter(lambda X: X<=limit,xmins)
+            if xmins==[]: xmins = [min(x)]
+        
+        if sample!=[]:
+            step = float(len(xmins))/(sample-1)
+            index_curr=0
+            new_xmins=[]
+            for i in range (0,sample):
+                if round(index_curr)==len(xmins): index_curr-=1
+                new_xmins.append(xmins[int(round(index_curr))])
+                index_curr+=step
+            xmins = unique(new_xmins)
+            xmins.sort()
+        
+        if xmins==[]:
+            print '(PLFIT) Error: x must contain at least two unique values.\n'
+            alpha = 'Not a Number'
+            xmin = x[0]
+            D = 'Not a Number'
+            return [alpha,xmin,D]
+        
+        xmax   = max(x)
+        
+        z      = x
+        z.sort()
+        datA=[]
+        datB=[]
+
+        for xm in range(0,len(xmins)):
+            xmin = xmins[xm]
+            z    = filter(lambda X:X>=xmin,z)
+            n    = len(z)
+            # estimate alpha via direct maximization of likelihood function
+
+            # force iterative calculation 
+            L       = []
+            slogz   = sum(map(log,z))
+            xminvec = map(float,range(1,xmin))
+            for k in range(0,len(vec)):
+                L.append(-vec[k]*float(slogz) - float(n)*log(float(zvec[k]) - sum(map(lambda X:pow(float(X),-vec[k]),xminvec))))
+            
+            
+            I = L.index(max(L))
+            # compute KS statistic
+            fit = reduce(lambda X,Y: X+[Y+X[-1]],\
+                         (map(lambda X: pow(X,-vec[I])/(float(zvec[I])-sum(map(lambda X: pow(X,-vec[I]),map(float,range(1,xmin))))),range(xmin,xmax+1))),[0])[1:]
+            cdi=[]
+            for XM in range(xmin,xmax+1):
+                cdi.append(len(filter(lambda X: floor(X)<=XM,z))/float(n))
+            
+            datA.append(max( map(lambda X: abs(fit[X] - cdi[X]),range(0,xmax-xmin+1))))
+            datB.append(vec[I])
+        # select the index for the minimum value of D
+        I = datA.index(min(datA))
+        xmin  = xmins[I]
+        z     = filter(lambda X:X>=xmin,x)
+        n     = len(z)
+        alpha = datB[I]
+        if finite: alpha = alpha*(n-1.)/n+1./n  # finite-size correction
+        if n < 50 and not finite and not nowarn:
+            print '(PLFIT) Warning: finite-size bias may be present.\n'
+        
+        L     = -alpha*sum(map(log,z)) - n*log(zvec[vec.index(max(filter(lambda X:X<=alpha,vec)))] - \
+                                              sum(map(lambda X: pow(X,-alpha),range(1,xmin))))
+    else:
+        print '(PLFIT) Error: x must contain only reals or only integers.\n'
+        alpha = []
+        xmin  = []
+        L     = []
+
+    return [alpha,xmin,L]
+
+
+# helper functions (unique and zeta)
+
+
+def unique(seq): 
+    # not order preserving 
+    set = {} 
+    map(set.__setitem__, seq, []) 
+    return set.keys()
+
+def _polyval(coeffs, x):
+    p = coeffs[0]
+    for c in coeffs[1:]:
+        p = c + x*p
+    return p
+
+_zeta_int = [\
+-0.5,
+0.0,
+1.6449340668482264365,1.2020569031595942854,1.0823232337111381915,
+1.0369277551433699263,1.0173430619844491397,1.0083492773819228268,
+1.0040773561979443394,1.0020083928260822144,1.0009945751278180853,
+1.0004941886041194646,1.0002460865533080483,1.0001227133475784891,
+1.0000612481350587048,1.0000305882363070205,1.0000152822594086519,
+1.0000076371976378998,1.0000038172932649998,1.0000019082127165539,
+1.0000009539620338728,1.0000004769329867878,1.0000002384505027277,
+1.0000001192199259653,1.0000000596081890513,1.0000000298035035147,
+1.0000000149015548284]
+
+_zeta_P = [-3.50000000087575873, -0.701274355654678147,
+-0.0672313458590012612, -0.00398731457954257841,
+-0.000160948723019303141, -4.67633010038383371e-6,
+-1.02078104417700585e-7, -1.68030037095896287e-9,
+-1.85231868742346722e-11][::-1]
+
+_zeta_Q = [1.00000000000000000, -0.936552848762465319,
+-0.0588835413263763741, -0.00441498861482948666,
+-0.000143416758067432622, -5.10691659585090782e-6,
+-9.58813053268913799e-8, -1.72963791443181972e-9,
+-1.83527919681474132e-11][::-1]
+
+_zeta_1 = [3.03768838606128127e-10, -1.21924525236601262e-8,
+2.01201845887608893e-7, -1.53917240683468381e-6,
+-5.09890411005967954e-7, 0.000122464707271619326,
+-0.000905721539353130232, -0.00239315326074843037,
+0.084239750013159168, 0.418938517907442414, 0.500000001921884009]
+
+_zeta_0 = [-3.46092485016748794e-10, -6.42610089468292485e-9,
+1.76409071536679773e-7, -1.47141263991560698e-6, -6.38880222546167613e-7,
+0.000122641099800668209, -0.000905894913516772796, -0.00239303348507992713,
+0.0842396947501199816, 0.418938533204660256, 0.500000000000000052]
+
+def zeta(s):
+    """
+    Riemann zeta function, real argument
+    """
+    if not isinstance(s, (float, int)):
+        try:
+            s = float(s)
+        except (ValueError, TypeError):
+            try:
+                s = complex(s)
+                if not s.imag:
+                    return complex(zeta(s.real))
+            except (ValueError, TypeError):
+                pass
+            raise NotImplementedError
+    if s == 1:
+        raise ValueError("zeta(1) pole")
+    if s >= 27:
+        return 1.0 + 2.0**(-s) + 3.0**(-s)
+    n = int(s)
+    if n == s:
+        if n >= 0:
+            return _zeta_int[n]
+        if not (n % 2):
+            return 0.0
+    if s <= 0.0:
+        return 0
+    if s <= 2.0:
+        if s <= 1.0:
+            return _polyval(_zeta_0,s)/(s-1)
+        return _polyval(_zeta_1,s)/(s-1)
+    z = _polyval(_zeta_P,s) / _polyval(_zeta_Q,s)
+    return 1.0 + 2.0**(-s) + 3.0**(-s) + 4.0**(-s)*z
+
+
+
+    
Index: trunk/UTIL/PYTHON/powerlaw/previous_tests/plot_granulo.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/previous_tests/plot_granulo.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/previous_tests/plot_granulo.py	(revision 1053)
@@ -0,0 +1,96 @@
+#! /usr/bin/env python
+
+import numpy as np
+from scipy import ndimage
+import matplotlib.pyplot as plt
+
+from myplot import getfield
+from netCDF4 import Dataset
+from scipy import ndimage
+
+def disk_structure(n):
+    struct = np.zeros((2 * n + 1, 2 * n + 1))
+    x, y = np.indices((2 * n + 1, 2 * n + 1))
+    mask = (x - n)**2 + (y - n)**2 <= n**2
+    struct[mask] = 1
+    return struct.astype(np.bool)
+
+
+def granulometry(data, sizes=None):
+    s = max(data.shape)
+    if sizes == None:
+        sizes = range(1, s/2, 2)
+    granulo = [ndimage.binary_opening(data, \
+            structure=disk_structure(n)).sum() for n in sizes]
+    return granulo
+
+
+#np.random.seed(1)
+#n = 10
+#l = 256
+#im = np.zeros((l, l))
+#points = l*np.random.random((2, n**2))
+#im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
+#im = ndimage.gaussian_filter(im, sigma=l/(4.*n))
+#mask = im > im.mean()
+
+namefile = "/home/aymeric/Big_Data/psfc_f18.nc"
+nc  = Dataset(namefile)
+psfc_full=getfield(nc,'PSFC')
+
+#mask = psfc_full[250,100:200,100:200]
+mask = psfc_full[250,:,:]
+mask = ndimage.laplace(mask)
+#mask = mask > 0.3
+
+print mask.std()
+print mask.mean()
+
+lim = mask.std() * 12.
+
+ind1 = np.where(mask > lim)
+ind2 = np.where(mask < lim)
+mask[ind1] = 1
+mask[ind2] = 0
+
+print mask
+
+
+labeled_array, num_features = ndimage.label(mask)
+
+print num_features
+
+plt.pcolor(labeled_array)
+plt.show()
+exit()
+
+
+
+
+mask = mask < mask.mean() - 0.3
+
+#mask = mask > 0.3
+
+
+#plt.pcolor(mask)
+#plt.show()
+
+granulo = granulometry(mask, sizes=[2,5,20])
+
+print granulo
+
+plt.figure(figsize=(6, 2.2))
+
+plt.subplot(121)
+plt.imshow(mask, cmap=plt.cm.gray)
+opened = ndimage.binary_opening(mask, structure=disk_structure(10))
+opened_more = ndimage.binary_opening(mask, structure=disk_structure(14))
+plt.contour(opened, [0.5], colors='b', linewidths=2)
+plt.contour(opened_more, [0.5], colors='r', linewidths=2)
+plt.axis('off')
+plt.subplot(122)
+#plt.plot(np.arange(2, 19, 4), granulo, 'ok', ms=8)
+
+
+plt.subplots_adjust(wspace=0.02, hspace=0.15, top=0.95, bottom=0.15, left=0, right=0.95)
+plt.show()
Index: trunk/UTIL/PYTHON/powerlaw/previous_tests/plplot.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/previous_tests/plplot.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/previous_tests/plplot.py	(revision 1053)
@@ -0,0 +1,201 @@
+import matplotlib.pyplot as plt
+from math import *
+
+# function h=plplot(x, xmin, alpha)
+# PLPLOT visualizes a power-law distributional model with empirical data.
+#    Source: http://www.santafe.edu/~aaronc/powerlaws/
+#
+#    PLPLOT  REQUIRES  the use of the free library: matplotlib
+#
+#    PLPLOT(x, xmin, alpha) plots (on log axes) the data contained in x 
+#    and a power-law distribution of the form p(x) ~ x^-alpha for 
+#    x >= xmin. For additional customization, PLPLOT returns a pair of 
+#    handles, one to the empirical and one to the fitted data series. By 
+#    default, the empirical data is plotted as 'bo' and the fitted form is
+#    plotted as 'k--'. PLPLOT automatically detects whether x is composed 
+#    of real or integer values, and applies the appropriate plotting 
+#    method. For discrete data, if min(x) > 50, PLFIT uses the continuous 
+#    approximation, which is a reliable in this regime.
+#
+#    Example:
+#       xmin  = 47;
+#       alpha = 2.71;
+#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
+#       h = plplot(x,xmin,alpha);
+#
+#    For more information, try 'type plplot'
+#
+#    See also PLFIT, PLVAR, PLPVA
+
+# Version 1.0   (2008 February)
+# Ported to python by Joel Ornstein(2011 July)
+# (joel_ornstein@hmc.edu)
+
+# Copyright (C) 2008-2011 Aaron Clauset (Santa Fe Institute)
+# Distributed under GPL 2.0
+# http://www.gnu.org/copyleft/gpl.html
+# PLFIT comes with ABSOLUTELY NO WARRANTY
+# 
+# The 'zeta' helper function is modified from the open-source library 'mpmath'
+#   mpmath: a Python library for arbitrary-precision floating-point arithmetic
+#   http://code.google.com/p/mpmath/
+#   version 0.17 (February 2011) by Fredrik Johansson and others
+# 
+# No Notes
+#
+
+def plplot(x,xmin,alpha):
+        # select method (discrete or continuous) for fitting
+    if     reduce(lambda X,Y:X==True and floor(Y)==float(Y),x,True): f_dattype = 'INTS'
+    elif reduce(lambda X,Y:X==True and (type(Y)==int or type(Y)==float or type(Y)==long),x,True):    f_dattype = 'REAL'
+    else:                 f_dattype = 'UNKN'
+    
+    if f_dattype=='INTS' and min(x) > 1000 and len(x)>100:
+        f_dattype = 'REAL'
+    plt.close()
+    plt.ion()
+    h=[[],[]]
+    # estimate xmin and alpha, accordingly
+    if f_dattype== 'REAL':
+        n = len(x)
+        c1 = sorted(x)
+        c2 = map(lambda X:X/float(n),range(n,0,-1))
+        q = sorted(filter(lambda X:X>=xmin,x))
+        cf = map(lambda X:pow(float(X)/xmin,1.-alpha),q)
+        cf = map(lambda X:X*float(c2[c1.index(q[0])]),cf)
+
+        h[0]=plt.loglog(c1, c2, 'bo',markersize=8,markerfacecolor=[1,1,1],markeredgecolor=[0,0,1])
+        h[1]=plt.loglog(q, cf, 'k--',linewidth=2)
+        
+        xr1 = pow(10,floor(log(min(x),10)))
+        xr2 = pow(10,ceil(log(min(x),10)))
+        yr1 = pow(10,floor(log(1./n,10)))
+        yr2 = 1
+        
+
+        plt.axhspan(ymin=yr1,ymax=yr2,xmin=xr1,xmax=xr2)
+        plt.ylabel('Pr(X >= x)',fontsize=16);
+        plt.xlabel('x',fontsize=16)
+        plt.draw()
+        
+    elif f_dattype== 'INTS':
+        n = len(x)
+        q = sorted(unique(x))
+        c=[]
+        for Q in q:
+            c.append(len(filter(lambda X: floor(X)==Q,x))/float(n))
+        c1 = q+[q[-1]+1]
+        c2 = map(lambda Z: 1.-Z,reduce(lambda X,Y: X+[Y+X[-1]],c,[0]))
+        c2 = filter(lambda X:float(X)>=pow(10,-10.),c2)
+        c1 = c1[0:len(c2)]
+        cf = map(lambda X:pow(X,-alpha)/(float(zeta(alpha)) - sum(map(lambda Y:pow(Y,-alpha),range(1,xmin)))),range(xmin,q[-1]+1))
+        cf1 = range(xmin,q[-1]+2)
+        cf2 = map(lambda Z: 1.-Z,reduce(lambda X,Y: X+[Y+X[-1]],cf,[0]))
+        cf2 = map(lambda X: X*float(c2[c1.index(xmin)]),cf2)
+
+        h[0]=plt.loglog(c1, c2, 'bo',markersize=8,markerfacecolor=[1,1,1],markeredgecolor=[0,0,1])
+        h[1]=plt.loglog(cf1, cf2, 'k--',linewidth=2)
+        
+        xr1 = pow(10,floor(log(min(x),10)))
+        xr2 = pow(10,ceil(log(min(x),10)))
+        yr1 = pow(10,floor(log(1./n,10)))
+        yr2 = 1
+
+
+        plt.axhspan(ymin=yr1,ymax=yr2,xmin=xr1,xmax=xr2)
+        plt.ylabel('Pr(X >= x)',fontsize=16);
+        plt.xlabel('x',fontsize=16)
+        plt.draw()
+                 
+          
+
+    return h
+
+# helper functions (unique and zeta)
+
+
+def unique(seq): 
+    # not order preserving 
+    set = {} 
+    map(set.__setitem__, seq, []) 
+    return set.keys()
+
+def _polyval(coeffs, x):
+    p = coeffs[0]
+    for c in coeffs[1:]:
+        p = c + x*p
+    return p
+
+_zeta_int = [\
+-0.5,
+0.0,
+1.6449340668482264365,1.2020569031595942854,1.0823232337111381915,
+1.0369277551433699263,1.0173430619844491397,1.0083492773819228268,
+1.0040773561979443394,1.0020083928260822144,1.0009945751278180853,
+1.0004941886041194646,1.0002460865533080483,1.0001227133475784891,
+1.0000612481350587048,1.0000305882363070205,1.0000152822594086519,
+1.0000076371976378998,1.0000038172932649998,1.0000019082127165539,
+1.0000009539620338728,1.0000004769329867878,1.0000002384505027277,
+1.0000001192199259653,1.0000000596081890513,1.0000000298035035147,
+1.0000000149015548284]
+
+_zeta_P = [-3.50000000087575873, -0.701274355654678147,
+-0.0672313458590012612, -0.00398731457954257841,
+-0.000160948723019303141, -4.67633010038383371e-6,
+-1.02078104417700585e-7, -1.68030037095896287e-9,
+-1.85231868742346722e-11][::-1]
+
+_zeta_Q = [1.00000000000000000, -0.936552848762465319,
+-0.0588835413263763741, -0.00441498861482948666,
+-0.000143416758067432622, -5.10691659585090782e-6,
+-9.58813053268913799e-8, -1.72963791443181972e-9,
+-1.83527919681474132e-11][::-1]
+
+_zeta_1 = [3.03768838606128127e-10, -1.21924525236601262e-8,
+2.01201845887608893e-7, -1.53917240683468381e-6,
+-5.09890411005967954e-7, 0.000122464707271619326,
+-0.000905721539353130232, -0.00239315326074843037,
+0.084239750013159168, 0.418938517907442414, 0.500000001921884009]
+
+_zeta_0 = [-3.46092485016748794e-10, -6.42610089468292485e-9,
+1.76409071536679773e-7, -1.47141263991560698e-6, -6.38880222546167613e-7,
+0.000122641099800668209, -0.000905894913516772796, -0.00239303348507992713,
+0.0842396947501199816, 0.418938533204660256, 0.500000000000000052]
+
+def zeta(s):
+    """
+    Riemann zeta function, real argument
+    """
+    if not isinstance(s, (float, int)):
+        try:
+            s = float(s)
+        except (ValueError, TypeError):
+            try:
+                s = complex(s)
+                if not s.imag:
+                    return complex(zeta(s.real))
+            except (ValueError, TypeError):
+                pass
+            raise NotImplementedError
+    if s == 1:
+        raise ValueError("zeta(1) pole")
+    if s >= 27:
+        return 1.0 + 2.0**(-s) + 3.0**(-s)
+    n = int(s)
+    if n == s:
+        if n >= 0:
+            return _zeta_int[n]
+        if not (n % 2):
+            return 0.0
+    if s <= 0.0:
+        return 0
+    if s <= 2.0:
+        if s <= 1.0:
+            return _polyval(_zeta_0,s)/(s-1)
+        return _polyval(_zeta_1,s)/(s-1)
+    z = _polyval(_zeta_P,s) / _polyval(_zeta_Q,s)
+    return 1.0 + 2.0**(-s) + 3.0**(-s) + 4.0**(-s)*z
+
+
+
+    
Index: trunk/UTIL/PYTHON/powerlaw/previous_tests/plpva.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/previous_tests/plpva.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/previous_tests/plpva.py	(revision 1053)
@@ -0,0 +1,476 @@
+from math import *
+from random import *
+
+# function [alpha, xmin, n]=plpva(x, xmin, varargin)
+# PLPVA calculates the p-value for the given power-law fit to some data.
+#    Source: http://www.santafe.edu/~aaronc/powerlaws/
+# 
+#    PLPVA(x, xmin) takes data x and given lower cutoff for the power-law
+#    behavior xmin and computes the corresponding p-value for the
+#    Kolmogorov-Smirnov test, according to the method described in 
+#    Clauset, Shalizi, Newman (2007).
+#    PLPVA automatically detects whether x is composed of real or integer
+#    values, and applies the appropriate method. For discrete data, if
+#    min(x) > 1000, PLPVA uses the continuous approximation, which is 
+#    a reliable in this regime.
+#   
+#    The fitting procedure works as follows:
+#    1) For each possible choice of x_min, we estimate alpha via the 
+#       method of maximum likelihood, and calculate the Kolmogorov-Smirnov
+#       goodness-of-fit statistic D.
+#    2) We then select as our estimate of x_min, the value that gives the
+#       minimum value D over all values of x_min.
+#
+#    Note that this procedure gives no estimate of the uncertainty of the 
+#    fitted parameters, nor of the validity of the fit.
+#
+#    Example:
+#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
+#       [p, gof] = plpva(x, 1);
+#
+#    For more information, try 'type plpva'
+#
+#    See also PLFIT, PLVAR
+
+# Version 1.0.7 (2009 October)
+# Copyright (C) 2008-2011 Aaron Clauset (Santa Fe Institute)
+
+# Ported to Python by Joel Ornstein (2011 July)
+#(joel_ornstein@hmc.edu)
+
+# Distributed under GPL 2.0
+# http://www.gnu.org/copyleft/gpl.html
+# PLPVA comes with ABSOLUTELY NO WARRANTY
+#
+# 
+# The 'zeta' helper function is modified from the open-source library 'mpmath'
+#   mpmath: a Python library for arbitrary-precision floating-point arithmetic
+#   http://code.google.com/p/mpmath/
+#   version 0.17 (February 2011) by Fredrik Johansson and others
+# 
+
+# Notes:
+# 
+# 1. In order to implement the integer-based methods in Matlab, the numeric
+#    maximization of the log-likelihood function was used. This requires
+#    that we specify the range of scaling parameters considered. We set
+#    this range to be 1.50 to 3.50 at 0.01 intervals by default. 
+#    This range can be set by the user like so,
+#
+#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
+#       a = plpva(x,1,'range',[1.50,3.50,0.01])
+#    
+# 2. PLPVA can be told to limit the range of values considered as estimates
+#    for xmin in two ways. First, it can be instructed to sample these
+#    possible values like so,
+#    
+#       a = plpva(x,1,'sample',100);
+#    
+#    which uses 100 uniformly distributed values on the sorted list of
+#    unique values in the data set. Second, it can simply omit all
+#    candidates above a hard limit, like so
+#    
+#       a = plpva(x,1,'limit',3.4);
+#    
+#    Finally, it can be forced to use a fixed value, like so
+#    
+#       a = plpva(x,1,'xmin',1);
+#    
+#    In the case of discrete data, it rounds the limit to the nearest
+#    integer.
+# 
+# 3. The default number of semiparametric repetitions of the fitting
+# procedure is 1000. This number can be changed like so
+#    
+#       p = plpva(x, 1,'reps',10000);
+# 
+# 4. To silence the textual output to the screen, do this
+#    
+#       p = plpva(x, 1,'reps',10000,'silent');
+# 
+
+def plpva(x,xmin, *varargin):
+    vec     = []
+    sample  = []
+    xminx   = []
+    limit   = []
+    Bt      = []
+    quiet   = False
+    
+
+    # parse command-line parameters trap for bad input
+    i=0 
+    while i<len(varargin): 
+        argok = 1 
+        if type(varargin[i])==str: 
+            if varargin[i]=='range':
+                Range = varargin[i+1]
+                if Range[1]>Range[0]:
+                    argok=0
+                    vec=[]
+                try:
+                    vec=map(lambda X:X*float(Range[2])+Range[0],\
+                            range(int((Range[1]-Range[0])/Range[2])))
+                    
+                    
+                except:
+                    argok=0
+                    vec=[]
+                    
+
+                if Range[0]>=Range[1]:
+                    argok=0
+                    vec=[]
+                    i-=1
+
+                i+=1
+                    
+                            
+            elif varargin[i]== 'sample':
+                sample  = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'limit':
+                limit   = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'xmin':
+                xminx   = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'reps':
+                Bt   = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'silent':       quiet  = True    
+ 
+            else: argok=0 
+        
+      
+        if not argok:
+            print '(PLPVA) Ignoring invalid argument #',i+1 
+      
+        i = i+1 
+
+    if vec!=[] and (type(vec)!=list or min(vec)<=1):
+        print '(PLPVA) Error: ''range'' argument must contain a vector or minimum <= 1. using default.\n'                        
+              
+        vec = []
+    
+    if sample!=[] and sample<2:
+        print'(PLPVA) Error: ''sample'' argument must be a positive integer > 1. using default.\n'
+        sample = []
+    
+    if limit!=[] and limit<min(x):
+        print'(PLPVA) Error: ''limit'' argument must be a positive value >= 1. using default.\n'
+        limit = []
+    
+    if xminx!=[] and xminx>=max(x):
+        print'(PLPVA) Error: ''xmin'' argument must be a positive value < max(x). using default behavior.\n'
+        xminx = []
+    if Bt!=[] and Bt<2:
+        print '(PLPVA) Error: ''reps'' argument must be a positive value > 1; using default.\n'
+        Bt = [];
+
+
+    # select method (discrete or continuous) for fitting
+    if     reduce(lambda X,Y:X==True and floor(Y)==float(Y),x,True): f_dattype = 'INTS'
+    elif reduce(lambda X,Y:X==True and (type(Y)==int or type(Y)==float or type(Y)==long),x,True):    f_dattype = 'REAL'
+    else:                 f_dattype = 'UNKN'
+    
+    if f_dattype=='INTS' and min(x) > 1000 and len(x)>100:
+        f_dattype = 'REAL'
+    
+    N=len(x)
+    
+    if Bt==[]: Bt=1000
+    nof = []
+    if not quiet:
+        print 'Power-law Distribution, p-value calculation\n'
+        print '   Copyright 2007-2010 Aaron Clauset\n'
+        print '   Warning: This can be a slow calculation; please be patient.\n'
+        print '   n    = ',len(x),'\n   xmin = ',xmin,'\n   reps = ',Bt
+    # estimate xmin and alpha, accordingly
+    if f_dattype== 'REAL':
+        # compute D for the empirical distribution
+        z = filter(lambda X:X>=xmin,x)
+        z.sort()
+        nz = len(z)
+        y = filter(lambda X:X<xmin,x)
+        ny = len(y)
+        alpha = 1 + float(nz) / sum(map(lambda X: log(float(X)/xmin),z))
+        cz    = map(lambda X:X/float(nz),range(0,nz))
+        cf    = map(lambda X:1-pow((xmin/float(X)),(alpha-1)),z)
+        gof   = max( map(lambda X:abs(cz[X] - cf[X]),range(0,nz)))
+        pz    = nz/float(N)
+        
+        # compute distribution of gofs from semi-parametric bootstrap
+        # of entire data set with fit
+        for B in range(0,Bt):
+            # semi-parametric bootstrap of data
+            n1=0
+            for i in range(0,N):
+                if random()>pz: n1+=1
+
+            q1=[]
+            for i in range(0,n1):
+                q1.append(y[int(floor(ny*random()))])
+            n2 = N-n1;
+            q2=[]
+            for i in range (0,n2):
+                q2.append(random())
+            q2 = map(lambda X:xmin*pow(1.-X,(-1./(alpha-1.))),q2)
+
+            q  = q1+q2
+            q.sort()
+
+            
+
+            # estimate xmin and alpha via GoF-method
+            qmins = unique(q);
+            qmins.sort()
+            qmins = qmins[0:-1];
+            
+            if xminx!=[]:
+                
+                qmins = [min(filter(lambda X: X>=xminx,qmins))]
+                
+            
+            if limit!=[]:
+                qmins=filter(lambda X: X<=limit,qmins)
+                if qmins==[]: qmins=[min(q)]
+                
+            if sample!=[]:
+                step = float(len(qmins))/(sample-1)
+                index_curr=0
+                new_qmins=[]
+                for i in range (0,sample):
+                    if round(index_curr)==len(qmins): index_curr-=1
+                    new_qmins.append(qmins[int(round(index_curr))])
+                    index_curr+=step
+                qmins = unique(new_qmins)
+                qmins.sort()
+            
+            
+        
+            dat   = []
+            
+            for qm in range(0,len(qmins)):
+                qmin = qmins[qm]
+                zq   = filter(lambda X:X>=qmin,q)
+                nq   = len(zq)
+               
+                a    = float(nq)/sum(map(lambda X: log(float(X)/qmin),zq))
+                cf   = map(lambda X:1-pow((float(qmin)/X),a),zq)
+                
+                dat.append( max( map(lambda X: abs(cf[X]-float(X)/nq),range(0,nq))))
+            if not quiet:
+                print '['+str(B+1)+']\tp = ',round(len(filter(lambda X:X>=gof,nof))/float(B+1),3)
+            nof.append(min(dat))
+        p=len(filter(lambda X:X>=gof,nof))/float(Bt+1) 
+
+    elif f_dattype== 'INTS':
+        x=map(int,x)
+        if vec==[]:
+            for X in range(150,351):
+                vec.append(X/100.)    # covers range of most practical 
+                                    # scaling parameters
+        zvec = map(zeta, vec)
+        z = filter (lambda X:X>=xmin,x)
+        nz = len(z)
+        y = filter (lambda X:X<xmin,x)
+        ny = len(y)
+        xmax=max(z)
+        L       = []
+        for k in range(0,len(vec)):
+            L.append(-vec[k]*float(sum(map(log,z))) - float(nz)*log(float(zvec[k]) - sum(map(lambda X:pow(float(X),-vec[k]),range(1,xmin)))))
+        I = L.index(max(L))
+        alpha=vec[I]
+        fit = reduce(lambda X,Y: X+[Y+X[-1]],\
+                     (map(lambda X: pow(X,-alpha)/(float(zvec[I])-sum(map(lambda X: pow(X,-alpha),map(float,range(1,xmin))))),range(xmin,xmax+1))),[0])[1:]
+        cdi=[]
+        for XM in range(xmin,xmax+1):
+            cdi.append(len(filter(lambda X: floor(X)<=XM,z))/float(nz))
+        
+        gof= max( map(lambda X: abs(fit[X] - cdi[X]),range(0,xmax-xmin+1)))
+        pz=nz/float(N)
+        mmax = 20*xmax
+        pdf=[]
+        for i in range(0,xmin):
+            pdf.append(0)
+        pdf += map(lambda X: pow(X,-alpha)/(float(zvec[I])-sum(map(lambda X: pow(X,-alpha),map(float,range(1,xmin))))),range(xmin,mmax+1))
+        cdf = reduce(lambda X,Y: X+[Y+X[-1]],pdf,[0])[1:]+[1]
+
+        # compute distribution of gofs from semi-parametric bootstrap
+        # of entire data set with fit
+        for B in range(0,Bt):
+            # semi-parametric bootstrap of data
+            n1=0
+            for i in range(0,N):
+                if random()>pz: n1+=1
+
+            q1=[]
+            for i in range(0,n1):
+                q1.append( y[int(floor(ny*random()))])
+            n2 = N-n1;
+            r2=[]
+            for i in range (0,n2):
+                r2.append(random())
+            r2.sort()
+            c=0
+            k=0
+            q2=[]
+            for i in range(xmin,mmax+2):
+                while c<n2 and r2[c]<=cdf[i]:c+=1
+                for k in range(k,c):q2.append(i)
+                k=c
+                if k>=n2: break
+
+            q  = q1+q2
+            qmins = unique(q)
+            qmins.sort()
+            qmins = qmins[0:-1]
+
+            if xminx!=[]:
+                
+                qmins = [min(filter(lambda X: X>=xminx,qmins))]
+                
+            
+            if limit!=[]:
+                qmins=filter(lambda X: X<=limit,qmins)
+                if qmins==[]: qmins=[min(q)]
+                
+            if sample!=[]:
+                step = float(len(qmins))/(sample-1)
+                index_curr=0
+                new_qmins=[]
+                for i in range (0,sample):
+                    if round(index_curr)==len(qmins): index_curr-=1
+                    new_qmins.append(qmins[int(round(index_curr))])
+                    index_curr+=step
+                qmins = unique(new_qmins)
+                qmins.sort()
+        
+            qmax   = max(q)
+            dat=[]
+            zq=q
+
+
+            for qm in range(0,len(qmins)):
+                qmin = qmins[qm]
+                zq    = filter(lambda X:X>=qmin,zq)
+                nq    = len(zq)
+                if nq>1:
+                    L=[]
+                    slogzq = sum(map(log,zq))
+                    qminvec = range (1,qmin)
+                    for k in range (1,len(vec)):
+                        L.append(-vec[k]*float(slogzq) - float(nq)*log(float(zvec[k]) - sum(map(lambda X:pow(float(X),-vec[k]),qminvec))))
+            
+                    I = L.index(max(L))
+                    fit = reduce(lambda X,Y: X+[Y+X[-1]],\
+                                 (map(lambda X: pow(X,-vec[I])/(float(zvec[I])-sum(map(lambda X: pow(X,-vec[I]),map(float,range(1,qmin))))),range(qmin,qmax+1))),[0])[1:]
+                    cdi=[]
+                    for QM in range(qmin,qmax+1):
+                        cdi.append(len(filter(lambda X: floor(X)<=QM,zq))/float(nq))
+                    
+                    dat.append(max( map(lambda X: abs(fit[X] - cdi[X]),range(0,qmax-qmin+1))))
+                else: dat[qm]=float('-inf')
+            if not quiet:
+                print '['+str(B+1)+']\tp = ',round(len(filter(lambda X:X>=gof,nof))/float(B+1),3)
+            nof.append(min(dat))
+        
+        p=len(filter(lambda X:X>=gof,nof))/float(Bt)
+    else:
+        print '(PLPVA) Error: x must contain only reals or only integers.\n'
+        p = []
+        gof  = []
+        
+
+    return [p,gof]
+
+
+# helper functions (unique and zeta)
+
+
+def unique(seq): 
+    # not order preserving 
+    set = {} 
+    map(set.__setitem__, seq, []) 
+    return set.keys()
+
+def _polyval(coeffs, x):
+    p = coeffs[0]
+    for c in coeffs[1:]:
+        p = c + x*p
+    return p
+
+_zeta_int = [\
+-0.5,
+0.0,
+1.6449340668482264365,1.2020569031595942854,1.0823232337111381915,
+1.0369277551433699263,1.0173430619844491397,1.0083492773819228268,
+1.0040773561979443394,1.0020083928260822144,1.0009945751278180853,
+1.0004941886041194646,1.0002460865533080483,1.0001227133475784891,
+1.0000612481350587048,1.0000305882363070205,1.0000152822594086519,
+1.0000076371976378998,1.0000038172932649998,1.0000019082127165539,
+1.0000009539620338728,1.0000004769329867878,1.0000002384505027277,
+1.0000001192199259653,1.0000000596081890513,1.0000000298035035147,
+1.0000000149015548284]
+
+_zeta_P = [-3.50000000087575873, -0.701274355654678147,
+-0.0672313458590012612, -0.00398731457954257841,
+-0.000160948723019303141, -4.67633010038383371e-6,
+-1.02078104417700585e-7, -1.68030037095896287e-9,
+-1.85231868742346722e-11][::-1]
+
+_zeta_Q = [1.00000000000000000, -0.936552848762465319,
+-0.0588835413263763741, -0.00441498861482948666,
+-0.000143416758067432622, -5.10691659585090782e-6,
+-9.58813053268913799e-8, -1.72963791443181972e-9,
+-1.83527919681474132e-11][::-1]
+
+_zeta_1 = [3.03768838606128127e-10, -1.21924525236601262e-8,
+2.01201845887608893e-7, -1.53917240683468381e-6,
+-5.09890411005967954e-7, 0.000122464707271619326,
+-0.000905721539353130232, -0.00239315326074843037,
+0.084239750013159168, 0.418938517907442414, 0.500000001921884009]
+
+_zeta_0 = [-3.46092485016748794e-10, -6.42610089468292485e-9,
+1.76409071536679773e-7, -1.47141263991560698e-6, -6.38880222546167613e-7,
+0.000122641099800668209, -0.000905894913516772796, -0.00239303348507992713,
+0.0842396947501199816, 0.418938533204660256, 0.500000000000000052]
+
+def zeta(s):
+    """
+    Riemann zeta function, real argument
+    """
+    if not isinstance(s, (float, int)):
+        try:
+            s = float(s)
+        except (ValueError, TypeError):
+            try:
+                s = complex(s)
+                if not s.imag:
+                    return complex(zeta(s.real))
+            except (ValueError, TypeError):
+                pass
+            raise NotImplementedError
+    if s == 1:
+        raise ValueError("zeta(1) pole")
+    if s >= 27:
+        return 1.0 + 2.0**(-s) + 3.0**(-s)
+    n = int(s)
+    if n == s:
+        if n >= 0:
+            return _zeta_int[n]
+        if not (n % 2):
+            return 0.0
+    if s <= 0.0:
+        return 0
+    if s <= 2.0:
+        if s <= 1.0:
+            return _polyval(_zeta_0,s)/(s-1)
+        return _polyval(_zeta_1,s)/(s-1)
+    z = _polyval(_zeta_P,s) / _polyval(_zeta_Q,s)
+    return 1.0 + 2.0**(-s) + 3.0**(-s) + 4.0**(-s)*z
+
+
+
+    
Index: trunk/UTIL/PYTHON/powerlaw/previous_tests/plvar.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/previous_tests/plvar.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/previous_tests/plvar.py	(revision 1053)
@@ -0,0 +1,444 @@
+from math import *
+from random import *
+
+# function [alpha, xmin, n]=plvar(x, varargin)
+# PLVAR estimates the uncertainty in the estimated power-law parameters.
+#    Source: http://www.santafe.edu/~aaronc/powerlaws/
+# 
+#    PLVAR(x) takes a vector of observations x and returns estimated
+#    uncertainties in the estimated power-law parameters, based on the
+#    nonparametric approach described in Clauset, Shalizi, Newman (2007). 
+#    PLVAR automatically detects whether x is composed of real or integer
+#    values, and applies the appropriate method. For discrete data, if
+#    min(x) > 1000, PLVAR uses the continuous approximation, which is 
+#    a reliable in this regime.
+#   
+#    The fitting procedure works as follows:
+#    1) For each possible choice of x_min, we estimate alpha via the 
+#       method of maximum likelihood, and calculate the Kolmogorov-Smirnov
+#       goodness-of-fit statistic D.
+#    2) We then select as our estimate of x_min, the value that gives the
+#       minimum value D over all values of x_min.
+#
+#    Note that this procedure gives no estimate of the validity of the fit.
+#
+#    Example:
+#
+#       x = [500,150,90,81,75,75,70,65,60,58,49,47,40]
+#       [alpha, xmin, ntail] = plvar(x);
+#
+#    For more information, try 'type plvar'
+#
+#    See also PLFIT, PLPVA
+
+
+# Version 1.0.8 (2010 April)
+# Copyright (C) 2008-2011 Aaron Clauset (Santa Fe Institute)
+
+# Ported to Python by Joel Ornstein (2011 July)
+#(joel_ornstein@hmc.edu)
+
+# Distributed under GPL 2.0
+# http://www.gnu.org/copyleft/gpl.html
+# PLVAR comes with ABSOLUTELY NO WARRANTY
+#
+# 
+# The 'zeta' helper function is modified from the open-source library 'mpmath'
+#   mpmath: a Python library for arbitrary-precision floating-point arithmetic
+#   http://code.google.com/p/mpmath/
+#   version 0.17 (February 2011) by Fredrik Johansson and others
+# 
+
+# Notes:
+# 
+# 1. In order to implement the integer-based methods in Matlab, the numeric
+#    maximization of the log-likelihood function was used. This requires
+#    that we specify the range of scaling parameters considered. We set
+#    this range to be 1.50 to 3.50 at 0.01 intervals by default. 
+#    This range can be set by the user like so,
+#    
+#       a = plvar(x,'range',[1.50,3.50,0.01])
+#    
+#    
+# 2. PLVAR can be told to limit the range of values considered as estimates
+#    for xmin in three ways. First, it can be instructed to sample these
+#    possible values like so,
+#    
+#       a = plvar(x,'sample',100);
+#    
+#    which uses 100 uniformly distributed values on the sorted list of
+#    unique values in the data set. Second, it can simply omit all
+#    candidates above a hard limit, like so
+#    
+#       a = plvar(x,'limit',3.4);
+#    
+#    Finally, it can be forced to use a fixed value, like so
+#    
+#       a = plvar(x,'xmin',3.4);
+#    
+#    In the case of discrete data, it rounds the limit to the nearest
+#    integer.
+# 
+# 3. The default number of nonparametric repetitions of the fitting
+# procedure is 1000. This number can be changed like so
+#    
+#       a = plvar(x,'reps',10000);
+#    
+# 4. To silence the textual output to the screen, do this
+#    
+#       p = plvar(x,'silent');
+# 
+#
+
+def plvar(x,*varargin):
+    vec     = []
+    sample  = []
+    xminx   = []
+    limit   = []
+    Bt      = []
+    quiet   = False
+    
+
+    # parse command-line parameters trap for bad input
+    i=0 
+    while i<len(varargin): 
+        argok = 1 
+        if type(varargin[i])==str: 
+            if varargin[i]=='range':
+                Range = varargin[i+1]
+                if Range[1]>Range[0]:
+                    argok=0
+                    vec=[]
+                try:
+                    vec=map(lambda X:X*float(Range[2])+Range[0],\
+                            range(int((Range[1]-Range[0])/Range[2])))
+                    
+                    
+                except:
+                    argok=0
+                    vec=[]
+                    
+
+                if Range[0]>=Range[1]:
+                    argok=0
+                    vec=[]
+                    i-=1
+
+                i+=1
+                    
+                            
+            elif varargin[i]== 'sample':
+                sample  = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'limit':
+                limit   = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'xmin':
+                xminx   = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'reps':
+                Bt   = varargin[i+1]
+                i = i + 1
+            elif varargin[i]==  'silent':       quiet  = True    
+ 
+            else: argok=0 
+        
+      
+        if not argok:
+            print '(PLVAR) Ignoring invalid argument #',i+1 
+      
+        i = i+1 
+
+    if vec!=[] and (type(vec)!=list or min(vec)<=1):
+        print '(PLVAR) Error: ''range'' argument must contain a vector or minimum <= 1. using default.\n'                        
+              
+        vec = []
+    
+    if sample!=[] and sample<2:
+        print'(PLVAR) Error: ''sample'' argument must be a positive integer > 1. using default.\n'
+        sample = []
+    
+    if limit!=[] and limit<min(x):
+        print'(PLVAR) Error: ''limit'' argument must be a positive value >= 1. using default.\n'
+        limit = []
+    
+    if xminx!=[] and xminx>=max(x):
+        print'(PLVAR) Error: ''xmin'' argument must be a positive value < max(x). using default behavior.\n'
+        xminx = []
+    if Bt!=[] and Bt<2:
+        print '(PLVAR) Error: ''reps'' argument must be a positive value > 1; using default.\n'
+        Bt = [];
+
+
+    # select method (discrete or continuous) for fitting
+    if     reduce(lambda X,Y:X==True and floor(Y)==float(Y),x,True): f_dattype = 'INTS'
+    elif reduce(lambda X,Y:X==True and (type(Y)==int or type(Y)==float or type(Y)==long),x,True):    f_dattype = 'REAL'
+    else:                 f_dattype = 'UNKN'
+    
+    if f_dattype=='INTS' and min(x) > 1000 and len(x)>100:
+        f_dattype = 'REAL'
+    
+    N=len(x)
+    
+    if Bt==[]: Bt=1000
+    bofA = []
+    bofB = []
+    bofC = []
+    if not quiet:
+        print 'Power-law Distribution, parameter error calculation\n'
+        print '   Copyright 2007-2009 Aaron Clauset\n'
+        print '   Warning: This can be a slow calculation; please be patient.\n'
+        print '   n    = ',len(x),'\n   reps = ',Bt
+    # estimate xmin and alpha, accordingly
+    if f_dattype== 'REAL':
+
+        for B in range(0,Bt):
+            #  bootstrap resample
+            y = []
+            for i in range(0,N):
+                y.append(x[int(floor(N*random()))])
+            ymins = unique(y)
+            ymins.sort()
+            ymins=ymins[0:-1]
+            
+            if xminx!=[]:
+                
+                ymins = [min(filter(lambda X: X>=xminx,ymins))]
+                
+            
+            if limit!=[]:
+                qmins=filter(lambda X: X<=limit,qmins)
+                if qmins==[]: qmins=[min(y)]
+                
+            if sample!=[]:
+                step = float(len(ymins))/(sample-1)
+                index_curr=0
+                new_ymins=[]
+                for i in range (0,sample):
+                    if round(index_curr)==len(ymins): index_curr-=1
+                    new_ymins.append(ymins[int(round(index_curr))])
+                    index_curr+=step
+                ymins = unique(new_ymins)
+                ymins.sort()
+
+            z = sorted(y)
+
+            dat   = []
+            
+            for xm in range(0,len(ymins)):
+                xmin = ymins[xm]
+                z   = filter(lambda X:X>=xmin,z)
+                n   = len(z)
+                # estimate alpha using direct MLE
+                a    = float(n)/sum(map(lambda X: log(float(X)/xmin),z))
+                # compute KS statistic
+                cf   = map(lambda X:1-pow((float(xmin)/X),a),z)
+                
+                dat.append( max( map(lambda X: abs(cf[X]-float(X)/n),range(0,n))))
+            ymin = ymins[dat.index(min(dat))]
+            z = filter(lambda X: X>=ymin,y)
+            n = len(z)
+            alpha    = 1+float(n)/sum(map(lambda X: log(float(X)/ymin),z))
+            bofA.append(n)
+            bofB.append(ymin)
+            bofC.append(alpha)
+            # store distribution of estimated parameter values 
+            if not quiet:
+                print '['+str(B+1)+']\tntail = ',round(mean(bofA),3),' (',round(std(bofA),3),')','\txmin = ',\
+                      round(mean(bofB),3),' (',round(std(bofB),3),')','\talpha = ',round(mean(bofC),3),' (',round(std(bofC),3),')'
+
+        n = std(bofA)
+        xmin = std(bofB)
+        alpha = std(bofC)
+
+    elif f_dattype== 'INTS':
+        x=map(int,x)
+        if vec==[]:
+            for X in range(150,351):
+                vec.append(X/100.)    # covers range of most practical 
+                                    # scaling parameters
+        zvec = map(zeta, vec)
+
+        for B in range(0,Bt):
+            #  bootstrap resample
+            y = []
+            for i in range(0,N):
+                y.append(x[int(floor(N*random()))])
+            ymins = unique(y)
+            ymins.sort()
+            ymins=ymins[0:-1]
+            
+            if xminx!=[]:
+                
+                ymins = [min(filter(lambda X: X>=xminx,ymins))]
+                
+            
+            if limit!=[]:
+                qmins=filter(lambda X: X<=limit,qmins)
+                if qmins==[]: qmins=[min(y)]
+                
+            if sample!=[]:
+                step = float(len(ymins))/(sample-1)
+                index_curr=0
+                new_ymins=[]
+                for i in range (0,sample):
+                    if round(index_curr)==len(ymins): index_curr-=1
+                    new_ymins.append(ymins[int(round(index_curr))])
+                    index_curr+=step
+                ymins = unique(new_ymins)
+                ymins.sort()
+                
+            ymax = max(y)
+            z = sorted(y)
+
+            datA = []
+            datB = []
+            for xm in range(0,len(ymins)):
+                xmin = ymins[xm]
+                z = filter(lambda X:X>=xmin,z)
+                n = len(z)
+                L = []
+                slogz = sum(map(log,z))
+                xminvec = range (1,xmin)
+                for k in range (1,len(vec)):
+                    L.append(-vec[k]*float(slogz) - float(n)*log(float(zvec[k]) - sum(map(lambda X:pow(float(X),-vec[k]),xminvec))))
+        
+                I = L.index(max(L))
+
+                # compute KS statistic
+                fit = reduce(lambda X,Y: X+[Y+X[-1]],\
+                             (map(lambda X: pow(X,-vec[I])/(float(zvec[I])-sum(map(lambda X: pow(X,-vec[I]),map(float,range(1,xmin))))),range(xmin,ymax+1))),[0])[1:]
+                cdi=[]
+                for XM in range(xmin,ymax+1):
+                    cdi.append(len(filter(lambda X: floor(X)<=XM,z))/float(n))
+                
+                datA.append(max( map(lambda X: abs(fit[X] - cdi[X]),range(0,ymax-xmin+1))))
+                datB.append(vec[I])
+                
+            I = datA.index(min(datA))
+            ymin = ymins[I]
+            n = len(filter(lambda X:X>=ymin,y))
+            alpha = datB[I]
+
+            bofA.append(n)
+            bofB.append(ymin)
+            bofC.append(alpha)
+            # store distribution of estimated parameter values 
+            if not quiet:
+                print '['+str(B+1)+']\tntail = ',round(mean(bofA),3),' (',round(std(bofA),3),')','\txmin = ',\
+                      round(mean(bofB),3),' (',round(std(bofB),3),')','\talpha = ',round(mean(bofC),3),' (',round(std(bofC),3),')'
+        n = std(bofA)
+        xmin = std(bofB)
+        alpha = std(bofC)
+        
+    else:
+        print '(PLVAR) Error: x must contain only reals or only integers.\n'
+        n = []
+        xmin = []
+        alpha = []
+        
+
+    return [alpha,xmin,n]
+
+
+# helper functions (unique and zeta)
+def mean(L):
+    try:
+        return float(sum(L))/len(L)
+    except:
+        
+        return 0
+
+def std(L):
+    try:
+        u = mean(L)
+        return sqrt((1./(len(L)-1))*sum(map(lambda X: pow(X-u,2),L)))
+    except:
+        return 0
+        
+
+def unique(seq): 
+    # not order preserving 
+    set = {} 
+    map(set.__setitem__, seq, []) 
+    return set.keys()
+
+def _polyval(coeffs, x):
+    p = coeffs[0]
+    for c in coeffs[1:]:
+        p = c + x*p
+    return p
+
+_zeta_int = [\
+-0.5,
+0.0,
+1.6449340668482264365,1.2020569031595942854,1.0823232337111381915,
+1.0369277551433699263,1.0173430619844491397,1.0083492773819228268,
+1.0040773561979443394,1.0020083928260822144,1.0009945751278180853,
+1.0004941886041194646,1.0002460865533080483,1.0001227133475784891,
+1.0000612481350587048,1.0000305882363070205,1.0000152822594086519,
+1.0000076371976378998,1.0000038172932649998,1.0000019082127165539,
+1.0000009539620338728,1.0000004769329867878,1.0000002384505027277,
+1.0000001192199259653,1.0000000596081890513,1.0000000298035035147,
+1.0000000149015548284]
+
+_zeta_P = [-3.50000000087575873, -0.701274355654678147,
+-0.0672313458590012612, -0.00398731457954257841,
+-0.000160948723019303141, -4.67633010038383371e-6,
+-1.02078104417700585e-7, -1.68030037095896287e-9,
+-1.85231868742346722e-11][::-1]
+
+_zeta_Q = [1.00000000000000000, -0.936552848762465319,
+-0.0588835413263763741, -0.00441498861482948666,
+-0.000143416758067432622, -5.10691659585090782e-6,
+-9.58813053268913799e-8, -1.72963791443181972e-9,
+-1.83527919681474132e-11][::-1]
+
+_zeta_1 = [3.03768838606128127e-10, -1.21924525236601262e-8,
+2.01201845887608893e-7, -1.53917240683468381e-6,
+-5.09890411005967954e-7, 0.000122464707271619326,
+-0.000905721539353130232, -0.00239315326074843037,
+0.084239750013159168, 0.418938517907442414, 0.500000001921884009]
+
+_zeta_0 = [-3.46092485016748794e-10, -6.42610089468292485e-9,
+1.76409071536679773e-7, -1.47141263991560698e-6, -6.38880222546167613e-7,
+0.000122641099800668209, -0.000905894913516772796, -0.00239303348507992713,
+0.0842396947501199816, 0.418938533204660256, 0.500000000000000052]
+
+def zeta(s):
+    """
+    Riemann zeta function, real argument
+    """
+    if not isinstance(s, (float, int)):
+        try:
+            s = float(s)
+        except (ValueError, TypeError):
+            try:
+                s = complex(s)
+                if not s.imag:
+                    return complex(zeta(s.real))
+            except (ValueError, TypeError):
+                pass
+            raise NotImplementedError
+    if s == 1:
+        raise ValueError("zeta(1) pole")
+    if s >= 27:
+        return 1.0 + 2.0**(-s) + 3.0**(-s)
+    n = int(s)
+    if n == s:
+        if n >= 0:
+            return _zeta_int[n]
+        if not (n % 2):
+            return 0.0
+    if s <= 0.0:
+        return 0
+    if s <= 2.0:
+        if s <= 1.0:
+            return _polyval(_zeta_0,s)/(s-1)
+        return _polyval(_zeta_1,s)/(s-1)
+    z = _polyval(_zeta_P,s) / _polyval(_zeta_Q,s)
+    return 1.0 + 2.0**(-s) + 3.0**(-s) + 4.0**(-s)*z
+
+
+
+    
Index: trunk/UTIL/PYTHON/powerlaw/previous_tests/randht.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/previous_tests/randht.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/previous_tests/randht.py	(revision 1053)
@@ -0,0 +1,169 @@
+from math import *
+from random import *
+
+# function x=randht(n, varargin)
+
+# RANDHT generates n observations distributed as some continous heavy-
+# tailed distribution. Options are power law, log-normal, stretched 
+# exponential, power law with cutoff, and exponential. Can specify lower 
+# cutoff, if desired.
+# 
+#    Example:
+#       x = randht(10000,'powerlaw',alpha);
+#       x = randht(10000,'xmin',xmin,'powerlaw',alpha);
+#       x = randht(10000,'cutoff',alpha, lambda);
+#       x = randht(10000,'exponential',lambda);
+#       x = randht(10000,'lognormal',mu,sigma);
+#       x = randht(10000,'stretched',lambda,beta);
+#
+#    See also PLFIT, PLVAR, PLPVA
+#
+#    Source: http://www.santafe.edu/~aaronc/powerlaws/
+
+
+# Version 1.0.2 (2008 April)
+# Copyright (C) 2007 Aaron Clauset (Santa Fe Institute)
+
+# Ported to python by Joel Ornstein (2011 August)
+# (joel_ornstein@hmc.edu)
+
+# Distributed under GPL 2.0
+# http://www.gnu.org/copyleft/gpl.html
+# RANDHT comes with ABSOLUTELY NO WARRANTY
+# 
+# Notes:
+# 
+def randht(n,*varargin):
+    Type   = '';
+    xmin   = 1;
+    alpha  = 2.5;
+    beta   = 1;
+    Lambda = 1;
+    mu     = 1;
+    sigma  = 1;
+
+
+    # parse command-line parameters; trap for bad input
+    i=0; 
+    while i<len(varargin): 
+        argok = 1; 
+        if type(varargin[i])==str: 
+            if varargin[i] == 'xmin':
+                xmin = varargin[i+1]
+                i = i + 1
+            elif varargin[i] == 'powerlaw':
+                Type = 'PL'
+                alpha  = varargin[i+1]
+                i = i + 1
+            elif varargin[i] == 'cutoff':
+                Type = 'PC';
+                alpha  = varargin[i+1]
+                Lambda = varargin[i+2]
+                i = i + 2
+            elif varargin[i] == 'exponential':
+                Type = 'EX'
+                Lambda = varargin[i+1]
+                i = i + 1
+            elif varargin[i] == 'lognormal':
+                Type = 'LN';
+                mu = varargin[i+1]
+                sigma = varargin[i+2]
+                i = i + 2
+            elif varargin[i] == 'stretched':
+                Type = 'ST'
+                Lambda = varargin[i+1]
+                beta = varargin[i+2]
+                i = i + 2
+            else: argok=0
+        
+      
+        if not argok: 
+            print '(RANDHT) Ignoring invalid argument #' ,i+1 
+      
+        i = i+1 
+
+    if n<1:
+        print '(RANDHT) Error: invalid ''n'' argument; using default.\n'
+        n = 10000;
+
+    if xmin < 1:
+        print '(RANDHT) Error: invalid ''xmin'' argument; using default.\n'
+        xmin = 1;
+
+
+
+
+    x=[]
+    if Type == 'EX':
+        x=[]
+        for i in range(n):
+            x.append(xmin - (1./Lambda)*log(1-random()))
+    elif Type == 'LN':
+        y=[]
+        for i in range(10*n):
+            y.append(exp(mu+sigma*normalvariate(0,1)))
+
+        while True:
+            y= filter(lambda X:X>=xmin,y)
+            q = len(y)-n;
+            if q==0: break
+
+            if q>0:
+                r = range(len(y));
+                shuffle(r)
+                ytemp = []
+                for j in range(len(y)):
+                    if j not in r[0:q]:
+                        ytemp.append(y[j])
+                y=ytemp
+                break
+            if (q<0):
+                for j in range(10*n):
+                    y.append(exp(mu+sigma*normalvariate(0,1)))
+            
+        x = y
+        
+    elif Type =='ST':
+        x=[]
+        for i in range(n):
+            x.append(pow(pow(xmin,beta) - (1./Lambda)*log(1.-random()),(1./beta)))
+    elif Type == 'PC':
+        
+        x = []
+        y=[]
+        for i in range(10*n):
+            y.append(xmin - (1./Lambda)*log(1.-random()))
+        while True:
+            ytemp=[]
+            for i in range(10*n):
+                if random()<pow(y[i]/float(xmin),-alpha):ytemp.append(y[i])
+            y = ytemp
+            x = x+y
+            q = len(x)-n
+            if q==0: break;
+
+            if (q>0):
+                r = range(len(x))
+                shuffle(r)
+
+                xtemp = []
+                for j in range(len(x)):
+                    if j not in r[0:q]:
+                        xtemp.append(x[j])
+                x=xtemp
+                break;
+            
+            if (q<0):
+                y=[]
+                for j in range(10*n):
+                    y.append(xmin - (1./Lambda)*log(1.-random()))
+
+
+    else:
+        x=[]
+        for i in range(n):
+            x.append(xmin*pow(1.-random(),-1./(alpha-1.)))
+
+    return x
+                              
+
Index: trunk/UTIL/PYTHON/powerlaw/randht.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/randht.py	(revision 946)
+++ 	(revision )
@@ -1,169 +1,0 @@
-from math import *
-from random import *
-
-# function x=randht(n, varargin)
-
-# RANDHT generates n observations distributed as some continous heavy-
-# tailed distribution. Options are power law, log-normal, stretched 
-# exponential, power law with cutoff, and exponential. Can specify lower 
-# cutoff, if desired.
-# 
-#    Example:
-#       x = randht(10000,'powerlaw',alpha);
-#       x = randht(10000,'xmin',xmin,'powerlaw',alpha);
-#       x = randht(10000,'cutoff',alpha, lambda);
-#       x = randht(10000,'exponential',lambda);
-#       x = randht(10000,'lognormal',mu,sigma);
-#       x = randht(10000,'stretched',lambda,beta);
-#
-#    See also PLFIT, PLVAR, PLPVA
-#
-#    Source: http://www.santafe.edu/~aaronc/powerlaws/
-
-
-# Version 1.0.2 (2008 April)
-# Copyright (C) 2007 Aaron Clauset (Santa Fe Institute)
-
-# Ported to python by Joel Ornstein (2011 August)
-# (joel_ornstein@hmc.edu)
-
-# Distributed under GPL 2.0
-# http://www.gnu.org/copyleft/gpl.html
-# RANDHT comes with ABSOLUTELY NO WARRANTY
-# 
-# Notes:
-# 
-def randht(n,*varargin):
-    Type   = '';
-    xmin   = 1;
-    alpha  = 2.5;
-    beta   = 1;
-    Lambda = 1;
-    mu     = 1;
-    sigma  = 1;
-
-
-    # parse command-line parameters; trap for bad input
-    i=0; 
-    while i<len(varargin): 
-        argok = 1; 
-        if type(varargin[i])==str: 
-            if varargin[i] == 'xmin':
-                xmin = varargin[i+1]
-                i = i + 1
-            elif varargin[i] == 'powerlaw':
-                Type = 'PL'
-                alpha  = varargin[i+1]
-                i = i + 1
-            elif varargin[i] == 'cutoff':
-                Type = 'PC';
-                alpha  = varargin[i+1]
-                Lambda = varargin[i+2]
-                i = i + 2
-            elif varargin[i] == 'exponential':
-                Type = 'EX'
-                Lambda = varargin[i+1]
-                i = i + 1
-            elif varargin[i] == 'lognormal':
-                Type = 'LN';
-                mu = varargin[i+1]
-                sigma = varargin[i+2]
-                i = i + 2
-            elif varargin[i] == 'stretched':
-                Type = 'ST'
-                Lambda = varargin[i+1]
-                beta = varargin[i+2]
-                i = i + 2
-            else: argok=0
-        
-      
-        if not argok: 
-            print '(RANDHT) Ignoring invalid argument #' ,i+1 
-      
-        i = i+1 
-
-    if n<1:
-        print '(RANDHT) Error: invalid ''n'' argument; using default.\n'
-        n = 10000;
-
-    if xmin < 1:
-        print '(RANDHT) Error: invalid ''xmin'' argument; using default.\n'
-        xmin = 1;
-
-
-
-
-    x=[]
-    if Type == 'EX':
-        x=[]
-        for i in range(n):
-            x.append(xmin - (1./Lambda)*log(1-random()))
-    elif Type == 'LN':
-        y=[]
-        for i in range(10*n):
-            y.append(exp(mu+sigma*normalvariate(0,1)))
-
-        while True:
-            y= filter(lambda X:X>=xmin,y)
-            q = len(y)-n;
-            if q==0: break
-
-            if q>0:
-                r = range(len(y));
-                shuffle(r)
-                ytemp = []
-                for j in range(len(y)):
-                    if j not in r[0:q]:
-                        ytemp.append(y[j])
-                y=ytemp
-                break
-            if (q<0):
-                for j in range(10*n):
-                    y.append(exp(mu+sigma*normalvariate(0,1)))
-            
-        x = y
-        
-    elif Type =='ST':
-        x=[]
-        for i in range(n):
-            x.append(pow(pow(xmin,beta) - (1./Lambda)*log(1.-random()),(1./beta)))
-    elif Type == 'PC':
-        
-        x = []
-        y=[]
-        for i in range(10*n):
-            y.append(xmin - (1./Lambda)*log(1.-random()))
-        while True:
-            ytemp=[]
-            for i in range(10*n):
-                if random()<pow(y[i]/float(xmin),-alpha):ytemp.append(y[i])
-            y = ytemp
-            x = x+y
-            q = len(x)-n
-            if q==0: break;
-
-            if (q>0):
-                r = range(len(x))
-                shuffle(r)
-
-                xtemp = []
-                for j in range(len(x)):
-                    if j not in r[0:q]:
-                        xtemp.append(x[j])
-                x=xtemp
-                break;
-            
-            if (q<0):
-                y=[]
-                for j in range(10*n):
-                    y.append(xmin - (1./Lambda)*log(1.-random()))
-
-
-    else:
-        x=[]
-        for i in range(n):
-            x.append(xmin*pow(1.-random(),-1./(alpha-1.)))
-
-    return x
-                              
-
Index: trunk/UTIL/PYTHON/powerlaw/showdd.py
===================================================================
--- trunk/UTIL/PYTHON/powerlaw/showdd.py	(revision 1053)
+++ trunk/UTIL/PYTHON/powerlaw/showdd.py	(revision 1053)
@@ -0,0 +1,92 @@
+#! /usr/bin/env python
+
+from ppclass import pp
+import numpy as np
+from scipy import ndimage
+
+dx = 50. ; filefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.nc"
+lim = 4.
+#lim = -3.0
+#lim = -3.2
+
+#dx = 10. ; filefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.nc"
+
+## getplot pressure field
+## ... generic settings
+psfc = pp()
+psfc.file = filefile
+psfc.var = "PSFC"
+psfc.xcoeff = dx/1000.
+psfc.ycoeff = dx/1000.
+psfc.xlabel = "x distance (km)"
+psfc.ylabel = "y distance (km)"
+psfc.title = "Surface pressure $P_s$"
+#psfc.vmin = 482.5
+#psfc.vmax = 483.5
+psfc.div = 20
+#psfc.out = "png" ; psfc.includedate = False ; psfc.folder = "dd/"
+
+## calculate std
+std = np.std(psfc.getf())
+
+## calculate mean pressure
+## ... same results whether or not time average
+psfcmean = pp()
+psfcmean << psfc
+psfcmean.x="0,10000"
+psfcmean.y="0,10000"
+psfcmean.t="0,10000"
+psfcmean.getplot()
+
+#for t in [300]:
+#for t in range(0,500,50):
+#for t in [390,395,400,405,410,415]:
+#for t in [178,180,182]:
+
+for t in range(265,275,1):
+
+   ## plot pressure field
+   psfc.t = t
+   psfc.filename = "dd"+str(int(t))
+   psfc.getplot(extraplot=1)
+
+   ## calculate anomaly normalized by std
+   anomaly = (psfc - psfcmean) / std
+      
+   ## plot anomaly
+   anomaly.vmax = lim
+   anomaly.vmin = -4.0
+   anomaly.div = 32
+   anomaly.colorb = "gist_earth_r"
+   anomaly.colorb = "gist_stern"
+   anomaly.plotin = psfc
+   anomaly.title = "Anomaly $P_s - <P_s>^{xy}$"
+   anomaly.units = "$\sigma$"
+   #anomaly.plot()
+
+   ## sobel transform
+   field = anomaly.request[0][0][0][0][0][0].field
+   sx = ndimage.sobel(field, axis=0, mode='constant') 
+   sy = ndimage.sobel(field, axis=1, mode='constant')
+   tmp = np.hypot(sx, sy)
+
+   ## more prominent circles if we mutiply by field
+   tmp = tmp*field
+   #lim = -12.
+   #tmp[tmp>lim]=0.
+   #tmp[tmp<lim]=1.
+
+   #tmp = np.exp(-tmp)   
+
+   anomaly.request[0][0][0][0][0][0].field = tmp
+   anomaly.title = "Sobel transform" 
+   anomaly.units = "dimless" ; anomaly.vmax = None ; anomaly.vmin = None ; anomaly.div = 10
+   anomaly.colorb = "binary_r"
+   #anomaly.vmax = np.max(tmp)
+   anomaly.plot()
+
+   ### laplace transform
+   #anomaly.request[0][0][0][0][0][0].field=ndimage.laplace(anomaly.request[0][0][0][0][0][0].field)
+   #anomaly.title = "Sobel + Laplace transform"
+   #anomaly.units = "dimless" ; anomaly.vmax = None ; anomaly.vmin = None ; anomaly.div = 10
+   #anomaly.plot()
