def plot_knn_boundary():
    ## Training dataset preparation
    # use sklearn iris dataset
    iris_dataset = datasets.load_iris()


    # first two dimensions as the features
    # it's easy to plot boundary in 2D
    train_data = iris_dataset.data[:,:2] 
    print "init:",train_data

    # get labels
    labels = iris_dataset.target # labels
    print "init2:",labels

    
    ## Test dataset preparation
    h = 0.1
    
    x0_min = train_data[:,0].min() - 0.5
    x0_max = train_data[:,0].max() + 0.5
    
    x1_min = train_data[:,1].min() - 0.5
    x1_max = train_data[:,1].max() + 0.5
    
    x0_features, x1_features = np.meshgrid(np.arange(x0_min, x0_max, h), 
                                           np.arange(x1_min, x1_max, h))
    
    # test dataset are samples from the whole regions of feature domains
    test_data = np.c_[x0_features.ravel(), x1_features.ravel()]
    
    ## KNN classification
    p_labels = []   # prediction labels
    for test_sample in test_data:
        # knn prediction
        p_label = knn_predict(train_data, labels, test_sample, n_neighbors = 6)
        p_labels.append(p_label)
    
    # list to array
    p_labels = np.array(p_labels)
    p_labels = p_labels.reshape(x0_features.shape)
    
    ## Boundary plotting  边界策划
    pl.figure(1)
    pl.set_cmap(pl.cm.Paired)
    pl.pcolormesh(x0_features, x1_features, p_labels)
    
    pl.scatter(train_data[:,0], train_data[:,1], c = labels)
    # x y轴的名称
    pl.xlabel('feature 0')
    pl.ylabel('feature 1')

    # 设置x,y轴的上下限
    pl.xlim(x0_features.min(), x0_features.max())
    pl.ylim(x1_features.min(), x1_features.max())
    # 设置x,y轴记号
    pl.xticks(())
    pl.yticks(())
    
    pl.show()

def saveBEVImageWithAxes(data, outputname, cmap = None, xlabel = 'x [m]', ylabel = 'z [m]', rangeX = [-10, 10], rangeXpx = None, numDeltaX = 5, rangeZ = [7, 62], rangeZpx = None, numDeltaZ = 5, fontSize = 16):
    '''
    
    :param data:
    :param outputname:
    :param cmap:
    '''
    aspect_ratio = float(data.shape[1])/data.shape[0]
    fig = pylab.figure()
    Scale = 8
    # add +1 to get axis text
    fig.set_size_inches(Scale*aspect_ratio+1,Scale*1)
    ax = pylab.gca()
    #ax.set_axis_off()
    #fig.add_axes(ax)
    if cmap != None:
        pylab.set_cmap(cmap)
    
    #ax.imshow(data, interpolation='nearest', aspect = 'normal')
    ax.imshow(data, interpolation='nearest')
    
    if rangeXpx == None:
        rangeXpx = (0, data.shape[1])
    
    if rangeZpx == None:
        rangeZpx = (0, data.shape[0])
        
    modBev_plot(ax, rangeX, rangeXpx, numDeltaX, rangeZ, rangeZpx, numDeltaZ, fontSize, xlabel = xlabel, ylabel = ylabel)
    #plt.savefig(outputname, bbox_inches='tight', dpi = dpi)
    pylab.savefig(outputname, dpi = data.shape[0]/Scale)
    pylab.close()
    fig.clear()

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 twoDimOrderPlot(outpath, base_name, title, obj_name, base_filename, order_num, data, x):
    pl.figure('2d order image', facecolor='white', figsize=(8, 5))
    pl.cla()
    pl.title(title + ', ' + base_name + ", order " + str(order_num), fontsize=14)
    pl.xlabel('wavelength($\AA$)', fontsize=12)
    pl.ylabel('row (pixel)', fontsize=12)
    #pl.imshow(img, aspect='auto')
    #pl.imshow(data, vmin=0, vmax=1024, aspect='auto')
    
    pl.imshow(exposure.equalize_hist(data), origin='lower', 
                  extent=[x[0], x[-1], 0, data.shape[0]], aspect='auto')      
#     from matplotlib import colors
#     norm = colors.LogNorm(data.mean() + 0.5 * data.std(), data.max(), clip='True')
#     pl.imshow(data, norm=norm, origin='lower',
#                   extent=[x[0], x[-1], 0, data.shape[0]], aspect='auto')               
    pl.colorbar()
    pl.set_cmap('jet')
#     pl.set_cmap('Blues_r')
    fn = constructFileName(outpath, base_name, order_num, base_filename)
    pl.savefig(fn)
    log_fn(fn)
    pl.close()
    
#     np.save(fn[:fn.rfind('.')], data)
    
    return

def plot_figs(fig_num, elev, azim):
    fig = pl.figure(fig_num, figsize=(4, 3))
    pl.clf()
    ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=elev, azim=azim)

    pl.set_cmap(pl.cm.hot_r)

    pts = ax.scatter(a[::10], b[::10], c[::10], c=density,
                     marker='+', alpha=.4)

    Y = np.c_[a, b, c]
    U, pca_score, V = np.linalg.svd(Y, full_matrices=False)
    x_pca_axis, y_pca_axis, z_pca_axis = V.T*pca_score/pca_score.min()

#ax.quiver(0.1*x_pca_axis, 0.1*y_pca_axis, 0.1*z_pca_axis,
#                x_pca_axis, y_pca_axis, z_pca_axis,
#                color=(0.6, 0, 0))

    x_pca_axis, y_pca_axis, z_pca_axis = 3*V.T
    x_pca_plane = np.r_[x_pca_axis[:2], - x_pca_axis[1::-1]]
    y_pca_plane = np.r_[y_pca_axis[:2], - y_pca_axis[1::-1]]
    z_pca_plane = np.r_[z_pca_axis[:2], - z_pca_axis[1::-1]]
    x_pca_plane.shape = (2, 2)
    y_pca_plane.shape = (2, 2)
    z_pca_plane.shape = (2, 2)
    ax.plot_surface(x_pca_plane, y_pca_plane, z_pca_plane)
    ax.w_xaxis.set_ticklabels([])
    ax.w_yaxis.set_ticklabels([])
    ax.w_zaxis.set_ticklabels([])

def sparect_plot(outpath, base_name, order_num, obj, flat):

    pl.figure('spatially rectified', facecolor='white', figsize=(8, 5))
    pl.cla()
    pl.suptitle('spatially rectified, {}, order {}'.format(base_name, order_num), fontsize=14)
    pl.set_cmap('Blues_r')

    obj_plot = pl.subplot(2, 1, 1)
    try:
        obj_plot.imshow(exposure.equalize_hist(obj))
    except:
        obj_plot.imshow(obj)
    obj_plot.set_title('object')
#     obj_plot.set_ylim([1023, 0])
    obj_plot.set_xlim([0, 1023])
    
    flat_plot = pl.subplot(2, 1, 2)
    try:
        flat_plot.imshow(exposure.equalize_hist(flat))
    except:
        flat_plot.imshow(flat)
    flat_plot.set_title('flat')
#     flat_plot.set_ylim([1023, 0])
    flat_plot.set_xlim([0, 1023])
 
    pl.tight_layout()
    pl.savefig(constructFileName(outpath, base_name, order_num, 'sparect.png'))
    pl.close()

def cutouts_plot(outpath, base_name, order_num, obj, flat, top_trace, bot_trace, trace):
    
    pl.figure('traces', facecolor='white', figsize=(8, 5))
    pl.cla()
    pl.suptitle('order cutouts, {}, order {}'.format(base_name, order_num), fontsize=14)
    pl.set_cmap('Blues_r')

    obj_plot = pl.subplot(2, 1, 1)
    try:
        obj_plot.imshow(exposure.equalize_hist(obj))
    except:
        obj_plot.imshow(obj)
    obj_plot.plot(np.arange(1024), top_trace, 'y-', linewidth=1.5)
    obj_plot.plot(np.arange(1024), bot_trace, 'y-', linewidth=1.5)
    obj_plot.plot(np.arange(1024), trace, 'y-', linewidth=1.5)
    obj_plot.set_title('object')
#     obj_plot.set_ylim([1023, 0])
    obj_plot.set_xlim([0, 1023])
    
    flat_plot = pl.subplot(2, 1, 2)
    try:
        flat_plot.imshow(exposure.equalize_hist(flat))
    except:
        flat_plot.imshow(flat)
    flat_plot.plot(np.arange(1024), top_trace, 'y-', linewidth=1.5)
    flat_plot.plot(np.arange(1024), bot_trace, 'y-', linewidth=1.5)    
    flat_plot.plot(np.arange(1024), trace, 'y-', linewidth=1.5)    
    flat_plot.set_title('flat')
#     flat_plot.set_ylim([1023, 0])
    flat_plot.set_xlim([0, 1023])
 
    pl.tight_layout()
    pl.savefig(constructFileName(outpath, base_name, order_num, 'cutouts.png'))
    pl.close()

def analysis():
    feature_importances_array = np.load("Default_AlexNet.npy")
    mean_feature_importances = np.mean(feature_importances_array, axis=0)

    for feature_importance, metric in zip(mean_feature_importances, METRIC_LIST):
        print("{}\t{}".format(metric, feature_importance))

    time_indexes = np.arange(1, feature_importances_array.shape[0] + 1)
    feature_importances_cumsum = np.cumsum(feature_importances_array, axis=0)
    feature_importances_mean = feature_importances_cumsum
    for column_index in range(feature_importances_mean.shape[1]):
        feature_importances_mean[:, column_index] = feature_importances_cumsum[:, column_index] / time_indexes

    index_ranks = np.flipud(np.argsort(mean_feature_importances))

    chosen_records = np.cumsum(mean_feature_importances[index_ranks]) <= 0.95
    chosen_index_ranks = index_ranks[chosen_records]

    sorted_mean_feature_importances = mean_feature_importances[chosen_index_ranks]
    sorted_metric_list = np.array(METRIC_LIST)[chosen_index_ranks]

    remaining = np.sum(mean_feature_importances[index_ranks[~chosen_records]])
    print("remaining is {:.4f}.".format(remaining))
    sorted_mean_feature_importances = np.hstack((sorted_mean_feature_importances, remaining))
    sorted_metric_list = np.hstack((sorted_metric_list, 'others'))

    pylab.pie(sorted_mean_feature_importances, labels=sorted_metric_list, autopct='%1.1f%%', startangle=0)
    pylab.axis('equal')
    pylab.set_cmap('plasma')
    pylab.show()

def traces_plot(outpath, base_name, order_num, obj, flat, top_trace, bot_trace):
    
    pl.figure('traces', facecolor='white', figsize=(8, 5))
    pl.cla()
    pl.suptitle('order edge traces, {}, order {}'.format(base_name, order_num), fontsize=14)
    pl.set_cmap('Blues_r')
    pl.rcParams['ytick.labelsize'] = 8

    obj_plot = pl.subplot(1, 2, 1)
    obj_plot.imshow(exposure.equalize_hist(obj))
    obj_plot.plot(np.arange(1024), top_trace, 'y-', linewidth=1.5)
    obj_plot.plot(np.arange(1024), bot_trace, 'y-', linewidth=1.5)

    obj_plot.set_title('object')
    obj_plot.set_ylim([1023, 0])
    obj_plot.set_xlim([0, 1023])

    
    flat_plot = pl.subplot(1, 2, 2)
    flat_plot.imshow(exposure.equalize_hist(flat))
    flat_plot.plot(np.arange(1024), top_trace, 'y-', linewidth=1.5)
    flat_plot.plot(np.arange(1024), bot_trace, 'y-', linewidth=1.5)    
    flat_plot.set_title('flat')
    flat_plot.set_ylim([1023, 0])
    flat_plot.set_xlim([0, 1023])
 
    pl.tight_layout()
    pl.savefig(constructFileName(outpath, base_name, order_num, 'traces.png'))
    pl.close()

def plot_sa_diff_figure(control_dataset, data, sa_mask):
    f = plt.figure('sa_diff')
    plt.set_cmap('RdGy_r')
    graph_settings = (
	    ((-4, 4), np.arange(-4, 4.1, 2)),
	    ((-6, 6), np.arange(-6, 6.1, 3)),
	    ((-0.7, 0.7), np.arange(-0.6, 0.61, 0.3)),
	    ((-4, 4), np.arange(-4, 4.1, 2)),
	    ((-0.2, 0.2), np.arange(-0.2, 0.21, 0.1)))

    variables = ['precip', 'surf_temp', 'q', 'field1389', 'field1385']
    nice_names = {'precip': '$\Delta$Precip (mm/day)', 
		  'surf_temp': '$\Delta$Surf temp (K)', 
		  'q':'$\Delta$Humidity (g/kg)', 
		  'field1389': '$\Delta$NPP (g/m$^2$/day)', 
		  'field1385': '$\Delta$Soil moisture'}

    f.subplots_adjust(hspace=0.2, wspace=0.1)
    for i in range(len(variables)):
	variable = variables[i]
	ax = plt.subplot(2, 3, i + 1)
	ax.set_title(nice_names[variable])
	variable_diff = data['data']['1pct'][variable] - data['data']['ctrl'][variable]
	if variable == 'field1389':
	    variable_diff *= 24*60*60*1000 # per s to per day, kg to g.
	lons, lats = get_vars_from_control_dataset(control_dataset)
	vmin, vmax = graph_settings[i][0]
	#general_plot(control_dataset, variable_diff.mean(axis=0), vmin=graph_settings[i][0][0], vmax=graph_settings[i][0][1], loc='sa', sa_mask=sa_mask)
	plot_data = variable_diff.mean(axis=0)
	#plot_south_america(lons, lats, sa_mask, plot_data, vmin, vmax)

	if variable in ('surf_temp', 'precip', 'q'):
	    # unmasked.
	    data_masked = plot_data
	    plot_lons, plot_data = extend_data(lons, lats, data_masked)
	else:
	    data_masked = np.ma.array(plot_data, mask=sa_mask)
	    plot_lons, plot_data = extend_data(lons, lats, data_masked)

	lons, lats = np.meshgrid(plot_lons, lats)

	m = Basemap(projection='cyl', resolution='c', llcrnrlat=-60, urcrnrlat=15, llcrnrlon=-85, urcrnrlon=-32)
	x, y = m(lons, lats)

	m.pcolormesh(x, y, plot_data, vmin=vmin, vmax=vmax)

	m.drawcoastlines()
	if i == 0 or i == 3:
	    m.drawparallels(np.arange(-60.,15.,10.), labels=[1, 0, 0, 0], fontsize=10)
	elif i == 2 or i == 4:
	    m.drawparallels(np.arange(-60.,15.,10.), labels=[0, 1, 0, 0], fontsize=10)
	else:
	    m.drawparallels(np.arange(-60.,15.,10.))

	m.drawmeridians(np.arange(-90.,-30.,10.), labels=[0, 0, 0, 1], fontsize=10)

	cbar = m.colorbar(location='bottom', pad='7%', ticks=graph_settings[i][1])

def plotting():
    conf = [[0 for x in range(L)] for y in range(L)]
    for k in range(N):
        x, y = x_y(k, L)
        conf[x][y] = S[k]
    pylab.imshow(conf, extent=[0, L, 0, L], interpolation='nearest')
    pylab.set_cmap('hot')
    pylab.title('Local_'+ str(T) + '_' + str(L))
    pylab.savefig('plot_A2_local_'+ str(T) + '_' + str(L)+ '.png')
    pylab.show()

def cmap_smooth(I=None,axh=None,Nlevels=256,cmap_lin=None):
  if cmap_lin is None:
    cmap_lin = pl.cm.jet
  if I is None:
    if not axh:
      axh = pl.gca()
    ihandles = axh.findobj(matplotlib.image.AxesImage)
    ih = ihandles[-1]
    I = ih.get_array()
  levels = np.percentile(I.ravel(),list(np.linspace(0,100,Nlevels)))
  cmap_nonlin = nlcmap(cmap_lin,levels)
  pl.set_cmap(cmap_nonlin)

def plot_classification(X, y, y_pred, keys, title, clf):
    print 'plot_classification(X=%s, y=%s, y_pred=%s, keys=%s, title=%s)' % (X.shape, 
        y.shape, y_pred.shape, keys, title)
    h = .02 # step size in the mesh
    
    n_plots = len(keys)*(len(keys)-1)//2
    n_side = int(math.sqrt(float(n_plots)))

    cnt = 1
    for i0 in range(len(keys)):
        for i1 in range(i0+1, len(keys)): 
            # create a mesh to plot in
            x_min, x_max = X[:,i0].min()-1, X[:,i0].max()+1
            y_min, y_max = X[:,i1].min()-1, X[:,i1].max()+1
            xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
            
            pl.set_cmap(pl.cm.Paired)

            # 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].
            print 'subplot(%d, %d, cnt=%d)' % (n_side, n_side, cnt)
            pl.subplot(n_side, n_side, cnt)
            print 'xx.size=%s, xx.shape=%s, X.shape=%s' % (xx.size, xx.shape, X.shape)
            points = np.zeros([xx.size, X.shape[1]])
            points[:,i0] = xx.ravel() 
            points[:,i1] = yy.ravel()
            Z = clf.predict(points)

            # Put the result into a color plot
            Z = Z.reshape(xx.shape)
            pl.set_cmap(pl.cm.Paired)
            pl.contourf(xx, yy, Z)
            pl.axis('tight')
            
            #pl.xlabel(keys[0])
            #pl.ylabel(keys[1])

            # Plot also the training points
            #pl.scatter(X[:,0], X[:,1], c=y)
            
            plot_2d_histo_raw(X[:,i0], X[:,i1], y, keys[i0], keys[i1], x_max-x_min, y_max-y_min)

            #pl.title('%s vs %s' % (keys[i1], keys[i0]))

            pl.axis('tight')
            cnt +=1 
            if cnt > n_side ** 2:
                break
        if cnt > n_side ** 2:
            break
            
    pl.savefig(os.path.join('results', '%s.png' % title)) 
    pl.show()

def plot(polylines, img=None):
    import pylab as plt
    plt.figure(1)

    for n, c in enumerate(polylines):
        x = c[:, 0]
        y = c[:, 1]
        plt.plot(x, y, linewidth=3)
        plt.text(x[-1], y[-1], str(n + 1))
    if img is not None:
        plt.imshow(img, interpolation='none')
        plt.set_cmap('gray')
    plt.show()

 def _heatmap_engine(self, suptitle, numIter=None, filter_count=None, saveas=None, show=True, cmap="Blues"):
     '''
     An alternative way to visualize the annotations per iterative via
     a heat-map like construct. The rows are the GO annotations and the columns
     are iterations. The color intensity indicates how much a given annotation
     was present in a given iteration
     '''
     res = self.annotation_dat if numIter is None else self.annotation_dat[0:numIter]
     depth = self.depth
     
     pl.figure(num=1,figsize=(20,8))
 
     #AS is complete Annotation Set, max_count is the maximum number
     # of genes that appears in any single annotation entry
     (AS, max_count) = self._common_Y()
     
     #map is a grid Y=annotation X=iteration M(X,Y) = scaled count (c/max c)
     M = sp.zeros( (len(AS), len(res) ) )
     for (col, dat, _) in res:
         if len(dat) < 1: continue
         for (l,d,c,_) in dat:
             row = AS.index((d,l))
             M[row,col] = c #(c*1.0)/max_count
         
     #filter rows / pathways which dont show up often
     if not filter_count is None:
         assert type(filter_count) == int
         vals = sp.sum(M, axis=1)
         M = M[ vals >= filter_count, :]  #only pathways with at least filter count over all iterations
         AS = [ x for i,x in enumerate(AS) if vals[i] >= filter_count ]
         
     pl.imshow(M, interpolation='nearest', aspect=len(res)*1.0/len(AS), origin='lower')
     #for (l,d) in AS:
     #    row = AS.index((d,l))
     #    pl.text(-0.5, row, d[0:40], color="white", verticalalignment='center', fontsize=9)
     
     (descs, _labels) = zip(*AS)
     descs = [ d[0:30] for d in descs]  #truncate long descriptions
     ylocations = sp.array(range(len(descs)))
     pl.yticks(ylocations, descs, fontsize=9, verticalalignment='center')
     
     pl.set_cmap(cmap)
     pl.xticks(range(len(res)), range(1, len(res)+1))
     pl.ylabel("GO Annotation at Depth %d"%depth, fontsize=16)
     pl.xlabel("Iteration", fontsize=16)
     pl.colorbar(ticks=range(1, max_count+1))       
     if not suptitle is None:
         pl.title("Ontology Annotations per Iteration: %s"%suptitle, fontsize=18)
     if not saveas is None:
         pl.savefig(saveas)   
     if show: pl.show()

def showIMG(IMG, extent = None, ticks = False):
	pylab.imshow(IMG, origin = 'lower', extent = extent)

	pylab.set_cmap(pylab.cm.gray)
	
	if extent == None:
		pylab.ylim((0, IMG.shape[0] - 1))
		pylab.xlim((0, IMG.shape[1] - 1))
	else:
		pylab.ylim((extent[2], extent[3]))
		pylab.xlim((extent[0], extent[1]))
	
	if not ticks:
		pylab.gca().get_xaxis().set_ticks([])
		pylab.gca().get_yaxis().set_ticks([])
	pylab.gca().invert_yaxis()

def plot_sa_seasonal_figure(control_dataset, data, sa_mask):
    f = plt.figure('sa_seasonal')
    plt.set_cmap('RdGy_r')
    graph_settings = (
	    ((-4, 4), np.arange(-4, 4.1, 2)),
	    ((-6, 6), np.arange(-6, 6.1, 3)))

    variables = ['precip', 'surf_temp']
    nice_names = {'precip': '$\Delta$Precip (mm/day)', 
		  'surf_temp': '$\Delta$Surf temp (K)'}
    f.subplots_adjust(hspace=0.2, wspace=0.1)

    for j in range(len(variables)):
	for i, roll in enumerate([1, 10, 7, 4]):
	    titles = ('DJF', 'MAM', 'JJA', 'SON')
	    variable = variables[j]
	    ax = plt.subplot(2, 4, i + j * 4 + 1)
	    ax.set_title(titles[i])
	    variable_diff = data['data']['1pct'][variable] - data['data']['ctrl'][variable]
	    lons, lats = get_vars_from_control_dataset(control_dataset)
	    vmin, vmax = graph_settings[j][0]
	    plot_data = np.roll(variable_diff, roll, axis=0)[:3].mean(axis=0)

	    # unmasked.
	    data_masked = plot_data
	    plot_lons, plot_data = extend_data(lons, lats, data_masked)

	    lons, lats = np.meshgrid(plot_lons, lats)

	    m = Basemap(projection='cyl', resolution='c', llcrnrlat=-60, urcrnrlat=15, llcrnrlon=-85, urcrnrlon=-32)
	    x, y = m(lons, lats)

	    m.pcolormesh(x, y, plot_data, vmin=vmin, vmax=vmax)

	    m.drawcoastlines()
	    if i == 0:
		m.drawparallels(np.arange(-60.,15.,10.), labels=[1, 0, 0, 0], fontsize=10)
		ax.set_ylabel(nice_names[variable])
		ax.get_yaxis().set_label_coords(-0.25, 0.5)
	    elif i == 3:
		m.drawparallels(np.arange(-60.,15.,10.), labels=[0, 1, 0, 0], fontsize=10)
	    else:
		m.drawparallels(np.arange(-60.,15.,10.))

	    m.drawmeridians(np.arange(-90.,-30.,10.), labels=[0, 0, 0, 1], fontsize=10)

	    cbar = m.colorbar(location='bottom', pad='7%', ticks=graph_settings[j][1])

def _show(S, L, T, nsteps, suffix, show):
    def x_y(k, L):
        y = k // L
        x = k - y * L
        return x, y

    conf = [[0 for x in range(L)] for y in range(L)]
    for k in range(N):
        x, y = x_y(k, L)
        conf[x][y] = S[k]

    pylab.imshow(conf, extent=[0, L, 0, L], interpolation='nearest')
    pylab.set_cmap('hot')
    pylab.title('Local_'+ str(T) + '_' + str(L))
    pylab.savefig('plot_A2_local_'+ str(T) + '_' + str(L)+ suffix + '_' + str(nsteps) + '.png')
    if show:
        pylab.show()

def plot_classifier(X, Y, models, classMap=None):
  """ Plots classifier or classifiers on 2d plot """
  #handle single model
  if not isinstance(models, types.ListType):
    models = [models]
    titles = ["classifier"]
  else:
    titles = []
    for model in models:
      titles.append("Classifier...")
  
  #colors for different decisions
  colors = ["red", "green", "blue", "yellow", "black"]

  # create a mesh to plot in
  h = 500  # step size in mesh
  x_min, x_max = X[:, 0].min() - .002, X[:, 0].max() + .002
  y_min, y_max = X[:, 1].min() - .002, X[:, 1].max() + .002
  xx, yy = np.meshgrid(np.arange(x_min, x_max, (x_max-x_min)/h),
                       np.arange(y_min, y_max, (y_max-y_min)/h))
  
  #set cmap
  pl.set_cmap(pl.cm.Paired)
  for mi, clf in enumerate( models ):
    # Plot the decision boundary. For that, we will asign a color to each
    # point in the mesh [x_min, m_max]x[y_min, y_max].
    pl.subplot(1, len(models), mi + 1)
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    pl.set_cmap(pl.cm.Paired)
    pl.contourf(xx, yy, Z)
    #pl.axis('off')

    # Plot also the training points
    if classMap:
      for c,i in classMap.iteritems():
        x = X[Y==i] 
        pl.plot(x[:,0], x[:,1], "o", c=colors[i], label=c) 

    #legend and title
    pl.legend()
    pl.title(titles[mi])
  pl.show()

def show_spins(S0, S1, L, label):
    pylab.set_cmap('hot')
    conf0 = [[0 for x in range(L)] for y in range(L)]
    conf1 = [[0 for x in range(L)] for y in range(L)]
    for k in range(N):
        y = k // L
        x = k - y * L
        conf0[x][y] = S0[k]
        conf1[x][y] = S1[k]
    pylab.subplot(1, 2, 1)
    pylab.imshow(conf0, extent=[0, L, 0, L], interpolation='nearest')
    pylab.title('S0 ' + label)
    pylab.subplot(1, 2, 2)
    pylab.imshow(conf1, extent=[0, L, 0, L], interpolation='nearest')
    pylab.title('S1 ' + label)
    pylab.tight_layout()
    pylab.savefig('plot_' + label + '.png')
    pylab.close()