def plotBode(fAg, lowerFreq, higherFreq = None):
    """'  plot Bode diagram using matplotlib
        Default frequency width is 3 decades
    '"""
    import pylab as py
    #import pdb; pdb.set_trace()
    if ( higherFreq == None):
        rangeAt = 1000.0  # 3 decade
    else:
        assert higherFreq > lowerFreq
        rangeAt = float(higherFreq)/lowerFreq

    N = 128
    lstScannAngFreqAt = [2*sc.pi*1j*lowerFreq*sc.exp(
                                                sc.log(rangeAt)*float(k)/N)
                                for k in range(N)]

    t              = [        lowerFreq*sc.exp(sc.log(rangeAt)*float(k)/N)
                                for k in range(N)]

    py.subplot(211)
    py.loglog( t, [(abs(fAg(x))) for x in lstScannAngFreqAt]  )
    py.ylabel('gain')
    py.grid(True)

    py.subplot(212)
    py.semilogx( t, [sc.arctan2(fAg(zz).imag, fAg(zz).real)
                        for zz in lstScannAngFreqAt])
    py.ylabel('phase')
    py.grid(True)
    py.gca().xaxis.grid(True, which='minor')  # minor grid on too
    py.show()

 def test_set_plot_limits(self):
     """set_limits() properly sets limits
     """
     args = MagicMock(savefig="", xlim=[-2, 2], ylim=[-3, 3])
     plot_lib.set_limits(args)
     self.assertEqual(pl.gca().get_xlim(), (-2.0, 2.0))
     self.assertEqual(pl.gca().get_ylim(), (-3.0, 3.0))

def trace(data, name, format='png', datarange=(None, None), suffix='', path='./', rows=1, columns=1, 
    num=1, last=True, fontmap = None, verbose=1):
    """
    Generates trace plot from an array of data.

    :Arguments:
        data: array or list
            Usually a trace from an MCMC sample.

        name: string
            The name of the trace.
            
        datarange: tuple or list
            Preferred y-range of trace (defaults to (None,None)).

        format (optional): string
            Graphic output format (defaults to png).

        suffix (optional): string
            Filename suffix.

        path (optional): string
            Specifies location for saving plots (defaults to local directory).
            
        fontmap (optional): dict
            Font map for plot.

    """

    if fontmap is None: fontmap = {1:10, 2:8, 3:6, 4:5, 5:4}

    # Stand-alone plot or subplot?
    standalone = rows==1 and columns==1 and num==1

    if standalone:
        if verbose>0:
            print_('Plotting', name)
        figure()

    subplot(rows, columns, num)
    pyplot(data.tolist())
    ylim(datarange)

    # Plot options
    title('\n\n   %s trace'%name, x=0., y=1., ha='left', va='top', fontsize='small')

    # Smaller tick labels
    tlabels = gca().get_xticklabels()
    setp(tlabels, 'fontsize', fontmap[rows/2])

    tlabels = gca().get_yticklabels()
    setp(tlabels, 'fontsize', fontmap[rows/2])

    if standalone:
        if not os.path.exists(path):
            os.mkdir(path)
        if not path.endswith('/'):
            path += '/'
        # Save to file
        savefig("%s%s%s.%s" % (path, name, suffix, format))

def plot_file_color(base, thin=True, start=0, size=14, save=False):
    conf, track, pegs = load(base)

    fig = pl.figure(figsize=(size,size*conf['top']/conf['wall']))

    track = track[start:]
    x = track[:,0];   y = track[:,1]
    t = np.linspace(0,1,x.shape[0])
    points = np.array([x,y]).transpose().reshape(-1,1,2)
    segs = np.concatenate([points[:-1],points[1:]],axis=1)
    lc = LineCollection(segs, linewidths=0.25, cmap=pl.cm.coolwarm)
    lc.set_array(t)
    pl.gca().add_collection(lc)

    #pl.scatter(x, y, c=np.arange(len(x)),linestyle='-',cmap=pl.cm.coolwarm)
    #pl.plot(track[-1000000:,0], track[-1000000:,1], '-', linewidth=0.0125, alpha=0.8)
    for peg in pegs:
        pl.gca().add_artist(pl.Circle(peg, conf['radius'], color='k', alpha=0.3))
    pl.xlim(0, conf['wall'])
    pl.ylim(0, conf['top'])
    pl.xticks([])
    pl.yticks([])
    pl.tight_layout()
    pl.show()
    if save:
        pl.savefig(base+".png", dpi=200)

def show_image(*imgs): 
    for idx, img in enumerate(imgs):
        subplot = 101 + len(imgs)*10 +idx
        pl.subplot(subplot)
        pl.imshow(img, cmap=pl.cm.gray)   
        pl.gca().set_axis_off()    
    pl.subplots_adjust(0.02, 0, 0.98, 1, 0.02, 0)        

def _plotResults(naDist1, naDist2, lfOneKnn, lf5Knn):
    plt.clf()

    plt.subplot(311)
    plt.scatter(naDist1[:, 0], naDist1[:, 1])

    plt.scatter(naDist2[:, 0], naDist2[:, 1], color="r")

    # plt.ylabel( 'Feature 2' )
    # plt.xlabel( 'Feature 1' )
    # gca().annotate( '', xy=( .8, 0 ), xytext=( -.3 , 0 ), arrowprops=dict(facecolor='red', shrink=0.05) )
    gca().annotate("", xy=(0.7, 0), xytext=(1.5, 0), arrowprops=dict(facecolor="black", shrink=0.05))
    plt.title("Data Distribution")

    plt.subplot(312)

    plt.plot(range(len(lfOneKnn)), lfOneKnn)

    plt.ylabel("1-KNN Value")
    # plt.xlabel( 'Distribution Merge' )

    plt.title("1-KNN Performance")

    plt.subplot(313)

    plt.plot(range(len(lf5Knn)), lf5Knn)

    plt.ylabel("% Correct Classification")
    # plt.xlabel( 'Distribution Merge' )

    plt.title("5-KNN Performance")

    plt.subplots_adjust()

    plt.show()

def auto_shift(offset):
    """
    Return a y-offset coordinate transform for the current axes.

    Each call to auto_shift increases the y-offset for the next line by
    the given number of points (with 72 points per inch).

    Example::

        from matplotlib import pyplot as plt
        from bumps.plotutil import auto_shift
        trans = auto_shift(plt.gca())
        plot(x, y, hold=True, trans=trans)
    """
    from matplotlib.transforms import ScaledTranslation
    import pylab
    ax = pylab.gca()
    if ax.lines and hasattr(ax, '_auto_shift'):
        ax._auto_shift += offset
    else:
        ax._auto_shift = 0
    trans = pylab.gca().transData
    if ax._auto_shift:
        trans += ScaledTranslation(0, ax._auto_shift/72.,
                                   pylab.gcf().dpi_scale_trans)
    return trans

def plotNumVariants(depth):
    plt.figure(figsize=(7, 5), dpi=300)
    depth.apply(lambda x: x.dropna().size / 1000.0).unstack("POP").plot.bar(ax=plt.gca())
    plt.gcf().subplots_adjust(bottom=0.30)
    plt.ylabel(r"Num of variants ($\times$1000)")
    renameLegend(ax=plt.gca())
    plt.savefig(pathPlots + "numVariants.pdf")

def plot(filename):
    """
    Read and print all command line arguments
    """
    import pylab

    canvas = pylab.gcf().canvas
    d = data(filename)
    if len(d.v.shape) > 2:
        pylab.gca().pcolormesh(d.v[0, :, :])
        pylab.xlabel(d.xlabel)
        pylab.ylabel(d.ylabel)
    elif len(d.v.shape) > 1:
        if filename.lower().endswith('bt4'):
            offset = 1
        else:
            offset = 0
        pylab.gca().pcolorfast(d.v[:, offset:])
        pylab.xlabel(d.xlabel)
        pylab.ylabel(d.ylabel)
    else:
        pylab.plot(d.x, d.v)
        pylab.xlabel(d.xlabel)
        pylab.ylabel(d.vlabel)
    pylab.show()

def _show_rates(rate, wo, wt, attenuator, tau_NP, tau_P):
    import pylab

    #pylab.figure()
    pylab.errorbar(rate, wt[0], yerr=wt[1], fmt='g.', label='attenuated')
    pylab.errorbar(rate, wo[0], yerr=wo[1], fmt='b.', label='unattenuated')

    pylab.xscale('log')
    pylab.yscale('log')
    pylab.xlabel('incident rate (counts/second)')
    pylab.ylabel('observed rate (counts/second)')
    pylab.legend(loc='best')
    pylab.grid(True)
    pylab.plot(rate, rate/attenuator, 'g-', label='target')
    pylab.plot(rate, rate, 'b-', label='target')

    Ipeak, Rpeak = peak_rate(tau_NP=tau_NP, tau_P=tau_P)
    if rate[0] <= Ipeak <= rate[-1]:
        pylab.axvline(x=Ipeak, ls='--', c='b')
        pylab.text(x=Ipeak, y=0.05, s=' %g'%Ipeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_xaxis_transform())
    if False:
        pylab.axhline(y=Rpeak, ls='--', c='b')
        pylab.text(y=Rpeak, x=0.05, s=' %g\n'%Rpeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_yaxis_transform())

def learning_curve_plot(
    predictors, valid_X, valid_Y, valid_W=None, train_X=None, train_Y=None, train_W=None, log_scale=False
):

    if hasattr(valid_W, "values"):
        valid_W = valid_W.values
    if train_W != None and hasattr(train_W, "values"):
        train_W = weights.values

    with_train = True if (train_X != None and valid_X != None) else False

    try:
        predictors = predictors.items()
    except:
        predictors = {"": predictors}.items()

    for name, predictor in predictors:
        iterations = np.arange(1, predictor.n_estimators + 1)
        p, = plt.plot(
            iterations, _step_wise_performance(predictor, valid_X, valid_Y, valid_W), "-", label=name + " (test)"
        )

        if with_train:
            plt.plot(
                iterations,
                _step_wise_performance(predictor, train_X, train_Y, train_W),
                "--",
                color=p.get_color(),
                label=name + " (train)",
            )

        plt.legend(loc="best")

    if log_scale:
        plt.gca().set_xscale("log")

def plot(D,title):
    im=plt.imshow(D,interpolation='nearest',cmap='Reds')
    plt.gca().xaxis.tick_top()
    x=np.arange(D.index.shape[0])
    plt.colorbar(im)
    plt.gca().tick_params(axis='both', which='major', labelsize=10)
    plt.title(title,y=1.03)

def plot_INT_footprint(center_ra,center_dec):
    #using full detector sizes for now because 
    detector_dra = 4100.*0.33/3600. # 2154 pixels * 0.33"/pix, /3600 to get deg
    detector_ddec = 2048.*0.33/3600. # 2154 pixels * 0.33"/pix, /3600 to get deg
    # draw footprint of chip 4
    rect= plt.Rectangle((center_ra-detector_dra/2.,center_dec-detector_ddec/2.), detector_dra, detector_ddec,fill=False, color='k')
    plt.gca().add_artist(rect)
    # draw footprint of chip 3
    # assuming chip 3 is NORTH and a smidge WEST of chip 4
    offset_dec = detector_ddec+17./3600. # 17 arcsec gap in RA between 
    offset_ra = -9.5/3600. # 9.5 arcsec offset toward N
    rect= plt.Rectangle((center_ra+offset_ra-detector_dra/2.,center_dec+offset_dec-detector_ddec/2.), detector_dra, detector_ddec,fill=False, color='k')
    plt.gca().add_artist(rect)

    # draw footprint of chip 1
    # assuming chip 1 is SOUTH and a smidge EAST of chip 4
    offset_dec = -1*detector_ddec-22.7/3600. # 17 arcsec gap in RA between 
    offset_ra = +3.18/3600. # 9.5 arcsec offset toward N
    rect= plt.Rectangle((center_ra+offset_ra-detector_dra/2.,center_dec+offset_dec-detector_ddec/2.), detector_dra, detector_ddec,fill=False, color='k')
    plt.gca().add_artist(rect)

    # draw footprint of chip 2
    # assuming chip 2 is WEST of chip 4
    offset_dec = detector_ddec/2.-detector_dra-19.2/3600. # hard to explain
    offset_ra =  -.5*detector_dra-23.8/3600.# hard to explain
    # this chip is rotated 90 deg, so detecter_dra and detector_ddec are interchanged
    rect= plt.Rectangle((center_ra+offset_ra,center_dec+offset_dec), -1.*detector_ddec, detector_dra,fill=False, color='k')
    plt.gca().add_artist(rect)

    # adding guide camera
    offset_dec = -2*detector_ddec-(22.7+98.1)/3600. # hard to explain
    offset_ra =  detector_dra/2-(3.18+649.9)/3600.# hard to explain
    # this chip is rotated 90 deg, so detecter_dra and detector_ddec are interchanged
    rect= plt.Rectangle((center_ra+offset_ra,center_dec+offset_dec), -7./60., 7./60,fill=False, color='k')
    plt.gca().add_artist(rect)

def validate_vs_jwpsf_nircam():

    models = [ ('NIRCam','F200W', 'f200w_perfect_offset', '/Users/mperrin/software/jwpsf_v3.0/data/NIRCam/OPD/perfect_opd.fits', 0.034,True),
            ('NIRCam','F200W', 'f200w_perfect', '/Users/mperrin/software/jwpsf_v3.0/data/NIRCam/OPD/perfect_opd.fits', 0.034,False),
            ('NIRCam','F200W', 'f200w', '/Users/mperrin/software/jwpsf_v3.0/data/NIRCam/OPD/nircam_obs_w_rsrv1.fits', 0.034,True),
                ('MIRI','F1000W', 'f1000w', '/Users/mperrin/software/jwpsf_v3.0/data/MIRI/OPD/MIRI_OPDisim1.fits', 0.11,True)]


    fig = P.figure(1, figsize=(13,8.5), dpi=80)
    oversamp=4
    for params in models:

        nc = webbpsf_core.Instrument(params[0])
        nc.filter = params[1]
        nc.pupilopd = params[3] #'/Users/mperrin/software/jwpsf_v3.0/data/NIRCam/OPD/nircam_obs_w_rsrv1.fits'
        nc.pixelscale = params[4] #0.034 # this is wrong, but compute this way to match JWPSF exactly
        if params[5]:
            # offset by half a pixel to match the JWPSF convention
            nc.options['source_offset_r'] = params[4]/2 * N.sqrt(2)/oversamp  # offset half a pixel each in X and Y
            nc.options['source_offset_theta'] = -45


        jw_fn = 'jwpsf_%s_%s.fits' % (params[0].lower(), params[2].lower())
        my_fn = 'test_vs_' + jw_fn

        if not os.path.exists( my_fn):
            my_psf = nc.calcPSF(my_fn, oversample=oversamp, fov_pixels=512./oversamp)
        else:
            my_psf = fits.open(my_fn)

        jw_psf = fits.open(jw_fn)
        jw_psf[0].header.update('PIXELSCL', jw_psf[0].header['CDELT1']*3600)


        P.clf()
        #P.subplots_adjust(top=0.95, bottom=0.05, left=0.01, right=0.99)
        P.subplot(231)
        titlestr = "%s %s, \n"%  (params[0], params[2])
        poppy.display_PSF(my_psf, title=titlestr+"computed with WebbPSF" , colorbar=False)
        P.subplot(232)
        poppy.display_PSF(jw_psf, title=titlestr+"computed with JWPSF" , colorbar=False)
        P.subplot(233)
        poppy.display_PSF_difference(my_psf,jw_psf, title=titlestr+'Difference Image', colorbar=False)

        imagecrop = 30*params[4]

        P.subplot(234)
        poppy.display_PSF(my_psf, title=titlestr+"computed with WebbPSF", colorbar=False, imagecrop=imagecrop)
        centroid = poppy.measure_centroid(my_psf)
        P.gca().set_xlabel("centroid = (%.3f,%.3f)" % centroid)

        P.subplot(235)
        poppy.display_PSF(jw_psf, title=titlestr+"computed with JWPSF", colorbar=False, imagecrop=imagecrop)
        centroid = poppy.measure_centroid(jw_psf)
        P.gca().set_xlabel("centroid = (%.3f,%.3f)" % centroid)

        P.subplot(236)
        poppy.display_PSF_difference(my_psf,jw_psf, title='Difference Image', colorbar=False, imagecrop=imagecrop)

        P.savefig("results_vs_jwpsf_%s_%s.pdf" % (params[0], params[2]))

def plot_descur(RBC, RBS, ZBC, ZBS, max_m, max_n, phi_plane=0., nu=150, debug=1, **plkwargs):
    """ Companion plotting routine to read_descur_data
    """
    theta_arr = linspace(0,2*pi,nu,endpoint=True)
    r_c = 0*theta_arr
    z_c = 0*theta_arr
    r_c2 = 0*theta_arr
    z_c2 = 0*theta_arr
    for (i,theta) in enumerate(theta_arr):
        for ixm in range(0, max_m):
            if ixm==0: start_n=0
            else: start_n = -max_n
            for ixn in range(start_n, 1+max_n):
                arg = ixn*phi_plane + ixm*theta
                sinarg = sin(arg)
                cosarg = cos(arg)
                r_c[i] += (RBC[ixn,ixm]*cosarg) + (RBS[ixn,ixm]*sinarg)
                z_c[i] += (ZBC[ixn,ixm]*cosarg) + (ZBS[ixn,ixm]*sinarg)
                r_c2[i] += (RBC[ixn,ixm]*cosarg) #Shaun modification to show the stellarator symetric version
                z_c2[i] += (ZBS[ixn,ixm]*sinarg) #Shaun modification to show the stellarator symetric version

    pl.plot(r_c, z_c, **plkwargs)
    pl.plot(r_c2, z_c2, **plkwargs)
    pl.gca().set_aspect('equal')
    pl.show()
    if debug>3:
        import pdb; pdb.set_trace()
        'debugging, c to continue'

    def plot_fixed_x(self, x_values, x_derivative=0, steps=1000, smooth=0, simple='auto', ymin="auto", ymax="auto", format=True, clear=1):
        """
        plots the data at fixed x-value, so z vs x
        """
        if simple=='auto': simple=self.simple
        
        # get the min and max
        if ymin=="auto": ymin = self.ymin
        if ymax=="auto": ymax = self.ymax
        if clear: _pylab.gca().clear()

        if not type(x_values) in [type([]), type(_pylab.array([]))]: x_values = [x_values]

        for x in x_values:
    
            # define a new simple function to plot, then plot it
            def f(y): return self.evaluate(x, y, x_derivative, smooth, simple)
            _pylab_help.plot_function(f, ymin, ymax, steps, 0, False)

            # label it
            a = _pylab.gca()
            a.set_xlabel(self.ylabel)
            if x_derivative: a.set_ylabel(str(x_derivative)+" "+str(self.xlabel)+" derivative of "+self.zlabel)
            else:            a.set_ylabel(self.zlabel)
            a.set_title(self._path+"\nSpline array plot at fixed x = "+self.xlabel)
            a.get_lines()[-1].set_label("x ("+self.xlabel+") = "+str(x))

        if format: _s.format_figure()
        return a

 def plot_surface(self, debug=debug, nu=150, phi=0, hold=True, linecolor='gray', marker=',',**plkwargs):
     (r,z)=self.make_curve(nu=nu, phi=phi)
     plkwargs.update({'color':linecolor, 'linewidth':0.3})
     if (linecolor != None) and (linecolor != ''): 
         pl.plot(r,z, hold=hold, **plkwargs)
     pl.gca().set_aspect('equal')
     if marker != '': pl.plot(r,z,marker, hold=1)        

    def plot_fixed_y(self, y_values, x_derivative=0, steps=1000, smooth=0, simple='auto', xmin="auto", xmax="auto", format=True, clear=1):
        """
        plots the data at a fixed y-value, so z vs y
        """
        if simple=='auto': simple=self.simple


        # get the min and max
        if xmin=="auto": xmin = self.xmin
        if xmax=="auto": xmax = self.xmax
        if clear: _pylab.gca().clear()

        if not type(y_values) in [type([]), type(_pylab.array([]))]: y_values = [y_values]

        for y in y_values:
            # define a new simple function to plot, then plot it
            def f(x): return self.evaluate(x, y, x_derivative, smooth, simple)
            _pylab_help.plot_function(f, xmin, xmax, steps, 0, True)

            # label it
            a = _pylab.gca()
            th = "th"
            if x_derivative == 1: th = "st"
            if x_derivative == 2: th = "nd"
            if x_derivative == 3: th = "rd"
            if x_derivative: a.set_ylabel(str(x_derivative)+th+" "+self.xlabel+" derivative of "+self.zlabel+" spline")
            else:            a.set_ylabel(self.zlabel)
            a.set_xlabel(self.xlabel)
            a.set_title(self._path+"\nSpline array plot at fixed y "+self.ylabel)
            a.get_lines()[-1].set_label("y ("+self.ylabel+") = "+str(y))
        if format: _s.format_figure()
        return a

def plot_features(im, features, num_to_plot=100, colors=["blue"]):

  plt.imshow(im)

  for i in range(min(features.shape[0], num_to_plot)):
    x = features[i,0]
    y = features[i,1]
    scale = features[i,2]
    rot = features[i,3]


    color = colors[i % len(colors)]

    box = patches.Rectangle((-scale/2,-scale/2), scale, scale, 
      edgecolor=color, facecolor="none", lw=1)
    arrow = patches.Arrow(0, -scale/2, 0, scale, 
      width=10, edgecolor=color, facecolor="none")
    t_start = plt.gca().transData
    transform = mpl.transforms.Affine2D().rotate(rot).translate(x,y) + t_start

    box.set_transform(transform)
    arrow.set_transform(transform)
    plt.gca().add_artist(box)
    plt.gca().add_artist(arrow)


  plt.axis('off')
  plt.set_cmap('gray')
  plt.show()

def plotLL(fname='out4.npy'):
    plt.figure()
    h= np.linspace(0,1,21)
    g= np.linspace(0,1,21)
    m=np.linspace(0,2,21)
    d=np.linspace(0,2,21)
    out=np.load(fname)
    print np.nanmax(out),np.nanmin(out)
    rang=np.nanmax(out)-np.nanmin(out)
    maxloc= np.squeeze(np.array((np.nanmax(out)==out).nonzero()))
    H,G=np.meshgrid(h,g)
    print maxloc
    for mm in range(m.size/2):
        for dd in range(d.size/2):
            plt.subplot(10,10,(9-mm)*10+dd+1)
            plt.pcolormesh(h,g,out[:,:,mm*2,dd*2].T,
                           vmax=np.nanmax(out),vmin=np.nanmax(out)-rang/4.)
            plt.gca().set_xticks([])
            plt.gca().set_yticks([])
            if mm==maxloc[2]/2 and dd==maxloc[3]/2:
                plt.plot(h[maxloc[0]],g[maxloc[1]],'ow',ms=8)
            if dd==0:
                print mm,dd
                plt.ylabel('%.1f'%m[mm*2])
            if mm==0: plt.xlabel('%.1f'%d[dd*2])
    plt.title(fname[:6])