def Contour(X,Y,Z, label='', levels=None, cmapidx=0, colors=None, fmt='%g', lwd=1, fsz=10,
        inline=0, wire=True, cbar=True, zorder=None, markZero='', clabels=True):
    """
    Plot contour
    ============
    """
    L = None
    if levels != None:
        if not hasattr(levels, "__iter__"): # not a list or array...
            levels = linspace(Z.min(), Z.max(), levels)
    if colors==None:
        c1 = contourf (X,Y,Z, cmap=Cmap(cmapidx), levels=levels, zorder=None)
    else:
        c1 = contourf (X,Y,Z, colors=colors, levels=levels, zorder=None)
    if wire:
        c2 = contour (X,Y,Z, colors=('k'), levels=levels, linewidths=[lwd], zorder=None)
        if clabels:
            clabel (c2, inline=inline, fontsize=fsz)
    if cbar:
        cb = colorbar (c1, format=fmt)
        cb.ax.set_ylabel (label)
    if markZero:
        c3 = contour(X,Y,Z, levels=[0], colors=[markZero], linewidths=[2])
        if clabels:
            clabel(c3, inline=inline, fontsize=fsz)

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 contour(*args, **kwargs):
    """
    Make an SPH image of the given simulation and render it as contours.
    nlevels and levels are passed to pyplot's contour command.

    Other arguments are as for *image*.
    """

    import copy
    kwargs_image = copy.copy(kwargs)
    nlevels = kwargs_image.pop('nlevels',None)
    levels = kwargs_image.pop('levels',None)
    width = kwargs_image.get('width','10 kpc')
    kwargs_image['noplot']=True
    im = image(*args, **kwargs_image)
    res = im.shape

    units = kwargs_image.get('units',None)

    if isinstance(width, str) or issubclass(width.__class__, _units.UnitBase):
        if isinstance(width, str):
            width = _units.Unit(width)
        sim = args[0]
        width = width.in_units(sim['pos'].units, **sim.conversion_context())

    width = float(width)
    x,y = np.meshgrid(np.linspace(-width/2,width/2,res[0]),np.linspace(-width/2,width/2,res[0]))

    p.contour(x,y,im,nlevels=nlevels,levels=levels)

def contour(list):
    xrange = numpy.arange(-5, 5, 0.05)
    yrange = numpy.arange(-5, 5, 0.05)

    grid = matrix([[indicator(x, y, list) for y in yrange] for x in xrange])
    pylab.contour(xrange, yrange, grid, (-1.0, 0.0, 1.0), colors=('red', 'black', 'blue'), linewidths=(1, 3, 1))
    pylab.show()

        def doSubplot(multiplier=1.0, layout=(-1, -1)):
            if time_sec < 16200:
                xs, ys = xs_1, ys_1
                domain_bounds = bounds_1sthalf
                grid = grid_1
            else:
                xs, ys = xs_2, ys_2
                domain_bounds = bounds_2ndhalf
                grid = grid_2

            try:
                mo = ARPSModelObsFile("%s/%s/KCYS%03dan%06d" % (base_path, exp, min_ens, time_sec))
            except AssertionError:
                mo = ARPSModelObsFile("%s/%s/KCYS%03dan%06d" % (base_path, exp, min_ens, time_sec), mpi_config=(2, 12))
            except:
                print "Can't load reflectivity ..."
                mo = {'Z':np.zeros((1, 255, 255), dtype=np.float32)}

            pylab.contour(xs, ys, wind[exp]['w'][wdt][domain_bounds], levels=np.arange(2, 102, 2), styles='-', colors='k')
            pylab.contour(xs, ys, wind[exp]['w'][wdt][domain_bounds], levels=np.arange(-100, 0, 2), styles='--', colors='k')

            pylab.quiver(xs[thin], ys[thin], wind[exp]['u'][wdt][domain_bounds][thin], wind[exp]['v'][wdt][domain_bounds][thin])

            pylab.contourf(xs, ys, mo['Z'][0][domain_bounds], levels=np.arange(10, 85, 5), cmap=NWSRef, zorder=-10)

            grid.drawPolitical(scale_len=10)

            row, col = layout
            if col == 1:
                pylab.text(-0.075, 0.5, exp_names[exp], transform=pylab.gca().transAxes, rotation=90, ha='center', va='center', size=12 * multiplier)

def plot_debug(result):
    import pylab
    from pylab import contour, grid, show
    contours = (70.9/22.4, 70.9/22.4*3, 70.9/22.4*20)
    contour(result, contours)
    grid(True)
    show()

def test_masking():
    xi, yi = 20*random.rand(2, 100)
    hi = exp( -(xi-10.)**2/10.0 -(yi-10.)**2/5.0)
    h = grid.extrap_xy_to_grid(xi, yi, hi)
    pl.contour(grid.x_rho, grid.y_rho, h)
    print h
    pl.show()

def visualize_fits(fitsfile, figname, colorbar = False) :
    '''Saves fitsfile contents as an image, <figname>.png'''
    pl.figure()
    pl.contour(get_fits_obj(fitsfile))
    if colorbar:
        pl.colorbar()
    pl.savefig(figname+'.png')

def plot_haxby(activation, title):
    z = 25

    fig = pl.figure(figsize=(4, 5.4))
    fig.subplots_adjust(bottom=0., top=1., left=0., right=1.)
    pl.axis('off')
    # pl.title('SVM vectors')
    pl.imshow(np.rot90(mean_img[:, 4:58, z]), cmap=pl.cm.gray,
              interpolation='nearest')
    pl.imshow(np.rot90(activation[:, 4:58, z]), cmap=pl.cm.hot,
              interpolation='nearest')

    mask_house = nibabel.load(h.mask_house[0]).get_data()
    mask_face = nibabel.load(h.mask_face[0]).get_data()

    pl.contour(np.rot90(mask_house[:, 4:58, z].astype(np.bool)), contours=1,
            antialiased=False, linewidths=4., levels=[0],
            interpolation='nearest', colors=['blue'])

    pl.contour(np.rot90(mask_face[:, 4:58, z].astype(np.bool)), contours=1,
            antialiased=False, linewidths=4., levels=[0],
            interpolation='nearest', colors=['limegreen'])

    p_h = Rectangle((0, 0), 1, 1, fc="blue")
    p_f = Rectangle((0, 0), 1, 1, fc="limegreen")
    pl.legend([p_h, p_f], ["house", "face"])
    pl.title(title, x=.05, ha='left', y=.90, color='w', size=28)

def main():
    print("Running main method.")
    print("")
    
    # Generate the datapoints and plot them.
    classA, classB = gen_datapoints_helper()

    kernels = [LINEAR, POLYNOMIAL, RADIAL_BASIS, SIGMOID]
    for k in kernels:
        pylab.figure()
        pylab.title(kernel_name(k))
        pylab.plot([p[0] for p in classA],
                    [p[1] for p in classA],
                    'bo')
        pylab.plot([p[0] for p in classB],
                    [p[1] for p in classB],
                    'ro')

        # Add the sets together and shuffle them around.
        datapoints = gen_datapoints_from_classes(classA, classB)

        t= train(datapoints, k)

        # Plot the decision boundaries.
        xr=numpy.arange(-4, 4, 0.05)
        yr=numpy.arange(-4, 4, 0.05)
        grid=matrix([[indicator(t, x, y, k) for y in yr] for x in xr])
        pylab.contour(xr, yr, grid, (-1.0, 0.0, 1.0), colors=('red', 'black', 'blue'), linewidths=(1, 3, 1))

    # Now that we are done we show the ploted datapoints and the decision
    # boundary.
    pylab.show()

def main():
    import pylab
    
    # Create the gaussian data
    Xin, Yin = pylab.mgrid[0:201, 0:201]
    data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)
    
    pylab.matshow(data, cmap='gist_rainbow')
    
    params = fitgaussian(data)
    fit = gaussian(*params)
    
    pylab.contour(fit(*pylab.indices(data.shape)), cmap='copper')
    ax = pylab.gca()
    (height, x, y, width_x, width_y) = params
    
    pylab.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=ax.transAxes)
    
    pylab.show()

def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, max_iters=300, optimize=True, plot=True):
    """
    Show an example of a multimodal error surface for Gaussian process
    regression. Gene 939 has bimodal behaviour where the noisy mode is
    higher.
    """

    # Contour over a range of length scales and signal/noise ratios.
    length_scales = np.linspace(0.1, 60., resolution)
    log_SNRs = np.linspace(-3., 4., resolution)

    try:import pods
    except ImportError:
        print 'pods unavailable, see https://github.com/sods/ods for example datasets'
        return
    data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number)
    # data['Y'] = data['Y'][0::2, :]
    # data['X'] = data['X'][0::2, :]

    data['Y'] = data['Y'] - np.mean(data['Y'])

    lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.RBF)
    if plot:
        pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet)
        ax = pb.gca()
        pb.xlabel('length scale')
        pb.ylabel('log_10 SNR')

        xlim = ax.get_xlim()
        ylim = ax.get_ylim()

    # Now run a few optimizations
    models = []
    optim_point_x = np.empty(2)
    optim_point_y = np.empty(2)
    np.random.seed(seed=seed)
    for i in range(0, model_restarts):
        # kern = GPy.kern.RBF(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.))
        kern = GPy.kern.RBF(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50))

        m = GPy.models.GPRegression(data['X'], data['Y'], kernel=kern)
        m.likelihood.variance = np.random.uniform(1e-3, 1)
        optim_point_x[0] = m.rbf.lengthscale
        optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);

        # optimize
        if optimize:
            m.optimize('scg', xtol=1e-6, ftol=1e-6, max_iters=max_iters)

        optim_point_x[1] = m.rbf.lengthscale
        optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);

        if plot:
            pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1] - optim_point_x[0], optim_point_y[1] - optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
        models.append(m)

    if plot:
        ax.set_xlim(xlim)
        ax.set_ylim(ylim)
    return m # (models, lls)

def plotRadTilt(plot_data, plot_cmaps, plot_titles, grid, file_name, base_ref=None):
    subplot_base = 220
    n_plots = 4

    xs, ys, gs_x, gs_y, map = grid

    pylab.figure(figsize=(12,12))
    pylab.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.02, wspace=0.04)

    for plot in range(n_plots):
        pylab.subplot(subplot_base + plot + 1)
        cmap, vmin, vmax = plot_cmaps[plot]
        cmap.set_under("#ffffff", alpha=0.0)

        pylab.pcolormesh(xs - gs_x / 2, ys - gs_y / 2, plot_data[plot] >= -90, vmin=0, vmax=1, cmap=mask_cmap)

        pylab.pcolormesh(xs - gs_x / 2, ys - gs_y / 2, plot_data[plot], vmin=vmin, vmax=vmax, cmap=cmap)
        pylab.colorbar()

        if base_ref is not None:
            pylab.contour(xs, ys, base_ref, levels=[10.], colors='k', lw=0.5)

#       if plot == 0:
#           pylab.contour(xs, ys, refl_88D[:, :, 0], levels=np.arange(10, 80, 10), colors='#808080', lw=0.5)

#       pylab.plot(radar_x, radar_y, 'ko')

        pylab.title(plot_titles[plot])

        drawPolitical(map)

    pylab.savefig(file_name)
    pylab.close()

    return

def contourFromFunction(XYfunction,plotPoints=100,\
			xrange=None,yrange=None,numContours=20,alpha=1.0, contourLines=None):
	"""
	Given a 2D function, plots constant contours over the given
	range.  If the range is not given, the current plotting
	window range is used.
	"""
	
	# set up x and y ranges
	currentAxis = pylab.axis()
	if xrange is not None:
		xvalues = pylab.linspace(xrange[0],xrange[1],plotPoints)
	else:
		xvalues = pylab.linspace(currentAxis[0],currentAxis[1],plotPoints)
	if yrange is not None:
		yvalues = pylab.linspace(yrange[0],yrange[1],plotPoints)
	else:
		yvalues = pylab.linspace(currentAxis[2],currentAxis[3],plotPoints)
	
	#coordArray = _coordinateArray2D(xvalues,yvalues)
	# add extra dimension to this to make iterable?
	# bug here!  need to fix for contour plots
	z = map( lambda y: map(lambda x: XYfunction(x,y), xvalues), yvalues)
	if contourLines:
		pylab.contour(xvalues,yvalues,z,contourLines,alpha=alpha)
	else:
		pylab.contour(xvalues,yvalues,z,numContours,alpha=alpha)

def polarmap_contours(polarmap):  # Example docstring
   """
   Identifies the real and imaginary contours in a polar map.  Returns
   the real and imaginary contours as 2D vertex arrays together with
   the pairs of contours known to intersect. The coordinate system is
   normalized so x and y coordinates range from zero to one.
   Contour plotting requires origin='upper' for consistency with image
   coordinate system.
   """
   # Convert to polar and normalise
   normalized_polar = normalize_polar_channel(polarmap)
   figure_handle = plt.figure()
   # Real component
   re_contours_plot = plt.contour(normalized_polar.real, levels=[0], origin='upper')
   re_path_collections = re_contours_plot.collections[0]
   re_contour_paths = re_path_collections.get_paths()
   # Imaginary component
   im_contours_plot = plt.contour(normalized_polar.imag, levels=[0], origin='upper')
   im_path_collections = im_contours_plot.collections[0]
   im_contour_paths = im_path_collections.get_paths()
   plt.close(figure_handle)
   intersections = [ (re_ind, im_ind)
                     for (re_ind, re_path) in enumerate(re_contour_paths)
                     for (im_ind, im_path) in enumerate(im_contour_paths)
                    if im_path.intersects_path(re_path)]

   (ydim, xdim) = polarmap.shape
   # Contour vertices  0.5 pixel inset. Eg. (0,0)-(48,48)=>(0.5, 0.5)-(47.5, 47.5)
   # Returned values will not therefore reach limits of 0.0 and 1.0
   re_contours =  [remove_path_duplicates(re_path.vertices) / [ydim, xdim] \
                       for re_path in re_contour_paths]
   im_contours =  [remove_path_duplicates(im_path.vertices) / [ydim, xdim] \
                       for im_path in im_contour_paths]
   return (re_contours,  im_contours, intersections)

def plot_decision_region(x1_min, x2_min, x1_max, x2_max, clf, n_points=1000):
    print "Plotting decision region"
    x1, x2 = numpy.meshgrid(numpy.linspace(x1_min, x1_max, n_points), numpy.linspace(x2_min, x2_max, n_points))
    z = clf.decision_function(numpy.c_[x1.ravel(), x2.ravel()]).reshape(x1.shape)

    pylab.contour(x1, x2, z, levels=[0.0], linestyles="solid", linewidths=2.0)
    pylab.contour(x1, x2, z, levels=[-1.0, 1.0], linestyles="dashed", linewidths=1.0)

def plot_model(gm, dim):
    X, Y, Z, V = gm.density_on_grid(dim = dim)
    h = gm.plot(dim = dim)
    [i.set_linestyle('-.') for i in h]
    P.contour(X, Y, Z, V)
    data = gm.sample(200)
    P.plot(data[:, dim[0]], data[:,dim[1]], '.')

def kappa_residual_grid_plot(env, model, base_model, obj_index, with_contours=False, only_contours=False, with_colorbar=True):
    obj0,data0 = base_model['obj,data'][obj_index]
    obj1,data1 = model['obj,data'][obj_index]

    kappa = data1['kappa'] - data0['kappa']
       
    grid = obj1.basis._to_grid(kappa, obj1.basis.subdivision)
    R = obj1.basis.mapextent

    kw = {'extent': [-R,R,-R,R],
          'interpolation': 'nearest',
          'aspect': 'equal',
          'origin': 'upper',
          'cmap': cm.gist_stern,
          'fignum': False}
          #'vmin': -1,
          #'vmax':  1}

    if not only_contours:
        pl.matshow(grid, **kw)
    if only_contours or with_contours:
        kw.update({'colors':'k', 'linewidths':1, 'cmap':None})
        pl.contour(grid, **kw)
        kw.update({'colors':'k', 'linewidths':2, 'cmap':None})
        pl.contour(grid, [0], **kw)

    if with_colorbar:
        glscolorbar()
    return

 def view_pcem(self,PCEM):
   nbins = np.sqrt(len(PCEM))
   levels=[0.0,0.01,0.02,0.03,0.04,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5]
   #levels2=[0.0,0.01,0.05,0.1,0.2,0.3,0.4,0.5,1.0]
 
   T1 = PCEM[:,0].reshape( (nbins,nbins))
   T2 = PCEM[:,1].reshape( (nbins,nbins))
   E  = PCEM[:,2].reshape( (nbins,nbins))
   EC  = PCEM[:,3].reshape( (nbins,nbins))
   M  = PCEM[:,4].reshape( (nbins,nbins))
   pp.figure()
   pp.subplot(1,2,1)
   pp.contour( T1, T2, E, levels,linewidths=0.5,colors="k",linestyles='solid'  )
   pp.contourf( T1, T2, E, levels, alpha=0.75  )
   pp.xlabel( "proposal")
   pp.ylabel( "current")
   pp.title( "Error" )
   pp.colorbar()
   # pp.subplot(1,3,2)
   # pp.contour( T1, T2, EC, levels,linewidths=0.5,colors="k",linestyles='solid'  )
   # pp.contourf( T1, T2, EC, levels, alpha=0.75  )
   # pp.xlabel( "proposal")
   # pp.ylabel( "current")
   # pp.colorbar()
   # pp.title( "Corrected Error" )
   pp.subplot(1,2,2)
   pp.contour( T1, T2, M, 20,linewidths=1.0,colors="k",linestyles='solid'  )
   pp.contourf( T1, T2, M, 20  )
   pp.xlabel( "proposal")
   pp.ylabel( "current")
   pp.colorbar()
   pp.title("MU Z")

def plot_function_contour(x, y, function, **kwargs):
    """Make a contour plot of a function of two variables.

    Parameters
    ----------
    x, y : array_like of float
        The positions of the nodes of a plotting grid.
    function : function
        The function to plot.
    filling : bool
        Fill contours if True (default).
    num_contours : int
        The number of contours to plot, 50 by default.
    xlabel, ylabel : str, optional
        The axes labels. Empty by default.
    title : str, optional
        The title. Empty by default.

    """
    X, Y = numpy.meshgrid(x, y)
    Z = []
    for y_value in y:
        Z.append([])
        for x_value in x:
            Z[-1].append(function(x_value, y_value))
    Z = numpy.array(Z)
    num_contours = kwargs.get('num_contours', 50)
    if kwargs.get('filling', True):
        pylab.contourf(X, Y, Z, num_contours, cmap=pylab.cm.jet)
    else:
        pylab.contour(X, Y, Z, num_contours, cmap=pylab.cm.jet)
    pylab.xlabel(kwargs.get('xlabel', ''))
    pylab.ylabel(kwargs.get('ylabel', ''))
    pylab.title(kwargs.get('title', ''))