def _decorate_histogram(vstats):
    import pylab
    from matplotlib.transforms import blended_transform_factory as blend
    # Shade things inside 1-sigma
    pylab.axvspan(vstats.p68[0],vstats.p68[1],
                  color='gold',alpha=0.5,zorder=-1)
    # build transform with x=data, y=axes(0,1)
    ax = pylab.gca()
    transform = blend(ax.transData, ax.transAxes)

    l95,h95 = vstats.p95
    l68,h68 = vstats.p68
    def marker(s,v):
        if v < l95: s,v,ha = '<'+s,l95,'left'
        elif v > h95: s,v,ha = '>'+s,h95,'right'
        else: ha='center'
        pylab.text(v, 0.95, s, va='top', ha=ha,
                   transform=transform, zorder=3, color='g')
        #pylab.axvline(v)
    marker('|',vstats.median)
    marker('E',vstats.mean)
    marker('*',vstats.best)

    pylab.text(0.01, 0.95, vstats.label, zorder=2,
        backgroundcolor=(1,1,0,0.2),
        verticalalignment='top',
        horizontalalignment='left',
        transform=pylab.gca().transAxes)
    pylab.setp([pylab.gca().get_yticklabels()],visible=False)
    ticks = (l95, l68, vstats.median, h68, h95)
    labels = [format_value(v,h95-l95) for v in ticks]
    if len(labels[2]) > 5:
        # Drop 68% values if too many digits
        ticks,labels= ticks[0::2],labels[0::2]
    pylab.xticks(ticks, labels)

def plot_exclude():
    global exclude
    #c = '#FFFFFF'
    c = '#DDDDDD'
    #c = 'g'
    for l,u in exclude:
        pl.axvspan(l,u, fc=c, ec=c)

def BinnedChroPlot(initial, final, binsize=100, CenShading=[CenStart,CenEnd], name='Chromosome'):
    ''' This function creates a graph of initial and final Chromosomes binned per (binsize) '''
    
    Title = '%s usage, binned per %i'%(name,binsize)
    plt.figure(Title)
    plt.title(Title)
    
    global initialHeights   
    
    initialSize = len(initial)/binsize
    finalSize = len(final)/binsize
    
    initialHeights = [sum(initial[x:x+binsize]) for x in xrange(0,len(initial),binsize)] # idem, but binned in bargraph[x:x+binsize]) for x in xrange(0,ChroL,binsize)] #create binned initial Chro
    finalHeights = [sum(final[x:x+binsize]) for x in xrange(0,len(final),binsize)]        #create binned Chro

    if CenShading:
        if CenShading == [CenStart,CenEnd]: CenShading = [CenStart/binsize,CenEnd/binsize]
        plt.axvspan(CenShading[0], CenShading[1], color='0.85')
    plt.bar(xrange(initialSize), initialHeights, lw=0, width=1, color='r', label='initial', alpha=0.4)
    plt.bar(xrange(finalSize), finalHeights, lw=0, width=1, color='b', label='final', alpha=0.4)
    
    plt.ylabel('sites occupied (of %i)'%(binsize))
    plt.xlabel('bin number (%i sites per bin)'%(binsize))
    plt.legend()
    plt.show()

    def plot(self):
        f = pylab.figure(figsize=(8,4))
        co = [] #colors container
        for zScore, r in itertools.izip(self.zScores, self.log2Ratio):
            if zScore < self.pCut:
                if r > 0:
                    co.append(Colors().greenColor)
                elif r < 0:
                    co.append(Colors().redColor)
                else:
                    raise Exception
            else:
                co.append(Colors().blueColor)

        #print "Probability this is from a normal distribution: %.3e" %stats.normaltest(self.log2Ratio)[1]
        ax = f.add_subplot(121)
        pylab.axvline(self.meanLog2Ratio, color=Colors().redColor)
        pylab.axvspan(self.meanLog2Ratio-(2*self.stdLog2Ratio), 
                      self.meanLog2Ratio+(2*self.stdLog2Ratio), color=Colors().blueColor, alpha=0.2)
        his = pylab.hist(self.log2Ratio, bins=50, color=Colors().blueColor)
        pylab.xlabel("log2 Ratio %s/%s" %(self.sampleNames[1], self.sampleNames[0]))
        pylab.ylabel("Frequency")
        
        ax = f.add_subplot(122, aspect='equal')
        pylab.scatter(self.genes1, self.genes2, c=co, alpha=0.5)        
        pylab.ylabel("%s RPKM" %self.sampleNames[1])
        pylab.xlabel("%s RPKM" %self.sampleNames[0])
        pylab.yscale('log')
        pylab.xscale('log')
        pylab.tight_layout()

def annotated_plot(v,normfiltfunc,a,e):#plots a the Hotspot along with the annotated regions from the posterior decoding
	import pylab
	v2 = normfiltfunc(v)
	(s,f,b,pdf_m) = posterior_step(v2,a,e)
	post = s*f*b
	maxstates = post.argmax(axis=0)
	#label the footprints
	foots = (maxstates == 3)*1
	bins = diff(foots)
	start = where(bins == 1)[0] + 1
	stop = where(bins == -1)[0]
	for p,q in zip(start,stop):
		foot = pylab.axvspan(p, q, facecolor='r', alpha=0.5)
	#label the HS1
	hs1s = (maxstates == 0)*1
	bins = diff(hs1s)
	start = concatenate(([0],where(bins == 1)[0] + 1),1)#the first state is hs1. this accounts for that
	stop = where(bins == -1)[0]
	for p,q in zip(start,stop):
		hs1 = pylab.axvspan(p, q, facecolor='g', alpha=0.5)
	#label the HS2
	hs2s = (maxstates == 4)*1
	bins = diff(hs2s)
	start = where(bins == 1)[0] + 1
	stop = concatenate((where(bins == -1)[0],[len(v)-1]),1)#the last state is hs2
	for p,q in zip(start,stop):
		hs2 = pylab.axvspan(p, q, facecolor='c', alpha=0.5)
	pylab.plot(v)
	pylab.legend((hs1,foot,hs2,),('HS1','Footprint','HS2',))
	pylab.xlabel('DHS Coordinates')
	pylab.ylabel('DNase I Cuts')

def plot_likelihoods( p, name, score = None ):
    from pylab import plot, show, cla, clf, legend, figure, xlabel, ylabel, \
            title, axvspan
    old_odds = PssmParameters.singleton().binding_background_odds_prior
    PssmParameters.singleton().binding_background_odds_prior = 1.0
    (
            bind,
            back,
            cum_bind,
            cum_back,
            odds_ratio,
            cum_odds_ratio,
            p_bind,
            cum_p_bind,
            p_value_p_bind
    ) = get_pssm_likelihoods( p )
    scores = [ float(i) / (len(bind) - 1.0) for i in range( len( bind ) ) ]
    # cla()
    # clf()
    figure()
    plot( scores, bind, 'g-', label='binding' )
    plot( scores, back, 'b-', label='background' )
    plot( scores, cum_bind, 'g--', label='binding (cumulative)' )
    plot( scores, cum_back, 'b--', label='background (cumulative)' )
    plot( scores, p_bind, 'r-', label='p(binding)' )
    plot( scores, cum_p_bind, 'r--', label='p(binding) (cumulative)' )
    plot( scores, p_value_p_bind, 'y--', label='p(binding) (p-value)' )
    # legend( loc='center left' )
    xlabel( 'score' )
    title( name )
    if score:
        idx = get_likelihood_index( len( bind ), score )
        axvspan( score, score )
    show()
    PssmParameters.singleton().binding_background_odds_prior = old_odds

 def search_clusters(self, data):
     """This function performs the search for significant clusters. 
        Uses a class from eegpy.stats.cluster
     """
     plotid=0
     for i_ch in range(data[0].shape[2]):
         for i_b in range(data[0].shape[3]):
             #print i_ch,i_b
             cs_data = [d[:,:,i_ch,i_b].T for d in data]
             #print len(cs_data), [csd.shape for csd in cs_data]
             fs,sct,scp,ct,cp = ClusterSearch1d(cs_data,num_surrogates=self._num_surrogates).search()
             #print "CPM: fs.shape=",fs.shape
             #print ct, cp
             for i_c in range(len(sct)):
                 self._clusters.append( (i_ch,i_b,sct[i_c],scp[i_c]) ) 
             #Do plot if significant cluster found
             if len(sct)>0:
                 #print "CPM: fs.shape=",fs.shape, fs.dtype
                 #print "CPM: fs[499:510]", fs[498:510]
                 p.plot(np.array(fs))
                 p.title("Channel %i, band %i"%(i_ch,i_b))
                 for cluster in sct:
                     p.axvspan(cluster[0],cluster[1],color="y",alpha=0.2)
                 p.savefig("/tmp/fbpm%03d.png"%plotid)
                 plotid+=1
                 p.clf()

def draw_search_graph(plots):
    import pylab as pl

    fg = pl.figure()
    ax = fg.add_subplot(111)

    ax.yaxis.set_major_formatter(pl.FuncFormatter(labelfmt))

    for (values, attrs) in plots:
        indexes, width = pl.arange(len(values)), 1.0 / len(plots)

        yvalues = [x.result for x in values]
        xoffset = width * plots.index((values, attrs))
        ax.plot(indexes + xoffset, yvalues, **attrs)

        legend = ax.legend(loc='best')
        legend.get_frame().set_alpha(0.6)
        fg.canvas.draw()

    pl.ylabel('tradeoff improvement -->')
    pl.xlabel('number of tested configurations -->')
    pl.title('Search Graph')
    pl.axhspan(0.0, 0.0)
    pl.axvspan(0.0, 0.0)
    pl.grid(True)
    pl.show()

def draw_onsets(onsets):
    if not onsets:
        return
    
    # Draw the onsets
    for onsetLeft, onsetCenter, onsetRight in onsets:
        pylab.axvspan( xmin = onsetLeft, xmax = onsetRight, facecolor = 'green', linewidth = 0, alpha = 0.25)
        pylab.axvline( x = onsetCenter, color = 'black', linewidth = 1.1)

def plotRes(data, errors, r):
    import pylab
    pylab.figure()

    nObs = len(data)

    n, bins, patches = pylab.hist(data, 2*np.sqrt(nObs), fc=[.7,.7,.7])

    binSize = bins[1] - bins[0]
    x = np.arange(bins[0], bins[-1])
    

    means, sigs, pis, mVars, weights = r

    inds = np.argmax(weights, 1)

    for i in range(means.size):
        #print i
        c = pylab.cm.hsv(float(i)/means.size)
        n, bin_s, patches = pylab.hist(data[inds == i], bins, alpha=0.3, facecolor=c)
    
    ys = np.zeros_like(x)

    i = 0
    for m, s, p in zip(means, sigs, pis):
        c = pylab.cm.hsv(float(i)/means.size)
        y = nObs*p*binSize*np.exp(-(x-m)**2/(2*s**2))/np.sqrt(2*np.pi*s**2)
        ys += y

        i+= 1

        pylab.plot(x,y, lw=2, color=c)

    #pylab.plot(x, ys, lw=3)

    pylab.figure()

    ci = (r[4]*np.arange(r[0].size)[None,:]).sum(1)

    I = np.argsort(ci)
    cis = ci[I]
    cil = 0

    for i in range(means.size):
        c = pylab.cm.hsv(float(i)/means.size)

        print(c)

        pylab.axvline(means[i], color=c)
        pylab.axvspan(means[i] - sigs[i], means[i] + sigs[i], alpha=0.5, facecolor=c)

        cin = cis.searchsorted(i+0.5)

        pylab.axhspan(cil, cin, alpha=0.3, facecolor=c)

        cil = cin

    pylab.errorbar(data[I], np.arange(data.size), xerr=errors[I], fmt='.')

def plot_tree(T, res=None, title=None, cmap_id="Pastel2"):
    """Plots a given tree, containing hierarchical segmentation.

    Parameters
    ----------
    T: mir_eval.segment.tree
        A tree object containing the hierarchical segmentation.
    res: float
        Frame-rate resolution of the tree (None to use seconds).
    title: str
        Title for the plot. `None` for no title.
    cmap_id: str
        Color Map ID
    """
    def round_time(t, res=0.1):
        v = int(t / float(res)) * res
        return v

    # Get color map
    cmap = plt.get_cmap(cmap_id)

    # Get segments by level
    level_bounds = []
    for level in T.levels:
        if level == "root":
            continue
        segments = T.get_segments_in_level(level)
        level_bounds.append(segments)

    # Plot axvspans for each segment
    B = float(len(level_bounds))
    #plt.figure(figsize=figsize)
    for i, segments in enumerate(level_bounds):
        labels = utils.segment_labels_to_floats(segments)
        for segment, label in zip(segments, labels):
            #print i, label, cmap(label)
            if res is None:
                start = segment.start
                end = segment.end
                xlabel = "Time (seconds)"
            else:
                start = int(round_time(segment.start, res=res) / res)
                end = int(round_time(segment.end, res=res) / res)
                xlabel = "Time (frames)"
            plt.axvspan(start, end,
                        ymax=(len(level_bounds) - i) / B,
                        ymin=(len(level_bounds) - i - 1) / B,
                        facecolor=cmap(label))

    # Plot labels
    L = float(len(T.levels) - 1)
    plt.yticks(np.linspace(0, (L - 1) / L, num=L) + 1 / L / 2.,
               T.levels[1:][::-1])
    plt.xlabel(xlabel)
    if title is not None:
        plt.title(title)
    plt.gca().set_xlim([0, end])

def plot_labels(all_labels, gt_times, est_file, algo_ids=None, title=None,
                output_file=None):
    """Plots all the labels.

    Parameters
    ----------
    all_labels: list
        A list of np.arrays containing the labels of the boundaries, one array
        for each algorithm.
    gt_times: np.array
        Array with the ground truth boundaries.
    est_file: str
        Path to the estimated file (JSON file)
    algo_ids : list
        List of algorithm ids to to read boundaries from.
        If None, all algorithm ids are read.
    title : str
        Title of the plot. If None, the name of the file is printed instead.
    """
    N = len(all_labels)  # Number of lists of labels
    if algo_ids is None:
        algo_ids = io.get_algo_ids(est_file)

    # Translate ids
    for i, algo_id in enumerate(algo_ids):
        algo_ids[i] = translate_ids[algo_id]
    algo_ids = ["GT"] + algo_ids

    # Index the labels to normalize them
    for i, labels in enumerate(all_labels):
        all_labels[i] = mir_eval.util.index_labels(labels)[0]

    # Get color map
    cm = plt.get_cmap('gist_rainbow')
    max_label = max(max(labels) for labels in all_labels)

    # To intervals
    gt_inters = utils.times_to_intervals(gt_times)

    # Plot labels
    figsize = (6, 4)
    plt.figure(1, figsize=figsize, dpi=120, facecolor='w', edgecolor='k')
    for i, labels in enumerate(all_labels):
        for label, inter in zip(labels, gt_inters):
            plt.axvspan(inter[0], inter[1], ymin=i / float(N),
                        ymax=(i + 1) / float(N), alpha=0.6,
                        color=cm(label / float(max_label)))
        plt.axhline(i / float(N), color="k", linewidth=1)

    # Draw the boundary lines
    for bound in gt_times:
        plt.axvline(bound, color="g")

    # Format plot
    _plot_formatting(title, est_file, algo_ids, gt_times[-1], N,
                     output_file)

def draw_correl_graph(getgraph, opts):
    # http://matplotlib.sourceforge.net/index.html
    fg = pl.figure()
    ax = fg.add_subplot(111)

    bars = getgraph()

    for (values, attrs) in bars:
        indexes, width = pl.arange(len(values)), 1.0 / len(bars)

        yvalues = [x.speedup for x in values]
        xoffset = width * bars.index((values, attrs))
        ax.bar(indexes + xoffset, yvalues, width, picker=4000, **attrs)

        ax.legend(loc='lower left')
        fg.canvas.draw()

    # dynamic annotations
    def on_pick(event):
        ind = int(event.mouseevent.xdata)
        point = bars[0][0][ind]
        tooltip.set_position(
            (event.mouseevent.xdata, event.mouseevent.ydata))
        tooltip.set_text(point_descr(point))
        tooltip.set_visible(True)
        fg.canvas.draw()

    tooltip = ax.text(
        0, 0, "undef", bbox=dict(facecolor='white', alpha=0.8),
        verticalalignment='bottom', visible=False)

    # graph title
    try:
        title = 'Correlation Graph for %s' % (
            opts.id or opts.targets or bars[0][0][0].target)
    except: title = 'Correlation Graph'

    # redraw axis, set labels, legend, grid, ...
    def labelfmt(x, pos=0): return '%.2f%%' % (100.0 * x)
    ax.yaxis.set_major_formatter(pl.FuncFormatter(labelfmt))
    pl.ylabel('speedup (higher is better) -->')

    pl.xlabel('Configurations (ordered by decreasing speedup of '
              + bars[0][1]['label'] + ') -->')
    pl.title(title)
    pl.axhspan(0.0, 0.0)
    pl.axvspan(0.0, 0.0)
    pl.grid(True)

    if opts.outfile:
        fg.savefig(opts.outfile)

    if opts.show:
        fg.canvas.mpl_connect('pick_event', on_pick)
        pl.show()

def plot_one_track(file_struct, est_times, est_labels, boundaries_id, labels_id,
                   title=None):
    """Plots the results of one track, with ground truth if it exists."""
    # Set up the boundaries id
    bid_lid = boundaries_id
    if labels_id is not None:
        bid_lid += " + " + labels_id
    try:
        # Read file
        jam = jams.load(file_struct.ref_file)
        ann = jam.search(namespace='segment_.*')[0]
        ref_inters, ref_labels = ann.data.to_interval_values()

        # To times
        ref_times = utils.intervals_to_times(ref_inters)
        all_boundaries = [ref_times, est_times]
        all_labels = [ref_labels, est_labels]
        algo_ids = ["GT", bid_lid]
    except:
        logging.warning("No references found in %s. Not plotting groundtruth"
                        % file_struct.ref_file)
        all_boundaries = [est_times]
        all_labels = [est_labels]
        algo_ids = [bid_lid]

    N = len(all_boundaries)

    # Index the labels to normalize them
    for i, labels in enumerate(all_labels):
        all_labels[i] = mir_eval.util.index_labels(labels)[0]

    # Get color map
    cm = plt.get_cmap('gist_rainbow')
    max_label = max(max(labels) for labels in all_labels)

    figsize = (8, 4)
    plt.figure(1, figsize=figsize, dpi=120, facecolor='w', edgecolor='k')
    for i, boundaries in enumerate(all_boundaries):
        color = "b"
        if i == 0:
            color = "g"
        for b in boundaries:
            plt.axvline(b, i / float(N), (i + 1) / float(N), color=color)
        if labels_id is not None:
            labels = all_labels[i]
            inters = utils.times_to_intervals(boundaries)
            for label, inter in zip(labels, inters):
                plt.axvspan(inter[0], inter[1], ymin=i / float(N),
                            ymax=(i + 1) / float(N), alpha=0.6,
                            color=cm(label / float(max_label)))
        plt.axhline(i / float(N), color="k", linewidth=1)

    # Format plot
    _plot_formatting(title, os.path.basename(file_struct.audio_file), algo_ids,
                     all_boundaries[0][-1], N, None)

def csv2png(p):
    print p
    title, axis, data = get_data(p)
    dates = data[0]

    release_title, release_axis, release_data = get_data( py.path.local("release_dates.dat") )
    release_dates, release_names = release_data
 
    sprint_title, sprint_axis, sprint_data = get_data( py.path.local("sprint_dates.dat") )
    sprint_locations, sprint_begin_dates, sprint_end_dates = sprint_data
 
    ax = pylab.subplot(111)
    for i, d in enumerate(data[1:]):
        args = [dates, d, colors[i]]
        pylab.plot_date(linewidth=0.8, *args)

    ymax = max(pylab.yticks()[0]) #just below the legend
    for i, release_date in enumerate(release_dates):
        release_name = release_names[i]
        if greyscale:
            color = 0.3
        else:
            color = "g"
        pylab.axvline(release_date, linewidth=0.8, color=color, alpha=0.5)
        ax.text(release_date, ymax * 0.4, release_name,
                fontsize=10,
                horizontalalignment='right',
                verticalalignment='top',
                rotation='vertical')
    for i, location in enumerate(sprint_locations):
        begin = sprint_begin_dates[i]
        end   = sprint_end_dates[i]
        if float(begin) >= float(min(dates[0],dates[-1])):
            if greyscale:
                color = 0.8
            else:
                color = "y"
            pylab.axvspan(begin, end, linewidth=0, facecolor=color, alpha=0.5)
            ax.text(begin, ymax * 0.85, location,
                    fontsize=10,
                    horizontalalignment='right',
                    verticalalignment='top',
                    rotation='vertical')
    pylab.legend(axis[1:], "upper left")
    pylab.ylabel(axis[0])
    pylab.xlabel("")
    ticklabels = ax.get_xticklabels()
    pylab.setp(ticklabels, 'rotation', 45, size=9)
#    ax.autoscale_view()
    ax.grid(True)
    pylab.title(title)

    pylab.savefig(p.purebasename + ".png")
    pylab.savefig(p.purebasename + ".eps")
    py.process.cmdexec("epstopdf %s" % (p.purebasename + ".eps", ))

def plot_stiff_area(bin_f0, B, bin_spread_below, bin_spread_above, stiff_partials, sample_rate):
    stiff_bins = numpy.array(
        [ partials.mode_B2freq( bin_f0, i+1, B)
            for i in xrange(len(stiff_partials)) ]
            )
    for i, est in enumerate(stiff_bins):
        low = stft.bin2hertz(est - bin_spread_below, sample_rate)
        high = stft.bin2hertz(est + bin_spread_above, sample_rate)
        if i == 0:
            pylab.axvspan(low, high, color='c', alpha=0.3,
                label="stiff")
        else:
            pylab.axvspan(low, high, color='c', alpha=0.3)

def format_data(name, data, offset=0):
    currently = 0
    changes = []
    for x in data:
        start = x.tick
        if x.hist.caught_at > 0:
            p.axvspan(
                x.hist.caught_at, x.hist.caught_at + x.duration, alpha=0.25, label="Processing " + name, color="red"
            )
        if currently != x.value.digital:
            changes.append((start, offset + currently))
            changes.append((start, offset + x.value.digital))
            currently = x.value.digital
    return changes

    def plotGait2(self, graphNumber=1, background='gray', lineWidth=10):
        graph = self._figure.add_subplot(self._numberRows, 1, graphNumber)

        self.setBackground(background)
        rect = graph.patch
        rect.set_facecolor(background)

        y_RH, y_RF, y_LF, y_LH = 0.2, 0.4, 0.6, 0.8

        stop  = self._cd.configs.get(GeneralConfigEnum.STOP_TIME)
        start = self._cd.configs.get(GeneralConfigEnum.START_TIME)
        steps = self._cd.configs.get(GeneralConfigEnum.STEPS)
        delta = (stop - start)/steps

        times  = plt.arange(self._plotStartTime, self._plotStopTime, delta)
        valuesRH = self.values(self._channel_RH, times)
        valuesRF = self.values(self._channel_RF, times)
        valuesLF = self.values(self._channel_LF, times)
        valuesLH = self.values(self._channel_LH, times)

        self.setColorMap('RdYlGn')

        for i in range(len(times)):
            time  = times[i]
            lh, lf, rf, rh = valuesLH[i], valuesLF[i], valuesRF[i], valuesRH[i]

            if   lh and lf and rf and rh:                  support = 1.0
            elif (lh and rf and rh) or (lh and lf and rh): support = 0.95
            elif (lh and lf and rf) or (lf and rf and rh): support = 0.85
            elif lh and rh:                                support = 0.75
            elif (lh and rf) or (lf and rh):               support = 0.7
            elif lf and rf:                                support = 0.2
            elif (lh and lf) or (rf and rh):               support = 0.1
            else:                                          support = 0.0
            if support > 0.0:
                plt.axvspan(time - 2.0*delta, time, fc=self.mapValueToColor(support), ec='none')

        self.setColorMap('gray')
        self.plotChannel(self._channel_LH, graph, y_LH, lineWidth)
        self.plotChannel(self._channel_LF, graph, y_LF, lineWidth)
        self.plotChannel(self._channel_RF, graph, y_RF, lineWidth)
        self.plotChannel(self._channel_RH, graph, y_RH, lineWidth)

        plt.xlim([self._plotStartTime, self._plotStopTime - 1])

        plt.ylim(0, 1)
        positions = (y_RH, y_RF, y_LF, y_LH)
        labels    = ('Right Hind', 'Right Fore', 'Left Fore','Left Hind')
        plt.yticks(positions, labels)
        return True

def plot_fit_and_tails(closeness, title):
    """Plot closeness centrality distribution, fit Gauss and mark outermosts and leaders areas in plot"""
    P.suptitle(title)
    m, sd = np.mean(closeness), np.std(closeness)
    outer = m - sd
    if outer > 0:
        P.axvspan(0, outer, alpha=0.2, color='c', label="Outermosts")
    leaders = m + sd
    if leaders < 1:
        P.axvspan(leaders, 0.2, alpha=0.2, color='y', label="Leaders")
    P.hist(closeness, color='#AB717A')
    x = np.linspace(0, 0.2)
    P.plot(x, mlab.normpdf(x, m, sd), '--', color='k')
    # P.xticks(np.arange(0, 0.2, 0.01))
    P.legend()

def ChroPlot(initial,final):
    ''' This function creates a graph of initial and final Chromosomes'''
    
    Title = 'Chromosome usage'
    plt.figure(Title)
    plt.title(Title)
    
    plt.plot(initial, 'r', lw=.05, label='initial')
    plt.plot(final, 'b', lw=.05, label='final')
    plt.axvspan(CenStart, CenEnd, color='none', lw=1.3, ec='k')
    
    plt.ylabel('CA')
    plt.xlabel('site')
    #plt.legend()       #labels invisible due to extremely thin lines   
    plt.show()