def test_fit_gaussian_2d(self):
        import matplotlib.pylab as plt
        # Create the gaussian data
        Xin, Yin = plt.mgrid[0:201, 0:201]
        data = self._gaussian(3, 100, 100, 20, 40)(Xin, Yin) + \
            np.random.random(Xin.shape)

        plt.clf()
        ax2 = plt.axes([0, 0.52, 1.0, 0.4])
        ax2.matshow(data, cmap=plt.cm.gist_earth_r)

        params = ap.fitGaussianImage(data)
        fit = self._gaussian(*params)

        ax2.contour(fit(*np.indices(data.shape)), cmap=plt.cm.cool)
        (height, x, y, width_x, width_y) = params

        plt.text(0.95, 0.05, """
        x : %.1f
        y : %.1f
        width_x : %.1f
        width_y : %.1f""" %(x, y, width_x, width_y),
             fontsize=16, horizontalalignment='right',
             verticalalignment='bottom', transform=ax2.transAxes)

        ax1 = plt.axes([0, 0.08, 1.0, 0.4])

def propagate_el_den_TB(Nx=50, Ny=40, mu=0.05, frame_num=100):
    
    ham = envtb.ldos.hamiltonian.HamiltonianTB(Ny, Nx)
    
    wf_final = wf_init_from_electron_density(ham, mu = mu)
    
    maxel =  1.2 * np.max(wf_final.wf1d)
    wf_final.plot_wave_function(maxel)
    plt.axes().set_aspect('equal')
    plt.savefig('../../../../Desktop/pics_2d/TB/0%d_2d.png' % 0)
    plt.close()
                    
    for i in xrange(frame_num):
        
        potential = envtb.ldos.potential.Potential1DFromFunction(
            lambda x: -5. * (Ny/2 - x) * 2. / Ny * np.sin(0.1 * i))
        
        ham2 = ham.apply_potential(potential)
        
        envtb.ldos.plotter.Plotter().plot_potential(
            ham2, ham, maxel = 5, minel = -5)
        
        plt.axes().set_aspect('equal')
        plt.savefig('../../../../Desktop/pics_2d/TB/pot%03d_2d.png' % i)
        plt.close()
        
        
        wf_init = wf_final
        wf_final = propagate_wave_function(wf_init, ham2, maxel = maxel, 
            file_out = '../../../../Desktop/pics_2d/TB/f%03d_2d.png' % i)
    
    return None

def display_grid(grid, **kwargs):
    fig = plt.figure()
    plt.axes().set_aspect('equal')

    if kwargs.get('mark_core_cells', True):
        core_cell_coords = grid._cell_nodes[1:-1, 1:-1]
        cellx, celly = core_cell_coords[:, :, 0], core_cell_coords[:, :, 1]
        plt.plot(cellx, celly, '-o', np.transpose(cellx), np.transpose(celly), '-o', color='red')

    if kwargs.get('mark_boundary_cells', True):
        boundary_cell_coords = grid._cell_nodes[0, :], \
                               grid._cell_nodes[-1, :], \
                               grid._cell_nodes[1:-1, 0], \
                               grid._cell_nodes[1:-1, -1]

        for coords in boundary_cell_coords:
            plt.plot(coords[:, 0], coords[:, 1], '-x', color='blue')

    if kwargs.get('show', False):
        plt.show()

    f = BytesIO()
    plt.savefig(f)

    return f

def propagate_gauss_graphene(Nx=30, Ny=30, frame_num=100):
    
    ham = envtb.ldos.hamiltonian.HamiltonianGraphene(Ny, Nx)
    
    ic = Nx/2 * Ny + Ny/2 + 2
    
    wf_final = wf_init_gaussian_wave_packet(
        ham.coords, ic, p0=[1.0, 0.0], sigma=7.)
    
    maxel = 0.8 * np.max(wf_final.wf1d)
    wf_final.plot_wave_function(maxel)
    plt.axes().set_aspect('equal')
    plt.savefig('../../../../Desktop/pics_2d/graphene/0%d_2d.png' % 0)
    plt.close()
    
    dt_new = 0.5
    NK_new = 12
    
    for i in xrange(frame_num):
        print 'frame %(i)d' % vars()
        wf_init = wf_final
        wf_final, dt_new, NK_new = propagate_wave_function(
            wf_init, ham, NK=NK_new, dt=dt_new, maxel = maxel, regime='SIL', 
            file_out = '../../../../Desktop/pics_2d/graphene/f%03d_2d.png'
            % i)
     
    return None

def plot(frame,dirname,clim=None,axis_limits=None):
    if not os.path.exists('./figures'):
        os.makedirs('./figures')
        
    try:
        sol=Solution(frame,file_format='petsc',read_aux=False,path='./saved_data/'+dirname+'/_p/',file_prefix='claw_p')
    except IOError:
        'Data file not found; please unzip the files in saved_data/.'
        return
    x=sol.state.grid.x.centers; y=sol.state.grid.y.centers
    mx=len(x); my=len(y)
    
    mp=sol.state.num_eqn    
    yy,xx = np.meshgrid(y,x)

    p=sol.state.q[0,:,:]
    if clim is not None:
        pl.pcolormesh(xx,yy,p,cmap=cm.RdBu_r)
    else:
        pl.pcolormesh(xx,yy,p,cmap=cm.Reds)
    pl.title("t= "+str(sol.state.t),fontsize=20)
    pl.xticks(size=20); pl.yticks(size=20)
    cb = pl.colorbar();

    if clim is not None:
        pl.clim(clim[0],clim[1]);
    imaxes = pl.gca(); pl.axes(cb.ax)
    pl.yticks(fontsize=20); pl.axes(imaxes)
    pl.axis('equal')
    if axis_limits is None:
        pl.axis([np.min(x),np.max(x),np.min(y),np.max(y)])
    else:
        pl.axis([axis_limits[0],axis_limits[1],axis_limits[2],axis_limits[3]])
    pl.savefig('./figures/'+dirname+'.png')
    pl.close()

def plot_p_leading_order(frame):
    mat = scipy.io.loadmat('sound-speed_2D-wave.mat')
    T=5; nt=T/0.5
    pp=mat['U'][nt,:,:]
    xx=mat['xx']
    yy=mat['yy']

    fig=pl.figure(figsize=(8, 3.5))
    #pl.title("t= "+str(sol.state.t),fontsize=20)
    pl.xticks(size=20); pl.yticks(size=20)
    pl.xlabel('x',fontsize=20); pl.ylabel('y',fontsize=20)
    #pl.pcolormesh(xx,yy,p_subxy,cmap=cm.OrRd)
    pl.pcolormesh(xx,yy,pp,cmap='RdBu_r')
    pl.autoscale(tight=True)
    cb = pl.colorbar(ticks=[0.5,1,1.5,2]);
    
    #pl.clim(ticks=[0.5,1,1.5,2])
    imaxes = pl.gca(); pl.axes(cb.ax)
    pl.yticks(fontsize=20); pl.axes(imaxes)
    #pl.xticks(fontsize=20); pl.axes(imaxes)
    #pl.axis('equal')
    pl.axis('tight')
    fig.tight_layout()
    pl.savefig('./_plots_to_paper/sound-speed_LO_t'+str(frame)+'_pcolor.png')
    pl.close()

def propagate_gauss_graphene_W90(Nx=30, Ny=20, magnetic_B=None,
                                 frame_num=100):
    
    ham_w90 = define_zigzag_ribbon_w90(
        "../../exampledata/02_graphene_3rdnn/graphene3rdnnlist.dat", 
        Ny, Nx, magnetic_B=magnetic_B)
    
    ham = envtb.ldos.hamiltonian.HamiltonianFromW90(ham_w90, Nx)
    
    ic = ham.Nx/2 * ham.Ny + ham.Ny/2 + 2 
    
    wf_final = wf_init_gaussian_wave_packet(ham.coords, ic, sigma = 10.)
    
    maxel = 0.8 * np.max(wf_final.wf1d)
    wf_final.plot_wave_function(maxel)
    plt.axes().set_aspect('equal')
    plt.savefig('../../../../Desktop/pics_2d/grapheneW90/0%d_2d.png' % 0)
    plt.close()

    for i in xrange(frame_num):
    
        wf_init = wf_final
        wf_final = propagate_wave_function(wf_init, ham, maxel=maxel,
            file_out='../../../../Desktop/pics_2d/grapheneW90/f%03d_2d.png' % i)[0]
    
    return None

def plotLFP ():
    print('Plotting LFP power spectral density...')

    colorspsd=array([[0.42,0.67,0.84],[0.42,0.83,0.59],[0.90,0.76,0.00],[0.90,0.32,0.00],[0.34,0.67,0.67],[0.42,0.82,0.83],[0.90,0.59,0.00],[0.33,0.67,0.47],[1.00,0.85,0.00],[0.71,0.82,0.41],[0.57,0.67,0.33],[1.00,0.38,0.60],[0.5,0.2,0.0],[0.0,0.2,0.5]]) 

    lfpv=[[] for c in range(len(sim.lfppops))]    
    # Get last modified .mat file if no input and plot
    for c in range(len(sim.lfppops)):
        lfpv[c] = sim.lfps[:,c]    
    lfptot = sum(lfpv)
        
    # plot pops separately
    plotPops = 0
    if plotPops:    
        figure() # Open a new figure
        for p in range(len(sim.lfppops)):
            psd(lfpv[p],Fs=200, linewidth= 2,color=colorspsd[p])
            xlabel('Frequency (Hz)')
            ylabel('Power')
            h=axes()
            h.set_yticklabels([])
        legend(['L2/3','L5A', 'L5B', 'L6'])

    # plot overall psd
    figure() # Open a new figure
    psd(lfptot,Fs=200, linewidth= 2)
    xlabel('Frequency (Hz)')
    ylabel('Power')
    h=axes()
    h.set_yticklabels([])

    show()

def plot_decision_surface(axes, clusters, X, Y=None):
    import matplotlib.pylab as pylab
    import numpy as np

    def kmeans_predict(clusters, X):
        from ..core import distance
        dist_m = distance.minkowski_dist(clusters, X)
        print 'dist_m:', dist_m.shape
        pred = np.argmin(dist_m, axis=0)
        print 'pred:', pred.shape
        return pred

    # step size in the mesh
    h = (np.max(X, axis=0) - np.min(X, axis=0)) / 100.0
    # create a mesh to plot in
    x_min = np.min(X, axis=0) - 1
    x_max = np.max(X, axis=0) + 1
    xx, yy = np.meshgrid(np.arange(x_min[0], x_max[0], h[0]),
                         np.arange(x_min[1], x_max[1], h[1]))
    Z = kmeans_predict(clusters, np.c_[xx.ravel(), yy.ravel()])
    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    pylab.set_cmap(pylab.cm.Paired)
    pylab.axes(axes)
    pylab.contourf(xx, yy, Z, cmap=pylab.cm.Paired)
    pylab.xlim(np.min(xx), np.max(xx))
    pylab.ylim(np.min(yy), np.max(yy))
    #pylab.axis('off')

    # Plot also the training points
    if Y is not None:
        pylab.scatter(X[:,0], X[:,1], c=Y)
    else:
        pylab.scatter(X[:,0], X[:,1])
    pylab.scatter(clusters[:,0], clusters[:,1], s=200, marker='x', color='white')

    def _stacking(self, star_list, x_list, y_list):
        """

        :param star_list:
        :return:
        """
        n_stars = len(star_list)

        shifteds = []
        for i in range(n_stars):
            xc, yc = x_list[i], y_list[i]
            data = star_list[i]
            x_shift = int(xc) - xc
            y_shift = int(yc) - yc
            shifted = interp.shift(data, [-y_shift, -x_shift], order=1)
            shifteds.append(shifted)
            print('=== object ===', i)
            import matplotlib.pylab as plt
            fig, ax1 = plt.subplots()
            im = ax1.matshow(np.log10(shifted), origin='lower')
            plt.axes(ax1)
            fig.colorbar(im)
            plt.show()

        combined = sum(shifteds)
        new=np.empty_like(combined)
        max_pix = np.max(combined)
        p = combined[combined>=max_pix/10**6]  #in the SIS regime
        new[combined < max_pix/10**6] = 0
        new[combined >= max_pix/10**6] = p
        kernel = util.kernel_norm(new)
        return kernel

    def __init__(self, ax, collection, mmc, img):
        self.colornormalizer = Normalize(vmin=0, vmax=1, clip=False)
        self.scat = plt.scatter(img[:, 0], img[:, 1], c=mmc.classvec)
        plt.gray()
        plt.setp(ax.get_yticklabels(), visible=False)
        ax.yaxis.set_tick_params(size=0)
        plt.setp(ax.get_xticklabels(), visible=False)
        ax.xaxis.set_tick_params(size=0)
        self.img = img
        self.canvas = ax.figure.canvas
        self.collection = collection
        #self.alpha_other = alpha_other
        self.mmc = mmc
        self.prevnewclazz = None

        self.xys = collection
        self.Npts = len(self.xys)
        
        self.lockedset = set([])

        self.lasso = LassoSelector(ax, onselect=self.onselect)#, lineprops = {:'prism'})
        self.lasso.disconnect_events()
        self.lasso.connect_event('button_press_event', self.lasso.onpress)
        self.lasso.connect_event('button_release_event', self.onrelease)
        self.lasso.connect_event('motion_notify_event', self.lasso.onmove)
        self.lasso.connect_event('draw_event', self.lasso.update_background)
        self.lasso.connect_event('key_press_event', self.onkeypressed)
        #self.lasso.connect_event('button_release_event', self.onrelease)
        self.ind = []
        self.slider_axis = plt.axes(slider_coords, visible = False)
        self.slider_axis2 = plt.axes(obj_fun_display_coords, visible = False)
        self.in_selection_slider = None
        newws = list(set(range(len(self.collection))) - self.lockedset)
        self.mmc.new_working_set(newws)
        self.lasso.line.set_visible(False)

def plot_fig(data):
    nullfmt = NullFormatter()
    x, y = data[:, 0], data[:, 1]
    left, width = 0.1, 0.65
    bottom, height = 0.1, 0.65
    bottom_h = left_h = left + width + 0.02

    rect_scatter = [left, bottom, width, height]
    rect_histx = [left, bottom_h, width, 0.2]
    rect_histy = [left_h, bottom, 0.2, height]

    plt.figure(1, figsize=(8, 8))

    axScatter = plt.axes(rect_scatter)
    plt.xlabel("eruptions")
    plt.ylabel("waiting")
    axHistx = plt.axes(rect_histx)
    axHisty = plt.axes(rect_histy)
    axHistx.xaxis.set_major_formatter(nullfmt)
    axHisty.yaxis.set_major_formatter(nullfmt)

    axScatter.scatter(x, y)

    # now determine nice limits by hand:
    axHistx.hist(x, bins=20)
    axHisty.hist(y, bins=20, orientation="horizontal")

    plt.savefig("Scatter_Plot.png")
    plt.show()

def plot_p(frame):
    sol=Solution(frame,file_format='petsc',read_aux=False,path='./_output/_p/',file_prefix='claw_p')
    x=sol.state.grid.x.centers; y=sol.state.grid.y.centers
    mx=len(x); my=len(y)
    
    mp=sol.state.num_eqn    
    yy,xx = np.meshgrid(y,x)

    p=sol.state.q[0,:,:]
    fig = pl.figure(figsize=(8, 3.5))
    #pl.title("t= "+str(sol.state.t),fontsize=20)
    pl.xticks(size=20); pl.yticks(size=20)
    pl.xlabel('x',fontsize=20); pl.ylabel('y',fontsize=20)
    #pl.pcolormesh(xx,yy,p_subxy,cmap=cm.OrRd)
    pl.pcolormesh(xx,yy,p,cmap='RdBu_r')
    pl.autoscale(tight=True)
    cb = pl.colorbar(ticks=[0.5,1,1.5,2]);
    
    #pl.clim(ticks=[0.5,1,1.5,2])
    imaxes = pl.gca(); pl.axes(cb.ax)
    pl.yticks(fontsize=20); pl.axes(imaxes)
    #pl.xticks(fontsize=20); pl.axes(imaxes)
    #pl.axis('equal')
    pl.axis('tight')
    fig.tight_layout()
    pl.savefig('./_plots_to_paper/sound-speed_FV_t'+str(frame)+'_pcolor.png')
    pl.close()

    def plot(self, data, topn):
        """ Plot data fit and residuals """
        from matplotlib import pylab as pl
        pl.axes([0.1, 0.4, 0.8, 0.4])   # leave room above the axes for the table
        self.fit_plot(data, topn=topn)

        pl.axes([0.1, 0.05, 0.8, 0.3]) 
        self.residual_plot(data, topn=topn)

def plot_scatter(x,y,show=1,nbins=100,xrange=None,yrange=None,title="",xtitle="",ytitle=""):

    from matplotlib.ticker import NullFormatter

    # the random data

    nullfmt   = NullFormatter()         # no labels

    # definitions for the axes
    left, width    = 0.1, 0.65
    bottom, height = 0.1, 0.65
    bottom_h = left_h = left+width+0.02

    rect_scatter = [left, bottom, width, height]
    rect_histx   = [left, bottom_h, width, 0.2]
    rect_histy   = [left_h, bottom, 0.2, height]

    # start with a rectangular Figure
    fig = plt.figure(figsize=(8,8))

    axScatter = plt.axes(rect_scatter)
    axHistx = plt.axes(rect_histx)
    axHisty = plt.axes(rect_histy)

    # no labels
    axHistx.xaxis.set_major_formatter(nullfmt)
    axHisty.yaxis.set_major_formatter(nullfmt)

    # now determine nice limits by hand:
    binwidth = np.array([x.max() - x.min(), y.max() - y.min()]).max() / nbins

    if xrange == None:
        xrange = np.array((x.min(), x.max()))
    if yrange == None:
        yrange = np.array((y.min(), y.max()))

    # the scatter plot:
    axScatter.scatter(x, y, marker='.', edgecolor='b', s=.1)

    axScatter.set_xlabel(xtitle)
    axScatter.set_ylabel(ytitle)

    axScatter.set_xlim( xrange )
    axScatter.set_ylim( yrange )

    bins_x = np.arange(xrange[0], xrange[1] + binwidth, binwidth)
    axHistx.hist(x, bins=nbins)
    axHisty.hist(y, bins=nbins, orientation='horizontal')

    axHistx.set_xlim( axScatter.get_xlim() )

    axHistx.set_title(title)
    axHisty.set_ylim( axScatter.get_ylim() )


    if show: plt.show()

    return fig

def newCreateAndSaveMultilineFig(xDataList, yDataList, xLabel="", yLabel="",
                                 figFileRoot="", fileExt='.png', xMin=0,
                                 xMax=0, yMin=0, yMax=0, legendFlag=1,
                                 legendFont=12, traceNameList=[],
                                 legLoc=(0, 0)):
    """This subroutine saves a figure with multiple lines."""

    figFileName = figFileRoot + fileExt
    colorDict = createColorDictWithDashes()

    if xMax == 0:
        curMax = 0
        for n in range(0, len(xDataList)):
            if type(xDataList[n]) == list:
                if max(xDataList[n]) > curMax:
                    curMax = max(xDataList[n])
            else:
                if xDataList[n].any() > curMax:
                    curMax = max(xDataList[n])
        xMax = curMax

    if yMax == 0:
        curMax = 0
        for n in range(0, len(yDataList)):
            if type(yDataList[n]) == list:
                if max(yDataList[n]) > curMax:
                    curMax = max(yDataList[n])
            else:
                if yDataList[n].any() > curMax:
                    curMax = max(yDataList[n])
        yMax = curMax

    plt.axes([0.1, 0.1, 0.71, 0.8])
    if traceNameList == []:
        for n in range(0, len(xDataList)):
            traceNameList.append("Trace_" + str(n))

    for n in range(0, len(xDataList)):
        xData = convert_list_to_array(xDataList[n])
        yData = convert_list_to_array(yDataList[n])
        tempPlot = plt.plot(
            xData, yData, colorDict[str(n + 1)], hold="True",
            label=traceNameList[n])

    plt.xlabel(xLabel)
    plt.ylabel(yLabel)
    plt.xlim(xMin, xMax)
    plt.ylim(yMin, yMax)
    plt.rc("legend", fontsize=legendFont)

    if legendFlag == 1:
        if legLoc != (0, 0):
            print(legLoc)
            plt.legend(loc=legLoc)
        else:
            plt.legend()
    plt.savefig(figFileName, dpi=300)
    plt.clf()

    def save_wave_function_pic(self, pic_out, maxel=None, figuresize=(20,10), **kwrds):
        if figuresize:
            plt.figure(figsize=figuresize) # in inches!

        self.plot_wave_function(maxel)
        plt.axes().set_aspect('equal')
        plt.savefig(pic_out)
        plt.close()
        return None

def image(Z,xnew,ynew,my_cmap=None,aspect='equal'):
	"""
Creates pretty image. You need to specify:
	"""
	imshow(log10(Z),extent=[xnew[0],xnew[-1],ynew[0],ynew[-1]], cmap=my_cmap)
	pylab.axes().set_aspect('equal')
	colorbar()
	circle2=Circle((0,0),1,color='k')
	gca().add_artist(circle2)
	savefig('tmp.png',transparent=True,dpi=150)

def ana():
    pl.figure()
    pl.clf()
    pl.axes([0.16,0.2,0.95-0.15,0.95-0.2])
    pl.plot(x,y1,'g:',label=r'$\psi$')
    pl.plot(x,y2,'-b',label=r'$\psi^*\psi$')
    pl.xlabel(r'$y$ [nm]')
    pl.ylabel(r'$\pi$ [arb. units]')
    ax = pl.gca()
    ax.yaxis.set_label_coords(xshift, 0.5)
    # pl.legend()
    pl.savefig('images/analytical3.pdf',transparent='true')

def plot_q(frame,file_prefix='claw',path='./_output/',plot_pcolor=True,plot_slices=True,slices_xlimits=None):
    import sys
    sys.path.append('.')
    import sw_eqns

    sol=Solution(frame,file_format='petsc',read_aux=False,path=path,file_prefix=file_prefix)
    x=sol.state.grid.x.centers; y=sol.state.grid.y.centers
    mx=len(x); my=len(y)

    h=sol.state.q[0,:,:]
    b=sw_eqns.bathymetry(x,y)[0,:,:]
    eta=h+b

    yy,xx = np.meshgrid(y,x)

    if frame < 10:
        str_frame = "000"+str(frame)
    elif frame < 100:
        str_frame = "00"+str(frame)
    elif frame < 1000:
        str_frame = "0"+str(frame)
    else:
        str_frame = str(frame)

    if plot_pcolor:
        pl.pcolormesh(xx,yy,eta)
        pl.title("t= "+str(sol.state.t),fontsize=20)
        pl.xlabel('x',fontsize=20); pl.ylabel('y',fontsize=20)
        pl.xticks(size=20); pl.yticks(size=20)
        cb = pl.colorbar(); 
        #pl.clim(colorbar_min,colorbar_max); 
        imaxes = pl.gca(); pl.axes(cb.ax)
        pl.yticks(fontsize=20); pl.axes(imaxes)
        pl.axis([np.min(x),np.max(x),np.min(y),np.max(y)])
        pl.savefig('./_plots/eta_'+str_frame+'.png')
        #pl.show()
        pl.close()
    if plot_slices:
        pl.figure(figsize=(8,3))
        pl.plot(x,eta[:,3*my/4.],'-r',lw=1)
        pl.plot(x,eta[:,my/4.],'--b',lw=1)
        #pl.plot(x,b[:,my/4],'-k',lw=1)
        #pl.plot(x,b[:,3*my/4],'--k',lw=1)
        pl.title("t= "+str(sol.state.t),fontsize=20)
        pl.xlabel('x',fontsize=20)
        pl.ylabel('Surface',fontsize=20)
        pl.xticks(size=20); pl.yticks(size=20)
        #pl.ylim([0.9998,1.0002])
        pl.xlim([0.0,20])
        if slices_xlimits is not None:
            pl.axis([slices_xlimits[0],slices_xlimits[1],np.min(p),np.max(p)])
        pl.savefig('./_plots/eta_'+str_frame+'_slices.png')
        pl.close()