def perform_interpolation(self,dataextrap,regridme,field_src,field_target,use_locstream):
		data = self.allocate()
		if self.geometry == 'surface':
			for kz in _np.arange(self.nz):
				field_src.data[:] = dataextrap[kz,:,:].transpose()
				field_target = regridme(field_src, field_target)
				if use_locstream:
					if self.nx == 1:
						data[kz,:,0] = field_target.data.copy()
					elif self.ny == 1:
						data[kz,0,:] = field_target.data.copy()
				else:
					data[kz,:,:] = field_target.data.transpose()[self.jmin:self.jmax+1, \
					                                             self.imin:self.imax+1]
					if self.debug and kz == 0:
						data_target_plt = _np.ma.masked_values(data[kz,:,:],self.xmsg)
						#data_target_plt = _np.ma.masked_values(field_target.data,self.xmsg)
						_plt.figure() ; _plt.contourf(data_target_plt[:,:],40) ; _plt.colorbar() ; 
						_plt.title('regridded') ; _plt.show()
		elif self.geometry == 'line':
			field_src.data[:] = dataextrap[:,:].transpose()
			field_target = regridme(field_src, field_target)
			if use_locstream:
				data[:,:] = _np.reshape(field_target.data.transpose(),(self.ny,self.nx))
			else:
				data[:,:] = field_target.data.transpose()[self.jmin:self.jmax+1,self.imin:self.imax+1]
		return data

def main_k_nearest_neighbour(k):
    X, y = make_blobs(n_samples=100,
                      n_features=2,
                      centers=2,
                      cluster_std=1.0,
                      center_box=(-10.0, 10.0))

    h = .4
    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 = np.c_[xx.ravel(), yy.ravel()]

    z_f = []
    for i_z in z:
        z_f.append(k_nearest_neighbour(X, y, i_z, k, False))

    zz = np.array(z_f).reshape(xx.shape)

    plt.figure()
    plt.contourf(xx, yy, zz, cmap=plt.cm.Paired)
    plt.axis('tight')
    plt.scatter(X[:, 0], X[:, 1], c=y)

    plt.show()

 def plot(x,y,field,filename,c=200):
     plt.figure()
     # define grid.
     xi = np.linspace(min(x),max(x),100)
     yi = np.linspace(min(y),max(y),100)
     # grid the data.
     si_lin = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     si_cub = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     print np.min(field)
     print np.max(field)
     plt.subplot(211)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_lin,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_lin,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Lineaarinen interpolointi')
     #plt.tight_layout()
     plt.subplot(212)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_cub,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_cub,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Kuubinen interpolointi')
     plt.savefig(filename)

def plot_decision_surface(axes, clusters, X, Y=None):
    import matplotlib.pylab as pylab
    import numpy as np

    def kmeans_predict(clusters, X):
        from ..core import distance
        dist_m = distance.minkowski_dist(clusters, X)
        print 'dist_m:', dist_m.shape
        pred = np.argmin(dist_m, axis=0)
        print 'pred:', pred.shape
        return pred

    # step size in the mesh
    h = (np.max(X, axis=0) - np.min(X, axis=0)) / 100.0
    # create a mesh to plot in
    x_min = np.min(X, axis=0) - 1
    x_max = np.max(X, axis=0) + 1
    xx, yy = np.meshgrid(np.arange(x_min[0], x_max[0], h[0]),
                         np.arange(x_min[1], x_max[1], h[1]))
    Z = kmeans_predict(clusters, np.c_[xx.ravel(), yy.ravel()])
    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    pylab.set_cmap(pylab.cm.Paired)
    pylab.axes(axes)
    pylab.contourf(xx, yy, Z, cmap=pylab.cm.Paired)
    pylab.xlim(np.min(xx), np.max(xx))
    pylab.ylim(np.min(yy), np.max(yy))
    #pylab.axis('off')

    # Plot also the training points
    if Y is not None:
        pylab.scatter(X[:,0], X[:,1], c=Y)
    else:
        pylab.scatter(X[:,0], X[:,1])
    pylab.scatter(clusters[:,0], clusters[:,1], s=200, marker='x', color='white')

    def plot_electric_field(self,sp=10,scales = 500000,cont_scale=90,savefigs = False):

        fig = plt.figure(figsize=(7.0, 7.0))
        xplot = self.x*1e6
        yplot = self.y*1e6
        X,Y = np.meshgrid(xplot,yplot)
        try:
            plt.contourf(X,Y,self.mode_field,cont_scale)
        except AttributeError:
             raise NotImplementedError("interpolate before plotting")

        plt.quiver(X[::sp,::sp], Y[::sp,::sp], np.real(self.E[::sp,::sp,0]), np.real(self.E[::sp,::sp,1]),scale = scales,headlength=7)
        plt.xlabel(r'$x(\mu m)$')
        plt.ylabel(r'$y(\mu m)$')
        #plt.title(r'mode$=$'+str(self.mode)+', '+'  $n_{eff}=$'+str(self.neff.real)+str(self.neff.imag)+'j')
        if savefigs == True:    
            plt.savefig('mode'+str(self.mode)+'.eps',bbox_inches ='tight')
        
            D = {}
            D['X'] = X
            D['Y'] = Y
            D['Z'] = self.mode_field
            D['u'] = np.real(self.E[::sp,::sp,0])
            D['v'] = np.real(self.E[::sp,::sp,1])
            D['scale'] = scales
            D['cont_scale'] = 90
            D['sp'] = sp
            savemat('mode'+str(self.mode)+'.mat',D)
        return None

	def drown_field(self,data,mask,drown):
		''' drown_field is a wrapper around the fortran code fill_msg_grid.
		depending on the output geometry, applies land extrapolation on 1 or N levels'''
		if self.geometry == 'surface':
			for kz in _np.arange(self.nz):
				tmpin = data[kz,:,:].transpose()
				if self.debug and kz == 0:
					tmpin_plt = _np.ma.masked_values(tmpin,self.xmsg)
					_plt.figure() ; _plt.contourf(tmpin_plt.transpose(),40) ; _plt.colorbar() ; 
					_plt.title('normalized before drown')
				if drown == 'ncl':
					tmpout = _fill.mod_poisson.poisxy1(tmpin,self.xmsg, self.guess, self.gtype, \
					self.nscan, self.epsx, self.relc)
				elif drown == 'sosie':
					tmpout = _mod_drown_sosie.mod_drown.drown(self.kew,tmpin,mask[kz,:,:].T,\
					nb_inc=200,nb_smooth=40)
				data[kz,:,:] = tmpout.transpose()
				if self.debug and kz == 0:
					_plt.figure() ; _plt.contourf(tmpout.transpose(),40) ; _plt.colorbar() ; 
					_plt.title('normalized after drown') 
					_plt.show()
		elif self.geometry == 'line':
			tmpin = data[:,:].transpose()
			if drown == 'ncl':
				tmpout = _fill.mod_poisson.poisxy1(tmpin,self.xmsg, self.guess, self.gtype, \
				self.nscan, self.epsx, self.relc)
			elif drown == 'sosie':
				tmpout = _mod_drown_sosie.mod_drown.drown(self.kew,tmpin,mask[:,:].T,\
				nb_inc=200,nb_smooth=40)
			data[:,:] = tmpout.transpose()
		return data

    def plot_field(self,x,y,u=None,v=None,F=None,contour=False,outdir=None,plot='quiver',figname='_field',format='eps'):
        outdir = self.set_dir(outdir)
        p = 64
        if F is None: F=self.calc_F(u,v)

        plt.close('all')
        plt.figure()
        #fig, axes = plt.subplots(nrows=1)
        if contour:
            plt.hold(True)
            plt.contourf(x,y,F)
        if plot=='quiver':
            plt.quiver(x[::p],y[::p],u[::p],v[::p],scale=0.1)
        
        if plot=='pcolor':
            plt.pcolormesh(x[::4],y[::4],F[::4],cmap=plt.cm.Pastel1)
            plt.colorbar()

        if plot=='stream':
            speed = F[::16]
            plt.streamplot(x[::16], y[::16], u[::16], v[::16], density=(1,1),color='k')

        plt.xlabel('$x$ (a.u.)')
        plt.ylabel('$y$ (a.u.)')
        
        plt.savefig(os.path.join(outdir,figname+'.'+format),format=format,dpi=320,bbox_inches='tight')
        plt.close()

def plot_electric(x,y,u,w,beta,a,V,Delta):
    ele = np.zeros([3,len(x),len(y)],dtype=np.complex128)
    for i,xx in enumerate(x):
        for j, yy in enumerate(y):
            ele[:,i,j] = electric_field(xx,yy,u,w,beta,a,V,0,0,Delta)
    abss = (np.abs(ele[0,:,:])**2 + np.abs(ele[1,:,:])**2 + np.abs(ele[2,:,:])**2)**0.5
    X,Y = np.meshgrid(x,y)
    fig = plt.figure()
    plt.contourf(X,Y,abss)
    plt.show()
    return 0

 def plot(self,key='AOD'):
     """
     Create a plot of a variable over the ORACLES study area. 
     
     Parameters
     ----------
     key : string
     See names for available datasets to plot.
     
     Modification history
     --------------------
     Written: Michael Diamond, 08/16/2016, Seattle, WA
     Modified: Michael Diamond, 08/21/2016
        -Added ORACLES routine flight plan, Walvis Bay (orange), and Ascension Island
     Modified: Michael Diamond, 09/02/2016, Swakopmund, Namibia
         -Updated flight track
     """
     plt.clf()
     size = 16
     font = 'Arial'
     m = Basemap(llcrnrlon=self.lon.min(),llcrnrlat=self.lat.min(),urcrnrlon=self.lon.max(),\
     urcrnrlat=self.lat.max(),projection='merc',resolution='i')
     m.drawparallels(np.arange(-180,180,5),labels=[1,0,0,0],fontsize=size,fontname=font)
     m.drawmeridians(np.arange(0,360,5),labels=[1,1,0,1],fontsize=size,fontname=font)
     m.drawmapboundary(linewidth=1.5)        
     m.drawcoastlines()
     m.drawcountries()
     m.drawmapboundary(fill_color='steelblue')
     m.fillcontinents(color='floralwhite',lake_color='steelblue',zorder=0)
     if key == 'AOD':
         m.pcolormesh(self.lon,self.lat,self.ds['%s' % key],cmap=self.colors['%s' % key],\
         latlon=True,vmin=self.v['%s' % key][0],vmax=self.v['%s' % key][1])
         cbar = m.colorbar()
         cbar.ax.tick_params(labelsize=size-2) 
         cbar.set_label('[%s]' % self.units['%s' % key],fontsize=size,fontname=font)
     elif key == 'ATYP':
         m.pcolormesh(self.lon,self.lat,self.ds['%s' % key],cmap=self.colors['%s' % key],\
         latlon=True,vmin=self.v['%s' % key][0],vmax=self.v['%s' % key][1])
         plt.contourf(np.array(([5,1],[3,2])),cmap=self.colors['%s' % key],levels=[0,1,2,3,4,5])
         cbar = m.colorbar(ticks=[0,1,2,3,4,5])
         cbar.ax.set_yticklabels(['Sea Salt','Sulphate','Organic C','Black C','Dust'])
         cbar.ax.tick_params(labelsize=size-2) 
     else:
         print('Error: Invalid key. Check names for available datasets.')
     m.scatter(14.5247,-22.9390,s=250,c='orange',marker='D',latlon=True)
     m.scatter(-14.3559,-7.9467,s=375,c='c',marker='*',latlon=True)
     m.scatter(-5.7089,-15.9650,s=375,c='chartreuse',marker='*',latlon=True)
     m.plot([14.5247,13,0],[-22.9390,-23,-10],c='w',linewidth=5,linestyle='dashed',latlon=True)
     m.plot([14.5247,13,0],[-22.9390,-23,-10],c='k',linewidth=3,linestyle='dashed',latlon=True)
     plt.title('%s from MSG SEVIRI on %s/%s/%s at %s UTC' % \
     (self.names['%s' % key],self.month,self.day,self.year,self.time),fontsize=size+4,fontname=font)
     plt.show()

def Decision_Surface(data, target, model, surface=True, probabilities=False, cell_size=.01):
    # Get bounds
    x_min, x_max = data[data.columns[0]].min(), data[data.columns[0]].max()
    y_min, y_max = data[data.columns[1]].min(), data[data.columns[1]].max()
    
    # Create a mesh
    xx, yy = np.meshgrid(np.arange(x_min, x_max, cell_size), np.arange(y_min, y_max, cell_size))
    meshed_data = pd.DataFrame(np.c_[xx.ravel(), yy.ravel()])
    
    # Add interactions
    for i in range(data.shape[1]):
        if i <= 1:
            continue
        
        meshed_data = np.c_[meshed_data, np.power(xx.ravel(), i)]

    if model != None:
        # Predict on the mesh
        if probabilities:
            Z = model.predict_proba(meshed_data)[:, 1].reshape(xx.shape)
        else:
            Z = model.predict(meshed_data).reshape(xx.shape)
    
    # Plot mesh and data
    if data.shape[1] > 2:
        plt.title("humor^(" + str(range(1,data.shape[1])) + ") and number_pets")
    else:
        plt.title("humor and number_pets")
    plt.xlabel("humor")
    plt.ylabel("number_pets")
    if surface and model != None:
        cs = plt.contourf(xx, yy, Z, cmap=plt.cm.coolwarm, alpha=0.4)
    color = ["blue" if t == 0 else "red" for t in target]
    plt.scatter(data[data.columns[0]], data[data.columns[1]], color=color)

def plot_contour(z,x,y,title="TITLE",xtitle="",ytitle="",xrange=None,yrange=None,plot_points=0,contour_levels=20,
                 cmap=None,cbar=True,fill=False,cbar_title="",show=1):

    fig = plt.figure()

    if fill:
        fig = plt.contourf(x, y, z.T, contour_levels, cmap=cmap, origin='lower')
    else:
        fig = plt.contour( x, y, z.T, contour_levels, cmap=cmap, origin='lower')

    if cbar:
        cbar = plt.colorbar(fig)
        cbar.ax.set_ylabel(cbar_title)


    plt.title(title)
    plt.xlabel(xtitle)
    plt.ylabel(ytitle)

    # the scatter plot:
    if plot_points:
        axScatter = plt.subplot(111)
        axScatter.scatter( np.outer(x,np.ones_like(y)), np.outer(np.ones_like(x),y))

    # set axes range
    plt.xlim(xrange)
    plt.ylim(yrange)

    if show:
        plt.show()

    return fig

def Plot2DwSVM(datax, datay):
    """
    plot data and corresponding SVM results.
    """
    clf = svm.SVC()
    clf.fit(datax, datay)
    step = 0.02
    x_min, x_max = datax[:,0].min(), datax[:,0].max()
    y_min, y_max = datax[:,1].min(), datax[:,1].max()
    xx, yy = numpy.meshgrid(numpy.arange(x_min, x_max, step), numpy.arange(y_min, y_max, step))
    Z = clf.decision_function(numpy.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    pylab.contourf(xx, yy, Z, 10, cmap=pylab.cm.Oranges)
    pylab.scatter(datax[datay == 1,0], datax[datay == 1,1], c='b', s=50)
    pylab.scatter(datax[datay != 1,0], datax[datay != 1,1], c='r', s=50)
    pylab.show()

def Decision_Surface(X, target, model, cell_size=.01, surface=True, points=True):
    # Get bounds
    x_min, x_max = X[:, 0].min(), X[:, 0].max()
    y_min, y_max = X[:, 1].min(), X[:, 1].max()
    
    # Create a mesh
    xx, yy = np.meshgrid(np.arange(x_min, x_max, cell_size), np.arange(y_min, y_max, cell_size))
    meshed_data = pd.DataFrame(np.c_[xx.ravel(), yy.ravel()])
    
    # Add interactions
    for i in range(X.shape[1]):
        if i <= 1:
            continue
        
        meshed_data = np.c_[meshed_data, np.power(xx.ravel(), i)]

    Z_flat = model.predict(meshed_data)
    Z = Z_flat.reshape(xx.shape)
    
    # Plot mesh and data
    plt.title("humor and number_pets")
    plt.xlabel("humor")
    plt.ylabel("number_pets")
    
    if surface:
        cs = plt.contourf(xx, yy, Z, color=colorizer(Z_flat))
    
    if points:
        plt.scatter(X[:, 0], X[:, 1], color=colorizer(target), linewidth=0, s=20)
    
    plt.xlabel("Feature 1")
    plt.xlabel("Feature 2")
    plt.show()

def plot(AB_final,lams,lamp2,title):
    X,Y = np.meshgrid(lams,lamp2)
    Z = AB_final
    
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    plt.gca().get_xaxis().get_major_formatter().set_useOffset(False)
    plt.gca().get_yaxis().get_major_formatter().set_useOffset(False)
    plt.title(title)
    ax.set_xlabel(r'$\lambda_p(\mu m)$')
    ax.set_ylabel(r'$\lambda_s(\mu m)$')
    ax.set_zlabel(r'$P_{idler}(dBm)$')

    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet, linewidth=0, antialiased=False)
    fig.colorbar(surf, shrink=0.5, aspect=5)
    plt.show()
    
    fig = plt.figure()
    axi = plt.contourf(X,Y,Z, cmap=cm.jet)
    plt.gca().get_xaxis().get_major_formatter().set_useOffset(False)
    plt.gca().get_yaxis().get_major_formatter().set_useOffset(False)
    plt.xlabel(r'$\lambda_p(\mu m)$')
    plt.ylabel(r'$\lambda_s(\mu m)$')
    plt.title(title)
    fig.colorbar(axi, shrink=0.5, aspect=5)
    plt.show()
    return 0

def contour(X,Y,Z, 
           extent=None, vrange=None, levels=None, extend='both', 
           inaxis=None,
           cmap=None, addcolorbar=True, clabel=None,
           smooth=True, smoothlen=None):
    '''Build a super fancy contour image'''
    
    # Build up some nice ranges and levels to be plotted 
    XX, YY = _getmesh(X,Y,Z)
    extent = _getextent(extent, XX, YY)
    vrange = _getvrange(vrange, XX,YY,Z, inaxis=inaxis)
    levels = _getlevels(levels, vrange)
    cmap   = _getcmap(cmap)
    
    # Smooth if needed
    if smooth:
        X,Y,Z = _smooth(X,Y,Z, smoothlen)
    
    cs = pylab.contourf(X, Y, Z, levels,
                       vmin=vrange[0], vmax=vrange[1],
                       extent=extent, extend='both',
                       cmap=cmap)
    ccs = pylab.contour(X, Y, Z, levels, vmin=vrange[0], vmax=vrange[1],
                        cmap=cmap)
    
    # setup a colorbar, add in the lines, and then return it all out.
    if addcolorbar:
        cb = colorbar(cs, ccs, levels=levels, clabel=clabel)
        return cs, ccs, cb
    
    return cs, ccs

def decision_surface(data, target, model, surface=True, probabilities=False, cell_size=0.01):
    plt.rcParams["figure.figsize"] = [10.0, 10.0]

    # Get bounds
    x_min, x_max = data[data.columns[0]].min(), data[data.columns[0]].max()
    y_min, y_max = data[data.columns[1]].min(), data[data.columns[1]].max()

    # Create a mesh
    xx, yy = np.meshgrid(np.arange(x_min, x_max, cell_size), np.arange(y_min, y_max, cell_size))
    meshed_data = pd.DataFrame(np.c_[xx.ravel(), yy.ravel()])

    if model != None:
        # Predict on the mesh
        if probabilities:
            Z = model.predict_proba(meshed_data)[:, 1].reshape(xx.shape)
        else:
            Z = model.predict(meshed_data).reshape(xx.shape)

    # Plot mesh and data
    plt.title("Decision Surface")

    plt.xlabel(data.columns[0])
    plt.ylabel(data.columns[1])

    if surface and model != None:
        cs = plt.contourf(xx, yy, Z, cmap=plt.cm.coolwarm_r, alpha=0.4)

    color = ["red" if t == 0 else "blue" for t in target]

    plt.scatter(data[data.columns[0]], data[data.columns[1]], color=color)
    plt.show()

def plot_solution(problem_data, solution):
    from numpy import linspace
    from matplotlib.pylab import contour, colorbar, contourf, xlabel, ylabel, title, show
    from matplotlib.mlab import griddata

    print(" * Preparing for plotting...")
    NN = problem_data["NN"]

    # Extract node coordinates seperately for plotting
    x, y = [0] * NN, [0] * NN
    for i, node in enumerate(problem_data["nodes"]):
        x[i] = node[0]
        y[i] = node[1]

    # Refine the contour plot mesh for a "smoother" image, by generating a
    # 200*200 grid
    xi = linspace(min(x), max(x), 200)
    yi = linspace(min(y), max(y), 200)

    # Approximate the mid values from neighbors
    zi = griddata(x, y, solution, xi, yi)

    print(" * Plotting...")
    # Plot the contour lines with black
    contour(xi, yi, zi, 15, linewidths=0.5, colors='k')
    # Plot the filled contour plot
    plot = contourf(xi, yi, zi, 15, antialiased=True)

    colorbar(plot, format="%.3f").set_label("T")
    xlabel('X')
    ylabel('Y')
    title("Contour plot of T values for {0}".format(problem_data["title"]))

    show()

def plot_svc(X, y, mysvc, bounds=None, grid=50):
    if bounds is None:
        xmin = np.min(X[:, 0], 0)
        xmax = np.max(X[:, 0], 0)
        ymin = np.min(X[:, 1], 0)
        ymax = np.max(X[:, 1], 0)
    else:
        xmin, ymin = bounds[0], bounds[0]
        xmax, ymax = bounds[1], bounds[1]
    aspect_ratio = (xmax - xmin) / (ymax - ymin)
    xgrid, ygrid = np.meshgrid(np.linspace(xmin, xmax, grid),
                              np.linspace(ymin, ymax, grid))
    plt.gca(aspect=aspect_ratio)
    plt.xlim(xmin, xmax)
    plt.ylim(ymin, ymax)
    plt.xticks([])
    plt.yticks([])
    plt.hold(True)
    plt.plot(X[y == 1, 0], X[y == 1, 1], 'bo')
    plt.plot(X[y == -1, 0], X[y == -1, 1], 'ro')
    
    box_xy = np.append(xgrid.reshape(xgrid.size, 1), ygrid.reshape(ygrid.size, 1), 1)
    if mysvc is not None:
        scores = mysvc.decision_function(box_xy)
    else:
        print 'You must have a valid SVC object.'
        return None;
    
    CS=plt.contourf(xgrid, ygrid, scores.reshape(xgrid.shape), alpha=0.5, cmap='jet_r')
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[0], colors='k', linestyles='solid', linewidths=1.5)
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[-1,1], colors='k', linestyles='dashed', linewidths=1)
    plt.plot(mysvc.support_vectors_[:,0], mysvc.support_vectors_[:,1], 'ko', markerfacecolor='none', markersize=10)
    CB = plt.colorbar(CS)

def plot(points, col):
    n = 256
    x = np.linspace(-100, 100, n)
    y = np.linspace(-100, 100, n)
    X, Y = np.meshgrid(x, y)
    
    xs = []
    ys = []
    
    pl.contourf(X, Y, fun([X, Y]), 8, alpha=.75, cmap='jet')
    pl.contour(X, Y, fun([X, Y]), 8, colors='black', linewidth=.5) 
    
    for i in range(len(points)):
        xs.append(points[i][0])
        ys.append(points[i][1])
    
    pl.plot(xs, ys, marker='o', linestyle='--', color=str(col), label='Square')

def plot_decision_surface(axes, clf, X, Y):
    import numpy as np
    h=.02 # step size in the mesh
    # create a mesh to plot in
    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)
    axes.set_cmap(pylab.cm.Paired)
    pylab.axes(axes)
    pylab.contourf(xx, yy, Z)
    pylab.axis('off')

    # Plot also the training points
    pylab.scatter(X[:,0], X[:,1], c=Y)