def __call__(self, **params):

       p = ParamOverrides(self, params)
       fig = plt.figure(figsize=(5, 5))

       # This one-liner works in Octave, but in matplotlib it
       # results in lines that are all connected across rows and columns,
       # so here we plot each line separately:
       #   plt.plot(x,y,"k-",transpose(x),transpose(y),"k-")
       # Here, the "k-" means plot in black using solid lines;
       # see matplotlib for more info.
       isint = plt.isinteractive() # Temporarily make non-interactive for
       # plotting
       plt.ioff()
       for r, c in zip(p.y[::p.skip], p.x[::p.skip]):
           plt.plot(c, r, "k-")
       for r, c in zip(np.transpose(p.y)[::p.skip],np.transpose(p.x)[::p.skip]):
           plt.plot(c, r, "k-")

       # Force last line avoid leaving cells open
       if p.skip != 1:
           plt.plot(p.x[-1], p.y[-1], "k-")
           plt.plot(np.transpose(p.x)[-1], np.transpose(p.y)[-1], "k-")

       plt.xlabel('x')
       plt.ylabel('y')
       # Currently sets the input range arbitrarily; should presumably figure out
       # what the actual possible range is for this simulation (which would presumably
       # be the maximum size of any GeneratorSheet?).
       plt.axis(p.axis)

       if isint: plt.ion()
       self._generate_figure(p)
       return fig

def plot_movie(x,t,u,NN,**kwargs):
	"""
	     x  | Spatial coordinate vector,  len (Nx)
	     t  | Temporal coordinate vector, len (Nt)
	  data  | The data in matrix form,    len (Nt,Nx)
	    NN  | integer giving interval of plotting.
	
	"""
	N,M = shape(u)
	clf()
	ion()
	line, = plot(x,u[0,:],'k-',label=parse_kwargs('label','$Wave$ $equation:$ $BW$',**kwargs),
		linewidth=2)
	plot(x,u[0,:],color='gray',alpha=0.75)
	line.axes.set_ylim(parse_kwargs('miny',-1,**kwargs),parse_kwargs('maxy',1,**kwargs))
	legend(loc=0)
	xlabel(parse_kwargs('xlabel','$x$',**kwargs))
	ylabel(parse_kwargs('ylabel','$u$',**kwargs))
	grid(True)
	for i in range(0,N,N/NN):
		title('$t={}$'.format(t[i]))
		line.set_ydata(u[i,:])
		plot(x,u[i,:],color='gray',alpha=0.2)
		xlim([min(x),max(x)])
		draw()
	line.set_ydata(u[-1,:])
	
	title('$t={}$'.format(t[-1]))

def kmr_test_plot(data, k, end_thresh):
    from matplotlib.pylab import ion, figure, draw, ioff, show, plot, cla
    ion()
    fig = figure()
    ax = fig.add_subplot(111)
    ax.grid(True)

    # get k centroids
    kmr = kmeans.kmeans_runner(k, end_thresh)
    kmr.init_data(data)
    print kmr.centroids

    plot(data[:,0], data[:,1], 'o')

    i = 0
    while kmr.stop_flag is False:
        kmr.iterate()
        #print kmr.centroids, kmr.itr_count
        plot(kmr.centroids[:, 0], kmr.centroids[:, 1], 'sr')
        time.sleep(.2)
        draw()
        i += 1

    print "N Iterations: %d" % (i)
    plot(kmr.centroids[:, 0], kmr.centroids[:, 1], 'g^', linewidth=3)

    ioff()
    show()
    print kmr.itr_count, kmr.centroids

    def matrix_plot(self, matrix, figure_name='matrix_plot.pdf'):
        import numpy
        from matplotlib import pylab
        def _blob(x,y,area,colour):
            hs = numpy.sqrt(area) / 2
            xcorners = numpy.array([x - hs, x + hs, x + hs, x - hs])
            ycorners = numpy.array([y - hs, y - hs, y + hs, y + hs])
            pylab.fill(xcorners, ycorners, colour, edgecolor=colour)
        reenable = False
        if pylab.isinteractive():
            pylab.ioff()
        pylab.clf()
        
        maxWeight = 2**numpy.ceil(numpy.log(numpy.max(numpy.abs(matrix)))/numpy.log(2))
        height, width = matrix.shape
        pylab.fill(numpy.array([0,width,width,0]),numpy.array([0,0,height,height]),'white')
        pylab.axis('off')
        pylab.axis('equal')
        for x in xrange(width):
            for y in xrange(height):
                _x = x+1
                _y = y+1
                w = matrix[y,x]
                if w > 0:
                    _blob(_x - 0.5, height - _y + 0.5, 0.2,'#0099CC')
                elif w < 0:
                    _blob(_x - 0.5, height - _y + 0.5, 0.2,'#660000')

        if reenable:
            pylab.ion()
        pylab.savefig(figure_name) 

    def do_fit_trans(self):
        f = self.fit_trans
        p0 = self.p0.copy()
        #        p0[:-3]=0.0
        #print p0,"call d0"
        #d0 = f(p0)
        if 0:
         for i in range(len(p0)):
            pt = p0.copy()
            pt[i] = p0[i]+0.001
            from matplotlib.pylab import clf, ion, title, plot, show
            print pt - p0, pt
            ion()
            clf()
            title("%d"%(i))
            plot(d0, f(pt) - d0, ",")
            show()
            if raw_input()[0] != " ": 
                break


        res = scipy.optimize.leastsq( f, p0, full_output=1)
        pfit, pcov, info, errmsg, ier = res
        if ier not in [1,2,3,4]:
            print s_sq, ier, errmsg
        else:
            residu = f(pfit) 
            s_sq = (residu**2).sum()/(len(residu)-len(p0))
        ubi = pfit[:9].reshape(3,3)
        print ("%.6f "*6)%(indexing.ubitocellpars(ubi))
        print pfit[9:12]
        self.g = grain( ubi, pfit[9:12].copy())

def example():
#    pl.ioff()
    pl.ion()

    import pandas
    from numpy.random import uniform

    n = 25
    m = pandas.DataFrame({
            'x': uniform(-1, 1, size=n),
            'y': uniform(-1, 1, size=n),
            'size': uniform(3, 10, size=n) ** 2,
            'color': uniform(0, 1, size=n),
    })

    # test using a custom index
    m['silly_index'] = ['%sth' % x for x in range(n)]
    m.set_index('silly_index', drop=True, inplace=True, verify_integrity=True)

    print m

    ax = pl.subplot(111)
    plt = ax.scatter(m['x'], m['y'])

    b = LassoBrowser(m, ax)
    print b.idxs

    #from viz.interact.pointbrowser import PointBrowser
    #pb = PointBrowser(m, plot=plt)

    pl.show()

    ip()

def count_barcodes(dataset, VERBOSE=0):
    '''Count the abundance of each barcode'''

    # Get the read filenames
    data_filenames = get_raw_read_files(dataset)
    datafile = data_filenames['adapter']

    # Count the abundance of each barcode
    bc_counts = defaultdict(int)
    rc = 0
    with open(datafile, 'r') as infile:
        for read in SeqIO.parse(infile, 'fastq'):
            bc_counts[read.seq.tostring()] += 1
            rc += 1
            if rc == maxreads:
                break
    
    print sorted(bc_counts.items(), key=lambda x:x[1], reverse=True)[:20]
    
    # Plot results
    plt.figure()
    ax=plt.subplot(111)
    plt.plot(range(1,len(bc_counts)+1), sorted(bc_counts.values(), reverse=True))
    ax.set_yscale('log')
    ax.set_xscale('log')
    plt.xlabel('barcode rank')
    plt.ylabel('abundance')

    plt.ion()
    plt.show()

def main():
    # u_t = u_xx
    dx = .1
    dt = .5
    timesteps = 100000

    x = np.arange(-10,10,dx)
    m = len(x)
    kappa = 50

    # u''(x) = (u(x + dx) - 2u(x) + u(x - dx)) / dx^2
    ones = lambda x: np.ones(x)
    A = np.diag(ones(m-1),k=-1) + -2*np.diag(ones(m)) + np.diag(ones(m-1),k=1)
    A *= kappa*(dx**2)

    U = 0*ones(m)
    for i in xrange(0,m):
        if x[i] > -2 and x[i] < 2:
            U[i] = 1

    p.ion()
    lines, = p.plot(x,U)
    for n in xrange(0,timesteps):
        U = U + dt*dudt(U,A)
        
        if n % 100 == 0:
            lines.set_ydata(U)
            p.draw()
    p.show()

    def plot_diagnostics(self):
        if 0:
            from matplotlib import pylab as plt
            plt.ion()
            plt.figure()

        _phot = phot.read_fits(self.lcfn,'optimum')
        with self.FigureManager('_0-aperture'):
            plotting.phot.aperture(_phot)

        with self.FigureManager('_1-background'):
            plotting.phot.background(_phot)

        if len(self.dfaper.npix.drop_duplicates()) > 1:
            with self.FigureManager('_2-noise_vs_aperture_size'):
                plotting.pipeline.noise_vs_aperture_size(self)

        with self.FigureManager("_3-fdt_t_roll_2D"):
            plotting.phot.detrend_t_roll_2D(_phot)

        with self.FigureManager("_4-fdt_t_roll_2D_zoom"):
            plotting.phot.detrend_t_roll_2D(_phot,zoom=True)

        with self.FigureManager("_5-fdt_t_rollmed"):
            plotting.phot.detrend_t_rollmed(_phot)

def rf_sub_size(Input, net_size, resolution):
    size = np.zeros((net_size*net_size,))
    coms = np.zeros((net_size*net_size, 2))
    R = np.zeros((net_size*resolution, net_size*resolution))
    Z = np.zeros((resolution, resolution))

    scale = 1.0/(resolution**2)

    X, Y = np.meshgrid(np.arange(Z.shape[0]), np.arange(Z.shape[1]))

    count_roi, count_nroi, count_tot = 0, 0, 0
    plt.ion()
    for i in range(net_size):
        for j in range(net_size):
            Z = np.abs(Input[i, j, ...] * (Input[i, j, ...] > 0) +
                       0.0 * (Input[i, j, ...] < 0))

            R[i*resolution:(i+1)*resolution, j*resolution:(j+1)*resolution] = Z
            size[i*net_size+j] = area_of_activation(Z) * scale

            d = np.unravel_index(Z.argmax(), Z.shape)
            Z = np.roll(Z, Z.shape[0]/2-d[0], axis=0)
            Z = np.roll(Z, Z.shape[1]/2-d[1], axis=1)

            xc = ((Z*Y).sum()/Z.sum() - Z.shape[0]/2 + d[0])/float(Z.shape[0])
            yc = ((Z*X).sum()/Z.sum() - Z.shape[1]/2 + d[1])/float(Z.shape[1])

            coms[i*net_size+j, 0] = (xc+1.0) % 1
            coms[i*net_size+j, 1] = (yc+1.0) % 1
    return coms, R, size

def plot(y, function):
    """ Show an animation of Poincare plot.

    --- arguments ---
    y: A list of initial values
    function: function which is argument of Runge-Kutta solver
    """
    h = dt
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.grid()
    time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
    plt.ion()

    for i in range(nmax + 1):
        for j in range(nstep):
            rk4 = RK.RK4(function)
            y = rk4.solve(y, j * h, h)
            # -pi <= theta <= pi
            while y[0] > pi:
                y[0] = y[0] - 2 * pi
            while y[0] < -pi:
                y[0] = y[0] + 2 * pi

        if ntransient <= i < nmax:          # <-- draw the poincare plots
            plt.scatter(y[0], y[1], s=2.0, marker='o', color='blue')
            time_text.set_text('n = %d' % i)
            plt.draw()

        if i == nmax:                       # <-- to stop the interactive mode
            plt.ioff()
            plt.scatter(y[0], y[1], s=2.0, marker='o', color='blue')
            time_text.set_text('n = %d' % i)
            plt.show()

    def eqDistribution(self, plot=True):
        """ Obtain and plot the equilibrium probabilities of each macrostate

        Parameters
        ----------
        plot : bool, optional, default=True
            Disable plotting of the probabilities by setting it to False

        Returns
        -------
        eq : ndarray
            An array of equilibrium probabilities of the macrostates

        Examples
        --------
        >>> model = Model(data)
        >>> model.markovModel(100, 5)
        >>> model.eqDistribution()
        """
        self._integrityCheck(postmsm=True)
        macroeq = np.ones(self.macronum) * -1
        for i in range(self.macronum):
            macroeq[i] = np.sum(self.msm.stationary_distribution[self.macro_ofmicro == i])

        if plot:
            from matplotlib import pylab as plt
            plt.ion()
            plt.figure()
            plt.bar(range(self.macronum), macroeq)
            plt.ylabel('Equilibrium probability')
            plt.xlabel('Macrostates')
            plt.xticks(np.arange(0.4, self.macronum+0.4, 1), range(self.macronum))
            plt.show()
        return macroeq

def demo():
    '''
    Load and plot a few CIB spectra.
    '''

    # define ell array.
    l = np.arange(100,4000)

    # get dictionary of CIBxCIB spectra.
    cl_cibcib = get_cl_cibcib(l)

    # plot
    import matplotlib.pylab as pl
    pl.ion()
    lw=2
    fs=18
    leg = []
    pl.clf()
    for band in ['857','545','353']:
        pl.semilogy(l, cl_cibcib['545',band],linewidth=lw)
        leg.append('545 x '+band)
    pl.xlabel(r'$\ell$',fontsize=fs)
    pl.ylabel(r'$C_\ell^{TT, CIB} [\mu K^2]$',fontsize=fs)
    pl.ylim(5e-2,6e3)
    pl.legend(leg, fontsize=fs)

def label_data(prefix, size=100, savename=None):
    from glob import glob
    from os.path import basename
    from PIL import Image
    from os.path import isfile
    if savename==None: savename=labelpath+'label_'+prefix+'.txt'
    # We want to avoid labeling an image twice, so keep track
    # of what we've labeled in previous labeling sessions.
    if isfile(savename):
        fileout = open(savename,'r')
        already_seen = [line.split(',')[0] for line in fileout]
        fileout.close()
    else: already_seen = []
    # Now reopen the file for appending.
    fileout = open(savename,'a')
    pl.ion()
    pl.figure(1,figsize=(9,9))
    files = glob(imgpath+prefix+'*.png')
    for file in np.random.choice(files, size=size, replace=False):
        if basename(file) in already_seen: continue
        pl.clf()
        pl.subplot(1,1,1)
        pl.imshow(np.array(Image.open(file)))
        pl.title(file)
        pl.axis('off')
        pl.draw()
        label = get_one_char()
        if label=='q': break
        fileout.write(basename(file)+','+label+'\n')
        print file,label
    fileout.close()
    return

def plot_average(filenames, save_plot=True, show_plot=False, dpi=100):

    ''' Plot Signal average from a list of averaged files. '''

    fname = get_files_from_list(filenames)

    # plot averages
    pl.ioff()  # switch off (interactive) plot visualisation
    factor = 1e15
    for fnavg in fname:
        name = fnavg[0:len(fnavg) - 4]
        basename = os.path.splitext(os.path.basename(name))[0]
        print fnavg
        # mne.read_evokeds provides a list or a single evoked based on condition.
        # here we assume only one evoked is returned (requires further handling)
        avg = mne.read_evokeds(fnavg)[0]
        ymin, ymax = avg.data.min(), avg.data.max()
        ymin *= factor * 1.1
        ymax *= factor * 1.1
        fig = pl.figure(basename, figsize=(10, 8), dpi=100)
        pl.clf()
        pl.ylim([ymin, ymax])
        pl.xlim([avg.times.min(), avg.times.max()])
        pl.plot(avg.times, avg.data.T * factor, color='black')
        pl.title(basename)

        # save figure
        fnfig = os.path.splitext(fnavg)[0] + '.png'
        pl.savefig(fnfig, dpi=dpi)

    pl.ion()  # switch on (interactive) plot visualisation

def integral_plots(vals, new_vals, n):
    
    """
    Parameters
    -------------------------------------------------------
        vals:   non-normalized values of the table
        
        new_vals: normalized values of the table
                  such that the maximum of each vertical
                  column is 1.
        
        n: number of points in theta and omega arrays
        (so the size of the 2D table is n x n)
        
        
    Return
    -------------------------------------------------------
        Plots a density plot of the table using imshow
        of the normalized table and also plots the integral 
        of the table along each axis (these are the integrals
        of the non-normalized values)
    """

    import matplotlib.gridspec as gridspec
    
    omega_par = np.linspace(0.1,0.4,n)
    theta_par = np.linspace(0.01,1.4,n)
    
    oo = []
    tt = []
    
    # integral in left side (sum of each row)
    for i in range(n):
        oo.append(sum(vals[i]))
    
    # integral in bottom side (sum of columns)
    for i in range(n):
        tt.append(sum(vals.T[i]))

    plt.ion()
    plt.figure(figsize=(11,9))

    gs  = gridspec.GridSpec(2, 2, width_ratios=[1,4],height_ratios=[4,1])
    ax1 = plt.subplot(gs[0])
    ax1.plot(oo,omega_par,  linewidth=2)
    plt.ylabel(r"$\Omega_{||}$", fontsize=20)

    ax2 = plt.subplot(gs[1])
    plt.imshow(new_vals, interpolation='nearest', origin="lower", extent=(theta_par.min(), theta_par.max(), omega_par.min(), omega_par.max()),aspect='auto')
    #plt.imshow(vals, interpolation='nearest', origin="lower", extent=(theta_par.min(), theta_par.max(), omega_par.min(), omega_par.max()),aspect='auto')
    plt.colorbar()

    plt.title("$P(t_s) \propto t_s $, maximum normalized to 1", fontsize=15)


    ax4 = plt.subplot(gs[3])
    ax4.plot(theta_par,tt, linewidth=2)
    plt.xlabel(r"$\theta_{||}$",fontsize=20)

def plot_pol_info_pat(problem_id, pol, times):
    fig, axes = plt.subplots(nrows=len(times), ncols=1)
    for idx, time in enumerate(times):
        df = pd.DataFrame(pol[:, :, time])
        df.plot(title='{0:s}: Whether to sell, depending on current price; time = {1:d}'.format(problem_id, time),
                     ax=axes[idx])
        axes[idx].set_xlabel('Number of stocks in inventory')
        axes[idx].set_ylabel('1->sell; 0->hold')
    plt.ion()

def assignpeaks( gr, pars, colfile, tol=None, omegatol=0.2 ):
    labels = np.zeros( colfile.nrows )
    cyf = cyfit.cyfit(pars)
    cyf.setscfc( colfile.sc, colfile.fc )
    hkl = np.zeros( cyf.XL.shape, np.float)
    drlv = np.zeros( colfile.nrows, np.float)
    kcalc = np.zeros( cyf.XL.shape, np.float)
    r_omega  = np.array(colfile.omega*np.pi/180.0,np.float)
    cyf.hkl( colfile.sc, colfile.fc, r_omega, gr, hkl, kcalc )
    # Make integer hkl
    wripaca.ih_drlv( hkl, drlv )
    # Now replace omegaobs by omegacalc where this is appropriate
    pre = np.eye(3).ravel()
    posti = np.dot(gv_general.wedgemat(pars.get('wedge')), 
                     gv_general.chimat(pars.get('chi'))).T.ravel()
    axis = np.array([0,0,-1],np.float)
    ub = np.linalg.inv(gr.ubi)
    gcalc = np.dot( ub, hkl.T ).T.copy()
    #print hkl
    romegacalc = np.zeros( r_omega.shape, np.float)
    romegaerr  = np.zeros( r_omega.shape, np.float)
    wripaca.omegacalcclose(gcalc,
                           pre,
                           posti,
                           axis,
                           r_omega,
                           romegacalc,
                           romegaerr,
                           pars.get('wavelength'),
                           )
    # OK, now we will accept anything within omegatol
    if 0:
        pylab.hist(romegaerr, bins=50)
        pylab.figure()
        pylab.plot(r_omega, romegaerr,",")
        pylab.show()
    r_omega = np.where( romegaerr < omegatol*np.pi/180.,  romegacalc, r_omega )
    cyf.hkl( colfile.sc, colfile.fc, r_omega, gr, hkl, kcalc )
    drlv_omegafree = drlv.copy()
    wripaca.ih_drlv( hkl, drlv_omegafree )
    m = drlv < tol
    sc = colfile.sc[m]
    fc = colfile.fc[m]
    omegaobs = colfile.omega[m].astype(np.float)
    omega = colfile.omega[m].astype(np.float)
    hkl = hkl[m]
    print "Got",m.sum(),"peaks",
    if 0:
        pylab.ion()
        pylab.title("tol = %f"%(tol))
        pylab.hist( drlv, bins=np.arange(0,0.05,0.001))
        pylab.hist( drlv_omegafree, bins=np.arange(0,0.05,0.001))
        pylab.show()
        raw_input("OK?")

    return dummycf( sc, fc, omega, omegaobs, hkl )

 def begin_draw():
   pl.ion()
   pl.figure( 1, figsize = ( 20, 20 ), dpi = 50 )
   pl.clf()
   pl.axis( "scaled" )
   pl.xlim( -4, 4 )
   pl.ylim( -2, 6 )
   pl.plot( [d[0] for d in data], [d[1] for d in data], "b." )
   for m in xrange( clusters ):
     plot_gaussian( mu[m], sigma[m], "r" )

def normalized_contour_max1(ptype, n):
    
    """
    Parameters
    -------------------------------------------------------
        ptype: type of the probability. It can be either
                "Gauss", "one_over_ts" or "ts"
                
        n: number of points in theta and omega arrays
    
    Return
    -------------------------------------------------------
        val:     non-normalized values of the 2D table.
        
        new_val: normalized values fo the 2D table
                 such that the maximum value of each
                 vertical column is 1.
                 
        It also plots the normalized density plot.
    """
    
    omega_par = np.linspace(0.1,0.4,n)
    theta_par = np.linspace(0.01,1.4,n)
    
    
    vals     = np.zeros((n,n))
    Nlist    = []
    new_vals = []


    # calculating the 2D table of distribution
    for i in range(n):
        for j in range(n):
            vals[i][j] = prob_Oparapar(5.,i, j, n, ptype)


    # calculating the normalization A
    for i in range(n):
        xx = (np.max(vals.T[i]))
        Nlist.append(1./xx)

    # calculating the normalized values
    for i in range(n):
        xx = vals.T[i] * Nlist[i]
        new_vals.append(xx)

    new_vals = np.array(new_vals)

    plt.ion()
    plt.imshow(new_vals.T, interpolation='nearest', origin="lower", extent=(theta_par.min(), theta_par.max(), omega_par.min(), omega_par.max()),aspect='auto')

    plt.xlabel(r"$\theta_{||}$",fontsize=20)
    plt.ylabel(r"$\Omega_{||}$", fontsize=20)
    plt.title("$p(t_s) = {0}$.formar(ptype), max of each column is 1", fontsize=15)
    return vals, new_vals