def plot_kde(opts, kde_data, centre = None):

    plt.figure()

    plt.contourf(kde_data[0], kde_data[1],
                 kde_data[2], 30)

    if centre:
        plt.plot(np.ones(2) * centre[0], (plt.ylim()), 'w:',
                label = 'KDE Centre')
        plt.plot((plt.xlim()), np.ones(2) * centre[1], 'w:')
    
    if opts.H0 == 100.0:
        plt.xlabel(r'X [Mpc h$^{-1}$]')
        plt.ylabel(r'Y [Mpc h$^{-1}$]')
    else:
        plt.xlabel('X [Mpc]')
        plt.ylabel('Y [Mpc]')
        
    plt.legend()
    plt.colorbar()

    output_file = opts.input_file + '.kde.pdf'
    plt.savefig(output_file)
    print ' KDE plot saved to:', output_file  
    
    plt.close()

 def contourf(self,vv=range(-10,0),**kwargs):
     fig= plt.figure(figsize=(9,8))
     #h.imshow(np.flipud(self.Z),extent=[bbox[0],bbox[1],bbox[3],bbox[2]])
     plt.contourf(self.grd.X,self.grd.Y,self.Z,vv,**kwargs)
     plt.colorbar()
     plt.axis('equal')
     return fig

def draw(data, classes, model, resolution=100):
    mycm = mpl.cm.get_cmap('Paired')
    
    one_min, one_max = data[:, 0].min()-0.1, data[:, 0].max()+0.1
    two_min, two_max = data[:, 1].min()-0.1, data[:, 1].max()+0.1
    xx1, xx2 = np.meshgrid(np.arange(one_min, one_max, (one_max-one_min)/resolution),
                     np.arange(two_min, two_max, (two_max-two_min)/resolution))
    
    inputs = np.c_[xx1.ravel(), xx2.ravel()]
    z = []
    for i in range(len(inputs)):
        z.append(predict(model, inputs[i])[0])
    result = np.array(z).reshape(xx1.shape)
    
    plt.contourf(xx1, xx2, result, cmap=mycm)
    plt.scatter(data[:, 0], data[:, 1], s=50, c=classes, cmap=mycm)
    
    t = np.zeros(15)
    for i in range(15):
        if i < 5:
            t[i] = 0
        elif i < 10:
            t[i] = 1
        else:
            t[i] = 2
    plt.scatter(model[:, 0], model[:, 1], s=150, c=t, cmap=mycm)
    
    plt.xlim([0, 10])
    plt.ylim([0, 10])
    
    plt.show()

def plot_qdens_log(r,z,dens, dens0):
    CS0 = plt.contourf(z,r,dens0, colors='k') #mark areas with data but zero density
    if dens.max() > 0.0:
        CS1 = plt.contourf(z, r, dens,norm=LogNorm(vmin=1e10, vmax=1e20))
        if drawCB:
            CB1 = plt.colorbar(CS1)
    else:
        CS1 = None
        if drawCB:
            CB1 = plt.colorbar(CS0, ticks=[])
    
    if dens.min() < 0.0:
        CS2 = plt.contourf(z, r, -dens,norm=LogNorm(vmin=1e10, vmax=1e20), hatches='x')
        if drawCB:
            CB2 = plt.colorbar(CS2)
    else:
        CS2 = None
        if drawCB:
            CB2 = plt.colorbar(CS0, ticks=[])
    
    if drawCB:
        CB1.set_label("Positive charge density [cm$^{-3}$]")
        CB2.set_label("Negative charge density [cm$^{-3}$]")
    plt.xlabel("z [um]")
    plt.ylabel("r [um]")
    if axisShape == "image":
        plt.axis('image')
    if cutR != None:
        plt.ylim(0,cutR)
    return (CS0, CS1, CS2)

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?')

def plot_decision_regions(X, y, classifier, resolution=0.02):

    # マーカーとカラーマップの準備
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # 決定領域のプロット　
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    # グリッドポイントの生成
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
      np.arange(x2_min, x2_max, resolution))
    # 各特徴量を1次元配列に変換して予測を実行
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    # 予測結果を元のグリッドポイントのデータサイズに変換
    Z = Z.reshape(xx1.shape)
    # グリッドポイントの等高線プロット
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    # 軸の範囲設定
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # クラスごとにサンプルをプロット
    for idx, cl in enumerate(np.unique(y)):
      plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1], alpha=0.8, c=cmap(idx),
        marker=markers[idx], label=cl)

def plot_3d(M, epoints, eloss, limits = [], xlim = [], ylim = [], 
            title = '', save_fig = False, outname = '3dmap.png'):
    """
    Plots 3d incident energy x energy loss graph. M is a matrix with the data
    as M[incident energy, energy loss]. epoints is the incident energies, and
    eloss the energy loss.
    """

    plt.figure()
    if len(limits[:]) == 0:
        fig = plt.contourf(epoints, eloss, M, 100)
        plt.colorbar()
    else:
        tks = []
        for i in range (6):
            tks.append(limits[0] + i*(limits[1]-limits[0])/5)
            
        z = np.linspace(limits[0],limits[1], 100, endpoint=True)
        fig = plt.contourf(epoints, eloss, M, z)
        plt.colorbar(ticks=tks)

    plt.xlabel('Incident Energy (eV)', fontsize=15)
    plt.ylabel('Energy Loss (eV)', fontsize=15)
    if len(xlim) != 0:    
        plt.xlim(xlim[0],xlim[1])
    if len(ylim) != 0:
        plt.ylim(ylim[0], ylim[1])
    plt.title(title)
    if save_fig is True:
        plt.savefig(outname, format='png', dpi=1000)

    def display(self, X, y, w):
        # create a mesh to plot in
        h = .02
        x_min, x_max = X[:, 1].min() - 1, X[:, 1].max() + 1
        y_min, y_max = X[:, 2].min() - 1, X[:, 2].max() + 1
        xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
                             np.arange(y_min, y_max, h))

        Z = np.inner(w, np.c_[np.ones(xx.ravel().shape),
                              xx.ravel(),
                              yy.ravel()])
        Z[Z >= 0] = 1
        Z[Z < 0] = -1
        Z = Z.reshape(xx.shape)
        # print 'Z', Z
        plt.figure()
        plt.contourf(xx, yy, Z, cmap=plt.cm.Paired, alpha=0.8)
        # Plot also the training points
        plt.scatter(X[:, 1], X[:, 2], c=y, 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 plot_landscape_2d(landscape, ds):
    """
    Plot the height of the landscape in 2 dimensions given a landscape object
    :param landscape: Landscape object with all data
    :param ds: Downsampling factor for only plotting every ds point
    :return: Plot the landscape in the x-y-coordinate system
    """

    # Construct the (x, y)-coordinate system
    x_grid = np.linspace(landscape.x_min, landscape.x_max, landscape.num_of_nodes_x)
    y_grid = np.linspace(landscape.y_max, landscape.y_min, landscape.num_of_nodes_y)
    x, y = np.meshgrid(x_grid[0::ds], y_grid[0::ds])
    z = landscape.arr[0::ds, 0::ds]

    # Decide color map and number of contour levels
    cmap = plt.get_cmap('terrain')
    v = np.linspace(min(landscape.coordinates[:, 2]), max(landscape.coordinates[:, 2]), 100, endpoint=True)
    plt.contourf(x, y, z, v, cmap=cmap)
    plt.colorbar(label='Height', spacing='uniform')

    # Title and labels
    plt.rcParams.update({'font.size': 14})
    plt.title('The landscape')
    plt.xlabel('x')
    plt.ylabel('y')

    plt.show()

def plot_ICON_clt(ICON_data_dict):
	"This function gets ICON data and plots corresponding satellite pictures."

	# Import and create basemap for plotting countries and coastlines
	from mpl_toolkits.basemap import Basemap

	# Create well formatted time_string
	time_string = datetime.fromtimestamp(int(unix_time_in)).strftime('%Y-%m-%d-%H-%M')

	# Plotting temperature data
	plt.figure()
	# Plot contourf plot with lat/lon regridded ICON data
	cmap_cloud   = plt.cm.gray
	levels_cloud = np.arange(0,101,10)
	plt.contourf(ICON_data_dict["ICON_X_mesh"], ICON_data_dict["ICON_Y_mesh"], ICON_data_dict["ICON_clt"], levels=levels_cloud, cmap=cmap_cloud)
	plt.colorbar()
	# Plot map data
	map = Basemap(llcrnrlon=4.0,llcrnrlat=47.0,urcrnrlon=15.0,urcrnrlat=55.0,
	             resolution='i')
	map.drawcoastlines()
	map.drawcountries()
	lat_ticks = [55.0,53.0,51.0,49.0,47.0]
	map.drawparallels(lat_ticks, labels=[1,0,0,0], linewidth=0.0)
	lon_ticks = [4.0,6.0,8.0,10.0,12.0,14.0]
	map.drawmeridians(lon_ticks, labels=[0,0,0,1], linewidth=0.0)
	# Save plot and show it
	plt.savefig(output_path + 'TotalCloudCover_' + time_string + '.png')

def plot2DInRealCrd( u, transverseProfile, stat ):
    R = stat["outerRadius"]
    v = np.linspace(0.0, 5E5, 1001)
    U, V = np.meshgrid(u,v)
    A = damping(U, V, stat["wavenumber"], stat["width"], stat["solver"]["eigenvalue"], R)
    energy = np.zeros(U.shape)
    for i in range(0, U.shape[1]):
        if (( U[0,i] > 0.0 ) or ( U[0,i] < -stat["width"])):
            energy[:,i] = transverseProfile[i]*A[:,i]
        else:
            energy[:,i] = transverseProfile[i]

    plt.contourf(U,V/1000.0, energy, 100, cmap="gist_heat", norm=mpl.colors.LogNorm())
    plt.xlabel("$u$ (nm)")
    plt.ylabel("$v \; (\mathrm{\mu m})$")
    fname = "Figures/profile2D_uvplane.jpeg"
    plt.savefig(fname, bbox_inches="tight", dpi=800)
    print ("Figure written to %s"%(fname))

    plt.clf()

    XYcompl = R*np.exp((U+1j*V)/R)
    transverse = XYcompl.real
    longitudinal = XYcompl.imag

    plt.contourf(longitudinal/1000.0,transverse,energy, 100, cmap="gist_heat", norm=mpl.colors.LogNorm())
    plt.gca().set_axis_bgcolor("#3E0000")
    plt.xlabel("$z \; (\mathrm{\mu m}$)")
    plt.ylabel("$x$ (nm)")
    fname = "Figures/profile2D_xyplane.jpeg"
    plt.savefig(fname, bbox_inches="tight", dpi=800)
    print ("Figure written to %s"%(fname))
    return longitudinal, transverse, energy

def plot_vert_vgradrho_rho_diff(show=True, save=False):
    plt.close('all')
    plt.figure()

    rho1 = np.array(g1a['Rho'])
    vgradrho1 = \
            g1a['V!Dn!N (up)']\
            *cr.calc_rusanov_alts_ausm(g1a['Altitude'],rho1)/g1a['Rho']

    rho2 = np.array(g2a['Rho'])
    vgradrho2 = \
            g2a['V!Dn!N (up)']\
            *cr.calc_rusanov_alts_ausm(g2a['Altitude'],rho2)/g2a['Rho']

    ilt = np.argmin(np.abs(g1a['LT'][2:-2, 0, 0]-whichlt))+2
    vgradrho1 = vgradrho1[ilt,2:-2,2:-2]
    vgradrho2 = vgradrho2[ilt,2:-2,2:-2]
    dvgradrho = np.array(vgradrho1-vgradrho2).T
    lat = np.array(g1a['dLat'][ilt, 2:-2, 2:-2]).T
    alt = np.array(g1a['Altitude'][ilt, 2:-2, 2:-2]/1000).T
    plt.contourf(lat, alt, dvgradrho, levels=np.linspace(-1, 1, 21)*1e-4,
                 cmap='seismic', extend='both')
    plt.xlim(-90, -40)
    plt.xlabel('Latitude')
    plt.ylabel('Altitude')
    plt.text(0.5, 0.95, 'Time: '+tstring, fontsize=15,
            horizontalalignment='center', transform=plt.gcf().transFigure)
    if show:
        plt.show()
    if save:
        plt.savefig(spath+'06_vert_vgradrho_rho_diff_%s.pdf' % tstring)
    return

def plot_vert_divv_diff(show=True, save=False):
    plt.close('all')
    plt.figure()

    velr = np.array(g1a['V!Dn!N (up)'])
    divv1 = cr.calc_div_vert(g1a['Altitude'], velr)

    velr = np.array(g2a['V!Dn!N (up)'])
    divv2 = cr.calc_div_vert(g2a['Altitude'], velr)

    ilt = np.argmin(np.abs(g1a['LT'][2:-2, 0, 0]-whichlt))+2
    divv1 = divv1[ilt,2:-2,2:-2]
    divv2 = divv2[ilt,2:-2,2:-2]
    ddivv = np.array(divv1-divv2).T
    lat = np.array(g1a['dLat'][ilt, 2:-2, 2:-2]).T
    alt = np.array(g1a['Altitude'][ilt, 2:-2, 2:-2]/1000).T
    plt.contourf(lat, alt, ddivv, levels=np.linspace(-1, 1, 21)*1e-4,
                 cmap='seismic', extend='both')
    plt.xlim(-90, -40)
    plt.xlabel('Latitude')
    plt.ylabel('Altitude')
    plt.text(0.5, 0.95, 'Time: '+tstring, fontsize=15,
            horizontalalignment='center', transform=plt.gcf().transFigure)
    if show:
        plt.show()
    if save:
        plt.savefig(spath+'04_vert_divv_diff%s.pdf' % tstring)
    return

 def plot(self):
     bounds = self.bounds
     x1 = np.linspace(bounds[0][0], bounds[0][1], 100)
     x2 = np.linspace(bounds[1][0], bounds[1][1], 100)
     X1, X2 = np.meshgrid(x1, x2)
     X = np.hstack((X1.reshape(100*100,1),X2.reshape(100*100,1)))
     Y = self.f(X)
     
     fig = plt.figure()
     ax = fig.gca(projection='3d')
     ax.plot_surface(X1, X2, Y.reshape((100,100)), rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
     ax.zaxis.set_major_locator(LinearLocator(10))
     ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
     ax.set_title(self.name)    
         
     plt.figure()
     size_color = np.linspace(-1.5, 6.0, 100, endpoint=True)
     plt.contourf(X1, X2, Y.reshape((100,100)),size_color)
     if (len(self.min)>1):    
         plt.plot(np.array(self.min)[:,0], np.array(self.min)[:,1], 'w.', markersize=20, label=u'Observations')
     else:
         plt.plot(self.min[0][0], self.min[0][1], 'w.', markersize=20, label=u'Observations')
     plt.colorbar()
     plt.xlabel('X1')
     plt.ylabel('X2')
     plt.title(self.name)
     # plt.show()
     savefig("object_true_2d.pdf")

def plot_decision_regions(
    X, y, classifier, xlab="X", ylab="y", legend_loc="upper left", test_idx=None, resolution=0.02
):
    # setup marker generator and color map
    markers = ("s", "x", "o", "^", "v")
    colors = ("red", "blue", "lightgreen", "gray", "cyan")
    cmap = ListedColormap(colors[: len(np.unique(y))])

    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution), np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    X_test, y_test = X[test_idx, :], y[test_idx]
    # plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1], alpha=0.8, c=cmap(idx), marker=markers[idx], label=cl)

    if test_idx:
        X_test, y_test = X[test_idx, :], y[test_idx]
        plt.scatter(X_test[:, 0], X_test[:, 1], c="", alpha=1.0, linewidth=1, marker="o", s=55, label="test set")

    plt.xlabel(xlab)
    plt.ylabel(ylab)
    plt.legend(loc=legend_loc)

def example1():
    """标准栅格情况下的画法
    当 x = [1,2,3,4], y = [1,2,3,4], 也就是一共是16个点，x和y相互对应的时候。这就叫标准栅格。
    """
    dx, dy = 0.01, 0.01
    x = np.arange(0.0, 10.0, 0.05)
    y = np.arange(0.0, 5.0, 0.05)
    x, y = np.meshgrid(x, y)
    z = np.sin(x) ** 10 + np.cos(10 + y * x) * np.cos(x)
    
    levels = MaxNLocator(nbins=15).tick_values(z.min(), z.max())
    cmap = plt.get_cmap('PiYG')
    norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
    
    plt.subplot(2, 1, 1)
    im = plt.pcolormesh(x, y, z, cmap=cmap, norm=norm)
    plt.colorbar()
    plt.axis([x.min(), x.max(), y.min(), y.max()])
    plt.title('pcolormesh with levels')
    
    plt.subplot(2, 1, 2)
    plt.contourf(x + dx / 2.,
                 y + dy / 2., z, levels=levels,
                 cmap=cmap)
    plt.colorbar()
    plt.title('contourf with levels')

    plt.show()

    def graph(self, x_axis, y_axis, image, params, result = None):
        #plot the sample data
        from matplotlib import pyplot
        pyplot.figure()
        pyplot.subplot(2,1,1)
        pyplot.contourf(x_axis, y_axis, image, alpha = 0.5)
        #plot the fit
        #sample the fit with more precision

        x_axis_fit = np.linspace(x_axis.min(), x_axis.max(), x_axis.size * 5)
        y_axis_fit = np.linspace(y_axis.min(), y_axis.max(), y_axis.size * 5)
        xx, yy = np.meshgrid(x_axis_fit, y_axis_fit)
        params['background_level'].value = 0                                                                                
        fit_1d = self.ion_model(params, xx, yy, use_1d=True)
        xx, yy = np.meshgrid(x_axis, y_axis)
        pyplot.subplot(2,1,2)
        ximage = image.sum(axis=0)
        ximage = ximage - ximage[0]
        a1d = np.ones_like(self.allions_1d)
        #low_values_indices = self.allions_1d - bg < (0.5)*params['amplitude'].value  # Where values are low
        low_values_indices = self.allions_1d < self.model_threshold  # Where values are low
        a1d[low_values_indices] = 0  # All low values set to 0
        pyplot.plot(x_axis, a1d * ximage.max())
        #pyplot.plot(x_axis, ximage * a1d, 'o', markersize=10)
        pyplot.plot(x_axis, ximage)
        pyplot.plot(x_axis_fit, fit_1d)
        pyplot.tight_layout()
        pyplot.show()

def example2():
    """散点栅格情况下的画法
    x, y = (1,7), (4,9), (11,2),... 这种没有规律的散点x,y坐标，我们称其为散点栅格
    """
    dx, dy = 0.05, 0.05
    x = np.array([[1, 2],
                  [1, 3]])
    y = np.array([[3, 3],
                  [1, 1]])
    z = np.array([[5, 10],[15, 20]])
                  
    levels = MaxNLocator(nbins=4).tick_values(z.min(), z.max())
    cmap = plt.get_cmap('PiYG')
    norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
    
    plt.subplot(2, 1, 1)
    im = plt.pcolormesh(x, y, z, cmap=cmap, norm=norm)
    plt.colorbar()
    plt.axis([x.min(), x.max(), y.min(), y.max()])
    plt.title('pcolormesh with levels')

    plt.subplot(2, 1, 2)
    plt.contourf(x + dx / 2.,
                 y + dy / 2., z, levels=levels,
                 cmap=cmap)
    plt.colorbar()
    plt.title('contourf with levels')

    plt.show()

def gridplot(z, npoints=1024, phi=None, theta=None, cbar=True, earthpt=None, vmin=None, vmax=None, aspect='auto', rasterized=False):
    importmpl()
    if phi == None:
        import gbm
        (phi, theta) = gbm.respgrid()
    p = np.linspace(0, 2*np.pi, npoints)
    t = np.linspace(0, np.pi, int(npoints/2))
    idx = np.nonzero(phi == 0)
    pwrap = phi[idx] + 2 * np.pi
    twrap = theta[idx]
    zwrap = z[idx]
    (pp, tt) = np.meshgrid(p, t)
    zgrid = mlab.griddata(np.hstack((phi, pwrap)), np.hstack((theta, twrap)), np.hstack((z, zwrap)), pp, tt, interp='linear')
    if earthpt != None:
        occlimit = np.cos(67 * np.pi/180.)
        eagrid = np.sum(pt2xyz((pp.ravel(), tt.ravel())).T * pt2xyz(earthpt), axis=1).reshape(pp.shape)
        zgrid[eagrid > occlimit] = 0
    h = plt.imshow(zgrid, extent=(0, 2*np.pi, np.pi, 0), vmin=vmin, vmax=vmax, aspect=aspect, rasterized=rasterized)
    if cbar:
        plt.colorbar(orientation='horizontal', aspect=30, shrink=1.0)
    if earthpt != None:
        plt.contourf(pp, tt, eagrid, [occlimit, 1.0], colors='0.75', label='Earth extent')
        plt.plot(earthpt[0], earthpt[1], marker=r'$\oplus$', ls='none', label='Earth')
    # plt.subplots_adjust(top=0.875)
    return (h, plt.gca())

def draw_graph(X, y, svc_linear, svc_poly2, svc_poly3, svc_rbf, h = 0.2):
    draw_init()

    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))

    titles = ['SVC with linear kernel',
              'SVC with polynomial degree 2 kernel',
              'SVC with polynomial degree 3 kernel',
              'SVC with RBF kernel']

    colors = np.array([x for x in 'bgrcmykbgrcmykbgrcmykbgrcmyk'])
    colors = np.hstack([colors] * 20)

    for i, clf in enumerate((svc_linear, svc_poly2, svc_poly3, svc_rbf)):
        plt.subplot(2, 2, i + 1)
        plt.subplots_adjust(wspace=0.4, hspace=0.4)

        Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

        Z = Z.reshape(xx.shape)
        plt.contourf(xx, yy, Z, cmap=plt.cm.Paired, alpha=0.8)

        plt.scatter(X[:, 0], X[:, 1], color=colors[y.tolist()].tolist())
        plt.xlim(xx.min(), xx.max())
        plt.ylim(yy.min(), yy.max())
        plt.xticks(())
        plt.yticks(())
        plt.title(titles[i])

    print(show_graph())