def plot_variable_importance(feature_importance, names_cols, save_name, save):
    """Show Variable importance graph."""    

    # scale by max importance first 20 variables in column names
    feature_importance = feature_importance / feature_importance.max()
    sorted_idx = np.argsort(feature_importance)[::-1][:20]
    barPos = np.arange(sorted_idx.shape[0]) + .8
    barPos = barPos[::-1]
    
    #plot.figure(num=None, facecolor='w', edgecolor='r') 
    plot.figure(num=None, facecolor='w') 
    plot.barh(barPos, feature_importance[sorted_idx]*100, align='center')
    plot.yticks(barPos, names_cols[sorted_idx])
    plot.xticks(np.arange(0, 120, 20), \
      ['0 %', '20 %', '40 %', '60 %', '80 %', '100 %'])    
    plot.margins(0.02)
    plot.subplots_adjust(bottom=0.15)
    
    plot.title('Variable Importance')
    
    if save:
	plot.savefig(save_name, bbox_inches='tight', dpi = 300)
	plot.close("all")
    else:
	plot.show()    

def command(args):
    from pylab import bar, yticks, subplots_adjust, show
    from numpy import arange

    import sr.tools.bom.bom as bom
    import sr.tools.bom.parts_db as parts_db

    db = parts_db.get_db()
    m = bom.MultiBoardBom(db)
    m.load_boards_args(args.arg)
    m.prime_cache()

    prices = []

    for srcode, pg in m.items():
        if srcode == "sr-nothing":
            continue

        prices.append((srcode, pg.get_price()))

    prices.sort(key=lambda x: x[1])

    bar(0, 0.8, bottom=range(0, len(prices)), width=[x[1] for x in prices],
        orientation='horizontal')

    yticks(arange(0, len(prices)) + 0.4, [x[0] for x in prices])

    subplots_adjust(left=0.35)

    show()

def savePlotSequence(fileBase,stride=1,imgRatio=(5,7),title=None,titleFrames=20,lastFrames=30):
	'''Save sequence of plots, each plot corresponding to one line in history. It is especially meant to be used for :yref:`yade.utils.makeVideo`.

	:param stride: only consider every stride-th line of history (default creates one frame per each line)
	:param title: Create title frame, where lines of title are separated with newlines (``\\n``) and optional subtitle is separated from title by double newline. 
	:param int titleFrames: Create this number of frames with title (by repeating its filename), determines how long the title will stand in the movie.
	:param int lastFrames: Repeat the last frame this number of times, so that the movie does not end abruptly.
	:return: List of filenames with consecutive frames.
	'''
	createPlots(subPlots=True,scatterSize=60,wider=True)
	sqrtFigs=math.sqrt(len(plots))
	pylab.gcf().set_size_inches(8*sqrtFigs,5*sqrtFigs) # better readable
	pylab.subplots_adjust(left=.05,right=.95,bottom=.05,top=.95) # make it more compact
	if len(plots)==1 and plots[plots.keys()[0]]==None: # only pure snapshot is there
		pylab.gcf().set_size_inches(5,5)
		pylab.subplots_adjust(left=0,right=1,bottom=0,top=1)
	#if not data.keys(): raise ValueError("plot.data is empty.")
	pltLen=max(len(data[data.keys()[0]]) if data else 0,len(imgData[imgData.keys()[0]]) if imgData else 0)
	if pltLen==0: raise ValueError("Both plot.data and plot.imgData are empty.")
	global current, currLineRefs
	ret=[]
	print 'Saving %d plot frames, it can take a while...'%(pltLen)
	for i,n in enumerate(range(0,pltLen,stride)):
		current=n
		for l in currLineRefs: l.update()
		out=fileBase+'-%03d.png'%i
		pylab.gcf().savefig(out)
		ret.append(out)
	if len(ret)==0: raise RuntimeError("No images created?!")
	if title:
		titleImgName=fileBase+'-title.png'
		createTitleFrame(titleImgName,Image.open(ret[-1]).size,title)
		ret=titleFrames*[titleImgName]+ret
	if lastFrames>1: ret+=(lastFrames-1)*[ret[-1]]
	return ret

def plot_vector_diff_F160W():
    """
    Using the xym2mat analysis on the F160W filter in the 2010 data set,
    plot up the positional offset vectors for each reference star in
    each star list relative to the average list. This should show
    us if there are systematic problems with the distortion solution
    or PSF variations.
    """
    x, y, m, xe, ye, me, cnt = load_catalog()

    dx = x - x[0]
    dy = y - y[0]
    dr = np.hypot(dx, dy)

    py.clf()
    py.subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95)

    for ii in range(1, x.shape[0]):
        idx = np.where((x[ii,:] != 0) & (y[ii,:] != 0) & (dr[ii,:] < 0.1) &
                       (xe < 0.05) & (ye < 0.05))[0]

        py.clf()
        q = py.quiver(x[0,idx], y[0,idx], dx[ii, idx], dy[ii, idx], scale=1.0)
        py.quiverkey(q, 0.9, 0.9, 0.03333, label='2 mas', color='red')
        py.title('{0} stars in list {1}'.format(len(idx), ii))
        py.xlim(0, 4500)
        py.ylim(0, 4500)

        foo = input('Continue?')
        if foo == 'q' or foo == 'Q':
            break

    def draw(self):

        print self.edgeno

        pos = 0
        dy = 8
        edgeno = self.edgeno
        edge = self.edges[edgeno]
        edgeprev = self.edges[edgeno-1]
        p = np.round(edge["top"](1024))
        top = min(p+2*dy, 2048)
        bot = min(p-2*dy, 2048)
        self.cutout = self.flat[1][bot:top,:].copy()

        pl.figure(1)
        pl.clf()
        start = 0
        dy = 512
        for i in xrange(2048/dy):
            pl.subplot(2048/dy,1,i+1)
            pl.xlim(start, start+dy)

            if i == 0: pl.title("edge %i] %s|%s" % (edgeno,
                edgeprev["Target_Name"], edge["Target_Name"]))


            pl.subplots_adjust(left=.07,right=.99,bottom=.05,top=.95)

            pl.imshow(self.flat[1][bot:top,start:start+dy], extent=(start,
                start+dy, bot, top), cmap='Greys', vmin=2000, vmax=6000)

            pix = np.arange(start, start+dy)
            pl.plot(pix, edge["top"](pix), 'r', linewidth=1)
            pl.plot(pix, edgeprev["bottom"](pix), 'r', linewidth=1)
            pl.plot(edge["xposs_top"], edge["yposs_top"], 'o')
            pl.plot(edgeprev["xposs_bot"], edgeprev["yposs_bot"], 'o')


            hpp = edge["hpps"]
            pl.axvline(hpp[0],ymax=.5, color='blue', linewidth=5)
            pl.axvline(hpp[1],ymax=.5, color='red', linewidth=5)

            hpp = edgeprev["hpps"]
            pl.axvline(hpp[0],ymin=.5,color='blue', linewidth=5)
            pl.axvline(hpp[1],ymin=.5,color='red', linewidth=5)


            if False:
                L = top-bot
                Lx = len(edge["xposs"])
                for i in xrange(Lx):
                    xp = edge["xposs"][i]
                    frac1 = (edge["top"](xp)-bot-1)/L
                    pl.axvline(xp,ymin=frac1)

                for xp in edgeprev["xposs"]: 
                    frac2 = (edgeprev["bottom"](xp)-bot)/L
                    pl.axvline(xp,ymax=frac2)

            start += dy

    def plot_bases(self, autoscale=True, stampsize=None):
        import pylab as plt

        N = len(self.psfbases)
        cols = int(np.ceil(np.sqrt(N)))
        rows = int(np.ceil(N / float(cols)))
        plt.clf()
        plt.subplots_adjust(hspace=0, wspace=0)

        cut = 0
        if stampsize is not None:
            H, W = self.shape
            assert H == W
            cut = max(0, (H - stampsize) / 2)

        ima = dict(interpolation="nearest", origin="lower")
        if autoscale:
            mx = self.psfbases.max()
            ima.update(vmin=-mx, vmax=mx)
        nil, xpows, ypows = self.polynomials(0.0, 0.0, powers=True)
        for i, (xp, yp, b) in enumerate(zip(xpows, ypows, self.psfbases)):
            plt.subplot(rows, cols, i + 1)

            if cut > 0:
                b = b[cut:-cut, cut:-cut]
            if autoscale:
                plt.imshow(b, **ima)
            else:
                mx = np.abs(b).max()
                plt.imshow(b, vmin=-mx, vmax=mx, **ima)
            plt.xticks([])
            plt.yticks([])
            plt.title("x^%i y^%i" % (xp, yp))
        plt.suptitle("PsfEx eigen-bases")

def walker_plot(samples, nwalkers, limit):
    s = samples.reshape(nwalkers, -1, 4)
    s = s[:, :limit, :]
    fig = P.figure(figsize=(8, 10))
    ax1 = P.subplot(4, 1, 1)
    ax2 = P.subplot(4, 1, 2)
    ax3 = P.subplot(4, 1, 3)
    ax4 = P.subplot(4, 1, 4)
    for n in range(len(s)):
        ax1.plot(s[n, :, 0], "k")
        ax2.plot(s[n, :, 1], "k")
        ax3.plot(s[n, :, 2], "k")
        ax4.plot(s[n, :, 3], "k")
    ax1.tick_params(axis="x", labelbottom="off")
    ax2.tick_params(axis="x", labelbottom="off")
    ax3.tick_params(axis="x", labelbottom="off")
    ax4.set_xlabel(r"step number")
    ax1.set_ylabel(r"$t_{smooth}$")
    ax2.set_ylabel(r"$\tau_{smooth}$")
    ax3.set_ylabel(r"$t_{disc}$")
    ax4.set_ylabel(r"$\tau_{disc}$")
    P.subplots_adjust(hspace=0.1)
    save_fig = (
        "/Users/becky/Projects/Green-Valley-Project/bayesian/find_t_tau/gv/not_clean/walkers_steps_all_"
        + str(time.strftime("%H_%M_%d_%m_%y"))
        + ".pdf"
    )
    fig.savefig(save_fig)
    return fig

def write_rgb():
    #g,r,z = [fitsio.read('detmap-%s.fits' % band) for band in 'grz']
    g,r,z = [fitsio.read('coadd-%s.fits' % band) for band in 'grz']

    plt.figure(figsize=(10,10))
    plt.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.95)

    plt.clf()
    for (im1,cc),scale in zip([(g,'b'),(r,'g'),(z,'r')],
                             [2.0, 1.2, 0.4]):
        im = im1 * scale
        im = im[im != 0]
        plt.hist(im.ravel(), histtype='step', color=cc,
                 range=[np.percentile(im, p) for p in (1,98)], bins=50)
    ps.savefig()
        
    #rgb = get_rgb_image(g,r,z, alpha=0.8, m=0.02)
    #rgb = get_rgb_image(g,r,z, alpha=16., m=0.005, m2=0.002,
    #rgb = get_rgb_image(g,r,z, alpha=32., m=0.01, m2=0.002,
    rgb = get_rgb_image(g,r,z, alpha=8., m=0.0, m2=0.0,
        scale_g = 2.,
        scale_r = 1.1,
        scale_z = 0.5,
        Q = 10)


    #for im in g,r,z:
    #    mn,mx = [np.percentile(im, p) for p in [20,99]]
    #    print 'mn,mx:', mn,mx
    
    plt.clf()
    plt.imshow(rgb, interpolation='nearest', origin='lower')
    ps.savefig()

    fitsio.write('rgb.fits', rgb)

def corner_plot(s, labels):
    x, y = s[:,0], s[:,1]
    fig = P.figure(figsize=(10,10))
    ax2 = P.subplot(223)
    ax2.set_xlabel(labels[0])
    ax2.set_ylabel(labels[1])
    im = triangle.histo2d(x, y, ax=ax2, extent=[[0, 13.807108309208775],[0, 3.0]])
    [l.set_rotation(45) for l in ax2.get_xticklabels()]
    [j.set_rotation(45) for j in ax2.get_yticklabels()]
    ax1 = P.subplot(221, xlim=[0, 13.807108309208775])
    ax1.tick_params(axis='x', labelbottom='off')
    ax1.tick_params(axis='y', labelleft='off')
    ax1.hist(x, bins=50, histtype='step', color='k', range=(0, 13.807108309208775))
    ax3 = P.subplot(224)
    ax3.tick_params(axis='x', labelbottom='off')
    ax3.tick_params(axis='y', labelleft='off')
    ax3.hist(y, bins=50, orientation='horizontal', histtype='step',color='k', range=(0,3))
    P.subplots_adjust(wspace=0.05)
    P.subplots_adjust(hspace=0.05)
    cbar_ax = fig.add_axes([0.55, 0.565, 0.02, 0.405])
    cb = fig.colorbar(im, cax = cbar_ax)
    cb.solids.set_edgecolor('face')
    cb.set_label(r'predicted SFR $[M_{\odot} yr^{-1}]$', labelpad = 20, fontsize=16)
    P.tight_layout()
    return fig

def plot_stable_features(X_train,y_train,featnames,**kwargs):
    from sklearn.linear_model import LassoLarsCV,RandomizedLasso

    n_resampling = kwargs.pop('n_resampling',200)
    n_jobs = kwargs.pop('n_jobs',-1)
    
    with warnings.catch_warnings():
        warnings.simplefilter('ignore', UserWarning)
        # estimate alphas via xvalidation 
        lars_cv = LassoLarsCV(cv=6,n_jobs=n_jobs).fit(X_train,y_train)        
        alphas = np.linspace(lars_cv.alphas_[0], .1 * lars_cv.alphas_[0], 6)

        clf = RandomizedLasso(alpha=alphas, random_state=42, n_jobs=n_jobs,
                              n_resampling=n_resampling)
        clf.fit(X_train,y_train)
        importances = clf.scores_ 
        indices = np.argsort(importances)[::-1]

        pl.bar(range(len(featnames)), importances[indices],
               color="r", align="center")
        pl.xticks(np.arange(len(featnames))+0.5,featnames[indices],
                  rotation=45,horizontalalignment='right')
        pl.xlim(-0.5,len(featnames)-0.5)
        pl.subplots_adjust(bottom=0.2)
        
        pl.ylim(0,np.max(importances)*1.01)
        pl.ylabel('Selection frequency (%) for %d resamplings '%n_resampling)
        pl.title("Stability Selection: Selection Frequencies")

def plot_importances(clf,featnames,outfile,**kwargs):

    pl.figure(figsize=(16,4))

    featnames = np.array(featnames)
    importances = clf.feature_importances_
    imp_std = np.std([tree.feature_importances_ for tree in clf.estimators_],
                     axis=0)
    indices = np.argsort(importances)[::-1]

    #for featname in featnames[indices]:
    #    print featname

    trunc_featnames = featnames[indices]
    trunc_featnames = trunc_featnames[0:24]
    trunc_importances = importances[indices]
    trunc_importances = trunc_importances[0:24]
    trunc_imp_std = imp_std[indices]
    trunc_imp_std = trunc_imp_std[0:24]

    pl.bar(range(len(trunc_featnames)), trunc_importances,
           color="r", yerr=trunc_imp_std, align="center")
    pl.xticks(np.arange(len(trunc_featnames))+0.5,trunc_featnames,rotation=45,
              horizontalalignment='right')
    pl.xlim(-0.5,len(trunc_featnames)-0.5)
    pl.ylim(0,np.max(trunc_importances+trunc_imp_std)*1.01)

#    pl.bar(range(len(featnames)), importances[indices],
#           color="r", yerr=imp_std[indices], align="center")
#    pl.xticks(np.arange(len(featnames))+0.5,featnames[indices],rotation=45,
#              horizontalalignment='right')
#    pl.xlim(-0.5,len(featnames)-0.5)
#    pl.ylim(0,np.max(importances+imp_std)*1.01)
    pl.subplots_adjust(bottom=0.2)
    pl.show()

def plot_cluster_context(sizes, densities, dir, name=None, k=None, suffix="png"):
    """
    so many conditionals!
    """
    print("plot_cluster_context(): plotting", name)
    if name is None:
        K = len(sizes)
        fn = "{}/clusters_{:04d}.{}".format(dir, K, suffix)
    else:
        fn = "{}/{}_context.{}".format(dir, name, suffix)
    if os.path.exists(fn):
        print("plot_cluster_context(): {} exists already".format(fn))
        return
    if k is None:
        fig = plt.figure(figsize=(6,6))
        plt.subplots_adjust(left=0.15, right=0.97, bottom=0.15, top=0.97)
        ms = 7.5
    else:
        fig = plt.figure(figsize=(4,4))
        plt.subplots_adjust(left=0.2, right=0.96, bottom=0.2, top=0.96)
        ms = 5.0
    plt.clf()
    if name is not None and k is None:
        plt.savefig(fn)
        print("plot_cluster_context(): wrote", fn)
        return
    _clusterplot(sizes, densities, k, ms=ms)
    _clusterlims(sizes, densities)
    plt.ylabel("cluster abundance-space density")
    plt.xlabel("number in abundance-space cluster")
    plt.loglog()
    [l.set_rotation(45) for l in plt.gca().get_xticklabels()]
    [l.set_rotation(45) for l in plt.gca().get_yticklabels()]
    plt.savefig(fn)
    print("plot_cluster_context(): wrote", fn)

 def create_fig(self):
     print "Creating fig..."
     self.fig_size = (14, 10)
     self.fig = pylab.figure(figsize=self.fig_size)
     pylab.subplots_adjust(hspace=0.4)
     pylab.subplots_adjust(wspace=0.35)
     return self.fig

    def drawCity(self):
        """
        作图
        :return:
        """
        pl.title("pm25 / time   " + str(self.numMonitors) + "_monitors")# give plot a title
        pl.xlabel('time')# make axis labels
        pl.ylabel('pm2.5')
        self.fill_cityPm25List()


        for monitorStr in self.cityPm25List:
            data = np.loadtxt(StringIO(monitorStr), dtype=np.dtype([("t", "S13"),("v", float)]))
            datestr = np.char.replace(data["t"], "T", " ")
            t = pl.datestr2num(datestr)
            v = data["v"]
            pl.plot_date(t, v, fmt="-o")



        pl.subplots_adjust(bottom=0.3)
        # pl.legend(loc=4)#指定legend的位置,读者可以自己help它的用法
        ax = pl.gca()
        ax.fmt_xdata = pl.DateFormatter('%Y-%m-%d %H:%M:%S')
        pl.xticks(rotation=70)
        # pl.xticks(t, datestr) # 如果以数据点为刻度，则注释掉这一行
        ax.xaxis.set_major_formatter(pl.DateFormatter('%Y-%m-%d %H:%M'))
        pl.grid()
        pl.show()# show the plot on the screen

    def gui_repr(self):
        """Generate a GUI to represent the sentence alignments
        """
        if __pylab_loaded__:
            fig_width = max(len(self.text_e), len(self.text_f)) + 1
            fig_height = 3
            pylab.figure(figsize=(fig_width*0.8, fig_height*0.8), facecolor='w')
            pylab.box(on=False)
            pylab.subplots_adjust(left=0, right=1, bottom=0, top=1)
            pylab.xlim(-1, fig_width - 1)
            pylab.ylim(0, fig_height)
            pylab.xticks([])
            pylab.yticks([])

            e = [0 for _ in xrange(len(self.text_e))]
            f = [0 for _ in xrange(len(self.text_f))]
            for (i, j) in self.align:
                e[i] = 1
                f[j] = 1
                # draw the middle line
                pylab.arrow(i, 2, j - i, -1, color='r')
            for i in xrange(len(e)):
                # draw e side line
                pylab.text(i, 2.5, self.text_e[i], ha = 'center', va = 'center',
                        rotation=30)
                if e[i] == 1:
                    pylab.arrow(i, 2.5, 0, -0.5, color='r', alpha=0.3, lw=2)
            for i in xrange(len(f)):
                # draw f side line
                 pylab.text(i, 0.5, self.text_f[i], ha = 'center', va = 'center',
                        rotation=30)
                 if f[i] == 1:
                    pylab.arrow(i, 0.5, 0, 0.5, color='r', alpha=0.3, lw=2)

            pylab.draw()

def clusterSample():
    clusters = [
        {"name": "M17", "distance": 2, "ra": "18:20:26", "dec": "-17:10:01"},
        {"name": "Wd2", "distance": 5, "ra": "10:23:58", "dec": "-57:45:49"},
        {"name": "Wd1", "distance": 5, "ra": "16:47:04", "dec": "-45:51:05"},
        {"name": "RSGC1", "distance": 6, "ra": "18:37:58", "dec": "-06:52:53"},
        {"name": "RSGC2", "distance": 6, "ra": "18:39:20", "dec": "-06:05:10"},
    ]

    py.close(1)
    py.figure(1, linewidth=2, figsize=(16, 10))
    py.subplots_adjust(left=0.05, right=0.97, bottom=0.1, top=0.95, wspace=0.2, hspace=0.25)

    for ii in range(len(clusters)):
        clust = clusters[ii]

        obj = ephem.FixedBody()
        obj._ra = ephem.hours(clust["ra"])
        obj._dec = ephem.degrees(clust["dec"])
        obj._epoch = 2000
        obj.compute()

        gal = ephem.Galactic(obj)

        longitude = math.degrees(float(gal.lon))
        print ""
        print "%-10s  %s %s  l = %.1f" % (clust["name"], clust["ra"], clust["dec"], longitude)

        py.subplot(2, 3, ii + 1)
        properMotions(longitude, clust["distance"], clust["name"])

def animation_compare(filename1, filename2, prefix, tend):
    num = 0
    for t in timestep(0,tend):
        t1,s1 = FieldInitializer.LoadState(filename1, t)
        t2,s2 = FieldInitializer.LoadState(filename2, t)
        rot1 = s1.CalculateRotationRodrigues()['z']
        rot2 = s2.CalculateRotationRodrigues()['z']

        pylab.figure(figsize=(5.12*2,6.12))
        rho = s2.CalculateRhoFourier().modulus()
        rot = s2.CalculateRotationRodrigues()['z']
        pylab.subplot(121)
        #pylab.imshow(rho*(rho>0.5), alpha=1, cmap=pylab.cm.bone_r)
        pylab.imshow(rot1)
        pylab.xticks([])
        pylab.yticks([])
        pylab.xlabel(r'$d_0$', fontsize=25)
        pylab.subplot(122)
        pylab.imshow(rot2)
        pylab.xticks([])
        pylab.yticks([])
        pylab.subplots_adjust(0.,0.,1.,1.,0.01,0.05)
        pylab.xlabel(r'$d_0/2$', fontsize=25)
        pylab.suptitle("Nabarro-Herring", fontsize=25)
        pylab.savefig("%s%04i.png" %(prefix, num))
        pylab.close('all')
        num = num + 1

def plot_data(yRange=None):
    '''
    Plots and saves the cell measurement data.  Returns nothing.
    '''
    fig = plt.figure(figsize=(18,12))
    ax = plt.subplot(111)
    plt.errorbar(range(len(avgCells.index)), avgCells[column], yerr=stdCells[column], fmt='o')
    ax = plt.gca()
    ax.set(xticks=range(len(avgCells.index)), xticklabels=avgCells.index)
    xlims = ax.get_xlim()
    ax.set_xlim([lim-1 for lim in xlims])
    # adjust yRange if it was specified
    if yRange!=None:
        ax.set_ylim(yRange)
        fileName = column + ' exlcuding outliers'
    else:
        fileName = column
    plt.subplots_adjust(bottom=0.2, right=0.98, left=0.05)
    plt.title(column)
    plt.ylabel('mm')
    locs, labels = plt.xticks()
    plt.setp(labels, rotation=90)
    mng = plt.get_current_fig_manager()
    mng.window.state('zoomed')
    #plt.show()
    path1 = 'Y:/Test data/ACT02/vision inspection/plot_100_cells/'
    path2 = 'Y:/Nate/git/nuvosun-python-lib/vision system/plot_100_cells/'
    fig.savefig(path1 + fileName, bbox_inches = 'tight')
    fig.savefig(path2 + fileName, bbox_inches = 'tight')
    plt.close()

def plot_C_gamma_grid_search(grid, C_range, gamma_range, score):
    '''
    Plots the scores computed on a grid. 
    
    Arguments: 
        grid - the grid search object created using GridSearchCV()
        C_range - the C parameter range 
        gamma_range - the gamma parameter range 
        score - the scoring function  
        
    
    '''

    # grid_scores_ contains parameter settings and scores
    # We extract just the scores
    scores = [x[1] for x in grid.grid_scores_]
    scores = np.array(scores).reshape(len(C_range), len(gamma_range))
    
    # draw heatmap of accuracy as a function of gamma and C
    pl.figure(figsize=(8, 6))
    pl.subplots_adjust(left=0.05, right=0.95, bottom=0.15, top=0.95)
    pl.imshow(scores, interpolation='nearest', cmap=pl.cm.spectral)
    pl.title("Grid search on C and gamma for best %s" % score)
    pl.xlabel('gamma')
    pl.ylabel('C')
    pl.colorbar()
    pl.xticks(np.arange(len(gamma_range)), gamma_range, rotation=45)
    pl.yticks(np.arange(len(C_range)), C_range)
    
    pl.show()

def walker_plot(samples, nwalkers, limit):
    s = samples.reshape(nwalkers, -1, 4)
    s = s[:,:limit, :]
    fig = P.figure(figsize=(8,10))
    ax1 = P.subplot(4,1,1)
    ax2 = P.subplot(4,1,2)
    ax3 = P.subplot(4,1,3)
    ax4 = P.subplot(4,1,4)
    for n in range(len(s)):
        ax1.plot(s[n,:,0], 'k')
        ax2.plot(s[n,:,1], 'k')
        ax3.plot(s[n,:,2], 'k')
        ax4.plot(s[n,:,3], 'k')
    ax1.tick_params(axis='x', labelbottom='off')
    ax2.tick_params(axis='x', labelbottom='off')
    ax3.tick_params(axis='x', labelbottom='off')
    ax4.set_xlabel(r'step number')
    ax1.set_ylabel(r'$t_{smooth}$')
    ax2.set_ylabel(r'$\tau_{smooth}$')
    ax3.set_ylabel(r'$t_{disc}$')
    ax4.set_ylabel(r'$\tau_{disc}$')
    P.subplots_adjust(hspace=0.1)
    save_fig = '/Users/becky/Projects/Green-Valley-Project/bayesian/find_t_tau/walkers_steps_red_s_'+str(time.strftime('%H_%M_%d_%m_%y'))+'.pdf'
    fig.savefig(save_fig)
    return fig