def albedo_parameter_plots(imarr, darr, params=None, plot_params=True,
                           ylabel='Reflectance', visible_only=True,
                           figsize=(12,7)):
    # from matplotlib import style
    # style.use('ggplot')
    if params is None:
        params = est_curve_params(darr, imarr)
    if visible_only:
        fig, axs = subplots(2, 3, figsize=figsize, sharey=False, sharex=True)
    else:
        fig, axs = subplots(2, 4, figsize=figsize, sharey=False, sharex=True)

    for i, ax in enumerate(axs.ravel()):
        if i >= imarr.shape[-1]:
            # This means I've got more axes than image bands so I'll skip plotting
            continue
        ax.scatter(darr.compressed(),imarr[...,i].compressed(), c='gold', alpha=0.2, edgecolor='none')
        cp = params[i]
        plotz = np.arange(darr.min(), darr.max(), 0.2)
        if plot_params:
            ax.plot(plotz, myR0(plotz, *cp), c='brown')
        ax.set_xlabel('Depth (m)')
        ax.set_ylabel(ylabel)
        btxt = "Band{b} $R_\infty = {R:.2f}$\n$A^{{toa}} = {A:.2f}$, $K_g = {Kg:.2f}$ "\
                .format(b=i+1, R=cp[0], A=cp[1], Kg=cp[2])
        ax.set_title(btxt)
    tight_layout()
    return fig

def iterTestRun(m, tspan):
    '''
    Method to perform iterations and generate plots
    '''
        
    # Test 16 random conditions - plot fft and phase diagrams
    fig, axes = subplots(4, 4)
    fig, axes2 = subplots(4, 4)
    for i in range(0, 4):
        for j in range(0, 4):
            # Run model
            yout = m.runModel(tspan, True)
            
            # Plot quiver diagram
            m.plotQuiver(yout, axes[i, j])
            
            # Plot Fourier transform
            freqs, aMags, rMags = m.calcFourier(yout, tspan)
            
            # Plot FFT
            axes2[i, j].plot(freqs, aMags, freqs, rMags)
            if i == 0:
                xlabel('Frequency (s^-1)')
            if j == 3:
                ylabel('Magnitude')
            
    legend(['Activator', 'Repressor'])

def plot_4_4_cdfs(df, param_name, teff_categories, age_categories, title='Somechart', figsize=(16, 4)):
    """

    """
    # sequential colors..
    colors = ['#fecc5c', '#fd8d3c', '#f03b20', '#bd0026']

    fig1, axes1 = plt.subplots(nrows=1, ncols=4, sharex=True, sharey=True, figsize=figsize)
    fig2, axes2 = plt.subplots(nrows=1, ncols=4, sharex=True, sharey=True, figsize=figsize)
    fig3, axes3 = plt.subplots(nrows=1, ncols=4, sharex=True, sharey=True, figsize=figsize)
    for teff_bin in xrange(len(teff_categories)):
        ax = axes1[teff_bin]
        ax.set_title("{} K".format(teff_categories[teff_bin]))

        serieses = []
        for age_bin in xrange(len(age_categories)):
            series = df[(df.teff_bins == teff_bin+1) & (df.age_bins == age_bin+1)][param_name]
            s = series.copy()
            s.sort()
            serieses.append(s)
            ax.plot(np.arange(s.size)/float(s.size), s, label=age_categories[age_bin], color=colors[age_bin])

        # plot the KS-score for each age-group (4choose2 graphs)
        plot_ks_scores(serieses, axes2[teff_bin])
        plot_ks_scores(serieses, axes3[teff_bin], log=True)

    axes2[0].axes.get_yaxis().set_visible(False)
    axes3[0].axes.get_yaxis().set_visible(False)
    add_colorful_yticks(axes2[0], colors, len(age_categories))
    fig1.suptitle(title, size='xx-large', y=1.08)
    fig3.suptitle("Distribution-Pairs KS scores", size='xx-large', y=1.00)
    axes1[-1].legend(loc='upper center', bbox_to_anchor=(-1.3, -0.05), fancybox=True, shadow=True, ncol=4)
    fig1.tight_layout()
    fig2.tight_layout()
    return fig1, fig2, fig3

def TablePlot():
    index_time=0
#    print [x[2] for x in RouteTablewithSeq]
        
    pl.ion()
    fig,ay=pl.subplots()
    fig,ax=pl.subplots()
    fig.set_tight_layout(True)
    
    idx_row = Index(np.arange(0,nNodes))
    idx_col = Index(np.arange(0,nPackets))
    df = DataFrame(cache_matrix[0,:,:], index=idx_row, columns=idx_col)
#    print df
    normal = pl.Normalize(0, 1)
    for index_time in range(len(RouteTablewithSeq_time)):
	vals=cache_matrix[index_time,:,:30]

#    fig = pl.figure(figsize=(15,8))
#    ax = fig.add_subplot(111, frameon=True, xticks=[], yticks=[])
#	print vals.shape
    	the_table=pl.table(cellText=vals, rowLabels=df.index, colLabels=df.columns, colWidths = [0.03]*vals.shape[1], loc='center', cellColours=pl.cm.hot(normal(vals)), fontsize=3)
	the_table.alpha=0
	for i in range(index_time+1):
		for j in range(vals.shape[0]):
			if (vals[j,i]==1):
				the_table._cells[(j+1, i)]._text.set_color('white')

	pl.title("Table at time: "+str(cache_time[index_time])+" Packet: "+str(index_time)+" Probability: "+str(p) )
    	pl.show()
	pl.pause(.0005)
	pl.clf()

def analyse_results(k,n,outpath=None):
    """Summarise multiple results"""

    if outpath != None:
        os.chdir(outpath)
    #add mirbase info
    k = k.merge(mirbase,left_on='name',right_on='mature1')
    ky1 = 'unique reads'
    ky2 = 'read count' #'RC'
    cols = ['name','freq','mean read count','mean_norm','total','perc','mirbase_id']
    print
    print ('found:')
    idcols,normcols = get_column_names(k)
    final = filter_expr_results(k,freq=.8,meanreads=200)
    print (final[cols])
    print ('-------------------------------')
    print ('%s total' %len(k))
    print ('%s with >=10 mean reads' %len(k[k['mean read count']>=10]))
    print ('%s found in 1 sample only' %len(k[k['freq']==1]))
    print ('top 10 account for %2.2f' %k['perc'][:10].sum())

    fig,ax = plt.subplots(nrows=1, ncols=1, figsize=(8,6))
    k.set_index('name')['total'][:10].plot(kind='barh',colormap='Spectral',ax=ax,log=True)
    plt.tight_layout()
    fig.savefig('srnabench_top_known.png')
    #fig = plot_read_count_dists(final)
    #fig.savefig('srnabench_known_counts.png')
    fig,ax = plt.subplots(figsize=(10,6))
    k[idcols].sum().plot(kind='bar',ax=ax)
    fig.savefig('srnabench_total_persample.png')
    print
    k.to_csv('srnabench_known_all.csv',index=False)
    return k

def plot_pis(indices,it):
    f,ax=pl.subplots(len(indices),sharex=True)
    for a in range(len(indices)):
        k=people[a]
        x=[]
        y=[]
        for j in range(settings.nlabels):
            for l in range(settings.nscores):
                x.append(l+j*settings.nlabels-0.4)
                y.append(pi(k,j,l))
        ax[a].bar(x,y)
        ax[a].get_yaxis().set_ticks([])
        ax[a].set_ylim(0,1)
    ax[0].set_xlim(-0.5,len(y)-0.5)
    pl.savefig('D:\Documents\images\pi_%03d.png' %it)
    
    f,ax=pl.subplots(len(indices),sharex=True)
    for a in range(len(indices)):
        i=indices[a]
        x=[]
        y=[]
        for j in range(settings.nlabels):
            x.append(j-0.4)
            y.append(results[i][j])
        ax[a].bar(x,y)
        ax[a].get_yaxis().set_ticks([])
        ax[a].set_ylim(0,1)
    ax[0].set_xlim(-0.5,len(y)-0.5)
    pl.savefig('D:\Documents\images\\result_%03d.png' %it)
    pl.close('all')

    def plot_network_representation(self):
        '''
            Plot the response of the network
        '''

        if self.W_type == 'dirichlet':
            # Sorting that emphasis balance
            balanced_indices_neurons = self.number_connections[:,0].argsort()[::-1]
        else:
            balanced_indices_neurons = np.arange(self.M)

        # Plot the population response
        plot_separation_y = 0.3*(np.max(self.network_representations) - np.min(self.network_representations))

        fig1, ax1 = plt.subplots(1)

        for r in xrange(self.R):
            ax1.plot(self.network_representations[r, :, balanced_indices_neurons] + np.arange(self.K)*plot_separation_y + r*(self.K+0.5)*plot_separation_y)

        ax1.autoscale(tight=True)

        # Plot Hinton graphs
        sf, ax = plt.subplots(self.R, 1)

        for r in xrange(self.R):
            hinton(self.W[r, balanced_indices_neurons].T, ax=ax[r])

def analyse_diff_res2(res1, res2, res3, res4):
    plt.figure(95)
    plt.clf()
    #plt.title('Energy vs time')
    f, (ax1, ax2) = plt.subplots(2, 1)
    f.subplots_adjust(hspace=0.25)
    ax1.plot(res3['sim'][0] / 86400, np.array(res3['sim'][4]) / 1e15, 
             label='Eulerian')
    ax1.plot(res4['sim'][0] / 86400, np.array(res4['sim'][4]) / 1e15, 
             label='semi-Lagrangian')
    ax1.set_xlabel('time (days)')
    ax1.set_ylabel('energy (PJ)')
    ax1.legend(loc='lower right')

    ax2.plot(res3['sim'][0] / 86400, np.array(res3['sim'][4]) / 1e15, 
             label='Eulerian')
    ax2.plot(res4['sim'][0] / 86400, np.array(res4['sim'][4]) / 1e15, 
             label='semi-Lagrangian')
    ax2.set_xlim((60, 100))
    ax2.set_ylim((2.65, 2.75))
    ax2.set_xlabel('time (days)')
    ax2.set_ylabel('energy (PJ)')

    print('{}: energy_diff={}'.format(res3['res'], res3['energy_diff']))
    print('{}: energy_diff={}'.format(res4['res'], res4['energy_diff']))

    pd1 = (np.array(res3['sim'][4]).max() - res3['sim'][4][-1] ) / np.array(res3['sim'][4]).max() * 100
    pd2 = (np.array(res4['sim'][4]).max() - res4['sim'][4][-1] ) / np.array(res4['sim'][4]).max() * 100
    print('% diff 1 {}'.format(pd1))
    print('% diff 2 {}'.format(pd2))

    plt.savefig('writeup/figures/task_d_energy.png')

    f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
    f.subplots_adjust(hspace=0.25, wspace=0.25, right=0.8)

    X, Y = res1['grid']
    ax1.set_title('(a) $\eta$ after 1 day - Eul.')
    ax1.contourf(X/1e3, Y/1e3, res1['sim'][1], 100, vmin=-0.01, vmax=0.01)
    ax2.set_title('(b) $\eta$ after 1 day - S-L')
    cf2 = ax2.contourf(X/1e3, Y/1e3, res2['sim'][1], 100, vmin=-0.01, vmax=0.01)

    cbar_ax1 = f.add_axes([0.85, 0.55, 0.05, 0.4])
    f.colorbar(cf2, cax=cbar_ax1)

    ax3.set_title('(c) $\eta$ after 100 days - Eul.')
    ax3.contourf(X/1e3, Y/1e3, res3['sim'][1], 100, vmin=-0.15, vmax=0.2)
    ax4.set_title('(d) $\eta$ after 100 days - S-L')
    cf4 = ax4.contourf(X/1e3, Y/1e3, res4['sim'][1], 100, vmin=-0.15, vmax=0.2)

    ax3.set_xlabel('x (km)') 
    ax4.set_xlabel('x (km)') 
    ax1.set_ylabel('y (km)') 
    ax3.set_ylabel('y (km)') 

    cbar_ax2 = f.add_axes([0.85, 0.05, 0.05, 0.4])
    f.colorbar(cf4, cax=cbar_ax2)
    plt.savefig('writeup/figures/task_d_eta.png')

def create_fig( file_name,
                small=True,
                marker=False,
                figsize=None,
                nrows=1,
                ncols=1,
                sharex=False,
                sharey=False ):
    if not isinstance( file_name, list ):
        file_name = [file_name]
    defaults = plot_defaults._get_defaults(small)
    params = plot_defaults._get_params( small )
    pylab.rcParams.update( params )
    pylab.rc('font', **defaults[ 'font' ])
    if figsize:
        fig, axs = pylab.subplots( nrows=nrows, ncols=ncols,
                                   sharex=sharex, sharey=sharey,
                                   figsize=figsize)
    else:
        fig, axs = pylab.subplots( nrows=nrows, ncols=ncols,
                                   sharex=sharex, sharey=sharey )
    if not isinstance(axs, list) and not isinstance(axs, ndarray):
        axs = [axs]
    for i, f_name in enumerate(file_name):
        curves = _input_from_txt(file_name = f_name, small=False)
        labels = curves.keys()
        labels.sort()
        for label in labels:
            curve = curves[ label ]
            number = label.split('_')[0]
            label  = '_'.join( label.split('_')[1:] )
            number = int(number)
            curr_defaults = plot_defaults._cycle_defaults( number, small=False )
            def current_default( key ):
                if key in curve.keys():
                    return curve[ key ]
                else:
                    if key in curr_defaults.keys():
                        return curr_defaults[key]
                    else:
                        if marker:
                            key += '_marker'
                            if key in curr_defaults.keys():
                                return curr_defaults[key]
                            else:
                                return None
            axs[i].plot( curve['x'], curve['y'],
                  color           = current_default( 'color'     ),
                  label           = label,
                  linewidth       = current_default( 'linewidth' ),
                  linestyle       = current_default( 'linestyle' ),
                  marker          = current_default( 'marker'    ),
                  markerfacecolor = current_default( 'markerfacecolor' ),
                  markeredgecolor = current_default( 'markeredgecolor' ),
                  markeredgewidth = current_default( 'markeredgewidth' ),
                  markersize      = current_default( 'markersize' )       )

    return fig

def events_per_sim_cycle_histograms():
    data = np.loadtxt("analysisData/eventsAvailableBySimCycle.csv", dtype=np.intc, delimiter = ",", skiprows=2)
    trimmedData = reject_first_last_outliers(data)
    fig, ax1 = pylab.subplots()
    outFile = outDir + 'eventsAvailableBySimCycle-histogram-std'
    pylab.title('Total LPs: %s; ' % "{:,}".format(total_lps) +
                'Total Sim Cycles: %s. '% "{:,}".format(len(trimmedData)))
    ax1.hist(trimmedData, bins=100, histtype='stepfilled')
    ax1.set_xlabel('Number of Events')
    ax1.set_ylabel('Number of Simulation Cycles')
    ax2=ax1.twinx()
    ax2.hist(trimmedData, bins=100, histtype='stepfilled')
    ax2.set_ylabel('Percent of Simulation Cycles')
    ax2.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(toPercentOfTotalSimCycles))
    display_graph(outFile)

    # ok, so now let's build a histogram of the % of LPs with active events.

    fig, ax1 = pylab.subplots()
    outFile = outDir + 'percentOfLPsWithAvailableEvents'
    pylab.title('Percent of LPs w/ Available Events as a Percentage of the Total Sim Cycles')
    
    ax1.hist(trimmedData.astype(float)/float(total_lps), bins=100, histtype='stepfilled')
    ax1.set_xlabel('Number of Events as a percentage of Total LPs')
    ax1.set_ylabel('Number of Sim Cycles said Percentage Occurs')
#    ax1 = pylab.gca()
    ax1.get_xaxis().set_major_formatter(mpl.ticker.FuncFormatter(toPercent))
#    ax.get_yaxis().set_major_formatter(mpl.ticker.FormatStrFormatter('%.1f%%'))
    ax2=ax1.twinx()
    ax2.hist(trimmedData, bins=100, histtype='stepfilled')
    ax2.set_ylabel('Percent of Simulation Cycles')
    ax2.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(toPercentOfTotalSimCycles))
    display_graph(outFile)

    # ok, let's present the histogram data using pandas series/value_counts.  much nicer plot.
    fig, ax1 = pylab.subplots()
    outFile = outDir + 'eventsAvailableBySimCycle-histogram-dual'
    pylab.title('Total LPs: %s; ' % "{:,}".format(total_lps) +
                'Total Sim Cycles: %s. '% "{:,}".format(len(trimmedData)))
    setNumOfSimCycles(len(data)+1)
    mean_events_available = np.mean(trimmedData)
    data = pd.Series(trimmedData).value_counts()
    data = data.sort_index()
    x_values = np.array(data.keys())
    y_values = np.array(data)
    ax1.plot(x_values, y_values)
    ax1.set_xlabel('Number of Events (Ave=%.2f)' % mean_events_available)
    ax1.set_ylabel('Number of Simulation Cycles')
    ax2=ax1.twinx()
    ax2.plot(x_values, y_values)
    ax2.set_ylabel('Percent of Simulation Cycles')
    ax2.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(toPercentOfTotalSimCycles))
    display_graph(outFile)

    return

def analyze(X, Y, data, classifier, label, outDir):
    #store the means for our ROC curve
    meanTPRate = 0.0
    meanFPRate = np.linspace(0, 1, 100)

    #start with the subplot
    pl.subplots()

    #now lets analyze the data
    for i, (train, test) in enumerate(data):
        #grab the data sets for training and testing
        xTrain, xTest, yTrain, yTest = X[train], X[test], Y[train], Y[test]
        #print xTrain
        #print yTrain

        #train the model
        classifier.fit(xTrain, yTrain)

        #now predict on the hold out
        predictions = classifier.predict(xTest)
        probabilities = classifier.predict_proba(xTest)

        #compute ROC and AUC
        fpr, tpr, thresholds = roc_curve(yTest, probabilities[:, 1])
        meanTPRate += interp(meanFPRate, fpr, tpr)
        meanTPRate[0] = 0.0
        rocAUC = auc(fpr,tpr)

        #now plot it out
        pl.plot(fpr, tpr, lw=1, label='ROC Iter %d (area = %0.2f)' % (i+1, rocAUC))

        #print "P: %s\nA: %s\n" % (predictions, yTest)

    #now plot the random chance line
    pl.plot([0,1], [0,1], '--', color=(0.6,0.6,0.6), label='Random Chance')

    #generate some stats for the mean plot
    meanTPRate /= len(data)
    meanTPRate[-1] = 1.0
    meanAUC = auc(meanFPRate, meanTPRate)

    #plot the average line
    pl.plot(meanFPRate, meanTPRate, 'k--', label='Mean ROC (area = %0.2f)' % (meanAUC), lw=2)

    #add some other plot parameters
    pl.xlim([-0.05, 1.05])
    pl.ylim([-0.05, 1.05])
    pl.xlabel('False Positive Rate')
    pl.ylabel('True Positive Rate')
    pl.title('ROC Curve (%s)' % (label))
    pl.legend(loc="lower right")

    pl.savefig('%s/%s.png' % (outDir, label))
    pl.close()
    logging.info("%s\t%f" % (label, meanAUC))

def analyse_isomirs(iso,outpath=None):
    """Analyse isomiR results in detail"""

    if iso is None:
        return
    if outpath != None:
        os.chdir(outpath)
    subcols = ['name','read','isoClass','NucVar','total','freq']
    iso = iso.sort_values('total', ascending=False)
    #filter very low abundance reads
    iso = iso[(iso.total>10) & (iso.freq>0.5)]
    top = get_top_isomirs(iso)
    top.to_csv('srnabench_isomirs_dominant.csv',index=False)
    print ('top isomiRs:')
    print (top[:20])
    print ('%s/%s with only 1 isomir' %(len(top[top.domisoperc==1]),len(top)))
    print ('different dominant isomir:', len(top[top.variant!='exact'])/float(len(top)))
    print ('mean dom isomir perc:', top.domisoperc.mean())
    print
    #stats
    fig,ax = plt.subplots(1,1)
    top.plot('isomirs','total',kind='scatter',logy=True,logx=True,alpha=0.8,s=50,ax=ax)
    ax.set_title('no. isomiRs per miRNA vs total adundance')
    ax.set_xlabel('no. isomiRs')
    ax.set_ylabel('total reads')
    fig.savefig('srnabench_isomir_counts.png')
    fig,ax = plt.subplots(1,1)
    #top.hist('domisoperc',bins=20,ax=ax)
    try:
        base.sns.distplot(top.domisoperc,bins=15,ax=ax,kde_kws={"lw": 2})
    except:
        pass
    fig.suptitle('distribution of dominant isomiR share of reads')
    fig.savefig('srnabench_isomir_domperc.png')

    x = iso[iso.name.isin(iso.name[:28])]
    bins=range(15,30,1)
    ax = x.hist('length',bins=bins,ax=ax,by='name',sharex=True,alpha=0.9)
    ax[-1,-1].set_xlabel('length')
    fig.suptitle('isomiR length distributions')
    fig.savefig('srnabench_isomir_lengths.png')
    plt.close('all')

    c=iso.variant.value_counts()
    #c=c[c>10]
    fig,ax = plt.subplots(1,1,figsize=(8,8))
    c.plot(kind='pie',colormap='Spectral',ax=ax, labels=None,legend=True,
             startangle=0,pctdistance=1.1,autopct='%.1f%%',fontsize=10)
    ax.set_title('isomiR class distribution')
    plt.tight_layout()
    fig.savefig('srnabench_isomir_classes.png')

    iso.to_csv('srnabench_isomirs_all.csv',index=False)
    return top

def GenomeChromosomewise(df, candSNPs=None, genes=None, axes=None,outliers=None):
    markerSize = 6
    fontsize = 6
    chrsize = df.reset_index().groupby('CHROM').POS.max()
    if axes is None:
        if chrsize.shape[0]>1:
            _, axes = plt.subplots(int(np.ceil(chrsize.shape[0] / 2.)), 2, sharey=True, dpi=200, figsize=(10, 6));
            ax = axes.reshape(-1)
        else:
            ax = [plt.subplots(1,1, sharey=True, dpi=200, figsize=(10, 6))[1]]

    for j, (chrom, a) in enumerate(df.groupby(level=0)):
        if candSNPs is not None:
            try:
                candSNPs.loc[chrom]
                for pos in candSNPs.loc[chrom].index.values:
                    ax[j].axvline(pos, color='r', linewidth=0.5, alpha=0.5)
                    ax[j].annotate(
                        '{:.0f},{:.2f}'.format(candSNPs['rank'].loc[(chrom, pos)], candSNPs.nu0.loc[(chrom, pos)]),
                        xy=(pos, a.max()), xytext=(pos, a.max()), fontsize=fontsize - 2)

            except:
                pass

        if genes is not None:
            try:
                X = genes.loc[chrom]
                if len(genes.loc[chrom].shape) == 1:
                    X = pd.DataFrame(X).T
                for _, row in X.iterrows():
                    ax[j].fill_between([row.start, row.end], a.min(), a.max(), color='r')
                    ax[j].annotate(row['name'], xy=(row.start, a.max()), xytext=(row.start, a.max()),
                                   fontsize=fontsize - 2)


            except:
                pass

        ax[j].scatter(a.loc[chrom].index, a.loc[chrom], s=markerSize, alpha=0.8, edgecolors='none')

        if outliers is not None:
            try:
                ax[j].scatter(outliers.loc[chrom].index, outliers.loc[chrom], s=markerSize, c='r', alpha=0.8, edgecolors='none')
            except:
                pass

        setSize(ax[j], fontsize)
        ax[j].set_xlim([-1000, chrsize[chrom] + 1000])
        # ax[j].set_title(chrom, fontsize=fontsize+2)
        annotate(chrom, ax=ax[j],fontsize=fontsize+4)
        ax[j].locator_params(axis='x', nbins=10)
    plt.tight_layout(pad=0.1)
    plt.gcf().subplots_adjust(bottom=0.1)

def printHeatMap(marginals, words, outFile):
    N = len(words)
    words_uni = [i.decode('UTF-8') for i in words]
    heatmap = np.zeros((N+1, N+1))
    for chart in marginals:
        heatmap[chart[0], chart[1]] = math.log(marginals[chart])
    fig, ax = plt.subplots()    
    mask = np.tri(heatmap.shape[0], k=0)
    heatmap = np.ma.array(heatmap, mask=mask)
    cmap = plt.cm.get_cmap('RdBu')
    cmap.set_bad('w')
    im = ax.pcolor(heatmap, cmap=cmap, alpha=0.8)
    font = mpl.font_manager.FontProperties(fname='/usr0/home/avneesh/spectral-scfg/data/wqy-microhei.ttf')
    ax.grid(True)
    ax.set_ylim([0,N])
    ax.invert_yaxis()
    ax.set_yticks(np.arange(heatmap.shape[1]-1)+0.5, minor=False)
    ax.set_yticklabels(words_uni, minor=False, fontproperties=font)
    ax.set_xticks(np.arange(heatmap.shape[0])+0.5, minor=True)
    ax.set_xticklabels(np.arange(heatmap.shape[0]), minor=True)
    ax.set_xticks([])
    cbar = fig.colorbar(im, use_gridspec=True)
    cbar.set_label('ln(sum)')
    ax.set_xlabel('Span End')
    ax.xaxis.set_label_position('top')
    ax.xaxis.tick_top()
    plt.ylabel('Span starting at word: ')
    plt.tight_layout()
    #ax.set_title('CKY Heat Map: Node Marginals')
    fig.savefig(outFile)    

def plotScalingFactor():
    r=2*1e-8
    l = 5e4
    dpi = 300
    j = 0
    for nu0 in [0.005, 0.1]:
        for s in [0.025, 0.1]:
            t = np.arange(0, 2 * (utl.logit(0.995) - utl.logit(nu0)) / s + 1., 1)
            fig, ax = plt.subplots(2, 1, figsize=(5.5, 2.5), dpi=dpi, sharex=True);
            nu(t, s=s, nu0=nu0).plot(color='k', legend=False, ax=ax[0])
            pplt.annotate(r'$s$={}, $\nu_0=${} ({} Sweep)'.format(s, nu0, ('Soft', 'Hard')[nu0 == 0.005]), fontsize=7,
                          ax=ax[0])
            pplt.setSize(ax=ax[0], fontsize=6)
            ax[0].set_ylabel(r'$\nu_t$')
            #
            H0 = H(t[0], s=s, nu0=nu0)
            Ht = H(t, s=s, nu0=nu0)
            df = pd.DataFrame([np.log(Ht / H0), -2 * r * t * l], columns=t, index=['log(Growth)', r'log(Decay)']).T
            df['log(Growth) + log(Decay)'] = df.sum(1)
            df.plot(ax=ax[1], grid=True, linewidth=2);
            ax[1].set_xlabel('Generations');
            ax[1].set_ylabel('Log(Scaling Factor)')
            ax[1].axvline(df.iloc[1:, 2].abs().idxmin(), color='k', linestyle='--', linewidth=0.5)
            # if j != 3:
            #     ax[1].legend_.remove()
            # else:
            ax[1].legend(['log(Growth)', r'log(Decay)', 'log(Growth) + log(Decay)'], bbox_to_anchor=(1.45, .75),
                         prop={'size': 6})
            pplt.setSize(ax[1], fontsize=6)

            plt.tight_layout(pad=0.1, rect=[0, 0, 0.7, 1])
            plt.gcf().subplots_adjust(bottom=0.15)
            pplt.savefig('decayFactors{}'.format(j), dpi=dpi)
            j += 1

def visData(dF="dannyFull.db"):
    f, ax = pl.subplots()
    m = visBg(ax)
    D = sq.connect(dF)
    cur = D.cursor()
    A = cur.execute(
        """select clg,clat,cluster
	from meta
	group by clg,clat,cluster"""
    )
    for w in A:
        clg = w[0]
        clat = w[1]
        x, y = m(clg, clat)
        cluster = w[2]
        ax.plot(x, y, "s", ms=10, color=pickColor(cluster))
    A = cur.execute(
        """select slg,slat
	from meta
	group by slg,slat"""
    )
    for w in A:
        slg = w[0]
        slat = w[1]
        x, y = m(slg, slat)
        ax.plot(x, y, "b*", ms=20)
    ax.plot(-float("inf"), -float("inf"), "b*", label="Servers")
    ax.plot(-float("inf"), -float("inf"), "bs", label="Clients")
    ax.legend(loc=4)
    pl.show()

def errorVis(e, t):
    f, ax = pl.subplots(1)
    ax.boxplot(e, widths=0.15)
    ax.set_xticks([1, 2])
    ax.set_xticklabels(t)
    ax.set_ylabel("Error in km", fontsize=20)
    pl.show()

 def plot(self, theta_a, theta_e, sigma_e, y, klic, N, T, D):
     # Set up an output figure
     fig, ax = pylab.subplots(3, 1)
     colours = ['b', 'g', 'r', 'c', 'm', 'y']
     # Plot smoothed densities
     for i in xrange(D):
         # Plot original theta values
         ax[0].plot(theta_a[:, i], ls='--', c=colours[i%len(colours)])
         # Plot data-estimated theta values
         ax[0].plot(theta_e[:,i], ls='-', c=colours[i%len(colours)])
         # Plot data-estimated confidence intervals
         ax[0].fill_between(numpy.arange(T),
             theta_e[:,i] - 2*numpy.sqrt(sigma_e[:,i,i]),
             theta_e[:,i] + 2*numpy.sqrt(sigma_e[:,i,i]),
             color=colours[i%len(colours)], alpha=.25)
         ax[1].plot(y[:,i], c=colours[i%len(colours)], ls='-')
     # Set axes labels and legends
     ax[0].set_ylabel('Theta')
     ax[0].set_title('Smoothed densities')
     ax[1].set_title('Observed pattern rates')
     ax[1].set_ylabel('Rate (patterns/second)')
     ax[2].set_title('KL Divergence')
     ax[2].set_ylabel('Bits')
     ax[2].plot(klic)
     ax[2].set_xlabel('Time (ms)')
     fig.tight_layout()
     pylab.show()

    def plot_calibration(self):
        """Plots the calibration function and data
        Arguments
        ------------
            cal: calibration instance

        Returns
        ------------
            figure
            axes
            calibration data graph
            calibration function graph
        """
        cal_function = self.calibrationFunction
        bin_no = np.logspace(np.log10(self.data.bin_no.min()), np.log10(self.data.bin_no.max()), 500)
        d = cal_function(bin_no)

        f, a = plt.subplots()

        cal_data, = a.plot(self.data.d,  self.data.bin_no, 'o', label='data',)
        cal_func, = a.plot(d, bin_no, label='function')

        a.loglog()

        a.set_xlim(0.9*self.data.d.min(), 1.1*self.data.d.max())
        a.set_xlabel('Diameter (nm)')

        a.set_ylim(0.9*self.data.bin_no.min(), 1.1*self.data.bin_no.max())
        a.set_ylabel('bin number')

        a.set_title('Calibration curve')
        a.legend(loc = 2)
        return f, a, cal_data, cal_func

def help_dcols():
	""" Show a plot demonstrating all the colours """
	from pylab import figure,clf,subplots,show
	from matplotlib.patches import Rectangle
	fig, axes = subplots(4,3)
	i=0
	print(axes)
	for axA in axes:
		for ax in axA:
			patch = Rectangle([0,0],1,1,color=cols[i])
			ax.add_patch(patch)
			ax.set_xticks([])
			ax.set_yticks([])
			ax.set_xticklabels([])
			ax.set_yticklabels([])			
			if colname[i] is not 'd_black':
				ax.text(0.2,0.4,colname[i],color='k')
			else:
				ax.text(0.2,0.4,colname[i],color='w')
			i+=1
	fig.subplots_adjust(bottom=0.25,left=0.03,right=0.98,top=0.98,wspace=0.05,hspace=0.05)
	fig.text(0.08,0.18,'To use: from durhamcolours import *')
	fig.text(0.08,0.12,'then e.g. ax.plot(x,y,colour=d_olive)')
	fig.text(0.08,0.06,"Note: no parentheses i.e. d_olive not 'd_olive'")
	
	fig.savefig('durhamcolors.png')
	fig.savefig('durhamcolors.pdf')
	
	show()