def plot_2d(x, y, mean, variance, ei, slice_at, v1_name, v2_name):
    h_fig = pplt.figure(figsize=(20, 8), dpi=100)
    pplt.subplot(131)
    h_mean = pplt.pcolormesh(x, y,
                             mean.reshape(x.shape[0], y.shape[0]))
    pplt.colorbar(h_mean)
    slice_at_list = np.squeeze(np.asarray(slice_at)).tolist()
    slice_at_string = str(["%.2f" % member for member in slice_at_list])
    pplt.xlabel(r'$' + v1_name + '$')
    pplt.ylabel(r'$' + v2_name + '$')
    pplt.title(r'Mean, slice along $( ' + v1_name + ',' + v2_name + ')$ at ' +
               slice_at_string)

    pplt.subplot(132)
    h_var = pplt.pcolormesh(x, y, 2*np.sqrt(variance.reshape(x.shape[0],
                                                   y.shape[0])))
    pplt.colorbar(h_var)
    pplt.xlabel(r'$' + v1_name + '$')
    pplt.ylabel(r'$' + v2_name + '$')
    pplt.title(r'2*Stdev, slice along $( ' + v1_name + ',' + v2_name + ')$' )

    pplt.subplot(133)
    h_ei = pplt.pcolormesh(x, y, ei.reshape(x.shape[0], y.shape[0]))
    pplt.colorbar(h_var)
    pplt.xlabel(r'$' + v1_name + '$')
    pplt.ylabel(r'$' + v2_name + '$')
    pplt.title(r'EI, slice along $( ' + v1_name + ',' + v2_name + ')$')

    pplt.draw()
    return (h_fig, h_mean, h_var, h_ei)

def Pcolor(xs, ys, zs, pcolor=True, contour=False, **options):
    """Makes a pseudocolor plot.
    
    xs:
    ys:
    zs:
    pcolor: boolean, whether to make a pseudocolor plot
    contour: boolean, whether to make a contour plot
    options: keyword args passed to pyplot.pcolor and/or pyplot.contour
    """
    Underride(options, linewidth=3, cmap=matplotlib.cm.Blues)

    X, Y = np.meshgrid(xs, ys)
    Z = zs

    x_formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)
    axes = pyplot.gca()
    axes.xaxis.set_major_formatter(x_formatter)

    if pcolor:
        pyplot.pcolormesh(X, Y, Z, **options)

    if contour:
        cs = pyplot.contour(X, Y, Z, **options)
        pyplot.clabel(cs, inline=1, fontsize=10)

def plt_data():
    t = [[0,1], [1,0], [1, 1], [0, 0]]
    t2 = [1, 1, 1, 0]
    X = np.array(t)
    Y = np.array(t2)

    h = .02  # step size in the mesh

    logreg = linear_model.LogisticRegression(C=1e5)

    # we create an instance of Neighbours Classifier and fit the data.
    logreg.fit(X, Y)
    # Plot the decision boundary. For that, we will assign a color to each
    # point in the mesh [x_min, m_max]x[y_min, y_max].
    x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
    y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
    Z = logreg.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    plt.figure(1, figsize=(4, 3))
    plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Paired)

    # Plot also the training points
    plt.scatter(X[:, 0], X[:, 1], c=Y, edgecolors='k', cmap=plt.cm.Paired)
    plt.xlabel('Sepal length')
    plt.ylabel('Sepal width')

    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())
    plt.xticks(())
    plt.yticks(())

    plt.show()

    def test_unimodality_of_GEV(self):
        x0 = 1500

        mu = 1000
        data = np.array([x0])

        ksi = np.arange(-2, 2, 0.01)
        sigma = np.arange(10, 8000, 10)

        n_ksi = len(ksi)
        n_sigma = len(sigma)

        z = np.zeros((n_ksi, n_sigma))

        for i, the_ksi in enumerate(ksi):
            for j, the_sigma in enumerate(sigma):
                z[i, j] = gevfit.objective_function_stationary_high([the_sigma, mu, the_ksi], data)


        sigma, ksi = np.meshgrid(sigma, ksi)
        z = np.ma.masked_where(z == gevfit.BIG_NUM, z)
        z = np.ma.masked_where(z > 9, z)

        plt.figure()
        plt.pcolormesh(ksi, sigma, z)
        plt.colorbar()
        plt.xlabel('$\\xi$')
        plt.ylabel('$\\sigma$')
        plt.title('$\\mu = %.1f, x = %.1f$' % (mu, x0))

        plt.show()


        pass

def prettyPicture(clf, X_test, y_test):
    x_min = 0.0;
    x_max = 1.0
    y_min = 0.0;
    y_max = 1.0

    # Plot the decision boundary. For that, we will assign a color to each
    # point in the mesh [x_min, m_max]x[y_min, y_max].
    h = .01  # step size in the mesh
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())

    plt.pcolormesh(xx, yy, Z)

    # Plot also the test points
    grade_sig = [X_test[ii][0] for ii in range(0, len(X_test)) if y_test[ii] == 0]
    bumpy_sig = [X_test[ii][1] for ii in range(0, len(X_test)) if y_test[ii] == 0]
    grade_bkg = [X_test[ii][0] for ii in range(0, len(X_test)) if y_test[ii] == 1]
    bumpy_bkg = [X_test[ii][1] for ii in range(0, len(X_test)) if y_test[ii] == 1]

    plt.scatter(grade_sig, bumpy_sig, color="b", label="fast")
    plt.scatter(grade_bkg, bumpy_bkg, color="r", label="slow")
    plt.legend()
    plt.xlabel("bumpiness")
    plt.ylabel("grade")

    plt.savefig("test.png")

    def visualize(self, output_file, width=2, show_charts=False):
        X = self.X

        # Create a grid of points
        x_min, x_max = min(X[:, 0] - width), max(X[:, 0] + width)
        y_min, y_max = min(X[:, 1] - width), max(X[:, 1] + width)
        xx,yy = np.meshgrid(np.arange(x_min, x_max, .05), np.arange(y_min,
            y_max, .05))

        # Flatten the grid so the values match spec for self.predict
        xx_flat = xx.flatten()
        yy_flat = yy.flatten()
        X_topredict = np.vstack((xx_flat,yy_flat)).T

        # Get the class predictions
        Y_hat = self.predict(X_topredict)
        Y_hat = Y_hat.reshape((xx.shape[0], xx.shape[1]))
        
        cMap = c.ListedColormap(['r','b','g'])

        # Visualize them.
        plt.figure()
        plt.pcolormesh(xx,yy,Y_hat, cmap=cMap)
        plt.scatter(X[:, 0], X[:, 1], c=self.C, cmap=cMap)
        plt.savefig(output_file)
        if show_charts:
            plt.show()

def plotmaptime():
    pcolormesh(yoko, time*1e6, absolute(Magcom))
    title("Reflection vs flux \n and time (1 us pulse) at 4.46 GHz")
    xlabel("Flux (V)")
    ylabel("Time (us)")
    #ylim(0, 1.5)
    colorbar()

def plot_spectrogram(raw_data, nfft, fs, channel_bottom, print_frequency_graph):
    data_shape = raw_data.shape

    print("Generating spectrogram...")
    plt_num = 1
    plt.clf()
    plt.figure(1)

    channel_data = []
    for i in range(0, data_shape[1]):
        plt.subplot(8, 2, plt_num)

        f, t, Sxx = signal.spectrogram(x=raw_data[:, i], nfft=nfft, fs=fs, noverlap=127, nperseg=128,
                                       scaling='density')  # returns PSD power per Hz
        plt.pcolormesh(t, f, Sxx)

        plt.xlabel('Time (sec)')
        plt.ylabel('Frequency (Hz)')
        plt.title('Channel %s' % (i + channel_bottom))
        plt_num += 1
        channel_data.append([f, t, Sxx])
        print("\tChannel %d spectrogram generated" % i)
    if print_frequency_graph:
        plt.show()
    return channel_data

def test_pcolormesh_global_with_wrap3():
    nx, ny = 33, 17
    xbnds = np.linspace(-1.875, 358.125, nx, endpoint=True)
    ybnds = np.linspace(91.25, -91.25, ny, endpoint=True)
    xbnds, ybnds = np.meshgrid(xbnds, ybnds)

    data = np.exp(np.sin(np.deg2rad(xbnds)) + np.cos(np.deg2rad(ybnds)))

    # this step is not necessary, but makes the plot even harder to do (i.e.
    # it really puts cartopy through its paces)
    ybnds = np.append(ybnds, ybnds[:, 1:2], axis=1)
    xbnds = np.append(xbnds, xbnds[:, 1:2] + 360, axis=1)
    data = np.ma.concatenate([data, data[:, 0:1]], axis=1)

    data = data[:-1, :-1]
    data = np.ma.masked_greater(data, 2.6)

    ax = plt.subplot(211, projection=ccrs.PlateCarree(-45))
    c = plt.pcolormesh(xbnds, ybnds, data, transform=ccrs.PlateCarree())
    assert c._wrapped_collection_fix is not None, \
        'No pcolormesh wrapping was done when it should have been.'

    ax.coastlines()
    ax.set_global()  # make sure everything is visible

    ax = plt.subplot(212, projection=ccrs.PlateCarree(-1.87499952))
    plt.pcolormesh(xbnds, ybnds, data, transform=ccrs.PlateCarree())
    ax.coastlines()
    ax.set_global()  # make sure everything is visible

def Picture(clf, X_test, y_test):
    x_min = 200.0
    x_max = 1000.0
    y_min = 600.0
    y_max = 2500.0

    # Plot the decision boundary. For that, we will assign a color to each
    # point in the mesh [x_min, m_max]x[y_min, y_max].
    h = 1  # step size in the mesh
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())

    plt.pcolormesh(xx, yy, Z, cmap=pl.cm.seismic)

    # Plot also the test points
    x1 = [X_test[ii][0] for ii in range(0, len(X_test)) if y_test[ii] == 0]
    y1 = [X_test[ii][1] for ii in range(0, len(X_test)) if y_test[ii] == 0]
    x2 = [X_test[ii][0] for ii in range(0, len(X_test)) if y_test[ii] == 1]
    y2 = [X_test[ii][1] for ii in range(0, len(X_test)) if y_test[ii] == 1]
    x3 = [X_test[ii][0] for ii in range(0, len(X_test)) if y_test[ii] == 2]
    y3 = [X_test[ii][1] for ii in range(0, len(X_test)) if y_test[ii] == 2]

    plt.scatter(x1, y1, color="b", label="class1")
    plt.scatter(x2, y2, color="r", label="class2")
    plt.scatter(x3, y3, color="g", label="class3")
    plt.legend()
    plt.xlabel("x")
    plt.ylabel("y")

    plt.savefig("testrf.png")

def train_Quasi_linear_SVM():
	from sklearn import svm
	clf = svm.SVC(kernel=get_KernelMatrix)
	X_train = np.r_[X1,X2]
	Y_train = np.r_[Y1,Y2]
	scatter(X[:,0],X[:,1],c='g')
	clf.fit(X_train, Y_train)

	y_pred = clf.predict(X_test)

	#scatter(X_test, y_pred)

	clf = svm.SVC(kernel=get_KernelMatrix)
	clf.fit(X, y)
	clf.predict(X_test1)

	# Plot the decision boundary. For that, we will assign a color to each
	# point in the mesh [x_min, m_max]x[y_min, y_max].
	h = 0.05
	x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
	y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
	xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
	Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

	# Put the result into a color plot
	Z = Z.reshape(xx.shape)
	plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Paired)

	# Plot also the training points
	plt.scatter(X_train[:, 0], X_train[:, 1], c=Y_train)
	plt.title('2-Class classification using Support Vector Machine with quasi-linear kernel')
	plt.axis('tight')
	plt.legend([Y_train[0], Y_train[-1]], ['negtive sample', 'postive sample'])
	plt.show()

def plotter(filename, xmin=0, xmax=300, Nx=2000,
            ymin=0, ymax=1, Ny=2000, sigma_x=3, sigma_y=0.01):
    root = '/home/cyneo/Work/Scans/Processed Data/Extracted CSV/'
    file1 = os.path.abspath(root + filename + '.csv')
    nx = linspace(xmin, xmax, Nx)
    ny = linspace(ymin, ymax, Ny)
    x, y = meshgrid(nx, ny)
    mastermesh = []
    with open(file1, 'r', encoding='utf8') as filein:
        file_reader = csv.reader(filein)
        next(file_reader)

        for word, frequency, inhubness, outhubness in file_reader:
            # want to feed the values into the center points
            if mastermesh == []:
                mastermesh = dgaussian(x, y, float(frequency),
                                       float(outhubness), sigma_x, sigma_y)
            else:
                mastermesh += dgaussian(x, y, float(frequency),
                                        float(outhubness), sigma_x, sigma_y)

    for x in range(len(mastermesh)):
        for y in range(len(mastermesh[x])):
            mastermesh[x, y] = np.log(mastermesh[x, y]+1)

    x, y = meshgrid(nx, ny)
    plt.pcolormesh(x, y, mastermesh)
    plt.show()
    outfile = os.path.abspath(root + filename + ' Array')
    np.save(outfile, mastermesh)

def plot_basins(f, Df, roots, xmin, xmax, ymin, ymax, numpoints=100, iters=15, colormap='brg'):
    '''Plot the basins of attraction of f.
    
    INPUTS:
    f       - A function handle. Should represent a function 
            from C to C.
    Df      - A function handle. Should be the derivative of f.
    roots   - An array of the zeros of f.
    xmin, xmax, ymin, ymax - Scalars that define the domain 
            for the plot.
    numpoints - A scalar that determines the resolution of 
            the plot. Defaults to 100.
    iters   - Number of times to iterate Newton's method. 
            Defaults to 15.
    colormap - A colormap to use in the plot. Defaults to 'brg'.    
    '''
    xreal = np.linspace(xmin, xmax, numpoints)
    ximag = np.linspace(ymin, ymax, numpoints)
    Xreal, Ximag = np.meshgrid(xreal, ximag)
    xold = Xreal+1j*Ximag
    n = 0
    while n <= iters:
        xnew = xold - f(xold)/Df(xold)
        xold = xnew
        n += 1 

    converged_to = np.empty_like(xnew)
    for i in xrange(xnew.shape[0]):
        for j in xrange(xnew.shape[1]):
            root = np.abs(roots-xnew[i,j]).argmin()
            converged_to[i,j] = root

    plt.pcolormesh(Xreal, Ximag, converged_to, cmap=colormap)

def pcolorRandom():
    "Makes a pcolormesh plot of randomly generated data pts."
    # make up some randomly distributed data
    npts = 100
    x = uniform(-3, 3, npts)
    y = uniform(-3, 3, npts)
    z = x * N.exp(-x ** 2 - y ** 2)

    # define grid.
    xi = N.arange(-3.1, 3.1, 0.05)
    yi = N.arange(-3.1, 3.1, 0.05)

    # grid the data.
    zi = griddata(x, y, z, xi, yi)

    # contour the gridded data, plotting dots at the randomly spaced data points.
    plt.pcolormesh(xi, yi, zi)	
    #CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
    #CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)
    plt.colorbar() # draw colorbar

    # plot data points.
    plt.scatter(x, y, marker='o', c='b', s=5)
    plt.xlim(-3, 3)
    plt.ylim(-3, 3)
    plt.title('griddata test (%d points)' % npts)
    plt.show()

def plot_eta(pT_lower_cut):

	properties_reco = [parse_file("/home/aashish/pythia_reco.dat", pT_lower_cut=pT_lower_cut), parse_file("/home/aashish/herwig_reco.dat", pT_lower_cut=pT_lower_cut), parse_file("/home/aashish/sherpa_reco.dat", pT_lower_cut=pT_lower_cut)]
	properties_truth = [parse_file("/home/aashish/pythia_truth.dat", pT_lower_cut=pT_lower_cut), parse_file("/home/aashish/herwig_truth.dat", pT_lower_cut=pT_lower_cut), parse_file("/home/aashish/sherpa_truth.dat", pT_lower_cut=pT_lower_cut)]
	labels = ["pythia", "herwig", "sherpa"]

	for prop_reco, prop_truth, label in zip(properties_reco, properties_truth, labels):
	
		x = prop_truth['hardest_eta']
		y = prop_reco['hardest_eta']

		H, xedges, yedges = np.histogram2d(x, y, bins=200, normed=1 )

		H = np.rot90(H)
		H = np.flipud(H)

		Hmasked = np.ma.masked_where(H == 0, H) # Mask pixels with a value of zero

		plt.pcolormesh(xedges,yedges, Hmasked)

		cbar = plt.colorbar()
		cbar.ax.set_ylabel('Counts')

		plt.xlim(0, 3)
		plt.ylim(0, 3)

		plt.xlabel('Truth $\eta$', fontsize=50, labelpad=75)
		plt.ylabel('Reco $\eta$', fontsize=50, labelpad=75)

		plt.gcf().set_size_inches(30, 30, forward=1)
		plt.gcf().set_snap(True)

		plt.savefig("plots/With MC/2D/eta_" + label + ".pdf")

		plt.clf()

    def test_tri_polar(self):
        # load data
        cubes = iris.load(tests.get_data_path(['NetCDF', 'ORCA2', 'votemper.nc']))
        cube = cubes[0]
        # The netCDF file has different data types for the points and
        # bounds of 'depth'. This wasn't previously supported, so we
        # emulate that old behaviour.
        cube.coord('depth').bounds = cube.coord('depth').bounds.astype(np.float32)

        # define a latitude trajectory (put coords in a different order to the cube, just to be awkward)
        latitudes = range(-90, 90, 2)
        longitudes = [-90]*len(latitudes)
        sample_points = [('longitude', longitudes), ('latitude', latitudes)]

        # extract
        sampled_cube = iris.analysis.trajectory.interpolate(cube, sample_points)
        self.assertCML(sampled_cube, ('trajectory', 'tri_polar_latitude_slice.cml'))

        # turn it upside down for the visualisation
        plot_cube = sampled_cube[0]
        plot_cube = plot_cube[::-1, :]

        plt.clf()
        plt.pcolormesh(plot_cube.data, vmin=cube.data.min(), vmax=cube.data.max())
        plt.colorbar()
        self.check_graphic()

        # Try to request linear interpolation.
        # Not allowed, as we have multi-dimensional coords.
        self.assertRaises(iris.exceptions.CoordinateMultiDimError, iris.analysis.trajectory.interpolate, cube, sample_points, method="linear")

        # Try to request unknown interpolation.
        self.assertRaises(ValueError, iris.analysis.trajectory.interpolate, cube, sample_points, method="linekar")

def write_potential(N=2.5, pphw=20, amplitude=1.0, sigmax=1e-1, sigmay=1e-1,
                    L=100., W=1.0, x_R0=0.05, y_R0=0.4, loop_type='Bell',
                    init_phase=0.0, shape='RAP', plot=True,
                    plot_dimensions=False, direction='right',
                    boundary_only=False, with_boundary=False, boundary_phase=0.0,
                    theta=0.0, smearing=False, verbose=True, linearized=False):

    p = Potential(N=N, pphw=pphw, amplitude=amplitude, sigmax=sigmax,
                  sigmay=sigmay, x_R0=x_R0, y_R0=y_R0, init_phase=init_phase,
                  shape=shape, L=L, W=W, loop_type=loop_type,
                  direction=direction, boundary_only=boundary_only,
                  with_boundary=with_boundary, theta=theta,
                  verbose=verbose, linearized=linearized)

    if not boundary_only:
        imag, imag_vector = p.imag, p.imag_vector
        real, real_vector = p.real, p.real_vector
    X, Y = p.X, p.Y

    if not boundary_only:
        if plot:
            import matplotlib.pyplot as plt
            if plot_dimensions:
                plt.figure(figsize=(L, W))
            plt.pcolormesh(X, Y, imag, cmap='RdBu_r')
            plt.savefig("imag.png")
            plt.pcolormesh(X, Y, real, cmap='RdBu_r')
            plt.savefig("real.png")

        np.savetxt("potential_imag.dat", zip(xrange(len(imag_vector)), imag_vector),
                   fmt=["%i", "%.12f"])
        np.savetxt("potential_real.dat", zip(xrange(len(real_vector)), real_vector),
                   fmt=["%i", "%.12f"])
        if shape != 'science':
            np.savez("potential_imag_xy.npz", X=X, Y=Y, P=imag_vector,
                     X_nodes=p.xnodes, Y_nodes=p.ynodes,
                     sigmax=sigmax, sigmay=sigmay)

    if shape == 'RAP':
        xi_lower, xi_upper = p.WG.get_boundary(theta=theta, smearing=smearing,
                                               boundary_phase=boundary_phase)
        # set last element to 0 (xi_lower) or W (xi_upper)
        print "WARNING: end of boundary not set zero!"
        # xi_lower[-1] = 0.0
        # xi_upper[-1] = W
        np.savetxt("upper.boundary", zip(xrange(p.nx), xi_upper))
        np.savetxt("lower.boundary", zip(xrange(p.nx), xi_lower))
        eps, delta = p.WG.get_cycle_parameters()
        np.savetxt("boundary.eps_delta", zip(eps, delta))
    if shape == 'RAP_TQD':
        eps_prime, delta_prime, theta_prime = p.WG.get_quantum_driving_parameters()
        xi_lower, xi_upper = p.WG.get_boundary(eps=eps_prime, delta=delta_prime,
                                               theta=theta_prime,
                                               smearing=smearing)
        # set last element to 0 (xi_lower) or W (xi_upper)
        xi_lower[-1] = 0.0
        xi_upper[-1] = W
        np.savetxt("upper.boundary", zip(xrange(p.nx), xi_upper))
        np.savetxt("lower.boundary", zip(xrange(p.nx), xi_lower))
        np.savetxt("boundary.eps_delta_theta", zip(eps_prime, delta_prime, theta_prime))

def iris_test(): # example code for LogisticRegression
    iris = datasets.load_iris()
    X = iris.data[:,:2] # first two features
    Y = iris.target
    h=.02 # step size in the mesh
    logreg = linear_model.LogisticRegression()
    logreg.fit(X,Y)

    x_min, x_max = X[:,0].min() - .5, X[:,0].max() + .5
    y_min, y_max = X[:,1].min() - .5, X[:,1].max() + .5
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
    Z = logreg.predict(np.c_[xx.ravel(), yy.ravel()])

    Z = Z.reshape(xx.shape)
    plt.figure(1, figsize=(4,3))
    plt.pcolormesh(xx,yy,Z,cmap=plt.cm.Paired)

    plt.scatter(X[:,0],X[:,1],c=Y,edgecolors='k',cmap=plt.cm.Paired)
    plt.xlabel('Sepal length')
    plt.ylabel('Sepal width')
    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())
    plt.xticks(())
    plt.yticks(())
    plt.show()

def train_plot(features, labels):
    y0, y1 = features[:, 2].min() * .9, features[:, 2].max() * 1.1
    x0, x1 = features[:, 0].min() * .9, features[:, 0].max() * 1.1
    X = np.linspace(x0, x1, 100)
    Y = np.linspace(y0, y1, 100)
    X, Y = np.meshgrid(X, Y)

    model = fit_model(1, features[:, (0, 2)], np.array(labels))
    C = predict(
        np.vstack([X.ravel(), Y.ravel()]).T, model).reshape(X.shape)
    if COLOUR_FIGURE:
        cmap = ListedColormap([(1., .6, .6), (.6, 1., .6), (.6, .6, 1.)])
    else:
        cmap = ListedColormap([(1., 1., 1.), (.2, .2, .2), (.6, .6, .6)])
    plt.xlim(x0, x1)
    plt.ylim(y0, y1)
    plt.xlabel(feature_names[0])
    plt.ylabel(feature_names[2])
    plt.pcolormesh(X, Y, C, cmap=cmap)
    if COLOUR_FIGURE:
        cmap = ListedColormap([(1., .0, .0), (.0, 1., .0), (.0, .0, 1.)])
        plt.scatter(features[:, 0], features[:, 2], c=labels, cmap=cmap)
    else:
        for lab, ma in zip(range(3), "Do^"):
            plt.plot(features[labels == lab, 0], features[
                     labels == lab, 2], ma, c=(1., 1., 1.))

def test_complete():
    fig = plt.figure('Figure with a label?', figsize=(10, 6))

    plt.suptitle('Can you fit any more in a figure?')

    # make some arbitrary data
    x, y = np.arange(8), np.arange(10)
    data = u = v = np.linspace(0, 10, 80).reshape(10, 8)
    v = np.sin(v * -0.6)

    plt.subplot(3, 3, 1)
    plt.plot(list(xrange(10)))

    plt.subplot(3, 3, 2)
    plt.contourf(data, hatches=['//', 'ooo'])
    plt.colorbar()

    plt.subplot(3, 3, 3)
    plt.pcolormesh(data)

    plt.subplot(3, 3, 4)
    plt.imshow(data)

    plt.subplot(3, 3, 5)
    plt.pcolor(data)

    plt.subplot(3, 3, 6)
    plt.streamplot(x, y, u, v)

    plt.subplot(3, 3, 7)
    plt.quiver(x, y, u, v)

    plt.subplot(3, 3, 8)
    plt.scatter(x, x**2, label='$x^2$')
    plt.legend(loc='upper left')

    plt.subplot(3, 3, 9)
    plt.errorbar(x, x * -0.5, xerr=0.2, yerr=0.4)

    ###### plotting is done, now test its pickle-ability #########

    # Uncomment to debug any unpicklable objects. This is slow (~200 seconds).
#    recursive_pickle(fig)

    result_fh = BytesIO()
    pickle.dump(fig, result_fh, pickle.HIGHEST_PROTOCOL)

    plt.close('all')

    # make doubly sure that there are no figures left
    assert_equal(plt._pylab_helpers.Gcf.figs, {})

    # wind back the fh and load in the figure
    result_fh.seek(0)
    fig = pickle.load(result_fh)

    # make sure there is now a figure manager
    assert_not_equal(plt._pylab_helpers.Gcf.figs, {})

    assert_equal(fig.get_label(), 'Figure with a label?')