Python pylab.subplots() Examples

def plot_pca(pX, palette='Spectral', labels=None, ax=None, colors=None):
    """Plot PCA result, input should be a dataframe"""

    if ax==None:
        fig,ax=plt.subplots(1,1,figsize=(6,6))
    cats = pX.index.unique()
    colors = sns.mpl_palette(palette, len(cats)+1)
    print (len(cats), len(colors))
    for c, i in zip(colors, cats):
        #print (i, len(pX.ix[i]))
        #if not i in pX.index: continue
        ax.scatter(pX.ix[i, 0], pX.ix[i, 1], color=c, s=90, label=i,
                   lw=.8, edgecolor='black', alpha=0.8)
    ax.set_xlabel('PC1')
    ax.set_ylabel('PC2')
    i=0
    if labels is not None:
        for n, point in pX.iterrows():
            l=labels[i]
            ax.text(point[0]+.1, point[1]+.1, str(l),fontsize=(9))
            i+=1
    ax.legend(fontsize=10,bbox_to_anchor=(1.5, 1.05))
    sns.despine()
    plt.tight_layout()
    return 

def contour_mult_mo(self,x,y,mo,xlabel='x',ylabel='y',title='',r0=0):
    '''Uses matplotlib to show slices of a molecular orbitals.'''
    import matplotlib.pyplot as plt
    
    # Plot slices
    f, pics = \
                plt.subplots(len(mo),1,sharex=True,sharey=True,figsize=(6,2+4*len(mo)))
    plt.suptitle(title)
    vmax = numpy.max(numpy.abs(mo))
    for i,pic in enumerate(pics):
      pic.contour(y,x,mo[i],50,linewidths=0.5,colors='k')
      pic.contourf(\
          y,x,mo[i],50,cmap=plt.cm.rainbow,vmax=vmax,vmin=-vmax)
      pic.set_ylabel(xlabel)  
      pic.set_xlabel(ylabel)  
      pic.set_title('Data Point %d' % (r0+i))
    
    f.subplots_adjust(left=0.15,bottom=0.05,top=0.95,right=0.95)
    f.show()
    return f,pics 

def __init__(self, project):
        self._ref_path = project.path
        self.project = project.copy()
        self.project._inspect_mode = True
        self.parent = None
        self.name = None
        self.fig, self.ax = None, None
        self.w_group = None
        self.w_node = None
        self.w_file = None
        self.w_text = None
        self.w_tab = None
        self.w_path = None
        self.w_type = None
        #         self.fig, self.ax = plt.subplots()

        self.create_widgets()
        self.connect_widgets()
        self.display() 

def plot_read_count_dists(counts, h=8, n=50):
    """Boxplots of read count distributions """

    scols,ncols = base.get_column_names(counts)
    df = counts.sort_values(by='mean_norm',ascending=False)[:n]
    df = df.set_index('name')[ncols]
    t = df.T
    w = int(h*(len(df)/60.0))+4
    fig, ax = plt.subplots(figsize=(w,h))
    if len(scols) > 1:
        sns.stripplot(data=t,linewidth=1.0,palette='coolwarm_r')
        ax.xaxis.grid(True)
    else:
        df.plot(kind='bar',ax=ax)
    sns.despine(offset=10,trim=True)
    ax.set_yscale('log')
    plt.setp(ax.xaxis.get_majorticklabels(), rotation=90)
    plt.ylabel('read count')
    #print (df.index)
    #plt.tight_layout()
    fig.subplots_adjust(bottom=0.2,top=0.9)
    return fig 

def plot_dos(self, ax=None, *args, **qwargs):
        """

        Args:
            *args:
            ax:
            **qwargs:

        Returns:

        """
        try:
            import pylab as plt
        except ImportError:
            import matplotlib.pyplot as plt
        if ax is None:
            fig, ax = plt.subplots(1, 1)
        ax.plot(self["output/dos_energies"], self["output/dos_total"], *args, **qwargs)
        ax.set_xlabel("Frequency [THz]")
        ax.set_ylabel("DOS")
        ax.set_title("Phonon DOS vs Energy")
        return ax 

def plot_contourf(self, ax=None, show_min_erg_path=False):
        """

        Args:
            ax:
            show_min_erg_path:

        Returns:

        """
        try:
            import pylab as plt
        except ImportError:
            import matplotlib.pyplot as plt
        x, y = self.meshgrid()
        if ax is None:
            fig, ax = plt.subplots(1, 1)
        ax.contourf(x, y, self.energies)
        if show_min_erg_path:
            plt.plot(self.get_minimum_energy_path(), self.temperatures, "w--")
        plt.xlabel("Volume [$\AA^3$]")
        plt.ylabel("Temperature [K]")
        return ax 

def plot_min_energy_path(self, *args, ax=None, **qwargs):
        """

        Args:
            *args:
            ax:
            **qwargs:

        Returns:

        """
        try:
            import pylab as plt
        except ImportError:
            import matplotlib.pyplot as plt
        if ax is None:
            fig, ax = plt.subplots(1, 1)
            ax.xlabel("Volume [$\AA^3$]")
            ax.ylabel("Temperature [K]")
        ax.plot(self.get_minimum_energy_path(), self.temperatures, *args, **qwargs)
        return ax 

def plot(self):
        """ Plot the layer data (for debugging)
        :return: The current figure
        """
        import pylab as pl
        aspect = self.nrows / float(self.ncols)
        figure_width = 6 #inches

        rows = max(1, int(np.sqrt(self.nlayers)))
        cols = int(np.ceil(self.nlayers/rows))
        # noinspection PyUnresolvedReferences
        pallette = {i:rgb for (i, rgb) in enumerate(pl.cm.jet(np.linspace(0, 1, 4), bytes=True))}
        f, a = pl.subplots(rows, cols)
        f.set_size_inches(6 * cols, 6 * rows)
        a = a.flatten()
        for i, label in enumerate(self.label_names):
            pl.sca(a[i])
            pl.title(label)
            pl.imshow(self.color_data)
            pl.imshow(colorize(self.label_data[:, :, i], pallette), alpha=0.5)
            # axis('off')
        return f 

def plot_fractions(df, label=None):
    """Process results of multiple mappings to get fractions
    of each annotations mapped
    label: plot this sample only"""

    fig,ax = plt.subplots(figsize=(8,8))
    df = df.set_index('label')
    df = df._get_numeric_data()
    if len(df) == 1:
        label = df.index[0]
    if label != None:
        ax = df.T.plot(y=label,kind='pie',colormap='Spectral',autopct='%.1f%%',
                      startangle=0, labels=None,legend=True,pctdistance=1.1,
                      fontsize=10, ax=ax)
    else:
        ax = df.plot(kind='barh',stacked=True,linewidth=0,cmap='Spectral',ax=ax)
        #ax.legend(ncol=2)
        ax.set_position([0.2,0.1,0.6,0.8])
        ax.legend(loc="best",bbox_to_anchor=(1.0, .9))
    plt.title('fractions mapped')
    #plt.tight_layout()
    return fig 

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() 

def get_subplot_axes(self):
        if not self.to_save and not self.to_show:
            return

        f, axes = pylab.subplots(2, 4, sharex='col', sharey='row')
        return [axes[0, 0], axes[0, 1], axes[0, 2], axes[0, 3], axes[1, 0], axes[1, 1], axes[1, 2], axes[1, 3]] 

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_results(k):
    fig,ax=plt.subplots(figsize=(8,6))
    ax.set_title('sRNAbench top 10')
    k.set_index('name')['read count'][:10].plot(kind='barh',colormap='Set2',ax=ax,log=True)
    plt.tight_layout()
    fig.savefig('srnabench_summary_known.png')
    return 

def plot_by_label(X, palette='Set1'):
    """Color scatter plot by dataframe index label"""

    import seaborn as sns
    cats = X.index.unique()
    colors = sns.mpl_palette(palette, len(cats))
    #sns.palplot(colors)
    f,ax = plt.subplots(figsize=(6,6))
    for c, i in zip(colors, cats):
        #print X.ix[i,0]
        ax.scatter(X.ix[i, 0], X.ix[i, 1], color=c, s=100, label=i,
                   lw=1, edgecolor='black')
    ax.legend(fontsize=10)
    sns.despine()
    return 

def plot_sample_variation(df):

    fig,axs=plt.subplots(2,1,figsize=(6,6))
    axs=axs.flat
    cols,ncols = mirdeep2.get_column_names(m)
    x = m.ix[2][cols]
    x.plot(kind='bar',ax=axs[0])
    x2 = m.ix[2][ncols]
    x2.plot(kind='bar',ax=axs[1])
    sns.despine(trim=True,offset=10)
    plt.tight_layout()
    return fig 

def plot_read_lengths(filename, df=None):
    """View read length distributions"""

    df = utils.fastq_to_dataframe(filename, size=5e5)
    x = analysis.read_length_dist(df)
    fig,ax=plt.subplots(1,1,figsize=(10,4))
    ax.bar(x[1][:-1],x[0], align='center')
    return fig 

def plot_array(self, val):
        """

        Args:
            val:

        Returns:

        """
        try:
            import pylab as plt
        except ImportError:
            import matplotlib.pyplot as plt

        plt.ioff()
        if self.fig is None:
            self.fig, self.ax = plt.subplots()
        else:
            self.ax.clear()

        # self.ax.set_title(self.name)
        if val.ndim == 1:
            self.ax.plot(val)
        elif val.ndim == 2:
            if len(val) == 1:
                self.ax.plot(val[0])
            else:
                self.ax.plot(val)
        elif val.ndim == 3:
            self.ax.plot(val[:, :, 0])
        # self.fig.canvas.draw()
        self.w_text.value = self.plot_to_html()
        plt.close() 

def plot_dense_stats(self):
        """
        Plot dense layers weight statistics

        :return: A plot
        :History: 2018-May-12 - Written - Henry Leung (University of Toronto)
        """
        self.has_model_check()
        dense_list = []
        for counter, layer in enumerate(self.keras_model.layers):
            if isinstance(layer, tfk.layers.Dense):
                dense_list.append(counter)

        denses = np.array(self.keras_model.layers)[dense_list]
        fig, ax = plt.subplots(1, figsize=(15, 10), dpi=100)
        for counter, dense in enumerate(denses):
            weight_temp = np.array(dense.get_weights())[0].flatten()
            ax.hist(weight_temp, 200, density=True, range=(-2., 2.), alpha=0.7,
                    label=f'Dense Layer {counter}, max: {weight_temp.max():.{2}f}, min: {weight_temp.min():.{2}f}, '
                          f'mean: {weight_temp.mean():.{2}f}, std: {weight_temp.std():.{2}f}')
        fig.suptitle(f'Dense Layers Weight Statistics of {self.folder_name}', fontsize=17)
        ax.set_xlabel('Weights', fontsize=17)
        ax.set_ylabel('Normalized Distribution', fontsize=17)
        ax.minorticks_on()
        ax.tick_params(labelsize=15, width=3, length=10, which='major')
        ax.tick_params(width=1.5, length=5, which='minor')
        ax.legend(loc='best', fontsize=15)
        fig.tight_layout(rect=[0, 0.00, 1, 0.96])
        fig.show()

        return fig 

def do_plot(data, a, b):
  """Create a plot."""
  fig, ax = plt.subplots()
  X, Y = zip(*data)
  plt.plot(X, Y, 'b-o', label="data")
  X = np.linspace(0, data[-1][0], 1000)
  # print a, b
  Y = a*np.power(X, b)
  plt.plot(X, Y, 'g-', label="a*X^b")
  ax.set_xlabel('Running Time')
  ax.set_ylabel('Size of Input')
  ax.set_title('Running Time Growth')
  ax.legend()
  plt.show()

# Algorithms, Part 1 by Kevin Wayne, Robert Sedgewick at Princeton University
# Lecture Week 1 Analysis of Algorithms Introduction (8:14)

# 07:19 SCIENTIFIC METHOD
#   * OBSERVE     some feature of the natural world (Computer itself)
#   * HYPOTHESIZE a model that is consistent with the observations
#   * PREDICT     events using the hypothesis
#   * VERIFY      the predictions by making further observations
#   * VALIDATE    by repeating until the hypothesis and observations agree.
#
# Principles:
#   * Experiments must be REPRODUCIBLE
#   * Hypothesis  must be FALSIFIABLE

#   Why is my program so slow?
#   Why does it run out of memory?
#
#   N^2 unacceptable because it does not scale
#   N*log(N) is almost linear
#
# QUESTION: Suppose that N == 1,000,000. Approximately how much faster is
# an algorithm that performs N*lg(N) operations versus one that performs
# N^2 operations? Recall that lg is log_2.
#   N^2     = (10^6)^2  =   10^12
#   N*lg(N) = 10^6*(20) = 2*10^7
#
# N^2        10^12
# ------ = ------- =  .5(10^5) = 5(10^-1)(10^5) = 5(10^4) = 5*10000 = 50,000
# N lg N   2*10^7

# Week 1 Observations (10:05)
# Power Law: a*N^b (06:22)
#   Math for straight-line in log-log plots with Power Law:
#     lg(T(N)) = b*lg(N) + c
#     T(N) = a*N^b, where a = 2^c 

def plot_rels(data, labels=None, colors=None, outfile="rels", latent=None, alpha=0.8):
    ns, n = data.shape
    if labels is None:
        labels = list(map(str, range(n)))
    ncol = 5
    # ncol = 4
    nrow = int(np.ceil(float(n * (n - 1) / 2) / ncol))
    #nrow=1
    #pylab.rcParams.update({'figure.autolayout': True})
    fig, axs = pylab.subplots(nrow, ncol)
    fig.set_size_inches(5 * ncol, 5 * nrow)
    #fig.set_canvas(pylab.gcf().canvas)
    pairs = list(combinations(range(n), 2))  #[:4]
    pairs = sorted(pairs, key=lambda q: q[0]**2+q[1]**2)  # Puts stronger relationships first
    if colors is not None:
        colors = (colors - np.min(colors)) / (np.max(colors) - np.min(colors)).clip(1e-7)

    for ax, pair in zip(axs.flat, pairs):
        if latent is None:
            ax.scatter(data[:, pair[0]], data[:, pair[1]], marker='.', edgecolors='none', alpha=alpha)
        else:
            # cs = 'rgbcmykrgbcmyk'
            markers = 'x+.o,<>^^<>,+x.'
            for j, ind in enumerate(np.unique(latent)):
                inds = (latent == ind)
                ax.scatter(data[inds, pair[0]], data[inds, pair[1]], c=colors[inds], cmap=pylab.get_cmap("jet"),
                           marker=markers[j], alpha=0.5, edgecolors='none', vmin=0, vmax=1)

        ax.set_xlabel(shorten(labels[pair[0]]))
        ax.set_ylabel(shorten(labels[pair[1]]))

    for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:
        ax.scatter(data[:, 0], data[:, 1], marker='.')

    pylab.rcParams['font.size'] = 12  #6
    pylab.draw()
    #fig.set_tight_layout(True)
    fig.tight_layout()
    for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:
        ax.set_visible(False)
    filename = outfile + '.png'
    if not os.path.exists(os.path.dirname(filename)):
        os.makedirs(os.path.dirname(filename))
    fig.savefig(outfile + '.png')  #df')
    pylab.close('all')
    return True 

def plot_rels(data, labels=None, colors=None, outfile="rels", latent=None, alpha=0.8, title=''):
    ns, n = data.shape
    if labels is None:
        labels = list(map(str, list(range(n))))
    ncol = 5
    nrow = int(np.ceil(float(n * (n - 1) / 2) / ncol))

    fig, axs = pylab.subplots(nrow, ncol)
    fig.set_size_inches(5 * ncol, 5 * nrow)
    pairs = list(combinations(list(range(n)), 2))
    if colors is not None:
        colors = (colors - np.min(colors)) / (np.max(colors) - np.min(colors))

    for ax, pair in zip(axs.flat, pairs):
        diff_x = max(data[:, pair[0]]) - min(data[:, pair[0]])
        diff_y = max(data[:, pair[1]]) - min(data[:, pair[1]])
        ax.set_xlim([min(data[:, pair[0]]) - 0.05 * diff_x, max(data[:, pair[0]]) + 0.05 * diff_x])
        ax.set_ylim([min(data[:, pair[1]]) - 0.05 * diff_y, max(data[:, pair[1]]) + 0.05 * diff_y])
        ax.scatter(data[:, pair[0]], data[:, pair[1]], c=colors, cmap=pylab.get_cmap("jet"),
                       marker='.', alpha=alpha, edgecolors='none', vmin=0, vmax=1)

        ax.set_xlabel(shorten(labels[pair[0]]))
        ax.set_ylabel(shorten(labels[pair[1]]))

    for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:
        ax.scatter(data[:, 0], data[:, 1], marker='.')

    fig.suptitle(title, fontsize=16)
    pylab.rcParams['font.size'] = 12  #6
    # pylab.draw()
    # fig.set_tight_layout(True)
    pylab.tight_layout()
    pylab.subplots_adjust(top=0.95)
    for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:
        ax.set_visible(False)
    filename = outfile + '.png'
    if not os.path.exists(os.path.dirname(filename)):
        os.makedirs(os.path.dirname(filename))
    fig.savefig(outfile + '.png')
    pylab.close('all')
    return True


# Hierarchical graph visualization utilities 

def generate_(self, inputs, mode='display', return_attend=False, return_all=False):
        # assert self.config['sample_stoch'], 'RNNLM sampling must be stochastic'
        # assert not self.config['sample_argmax'], 'RNNLM sampling cannot use argmax'

        args = dict(k=self.config['sample_beam'],
                    maxlen=self.config['max_len'],
                    stochastic=self.config['sample_stoch'] if mode == 'display' else None,
                    argmax=self.config['sample_argmax'] if mode == 'display' else None,
                    return_attend=return_attend)
        context, _, c_mask, _, Z, R = self.encoder.gtenc(inputs)
        # c_mask[0, 3] = c_mask[0, 3] * 0
        # L   = context.shape[1]
        # izz = np.concatenate([np.arange(3), np.asarray([1,2]), np.arange(3, L)])
        # context = context[:, izz, :]
        # c_mask  = c_mask[:, izz]
        # inputs  = inputs[:, izz]
        # context, _, c_mask, _ = self.encoder.encode(inputs)
        # import pylab as plt
        # # visualize_(plt.subplots(), Z[0][:, 300:], normal=False)
        # visualize_(plt.subplots(), context[0], normal=False)

        if 'explicit_loc' in self.config:
            if self.config['explicit_loc']:
                max_len = context.shape[1]
                expLoc  = np.eye(max_len, self.config['encode_max_len'], dtype='float32')[None, :, :]
                expLoc  = np.repeat(expLoc, context.shape[0], axis=0)
                context = np.concatenate([context, expLoc], axis=2)

        sample, score, ppp    = self.decoder.get_sample(context, c_mask, inputs, **args)
        if return_all:
            return sample, score, ppp

        if not args['stochastic']:
            score  = score / np.array([len(s) for s in sample])
            idz    = score.argmin()
            sample = sample[idz]
            score  = score.min()
            ppp    = ppp[idz]
        else:
            score /= float(len(sample))

        return sample, np.exp(score), ppp 

def analyse_(self, inputs, outputs, idx2word, inputs_unk=None, return_attend=False, name=None, display=False):
        def cut_zero(sample, idx2word, ppp=None, Lmax=None):
            if Lmax is None:
                Lmax = self.config['dec_voc_size']
            if ppp is None:
                if 0 not in sample:
                    return ['{}'.format(idx2word[w].encode('utf-8'))
                            if w < Lmax else '{}'.format(idx2word[inputs[w - Lmax]].encode('utf-8'))
                            for w in sample]

                return ['{}'.format(idx2word[w].encode('utf-8'))
                        if w < Lmax else '{}'.format(idx2word[inputs[w - Lmax]].encode('utf-8'))
                        for w in sample[:sample.index(0)]]
            else:
                if 0 not in sample:
                    return ['{0} ({1:1.1f})'.format(
                            idx2word[w].encode('utf-8'), p)
                            if w < Lmax
                            else '{0} ({1:1.1f})'.format(
                            idx2word[inputs[w - Lmax]].encode('utf-8'), p)
                            for w, p in zip(sample, ppp)]
                idz = sample.index(0)
                return ['{0} ({1:1.1f})'.format(
                        idx2word[w].encode('utf-8'), p)
                        if w < Lmax
                        else '{0} ({1:1.1f})'.format(
                        idx2word[inputs[w - Lmax]].encode('utf-8'), p)
                        for w, p in zip(sample[:idz], ppp[:idz])]

        if inputs_unk is None:
            result, _, ppp = self.generate_(inputs[None, :],
                                            return_attend=return_attend)
        else:
            result, _, ppp = self.generate_(inputs_unk[None, :],
                                            return_attend=return_attend)

        source = '{}'.format(' '.join(cut_zero(inputs.tolist(),  idx2word, Lmax=len(idx2word))))
        target = '{}'.format(' '.join(cut_zero(outputs.tolist(), idx2word, Lmax=len(idx2word))))
        decode = '{}'.format(' '.join(cut_zero(result, idx2word)))

        if display:
            print source
            print target
            print decode

            idz    = result.index(0)
            p1, p2 = [np.asarray(p) for p in zip(*ppp)]
            print p1.shape
            import pylab as plt
            # plt.rc('text', usetex=True)
            # plt.rc('font', family='serif')
            visualize_(plt.subplots(), 1 - p1[:idz, :].T, grid=True, name=name)
            visualize_(plt.subplots(), 1 - p2[:idz, :].T, name=name)

            # visualize_(plt.subplots(), 1 - np.mean(p2[:idz, :], axis=1, keepdims=True).T)
        return target == decode 

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 plot_read_stack(reads, refseq=None, by=None, cutoff=0, ax=None):
    """Plot read stack using coverage at each position"""

    reads = reads.copy()
    if by != None:
        reads = reads.sort_values(by)
    else:
        reads = reads.sort_values(by=['start','reads'],ascending=[1,0])
    if refseq != None:
        seqlen = len(refseq)
    else:
        offset = reads.start.min()
        reads.start = reads.start-offset
        reads.end = reads.end-offset
        seqlen = reads.end.max()
    reads = reads[reads.reads>cutoff]
    if len(reads)==0:
        return

    def pos_coverage(r, p):
        x = [r.reads if (i>=r.start and i<=r.end) else 0 for i in p]
        return pd.Series(x,index=p)

    reads = reads.set_index('reads',drop=False)
    p = range(1,seqlen+1)
    m = reads.apply( lambda x: pos_coverage(x,p), 1 )
    m = m.replace(0,1)
    import matplotlib
    matplotlib.use('agg')
    import pylab as plt
    import seaborn as sns
    from matplotlib.colors import LogNorm
    if ax == None:
        h = 12*len(m)/80+1
        fig,ax = plt.subplots(1,1,figsize=(12,h))
    ax = sns.heatmap(m,ax=ax,norm=LogNorm(vmin=m.min(), vmax=m.max()),
                cmap='Blues',cbar=False)
    plt.gca().xaxis.grid(True)
    start, end = ax.get_xlim()
    xlbls = np.arange(start, end, 5).astype(int)
    ax.xaxis.set_ticks(xlbls)
    ax.set_xticklabels(xlbls)
    ax.tick_params(axis='y', which='major', labelsize=9)
    for xmin in ax.xaxis.get_majorticklocs():
        ax.axvline(x=xmin,ls='--',lw=0.5,color='gray')
    #m.plot(kind='bar',width=.8,figsize=(12,2))
    return ax 

def analyse_(self, inputs, outputs, idx2word, inputs_unk=None, return_attend=False, name=None, display=False):
        def cut_zero(sample, idx2word, ppp=None, Lmax=None):
            if Lmax is None:
                Lmax = self.config['dec_voc_size']
            if ppp is None:
                if 0 not in sample:
                    return ['{}'.format(idx2word[w].encode('utf-8'))
                            if w < Lmax else '{}'.format(idx2word[inputs[w - Lmax]].encode('utf-8'))
                            for w in sample]

                return ['{}'.format(idx2word[w].encode('utf-8'))
                        if w < Lmax else '{}'.format(idx2word[inputs[w - Lmax]].encode('utf-8'))
                        for w in sample[:sample.index(0)]]
            else:
                if 0 not in sample:
                    return ['{0} ({1:1.1f})'.format(
                            idx2word[w].encode('utf-8'), p)
                            if w < Lmax
                            else '{0} ({1:1.1f})'.format(
                            idx2word[inputs[w - Lmax]].encode('utf-8'), p)
                            for w, p in zip(sample, ppp)]
                idz = sample.index(0)
                return ['{0} ({1:1.1f})'.format(
                        idx2word[w].encode('utf-8'), p)
                        if w < Lmax
                        else '{0} ({1:1.1f})'.format(
                        idx2word[inputs[w - Lmax]].encode('utf-8'), p)
                        for w, p in zip(sample[:idz], ppp[:idz])]

        if inputs_unk is None:
            result, _, ppp = self.generate_(inputs[None, :],
                                            return_attend=return_attend)
        else:
            result, _, ppp = self.generate_(inputs_unk[None, :],
                                            return_attend=return_attend)

        source = '{}'.format(' '.join(cut_zero(inputs.tolist(),  idx2word, Lmax=len(idx2word))))
        target = '{}'.format(' '.join(cut_zero(outputs.tolist(), idx2word, Lmax=len(idx2word))))
        decode = '{}'.format(' '.join(cut_zero(result, idx2word)))

        if display:
            print(source)
            print(target)
            print(decode)

            idz    = result.index(0)
            p1, p2 = [np.asarray(p) for p in zip(*ppp)]
            print(p1.shape)
            import pylab as plt
            # plt.rc('text', usetex=True)
            # plt.rc('font', family='serif')
            visualize_(plt.subplots(), 1 - p1[:idz, :].T, grid=True, name=name)
            visualize_(plt.subplots(), 1 - p2[:idz, :].T, name=name)

            # visualize_(plt.subplots(), 1 - np.mean(p2[:idz, :], axis=1, keepdims=True).T)
        return target == decode 

def generate_multiple(self, inputs, mode='display', return_attend=False, return_all=True, return_encoding=False):
        # assert self.config['sample_stoch'], 'RNNLM sampling must be stochastic'
        # assert not self.config['sample_argmax'], 'RNNLM sampling cannot use argmax'
        args = dict(k=self.config['sample_beam'],
                    maxlen=self.config['max_len'],
                    stochastic=self.config['sample_stoch'] if mode == 'display' else None,
                    argmax=self.config['sample_argmax'] if mode == 'display' else None,
                    return_attend=return_attend,
                    type=self.config['predict_type']
                    )
        '''
        Return the encoding of input.
            Similar to encoder.encode(), but gate values are returned as well
            I think only gtenc with attention
            default: with_context=False, return_sequence=True, return_embed=True
        '''

        """
        return
            context:  a list of vectors [nb_sample, max_len, 2*enc_hidden_dim], encoding of each time state (concatenate both forward and backward RNN)
            _:      embedding of text X [nb_sample, max_len, enc_embedd_dim]
            c_mask: mask, an array showing which elements in context are not 0 [nb_sample, max_len]
            _: encoding of end of X, seems not make sense for bidirectional model (head+tail) [nb_sample, 2*enc_hidden_dim]
            Z:  value of update gate, shape=(nb_sample, 1)
            R:  value of update gate, shape=(nb_sample, 1)
        but.. Z and R are not used here
        """
        context, _, c_mask, _, Z, R = self.encoder.gtenc(inputs)
        # c_mask[0, 3] = c_mask[0, 3] * 0
        # L   = context.shape[1]
        # izz = np.concatenate([np.arange(3), np.asarray([1,2]), np.arange(3, L)])
        # context = context[:, izz, :]
        # c_mask  = c_mask[:, izz]
        # inputs  = inputs[:, izz]
        # context, _, c_mask, _ = self.encoder.encode(inputs)
        # import pylab as plt
        # # visualize_(plt.subplots(), Z[0][:, 300:], normal=False)
        # visualize_(plt.subplots(), context[0], normal=False)

        if 'explicit_loc' in self.config: # no
            if self.config['explicit_loc']:
                max_len = context.shape[1]
                expLoc  = np.eye(max_len, self.config['encode_max_len'], dtype='float32')[None, :, :]
                expLoc  = np.repeat(expLoc, context.shape[0], axis=0)
                context = np.concatenate([context, expLoc], axis=2)

        sample, score, ppp, output_encoding    = self.decoder.get_sample(context, c_mask, inputs, **args)
        if return_all:
            if return_encoding:
                return context, sample, score, output_encoding
            else:
                return sample, score
        return sample, score 

def generate_(self, inputs, mode='display', return_attend=False, return_all=False):
        # assert self.config['sample_stoch'], 'RNNLM sampling must be stochastic'
        # assert not self.config['sample_argmax'], 'RNNLM sampling cannot use argmax'

        args = dict(k=self.config['sample_beam'],
                    maxlen=self.config['max_len'],
                    stochastic=self.config['sample_stoch'] if mode == 'display' else None,
                    argmax=self.config['sample_argmax'] if mode == 'display' else None,
                    return_attend=return_attend)
        context, _, c_mask, _, Z, R = self.encoder.gtenc(inputs)
        # c_mask[0, 3] = c_mask[0, 3] * 0
        # L   = context.shape[1]
        # izz = np.concatenate([np.arange(3), np.asarray([1,2]), np.arange(3, L)])
        # context = context[:, izz, :]
        # c_mask  = c_mask[:, izz]
        # inputs  = inputs[:, izz]
        # context, _, c_mask, _ = self.encoder.encode(inputs)
        # import pylab as plt
        # # visualize_(plt.subplots(), Z[0][:, 300:], normal=False)
        # visualize_(plt.subplots(), context[0], normal=False)

        if 'explicit_loc' in self.config:
            if self.config['explicit_loc']:
                max_len = context.shape[1]
                expLoc  = np.eye(max_len, self.config['encode_max_len'], dtype='float32')[None, :, :]
                expLoc  = np.repeat(expLoc, context.shape[0], axis=0)
                context = np.concatenate([context, expLoc], axis=2)

        sample, score, ppp, _    = self.decoder.get_sample(context, c_mask, inputs, **args)
        if return_all:
            return sample, score, ppp

        if not args['stochastic']:
            score  = score / np.array([len(s) for s in sample])
            idz    = score.argmin()
            sample = sample[idz]
            score  = score.min()
            ppp    = ppp[idz]
        else:
            score /= float(len(sample))

        return sample, np.exp(score), ppp 

def plot_one_dir(args, directory):
    """ The actual plotting code.

    Assumes that we'll be plotting from one directory, which usually means
    considering one random seed only, however it's better to have multiple
    random seeds so this code generalizes. For ES, we should store the output at
    *every* timestep, so A['TotalIterations'] should be like np.arange(...), but
    this generalizes in case Ray can help me run for many more iterations.
    """
    print("Now plotting based on directory {} ...".format(directory))

    ### Figure 1: The log.txt file.
    num = len(ATTRIBUTES)
    fig, axes = subplots(num, figsize=(12,3*num))
    for (dd, cc) in zip(directory, COLORS):
        A = np.genfromtxt(join(args.expdir, dd, 'log.txt'),
                          delimiter='\t', dtype=None, names=True)
        x = A['TotalIterations']
        for (i,attr) in enumerate(ATTRIBUTES):
            axes[i].plot(x, A[attr], '-', lw=lw, color=cc, label=dd)
            axes[i].set_ylabel(attr, fontsize=ysize)
            axes[i].tick_params(axis='x', labelsize=tick_size)
            axes[i].tick_params(axis='y', labelsize=tick_size)
            axes[i].legend(loc='best', ncol=1, prop={'size':legend_size})
    plt.tight_layout()
    plt.savefig(args.out+'_log.png')

    ### Figure 2: Error regions.
    num = len(directory)
    if num == 1: 
        num+= 1
    fig, axes = subplots(1,num, figsize=(12*num,10))
    for (i, (dd, cc)) in enumerate(zip(directory, COLORS)):
        A = np.genfromtxt(join(args.expdir, dd, 'log.txt'),
                          delimiter='\t', dtype=None, names=True)
        axes[i].plot(A['TotalIterations'], A["FinalAvgReturns"], 
                     color=cc, marker='x', ms=ms, lw=lw)
        axes[i].fill_between(A['TotalIterations'],
                             A["FinalAvgReturns"] - A["FinalStdReturns"],
                             A["FinalAvgReturns"] + A["FinalStdReturns"],
                             alpha = error_region_alpha,
                             facecolor='y')
        axes[i].set_ylim(ENV_TO_YLABELS[args.envname])
        axes[i].tick_params(axis='x', labelsize=tick_size)
        axes[i].tick_params(axis='y', labelsize=tick_size)
        axes[i].set_title("Mean Episode Rewards ({})".format(dd), fontsize=title_size)
        axes[i].set_xlabel("ES Iterations", fontsize=xsize)
        axes[i].set_ylabel("Rewards", fontsize=ysize)
    plt.tight_layout()
    plt.savefig(args.out+'_rewards_std.png') 

def galaxy10_confusion(confusion_mat):
    """
    NAME:
        galaxy10_confusion
    PURPOSE:
        to plot confusion matrix
    INPUT:
        confusion_mat (ndarray): An integer 0-9
    OUTPUT:
        (string): Name of the class
    HISTORY:
        2018-Feb-11 - Written - Henry Leung (University of Toronto)
    """
    import pylab as plt

    conf_arr = confusion_mat.astype(int)

    norm_conf = []
    a = np.max(conf_arr)
    for i in conf_arr:
        tmp_arr = []
        for j in i:
            tmp_arr.append(float(j) / float(a))
        norm_conf.append(tmp_arr)

    fig, ax = plt.subplots(1, figsize=(10, 10.5), dpi=100)
    fig.suptitle("Confusion Matrix for Galaxy10 trained by astroNN", fontsize=18)
    ax.set_aspect(1)
    ax.imshow(np.array(norm_conf), cmap=plt.get_cmap('Blues'), interpolation='nearest')

    width, height = conf_arr.shape

    for x in range(width):
        for y in range(height):
            ax.annotate(str(conf_arr[x][y]), xy=(y, x),
                        horizontalalignment='center',
                        verticalalignment='center')

    alphabet = '0123456789'
    plt.xticks(range(width), alphabet[:width], fontsize=20)
    plt.yticks(range(height), alphabet[:height], fontsize=20)
    ax.set_ylabel('Prediction class by astroNN', fontsize=18)
    ax.set_xlabel('True class', fontsize=18)
    fig.tight_layout(rect=[0, 0.00, 0.8, 0.96])
    fig.show()

    return None