def VisualizeAlm(alm,figno=1,max_l=None):
    """ Visualize a healpy a_lm vector """
    lmax = hp.Alm.getlmax(f_lm.size)
    l,m = hp.Alm.getlm(lmax)
    mag = np.zeros([lmax+1,lmax+1])
    phs = np.zeros([lmax+1,lmax+1])
    mag[m,l] = np.abs(alm)
    phs[m,l] = np.angle(alm)
    cl = hp.alm2cl(alm)
    # Decide the range of l to plot
    if max_l != None:
        max_l = (max_l if (max_l <= lmax) else lmax)
    else:
        max_l = lmax 
    print max_l
    plt.figure(figno)
    plt.clf()
    plt.subplot(211)
    plt.imshow(mag[0:max_l,0:max_l],interpolation='nearest',origin='lower')
    plt.colorbar()
    plt.subplot(212)
    plt.imshow(phs[0:max_l,0:max_l],interpolation='nearest',origin='lower')
    plt.colorbar()
    # plt.subplot(313)
    #plt.semilogy(cl[0:max_l])
    return {'mag':mag,'phs':phs,'cl':cl}

    def test_draw_xy_plot(self):
        """draw_xy_plot() properly produces an output html file
        """
        out_file = os.path.join(self.dir_name, "test.html")
        argv = (
            "p.plot -x x -y btrace ctrace -s o- --xlabel myxlabel "
            "--ylabel myylabel --title mytitle --theme darkgrid "
            "--context talk --palette muted -a .5 --nogrid "
            "--legend best --xlim 0 10 --ylim -10 10 "
            "--savefig {}".format(out_file)
        ).split()
        with patch("pandashells.lib.plot_lib.sys.argv", argv):
            pl.clf()
            df = pd.DataFrame({"x": range(10), "btrace": [-x for x in range(10)], "ctrace": [x for x in range(10)]})
            parser = argparse.ArgumentParser()
            arg_lib.add_args(parser, "io_in", "xy_plotting", "decorating", "example")

            parser.add_argument("-a", "--alpha", help="Set opacity", nargs=1, default=[1.0], type=float)
            args = parser.parse_args()
            plot_lib.draw_xy_plot(args, df)
            with open(out_file) as f:
                html = f.read()
                self.assertTrue("myxlabel" in html)
                self.assertTrue("myylabel" in html)
                self.assertTrue("mytitle" in html)
                self.assertTrue("btrace" in html)
                self.assertTrue("ctrace" in html)
                self.assertTrue("1" in html)
                self.assertTrue("10" in html)

def check_vpd_ks2_astrometry():
    """
    Check the VPD and quiver plots for our KS2-extracted, re-transformed astrometry.
    """
    catFile = workDir + '20.KS2_PMA/wd1_catalog.fits'
    tab = atpy.Table(catFile)

    good = (tab.xe_160 < 0.05) & (tab.ye_160 < 0.05) & \
        (tab.xe_814 < 0.05) & (tab.ye_814 < 0.05) & \
        (tab.me_814 < 0.05) & (tab.me_160 < 0.05)

    tab2 = tab.where(good)

    dx = (tab2.x_160 - tab2.x_814) * ast.scale['WFC'] * 1e3
    dy = (tab2.y_160 - tab2.y_814) * ast.scale['WFC'] * 1e3

    py.clf()
    q = py.quiver(tab2.x_814, tab2.y_814, dx, dy, scale=5e2)
    py.quiverkey(q, 0.95, 0.85, 5, '5 mas', color='red', labelcolor='red')
    py.savefig(workDir + '20.KS2_PMA/vec_diffs_ks2_all.png')

    py.clf()
    py.plot(dy, dx, 'k.', ms=2)
    lim = 30
    py.axis([-lim, lim, -lim, lim])
    py.xlabel('Y Proper Motion (mas)')
    py.ylabel('X Proper Motion (mas)')
    py.savefig(workDir + '20.KS2_PMA/vpd_ks2_all.png')

    idx = np.where((np.abs(dx) < 10) & (np.abs(dy) < 10))[0]
    print('Cluster Members (within dx < 10 mas and dy < 10 mas)')
    print(('   dx = {dx:6.2f} +/- {dxe:6.2f} mas'.format(dx=dx[idx].mean(),
                                                        dxe=dx[idx].std())))
    print(('   dy = {dy:6.2f} +/- {dye:6.2f} mas'.format(dy=dy[idx].mean(),
                                                        dye=dy[idx].std())))

def test_path():
  "generate and draw a random path"
  path = genpath()
  P.ion()
  P.clf()
  draw_path(P.gca(), path)
  P.draw()

def plotslice(pos,filename='',boxsize=100.):
    ng = pos.shape[0]
    M.clf()
    M.scatter(pos[ng/4,:,:,1].flatten(),pos[ng/4,:,:,2].flatten(),s=1.,lw=0.)
    M.axis('tight')
    if filename != '':
        M.savefig(filename)

def plot_ch():
    for job in jobs_orig:
        print "plane of", job.path
        pylab.clf()
        x_center = int((job.position(0)[0] + job.position(1)[0])/2)
        x_final = 50 + x_center
        #plane = np.concatenate((job.plane(y=50)[:, x_final:], 
        #                        job.plane(y=50)[:, :x_final]), axis=1)
        plane = job.plane(y=50)
        myplane = plane[plane < 0.0]
        p0 = myplane.min()
        p12 = np.median(myplane)
        p14 = np.median(myplane[myplane<p12])
        p34 = np.median(myplane[myplane>p12])
        p1 = myplane.max()
        contour_values = (p0, p14, p12, p34, p1)
        pylab.title(r'$u_x=%.4f,\  D_{-}=%.4f,\  D_{+}=%.4f,\ ch=%i$ ' %
                    (job.u_x, job.D_minus, job.D_plus, job.ch_objects))
        car = pylab.imshow(plane, vmin=-0.001, vmax=0.0, 
                           interpolation='nearest')
        pylab.contour(plane, contour_values, linestyles='dashed', 
                                             colors='white')
        pylab.grid(True)
        pylab.colorbar(car)
        #imgfilename = 'plane_r20-y50-u_x%.4fD%.4fch%03i.png' % \
        #              (job.u_x, job.D_minus, job.ch_objects)
        imgfilename = 'plane_%s.png' % job.job_id
        pylab.savefig(imgfilename)

    def __call__(self, n):
        if len(self.f.shape) == 3:
            # f = f[x,v,t], 2 dim in phase space
            ft = self.f[n,:,:]
            pylab.pcolormesh(self.X, self.V, ft.T, cmap = 'jet')
            pylab.colorbar()
            pylab.clim(0,0.38) # for Landau test case
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$v$', fontsize = 18)
            pylab.title('$N_x$ = %d, $N_v$ = %d, $t$ = %2.1f' % (self.x.N, self.v.N, self.it*self.t.width))
            pylab.savefig(self.path + self.filename)
            pylab.clf()
            return None

        if len(self.f.shape) == 2:
            # f = f[x], 1 dim in phase space
            ft = self.f[n,:]
            pylab.plot(self.x.gridvalues,ft,'ob')
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$f(x)$', fontsize = 18)
            pylab.savefig(self.path + self.filename)
            return None

def EBTransitPhase(tset, kid_x):
    ebp = atpy.Table('%s/eb_pars.txt' %dir, type = 'ascii')
    
    lc = tset.tables[1]   
    time = lc.TIME
    flux = lc.PDCSAP_FLUX
    nobs = len(time)
    lg = scipy.isfinite(time)
    pylab.figure(52)
    pylab.clf()
    pylab.plot(time[lg], flux[lg])
    npl = 2
    phase = scipy.zeros((npl, nobs))
    inTr = scipy.zeros((npl, nobs), 'int')
    period = ebp.P[ebp.KID == kid_x]
    for ipl in scipy.arange(npl):
        if ipl == 0: t0 = ebp.Ep1[ebp.KID == kid_x]
        if ipl == 1: t0 = ebp.Ep2[ebp.KID == kid_x]
        if ipl == 0: dur = ebp.Dur1[ebp.KID == kid_x]
        if ipl == 1: dur = ebp.Dur2[ebp.KID == kid_x]
        dur /= period
        counter = 0
        while (time[lg] - t0).min() < 0:
            t0 -= period
            counter += 1
            if counter > 1000: break
        ph = ((time - t0) % period) / period
        ph[ph < -0.5] += 1
        ph[ph > 0.5] -= 1
        phase[ipl,:] = ph
        inTr[ipl,:] = (abs(ph) <= dur/1.5)
    return phase, inTr

 def newFigLayer():
     pylab.clf()
     pylab.figure(figsize=(8, 8))
     pylab.axes([0.15, 0.15, 0.8, 0.81])
     pylab.axis([0.6, -0.4, -0.4, 0.6])
     pylab.xlabel(r"$\Delta$\textsf{RA from Sgr A* (arcsec)}")
     pylab.ylabel(r"$\Delta$\textsf{Dec. from Sgr A* (arcsec)}")

    def dovis(self):
        """
        Do runtime visualization. 
        """

        pylab.clf()

        phi = self.cc_data.get_var("phi")

        myg = self.cc_data.grid

        pylab.imshow(numpy.transpose(phi[myg.ilo:myg.ihi+1,
                                         myg.jlo:myg.jhi+1]), 
                     interpolation="nearest", origin="lower",
                     extent=[myg.xmin, myg.xmax, myg.ymin, myg.ymax])

        pylab.xlabel("x")
        pylab.ylabel("y")
        pylab.title("phi")

        pylab.colorbar()
        
        pylab.figtext(0.05,0.0125, "t = %10.5f" % self.cc_data.t)

        pylab.draw()

def field_map_ndar(ndar_field,t,ar_coorx,ar_coory,X,image_out,variable):

    ar_field=ndar_field[t,:]
    max_val=int(np.max(ndar_field))
    if variable==4:
        max_val=100.
    xmin=min(ar_coorx);xmax=max(ar_coorx)
    ymin=min(ar_coory);ymax=max(ar_coory)
    step=X
    nx=(xmax-xmin)/step+1
    ny=(ymax-ymin)/step+1

    ar_indx=np.array((ar_coorx-xmin)/step,int)
    ar_indy=np.array((ar_coory-ymin)/step,int)

    ar_map=np.ones((ny,nx))*-99.9
    ar_map[ar_indy,ar_indx]=ar_field

    ar_map2 = M.masked_where(ar_map <0, ar_map)
    ut.check_file_exist(image_out)

    pl.clf()
    pl.axes(axisbg='gray')
    pl.imshow(ar_map2, cmap=pl.cm.RdBu,
              interpolation='Nearest', origin='lower', vmax=max_val, vmin=0)
    pl.title('time step= '+ut.string(t,len(str(t))))
    pl.colorbar()
    pl.savefig(image_out)

def plot_mcmc_results(chain):
    # Pull m and b arrays out of the Markov chain.
    mm = [m for b,m in chain]
    bb = [b for b,m in chain]
    # Scatterplot of m,b posterior samples
    plt.clf()
    plt.contour(bgrid, mgrid, posterior, pdf_contour_levels(posterior))
    plt.plot(bb, mm, 'b.', alpha=0.1)
    plot_mb_setup()
    plt.show()
    
    # Histograms
    import triangle
    triangle.corner(chain, labels=['b','m'], extents=[0.99]*2)
    plt.show()
    
    # Traces
    plt.clf()
    plt.subplot(2,1,1)
    plt.plot(mm, 'k-')
    plt.ylim(mlo,mhi)
    plt.ylabel('m')
    plt.subplot(2,1,2)
    plt.plot(bb, 'k-')
    plt.ylabel('b')
    plt.ylim(blo,bhi)
    plt.show()

def plot_vector_diff_F160W():
    """
    Using the xym2mat analysis on the F160W filter in the 2010 data set,
    plot up the positional offset vectors for each reference star in
    each star list relative to the average list. This should show
    us if there are systematic problems with the distortion solution
    or PSF variations.
    """
    x, y, m, xe, ye, me, cnt = load_catalog()

    dx = x - x[0]
    dy = y - y[0]
    dr = np.hypot(dx, dy)

    py.clf()
    py.subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95)

    for ii in range(1, x.shape[0]):
        idx = np.where((x[ii,:] != 0) & (y[ii,:] != 0) & (dr[ii,:] < 0.1) &
                       (xe < 0.05) & (ye < 0.05))[0]

        py.clf()
        q = py.quiver(x[0,idx], y[0,idx], dx[ii, idx], dy[ii, idx], scale=1.0)
        py.quiverkey(q, 0.9, 0.9, 0.03333, label='2 mas', color='red')
        py.title('{0} stars in list {1}'.format(len(idx), ii))
        py.xlim(0, 4500)
        py.ylim(0, 4500)

        foo = input('Continue?')
        if foo == 'q' or foo == 'Q':
            break

    def test_pixpsf(self):
        tim,tinf = get_tractor_image_dr8(94, 2, 520, 'i', psf='kl-pix',
                                         roi=[500,600,500,600], nanomaggies=True)
        psf = tim.getPsf()
        print 'PSF', psf

        for i,(dx,dy) in enumerate([
                (0.,0.), (0.2,0.), (0.4,0), (0.6,0),
                (0., -0.2), (0., -0.4), (0., -0.6)]):
            px,py = 50.+dx, 50.+dy
            patch = psf.getPointSourcePatch(px, py)
            print 'Patch size:', patch.shape
            print 'x0,y0', patch.x0, patch.y0
            H,W = patch.shape
            XX,YY = np.meshgrid(np.arange(W), np.arange(H))
            im = patch.getImage()
            cx = patch.x0 + (XX * im).sum() / im.sum()
            cy = patch.y0 + (YY * im).sum() / im.sum()
            print 'cx,cy', cx,cy
            print 'px,py', px,py

            self.assertLess(np.abs(cx - px), 0.1)
            self.assertLess(np.abs(cy - py), 0.1)
            
            plt.clf()
            plt.imshow(patch.getImage(), interpolation='nearest', origin='lower')
            plt.title('dx,dy %f, %f' % (dx,dy))
            plt.savefig('pixpsf-%i.png' % i)

 def savepng(pre, img, title=None, **kwargs):
     fn = '%s-%s.png' % (pre, idstr)
     print 'Saving', fn
     plt.clf()
     plt.imshow(img, **kwargs)
     ax = plt.axis()
     if debug:
         print len(xplotx),len(allobjx)
         for i,(objx,objy,objc) in enumerate(zip(allobjx,allobjy,allobjc)):
             plt.plot(objx,objy,'-',c=objc)
             tempx = []
             tempx.append(xplotx[i])
             tempx.append(objx[0])
             tempy = []
             tempy.append(xploty[i])
             tempy.append(objy[0])
             plt.plot(tempx,tempy,'-',c='purple')
         plt.plot(pointx,pointy,'y.')
         plt.plot(xplotx,xploty,'xg')
     plt.axis(ax)
     if title is not None:
         plt.title(title)
     plt.colorbar()
     plt.gray()
     plt.savefig(fn)

    def compute_parameters(input_data):
        """
        Build the required HTML response to be displayed.
        """
        figure = pylab.figure(figsize=(7, 8))
        pylab.clf()

        matrix = input_data.ordered_weights
        matrix_size = input_data.number_of_regions
        labels = input_data.ordered_labels
        plot_title = "Connection Strength"

        order = numpy.arange(matrix_size)
        # Assumes order is shape (number_of_regions, )
        order_rows = order[:, numpy.newaxis]
        order_columns = order_rows.T

        axes = figure.gca()
        img = axes.matshow(matrix[order_rows, order_columns])
        axes.set_title(plot_title)
        figure.colorbar(img)

        if labels is None:
            return

        axes.set_yticks(numpy.arange(matrix_size))
        axes.set_yticklabels(list(labels[order]), fontsize=8)
        axes.set_xticks(numpy.arange(matrix_size))
        axes.set_xticklabels(list(labels[order]), fontsize=8, rotation=90)

        figure.canvas.draw()
        parameters = dict(mplh5ServerURL=TvbProfile.current.web.MPLH5_SERVER_URL,
                          figureNumber=figure.number, showFullToolbar=False)
        return parameters, {}

    def show(self,filename=None):
        self.fig = fig = P.figure(self.fignum,self.figsize)
        P.clf()

        header = self.maps[0].header

        self.grid = grid = ImageGrid(fig, (1, 1, 1), 
                                     nrows_ncols = (self.nrows, self.ncols),
                                     share_all=True,
                                     axes_class=(pywcsgrid2.Axes,
                                                 dict(header=header)),
                                    **self.grid_kwargs)

        for i,(map,lower,upper) in enumerate(zip(self.maps,self.lower_energies,self.upper_energies)):
            map.show(axes=grid[i], cax=grid[i].cax)
            format_energy=lambda x: '%.1f' % (x/1000.) if x < 1000 else '%.0f' % (x/1000.)
            lower_string=format_energy(lower)
            upper_string=format_energy(upper)
            grid[i].add_inner_title("%s-%s GeV" % (lower_string,upper_string), loc=2)

        if self.title is not None:
            fig.suptitle(self.title)

        tight_layout(fig)
        if filename is not None: P.savefig(filename)

def graphTSPResults(resultsDir,numberOfTours):
    x=[]
    y=[]
    files=[open("%s/objectiveFunctionReport.txt" % resultsDir),
       open("%s/fitnessReport.txt" % resultsDir)]
    for f in files:
        x.append([])
        y.append([])
        i=len(x)-1
        for line in f:
            line=line.split(',')
            if line[0] != "gen":
                x[i].append(int(line[0]))
                y[i].append(float(line[1] if i==1 else line[2]))
            
    ylen=len(y[0])
    pl.subplot(2,1,1)
    pl.plot(x[0],y[0],'bo')
    pl.ylabel('Minimum Distance')
    pl.title("TSP with a %s City Tour" % numberOfTours)
    pl.annotate("{0:,}".format(y[0][0]),xy=(x[0][0],y[0][0]),  xycoords='data',
             xytext=(30, -30), textcoords='offset points',
             arrowprops=dict(arrowstyle="->") )
    pl.annotate("{0:,}".format(y[0][ylen-1]),xy=(x[0][ylen-1],y[0][ylen-1]),  xycoords='data',
             xytext=(-30,30), textcoords='offset points',
             arrowprops=dict(arrowstyle="->") )

    pl.subplot(2,1,2)
    pl.plot(x[1],y[1],'go')
    pl.xlabel('Generation')
    pl.ylabel('Fitness')
    pl.savefig("%s/tsp_result.png" % resultsDir)
    pl.clf()

def do_Plot_item_a(Ustar, U, grad_phi):
    """ plotting stuff: here we plot the first row with Ustar and
        the second rows with Ud and grad phi, to see how they sum up """
    pylab.clf()
    pylab.cla()
    
    f = pylab.figure() 
    f.text(.5, .95, r"Top: $U = U_d + \nabla \phi$, where "  r"$U_d$ =  - $\sin^2 ( \pi x)\sin(2 \pi y)\hat x$ + $\sin^2(\pi y)\sin(2\pi x)\hat y$ and $\phi = \frac{1}{10} \cos(2\pi y) \cos(2\pi x$) ",  horizontalalignment='center', size=9)

    pylab.subplot(221) 
    pylab.imshow(Ustar[0])
    pylab.xlabel("Vector "r"$U_d$ in $\hat x$:" ,size=8)
    pylab.ylabel("U in terms of # of cells in " r"$\hat x$", size=8)

    pylab.subplot(222)
    pylab.imshow(Ustar[1])
    pylab.xlabel(r"Vector $\nabla \phi$ in $ \hat x $:", size=8)
    pylab.ylabel("U in terms of # of cells in " r"$\hat y$",size=8)

    pylab.subplot(223)
    pylab.imshow(U [0])
    pylab.ylabel("# of cells",size=8)
    pylab.xlabel("# of cells",size=8)
    
    
    pylab.subplot(224)
    pylab.imshow(grad_phi [0])
    pylab.xlabel("# of cells",size=8)
    pylab.ylabel("# of cells",size=8)
    pylab.savefig("plots/item_a_Ustar.png")
    
    return 0

def displayRetirementWithMonthlies(monthlies, rate, terms):
    plt.figure('retireMonth')
    plt.clf()
    for monthly in monthlies:
        xvals, yvals = retire(monthly, rate, terms)
        plt.plot(xvals, yvals, label = 'retire:'+str(monthly))
        plt.legend(loc = 'upper left')