def logskew(d,floors = [],col='',divvypk=[]):
    maxd = max(d.flatten())
    if len(floors) == 0:
        floors = 1./maxd * 2.**N.arange(-5,5.01,1.)
    print floors

    for fl in floors:
        dcopy = 1.*d
        wltf = N.where(dcopy < fl)
        dcopy[wltf] = 0.*dcopy[wltf] + fl

        logrho = N.log(dcopy)
        var  = N.var(decreaseres(logrho,f=16).flatten())
        #skew = SS.skew(logrho.flatten())
        #print fl,'log:',var,'lin:',N.var(dcopy.flatten())
        print fl,'log:',var#,SS.skew(dcopy.flatten())

        k,pk = power.pk(logrho)

        if (len(divvypk) > 0):
            M.loglog(k,pk/divvypk/var,col)

    if (len(divvypk) == 0):
        return pk,pk/var
    else:
        return

def extrapolate(n, y, tolerance=1e-15, plot=False, call_show=True):
    "Extrapolate functional value Y from sequence of values (n, y)."

    # Make sure we have NumPy arrays
    n = array(n)
    y = array(y)

    # Create initial "bound"
    Y0 = 0.99*y[-1]
    Y1 = 1.01*y[-1]

    # Compute initial interior points
    phi = (sqrt(5.0) + 1.0) / 2.0
    Y2 = Y1 - (Y1 - Y0) / phi
    Y3 = Y0 + (Y1 - Y0) / phi

    # Compute initial values
    F0, e, nn, ee, yy = _eval(n, y, Y0)
    F1, e, nn, ee, yy = _eval(n, y, Y1)
    F2, e, nn, ee, yy = _eval(n, y, Y2)
    F3, e, nn, ee, yy = _eval(n, y, Y3)

    # Solve using direct search (golden ratio fraction)
    while Y1 - Y0 > tolerance:

        if F2 < F3:
            Y1, F1 = Y3, F3
            Y3, F3 = Y2, F2
            Y2 = Y1 - (Y1 - Y0) / phi
            F2, e, nn, ee, yy = _eval(n, y, Y2)
        else:
            Y0, F0 = Y2, F2
            Y2, F2 = Y3, F3
            Y3 = Y0 + (Y1 - Y0) / phi
            F3, e, nn, ee, yy = _eval(n, y, Y3)

        print Y0, Y1

    # Compute reference value
    Y = 0.5*(Y0 + Y1)

    # Print results
    print
    print "Reference value:", Y

    # Plot result
    if plot:
        pylab.figure()
        pylab.subplot(2, 1, 1)
        pylab.title("Reference value: %g" % Y)
        pylab.semilogx(n, y, 'b-o')
        pylab.semilogx(nn, yy, 'g--')
        pylab.subplot(2, 1, 2)
        pylab.loglog(n, e, 'b-o')
        pylab.loglog(nn, ee, 'g--')
        pylab.grid(True)
        if call_show:
            pylab.show()

    return Y

    def romanzuniga07(wavelength, AKs, makePlot=False):
        filters = ['J', 'H', 'Ks', '[3.6]', '[4.5]', '[5.8]', '[8.0]']
        wave =      np.array([1.240, 1.664, 2.164, 3.545, 4.442, 5.675, 7.760])
        A_AKs =     np.array([2.299, 1.550, 1.000, 0.618, 0.525, 0.462, 0.455])
        A_AKs_err = np.array([0.530, 0.080, 0.000, 0.077, 0.063, 0.055, 0.059])
        
        # Interpolate over the curve
        spline_interp = interpolate.splrep(wave, A_AKs, k=3, s=0)

        A_AKs_at_wave = interpolate.splev(wavelength, spline_interp)
        A_at_wave = AKs * A_AKs_at_wave

        if makePlot:
            py.clf()
            py.errorbar(wave, A_AKs, yerr=A_AKs_err, fmt='bo', 
                        markerfacecolor='none', markeredgecolor='blue',
                        markeredgewidth=2)

            # Make an interpolated curve.
            wavePlot = np.arange(wave.min(), wave.max(), 0.1)
            extPlot = interpolate.splev(wavePlot, spline_interp)
            py.loglog(wavePlot, extPlot, 'k-')

            # Plot a marker for the computed value.
            py.plot(wavelength, A_AKs_at_wave, 'rs',
                    markerfacecolor='none', markeredgecolor='red',
                    markeredgewidth=2)
            py.xlabel('Wavelength (microns)')
            py.ylabel('Extinction (magnitudes)')
            py.title('Roman Zuniga et al. 2007')


        return A_at_wave

    def plotppf(self,x=None,xmin=None,alpha=None,dolog=True,**kwargs):
        """
        Plots the power-law-predicted value on the Y-axis against the real
        values along the X-axis.  Can be used as a diagnostic of the fit 
        quality.
        """
        if not(xmin): xmin=self._xmin
        if not(alpha): alpha=self._alpha
        if not(x): x=numpy.sort(self.data[self.data>xmin])
        else: x=numpy.sort(x[x>xmin])

        # N = M^(-alpha+1)
        # M = N^(1/(-alpha+1))
        
        m0 = min(x)
        N = (1.0+numpy.arange(len(x)))[::-1]
        xmodel = m0 * N**(1/(1-alpha)) / max(N)**(1/(1-alpha))
        
        if dolog:
            pylab.loglog(x,xmodel,'.',**kwargs)
            pylab.gca().set_xlim(min(x),max(x))
            pylab.gca().set_ylim(min(x),max(x))
        else:
            pylab.plot(x,xmodel,'.',**kwargs)
        pylab.plot([min(x),max(x)],[min(x),max(x)],'k--')
        pylab.xlabel("Real Value")
        pylab.ylabel("Power-Law Model Value")

def plotData(filename):
    gatetime, allanvar = np.loadtxt(filename, comments="#", delimiter=",", unpack=True)
    pl.loglog(gatetime, allanvar, '.-')
    pl.xlabel("Gate time (s)")
    pl.ylabel("Allan deviation (Hz)")
    pl.xlim([min(gatetime), max(gatetime)])
    pl.show()

def plot_calibrated(time_len=128, chan=1):
    """ Plot calibrated as a function of time """
    nar.load_data() if plot_old else nar.take_data(time_len)
        
    cal   = (nar.ts_x_on + nar.ts_x_off) * Td / (nar.ts_x_on/nar.ts_x_off - 1)
    p_tot = (nar.ts_x_on + nar.ts_x_off)
    time_series  = cal[:, chan]
    p_avg        = np.average(p_tot[:, chan])
    time_series  = time_series / p_avg
    spec_series = np.abs(np.fft.rfft(time_series))
    
    #uncal = nar.ts_x_on
    #uncal = uncal[:, chan]
    #spec_uncal = np.abs(np.fft.rfft(uncal))
    
    t = np.cumsum(np.ones([time_len]))
    t = t /np.max(t) * total_time
    tu = np.cumsum(np.ones([time_len/2+1]))/total_time
    
    plt.subplot(211)
    #plt.plot(t, uncal, c=c[1])
    plt.plot(t, time_series, c=c[0])
    plt.xlabel("Time (s)")
    plt.ylabel("Calibrated signal (K)")
    plt.subplot(212)
    #plt.loglog(tu, spec_uncal, c=c[1])
    plt.loglog(tu, spec_series, c=c[0])

    plt.xlabel("Spectrum (Hz)")
    plt.ylabel("")
    plt.show()

def snr_mat_f(mchvec, reds, lum_dist, fmin, fmax, fvec, finteg, tobs, sn_f):
	''''''
	mch_fmat=np.transpose(np.tile(mchvec, (len(reds), len(fvec), finteg, 1) ), axes=(0,3,1,2))
	z_fmat=np.transpose(np.tile(reds, (len(mchvec), len(fvec), finteg, 1) ),axes=(3,0,1,2))
	f_fmat=np.transpose(np.tile(fvec, (len(reds), len(mchvec), finteg, 1) ), axes=(0,1,3,2))
	finteg_fmat=np.transpose(np.tile(np.arange(finteg), (len(reds), len(mchvec), len(fvec), 1) ), axes=(0,1,2,3))
	stshape=np.shape(z_fmat) #Standard shape of all matrices that I will use.
	DL_fmat=np.transpose(np.tile(lum_dist, (len(mchvec), len(fvec), finteg, 1) ),axes=(3,0,1,2)) #Luminosity distance in Mpc.
	flim_fmat=A8.f_cut(1./4., 2.*mch_fmat*2.**(1./5.))*1./(1.+z_fmat) #The symmetric mass ratio is 1/4, since I assume equal masses.
	flim_det=np.maximum(np.minimum(fmax, flim_fmat), fmin) #The isco frequency limited to the detector window.
	tlim_fmat=CM.tafter(mch_fmat, f_fmat, flim_fmat, z_fmat)
	#By construction, f_mat cannot be smaller than fmin or larger than fmax (which are the limits imposed by the detector).
	fmin_fmat=np.minimum(f_fmat, flim_det) #I impose that the minimum frequency cannot be larger than the fisco.
	fmaxobs_fmat=flim_det.copy()
	#fmaxobs_fmat=fmin_fmat.copy()
	fmaxobs_fmat[tobs<tlim_fmat]=CM.fafter(mch_fmat[tobs<tlim_fmat], z_fmat[tobs<tlim_fmat], f_fmat[tobs<tlim_fmat], tobs)
	fmax_fmat=np.minimum(fmaxobs_fmat, flim_det) #The maximum frequency (after an observation tobs) cannot exceed fisco or the maximum frequency of the detector.
	integconst=(np.log10(fmax_fmat)-np.log10(fmin_fmat))*1./(finteg-1)
	finteg_fmat=fmin_fmat*10**(integconst*finteg_fmat)
	sn_vec=sn_f(fvec)##########
	sn_fmat=sn_f(finteg_fmat) #Noise spectral density.
	#htilde_fmat=A8.htilde_f(1./4., 2.*mch_fmat*2**(1./5.), z_fmat, DL_fmat, f_fmat)
	htilde_fmat=A8.htilde_f(1./4., 2.*mch_fmat*2**(1./5.), z_fmat, DL_fmat, finteg_fmat)
	py.loglog(finteg_fmat[0,0,:,0],htilde_fmat[0,0,:,0]**2.)
	py.loglog(finteg_fmat[0,0,:,0],sn_fmat[0,0,:,0])
	snrsq_int_fmat=4.*htilde_fmat**2./sn_fmat #Integrand of the S/N square.
	snrsq_int_m_fmat=0.5*(snrsq_int_fmat[:,:,:,1:]+snrsq_int_fmat[:,:,:,:-1]) #Integrand at the arithmetic mean of the infinitesimal intervals.
	df_fmat=np.diff(finteg_fmat, axis=3) #Infinitesimal intervals.
	snr_full_fmat=np.sqrt(np.sum(snrsq_int_m_fmat*df_fmat,axis=3)) #S/N as a function of redshift, mass and frequency.
	fopt=fvec[np.argmax(snr_full_fmat, axis=2)] #Frequency at which the S/N is maximum, for each pixel of redshift and mass.
	snr_opt=np.amax(snr_full_fmat, axis=2) #Maximum S/N at each pixel of redshift and mass.
	snr_min=snr_full_fmat[:,:,0]
	return snr_opt

def nishiyama09(wavelength, AKs, makePlot=False):
    # Data pulled from Nishiyama et al. 2009, Table 1

    filters = ['V', 'J', 'H', 'Ks', '[3.6]', '[4.5]', '[5.8]', '[8.0]']
    wave =      np.array([0.551, 1.25, 1.63, 2.14, 3.545, 4.442, 5.675, 7.760])
    A_AKs =     np.array([16.13, 3.02, 1.73, 1.00, 0.500, 0.390, 0.360, 0.430])
    A_AKs_err = np.array([0.04,  0.04, 0.03, 0.00, 0.010, 0.010, 0.010, 0.010])

    # Interpolate over the curve
    spline_interp = interpolate.splrep(wave, A_AKs, k=3, s=0)

    A_AKs_at_wave = interpolate.splev(wavelength, spline_interp)
    A_at_wave = AKs * A_AKs_at_wave

    if makePlot:
        py.clf()
        py.errorbar(wave, A_AKs, yerr=A_AKs_err, fmt='bo', 
                    markerfacecolor='none', markeredgecolor='blue',
                    markeredgewidth=2)
        
        # Make an interpolated curve.
        wavePlot = np.arange(wave.min(), wave.max(), 0.1)
        extPlot = interpolate.splev(wavePlot, spline_interp)
        py.loglog(wavePlot, extPlot, 'k-')

        # Plot a marker for the computed value.
        py.plot(wavelength, A_AKs_at_wave, 'rs',
                markerfacecolor='none', markeredgecolor='red',
                markeredgewidth=2)
        py.xlabel('Wavelength (microns)')
        py.ylabel('Extinction (magnitudes)')
        py.title('Nishiyama et al. 2009')

    
    return A_at_wave

def ConvIndicator(X, Y, pct=0.1, fs=14, eqaxis=False):
    """
    Convergence indicator icon
    ==========================
    """
    if len(X)<2: raise Exception('at least 2 points are required')
    xx, yy   = log10(X), log10(Y)
    p        = polyfit(xx, yy, 1)
    m        = round(p[0])
    xx0, xx1 = min(xx), max(xx)
    yy0, yy1 = min(yy), max(yy)
    dxx, dyy = xx1-xx0, yy1-yy0
    xxm, yym = (xx0+xx1)/2.0, (yy0+yy1)/2.0
    xxl, xxr = xxm-pct*dxx, xxm+pct*dxx
    shift    = 0.5*pct*dxx*m
    xm,  ym  = 10.0**xxm, 10.0**(yym-shift)
    xl,  xr  = 10.0**xxl, 10.0**xxr
    yl,  yr  = 10.0**(yym+m*(xxl-xxm)-shift),10.0**(yym+m*(xxr-xxm)-shift)
    loglog(X, Y)
    #plot(xm, ym, 'ro')
    #plot(xl, yl, 'go')
    #plot(xr, yr, 'mo')
    points = array([[xl,yl],[xr,yl],[xr,yr]])
    gca().add_patch(Polygon(points, ec='k', fc='None'))
    xxR = xxm+1.2*pct*dxx
    xR  = 10.0**xxR
    text(xR, ym, '%g'%m, ha='left', va='center', fontsize=fs)
    if eqaxis: axis('equal')
    return m

def plot_point(rank1, style=None):
    if not style:
        style = "r--"
    # plot lines to point (rank1,pop1)
    pop1 = pop[rank1]
    P.loglog([rank1, rank1], [0.1, pop1], style)
    P.loglog([1, rank1], [pop1, pop1], style)

def main():
    plt.ion()

    fil = FletcherFilter()
    Niter = 12
    logp = plt.zeros((Niter,2))
    for k in range(Niter):
        while True:
            #print k
            p = plt.rand(2)
            if not fil.dominated(p):
                break
        logp[k] = p
        fil.add(p, 0.0, 0.0)
        ff = fil.values[fil.valid]
        ff = plt.r_[[[1e-6,1]], ff[plt.argsort(ff[:,0])], [[1,1e-6]]]
        ww = plt.zeros((ff.shape[0] * 2 - 1, 2))
        ww[::2] = ff
        ww[1::2,0] = ff[1:,0]
        ww[1::2,1] = ff[:-1,1]
        plt.loglog(ww[:,0], ww[:,1], '-')
    plt.loglog(logp[:,0], logp[:,1], 'ys-', lw=2)
    plt.axis([0,1,0,1])
    plt.axis('equal')
    plt.grid()
        
    code.interact()

def test_integrators(tmax=100,
                     Nsteps=1000,
                     x0=[0.0,1.0],
                     v0=[1.0,0.0]):
    mag = lambda x: numpy.sqrt(numpy.dot(x,x))
    fprime = lambda x: -x/mag(x)**3
    energy = lambda x,v: 0.5*(v**2).sum(1) - 1./numpy.sqrt((x**2).sum(1))

    pylab.figure(1)
    
    for (method,N) in ((Euler,Nsteps),
                       (DriftKick,Nsteps),
                       (LeapFrog,Nsteps),
                       (VelocityVerlet,Nsteps),
                       (RungeKutta,Nsteps/4)):
        t = numpy.linspace(0,tmax,N+1)
        x_t,v_t = method(fprime,x0,v0,tmax,N)
        E = energy(x_t,v_t)

        delta_E = abs((E-E[0])/E[0])
        
        pylab.loglog(t[1:],delta_E[1:],
                     label=method.__name__)
    
    pylab.ylim(1E-8,1E4)
    pylab.legend(loc=2)
    pylab.xlabel(r'$\mathdefault{t}$')
    pylab.ylabel(r'$\mathdefault{|\Delta E/E|}$')
    
    pylab.show()

def plot_convergence_data_from_file( labels ):
    import pickle
    l2_error = pickle.load( open( 'l2-error.p', 'rb' ) )
    inf_error = pickle.load( open( 'inf-error.p', 'rb' ) )
    num_pts = pickle.load( open( 'num-grid-points.p', 'rb' ) )


    max_error = 0.
    min_error = 1.
    for i in xrange( len( l2_error ) ):
        pylab.figure( 1 )
        pylab.loglog( num_pts[i], l2_error[i], '-o', label = labels[i] )


        pylab.figure( 2 )
        pylab.loglog( num_pts[i], inf_error[i], '-o', label = labels[i] )

        max_error = max( max_error, numpy.array( l2_error[i] ).max() )
        min_error = min( min_error, numpy.array( l2_error[i] ).min() )

    #matplotlib.rcParams.update({'font.size': 16})
    pylab.figure( 1 )
    pylab.legend()
    pylab.xlabel(r'Number of grid points')
    pylab.ylabel(r'$\lVert f - \hat{f}\rVert_{\ell_2}$')#,fontsize=16)
    pylab.ylim(min_error/5.,5.*max_error)
    pylab.savefig('genz-corner-peak-10d-5e-1c-quartic-decay-l2-convergence.eps',dpi=1200)

    pylab.figure( 2 )
    #pylab.title(r'$\int_{-1}^1 f(x)$', fontsize=10)
    pylab.legend()
    pylab.xlabel(r'Number of grid points')
    pylab.ylabel(r'$\lVert I[f] - Q[f]\rVert_{\ell_2}$')
    pylab.ylim(1e-16,5*max_error)

def expy(d):
    colors = ['b','g','r','c','m','y','k','b--','g--','r--','c--','m--','y--']
    multi = 2.**N.arange(2,3.1,1.)
    #multi = [77.]
    for i in range(len(multi)):
        dinflate = N.exp((d)*multi[i])-1.
        k,p = power.pk(dinflate)

        var = N.var((d)*multi[i])
        lognocor = ((N.exp(var)-1)*N.exp(var)/var)
        lognovar = (N.exp(var)-1)*N.exp(var)
        print 1+2*var,'lognocor:',lognocor
        #varlog = N.var(N.log(dinflate+1.).flatten())
        k,plog = power.pk(N.log(dinflate+1))
        print p[0]/plog[0], p[-1]/plog[-1]

        #M.subplot(121)
        M.loglog(k,p/plog/lognocor,colors[i])
        #M.subplot(122)

        #M.loglog(k,plog,colors[i])
        #M.loglog([k[0],k[-1]],[lognovar,lognovar],colors[i])
        #M.loglog(k,1./plog**mul,colors[i])
        

    
    #preal = M.load('mill/s63/pm.pnl.dat')
    #plogreal = M.load('mill/s63/plogm.pnl.dat')    
    #M.loglog(preal[:,0],preal[:,1]/plogreal[:,1],'b')

    #M.xlabel(r'$k\ [\rm{Mpc}/h]$',fontsize=20)
    #M.ylabel(r'$P_\delta(k)/P_{\log (1+\delta)}(k)$',fontsize=20)

    M.show()

def MisorientationScalingCollapseCompareInset(misses, bdlengths, labels, alpha=2.5):
    """ good values for alpha seem to be 4, but 2.5 for experiment """
    colors = ['b', 'r', 'g', 'y']

    pl.rcParams.update({'legend.fontsize': 14,
            'legend.columnspacing':1.2,
        })
    for i, (mis, label, bdlength) in enumerate(zip(misses, labels, bdlengths)):
        t = mis*bdlength/(mis*bdlength).mean()
        dx = 5./100.
        y,x = np.histogram(t, bins=np.linspace(0, 5, 100))
        x = (x[:-1]+x[1:])/2
        y = y.astype('float')/y.sum() / dx
        pl.plot(x, y, colors[i]+'o--', label=label)

        alpha = fit_alpha(x, y)
        xt = np.linspace(0,5, 1000)
        pl.plot(xt, scaling(alpha, xt), colors[i]+'-', label=r"Fit, $\alpha$ = %0.1f" % alpha)

    pl.xlabel(r"$\theta / \theta_{av}$")
    pl.ylabel(r"$\theta_{av}\,P(\theta, \theta_{av})$")
    pl.legend(loc='lower right')

    ax = pl.axes([0.52, 0.52, 0.35, 0.35])
    for i, (mis, label, bdlength) in enumerate(zip(misses, labels, bdlengths)):
        t = mis*bdlength/(mis*bdlength).mean()
        dx = 5./100.
        y,x = np.histogram(t, bins=np.linspace(0, 5, 100))
        x = (x[:-1]+x[1:])/2
        y = y.astype('float')/y.sum() / dx
        pl.plot(x, y, colors[i]+'o-', label=label)
    pl.loglog()

    return x, y

def plotDependencyComponents():

    """Plot thoretical dependency between n_samples and n_components"""

    # range of admissible distortions
    eps_range = np.linspace(0.1, 0.99, 5)
    colors = pl.cm.Blues(np.linspace(0.3, 1.0, len(eps_range)))

    # range of number of samples to embed
    n_samples_range = np.logspace(1, 9, 9)

    
    pl.figure()

    for eps, color in zip(eps_range, colors):
        min_n_components = johnson_lindenstrauss_min_dim(n_samples_range, \
                                                         eps=eps)
        pl.loglog(n_samples_range, min_n_components, color=color)

    pl.legend(["eps = %.1f" % eps for eps in eps_range], \
              loc="lower right")

    pl.xlabel("Number of observations to eps-embed")
    pl.ylabel("Minimum number of dimensions")
    pl.title("Johnson-Lindenstrauss bounds:\nn_samples vs n_components")
    pl.show()

def plot_gain_drift(time_len=128, chan=1):
    """ Have a look at the gain drift """    
    nar.load_data() if plot_old else nar.take_data(time_len)

    time_series_on = nar.ts_x_on[:, chan]
    spec_series_on = np.abs(np.fft.fft(time_series_on))

    time_series_off = nar.ts_x_off[:, chan]
    spec_series_off = np.abs(np.fft.fft(time_series_off))
    
    #t  = np.cumsum(np.ones([time_len]))
    #tt = t/ np.max(t) * total_time
    #tu = np.cumsum(np.ones([time_len]))/total_time
    
    #print spec_series_off.shape
    #print tu.shape
    
    plt.subplot(211)
    plt.plot(time_series_on, label="X-on", c=c[1])
    #plt.plot(tt, 10*np.log10(time_series_off)), label="X-off", c=c[2])
    plt.xlabel("Time")
    plt.legend()
    plt.subplot(212)
    
    plt.loglog(spec_series_on, c=c[1])
    #plt.loglog(tu, spec_series_off, c=c[2])
    plt.xlabel("Spectrum")
    plt.show()

def Run(inputfile,input,shape,fieldtype,corrfunctype,symmetrytype,logbinning,binnum,outputfile,plot):
    """
    Performs RadialCorrelationFunctions on an inputfile that holds field data
    in the space separated values of arbitrary arrangement (i.e. rows / columns don't matter).
    
    Input is either by filename or by numpy.array.  If input is specified it is taken to be the 
    data array.  If not, the inputfile is opened and read for data.

    data structure: array.shape should be like (dimension, dimension[0], dimension[1], ...).
        examples:
            * one dimensional 3-component vector field with 256 values - (3, 256)
            * three dimensional 3-component vector field - (3, 256, 256, 256)
            * 3x3 3d tensor field - (3, 3, 256, 256, 256)
    """
    if input is None:
        import Data_IO
        data = Data_IO.ReadInScalarField(inputfile,shape) 
    else:
        data = input
    corr_func = RadialCorrelationFunctions(data,fieldtype,corrfunctype,symmetrytype,logbinning,binnum,fromState=False)
    if outputfile is None:
        outputfile = inputfile+'_correlationfunction.dat'
    Data_IO.OutputXY(corr_func[0],corr_func[1],outputfile)
    if plot:
        import pylab
        pylab.figure()
        pylab.loglog(corr_func[0],corr_func[1],'.-')
        pylab.xlabel(r'$C(r)$',fontsize=20)
        pylab.ylabel(r'$r$',fontsize=20)
        pylab.show()

def plot_cluster_context(sizes, densities, dir, name=None, k=None, suffix="png"):
    """
    so many conditionals!
    """
    print("plot_cluster_context(): plotting", name)
    if name is None:
        K = len(sizes)
        fn = "{}/clusters_{:04d}.{}".format(dir, K, suffix)
    else:
        fn = "{}/{}_context.{}".format(dir, name, suffix)
    if os.path.exists(fn):
        print("plot_cluster_context(): {} exists already".format(fn))
        return
    if k is None:
        fig = plt.figure(figsize=(6,6))
        plt.subplots_adjust(left=0.15, right=0.97, bottom=0.15, top=0.97)
        ms = 7.5
    else:
        fig = plt.figure(figsize=(4,4))
        plt.subplots_adjust(left=0.2, right=0.96, bottom=0.2, top=0.96)
        ms = 5.0
    plt.clf()
    if name is not None and k is None:
        plt.savefig(fn)
        print("plot_cluster_context(): wrote", fn)
        return
    _clusterplot(sizes, densities, k, ms=ms)
    _clusterlims(sizes, densities)
    plt.ylabel("cluster abundance-space density")
    plt.xlabel("number in abundance-space cluster")
    plt.loglog()
    [l.set_rotation(45) for l in plt.gca().get_xticklabels()]
    [l.set_rotation(45) for l in plt.gca().get_yticklabels()]
    plt.savefig(fn)
    print("plot_cluster_context(): wrote", fn)

def plot_session_figs(in_file):
    """Plot some graphs on the generated session file:
    - thp vs duration
    - cdf of duration
    """
    sessions = np.load(in_file)
    # duration vs min thp
    pylab.clf()
    pylab.plot(sessions[:, 3], sessions[:, 1], 'k, ')
    pylab.loglog()
    axes = pylab.gca()
    pylab.grid()
    axes = pylab.gca()
    for tick in axes.xaxis.get_major_ticks():
        tick.label1.set_fontsize(16)
    pylab.xlabel("Minimum throughput in kbps", size=16)
    for tick in axes.yaxis.get_major_ticks():
        tick.label1.set_fontsize(16)
    pylab.ylabel("Session duration per client (sec)", size=16)
    pylab.savefig('%s_session_duration_vs_min_thp.pdf' % in_file)
    # nb flows vs avg thp
    pylab.clf()
    pylab.plot(sessions[:, 4], sessions[:, 2], 'k, ')
    pylab.loglog()
    pylab.grid()
#    axes = pylab.gca()
    pylab.xlabel("Average throughput in kbps", size=16)
    pylab.ylabel("Nb of flows", size=16)
    pylab.savefig('%s_session_nb_fl_vs_avg_thp.pdf' % in_file)
    pylab.clf()
    import cdfplot
    cdfplot.cdfplotdata(sessions[:, 1], _xlabel='Duration in seconds',
            _title='Session durations', _fs_legend='x-large')