def plot_compare(t, u, v, fig_title='', file_name=''):
    """
    Compare two signals and plot the difference between them.

    Parameters
    ----------
    t : ndarray of floats
        Times (s) at which the signal is defined.
    u, v : ndarrays of floats
        Signal samples.
    fig_title : string
        Plot title.
    file_name : string
        File in which to save the plot.

    """

    p.clf()
    p.gcf().canvas.set_window_title(fig_title)
    p.subplot(211)
    p.plot(t, u, 'b', t, v, 'r')
    p.xlabel('t (s)')
    p.ylabel('u(t)')
    p.title(fig_title)
    p.gca().set_xlim(min(t), max(t))
    p.subplot(212)
    p.plot(t, 20*np.log10(abs(u-v)))
    p.xlabel('t (s)')
    p.ylabel('error (dB)')
    p.gca().set_xlim(min(t), max(t))
    p.draw_if_interactive()
    if file_name:
        p.savefig(file_name)

def plotNumVariants(depth):
    plt.figure(figsize=(7, 5), dpi=300)
    depth.apply(lambda x: x.dropna().size / 1000.0).unstack("POP").plot.bar(ax=plt.gca())
    plt.gcf().subplots_adjust(bottom=0.30)
    plt.ylabel(r"Num of variants ($\times$1000)")
    renameLegend(ax=plt.gca())
    plt.savefig(pathPlots + "numVariants.pdf")

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 compare_expvaluecollections(coll1, coll2, show=True, **kwargs):
    """
    Plot all subsystems of two ExpectationValueCollections.

    *Arguments*
        * *coll1*
            First :class:`pycppqed.expvalues.ExpectationValue.Collection`.

        * *coll2*
            Second :class:`pycppqed.expvalues.ExpectationValue.Collection`.

        * *show* (optional):
            If True pylab.show() is called finally. This means a plotting
            window will pop up automatically. (Default is True)

        * Any other arguments that the pylab plotting command can use.

    This function allows a fast comparison between two sets of expectation
    values that were obtained by different calculations.
    """
    import pylab
    s1 = coll1.subsystems
    s2 = coll2.subsystems
    assert len(s1) == len(s2)
    for i in range(len(s1)):
        pylab.figure()
        _compare_expvaluesubsystems(s1.values()[i], s2.values()[i],
                                        show=False, **kwargs)
        title = "%s vs. %s" % (s1.keys()[i], s2.keys()[i])
        if hasattr(pylab, "suptitle"): # For old versions not available.
            pylab.suptitle(title)
            pylab.gcf().canvas.set_window_title(title)
    if show:
        pylab.show()

def plot_shaded_lines(my_xticks, y1, y2, error1, error2, ylab, xlab, filename):
    plt.figure(figsize=(6,6))
    from matplotlib import rcParams
    rcParams.update({'figure.autolayout': True})

    x = range(0, len(y1))
    plt.plot(x, y1, 'k-', color="blue",  label='men')
    plt.fill_between(x, y1-error1, y1+error1, facecolor='blue', alpha=.2)

    plt.plot(x, y2, 'k-', color="red",  label='women')
    plt.fill_between(x, y2-error2, y2+error2, facecolor='red', alpha=.2)

    #if isinstance(x, (int, long, float, complex)):
    #    plt.xlim(np.min(x), np.max(x))
    plt.gcf().subplots_adjust(bottom=0.3)

    plt.xticks(x, my_xticks)
    plt.xticks(rotation=70, fontsize=14)
    plt.yticks(fontsize=14)
    #plt.setp(ax.get_xticklabels(), rotation='vertical', fontsize=14)
    plt.ylabel(ylab, fontsize=14)
    plt.xlabel(xlab, fontsize=14)
    plt.legend()

    plt.savefig(filename)

 def plot(self, key, graph_data, num, row, col, ax):
     pylab.subplots_adjust(wspace=.5, hspace=.5)
     pylab.title("Region %s" % str(key))
     pylab.setp(ax.get_yticklabels(), visible=False)
     pylab.setp(ax.get_xticklabels(), visible=False)
     pylab.gcf().set_size_inches((row+col, row+col))
     h = pylab.hist(graph_data, bins=self.bins)

def cycles_vs_align_len(classData, blastPath, floworder):
    pylab.clf()         
    l = []
    c = []
    for rowcol,leng,dif,seq in classData:
        t = template.deseqify(remove_spaces(seq), floworder)
        l.append(dif)
        c.append(len(t)/4)
    pylab.scatter(l,c)
    if "/" in blastPath:
        p = blastPath.strip().split("/")[-1]
        path = p.strip().split(".")[0]
    else:
        path = blastPath.strip().split(",")[0]
    if 'blastn' in path:
        pretitle = 'Relaxed Blast'
    if 'megablast' in path:
        pretitle = 'Strict Blast'
    pylab.title("%s Number of Flows vs. Alignment Length" % pretitle)
    pylab.xlabel("Length")
    pylab.ylabel("Flows")
    pylab.gcf().set_size_inches((8,4))
    pylab.axis([0,max(l)+10,0,max(c)+10])
    #x = [2.2*i for i in range(max(l))]
    #pylab.plot(x)
    pylab.savefig("alignment_vs_cycles_%s.png" % path)
    pylab.clf()

    def after_step(self, rbm, trainer, i):
        it = i + 1

        save = it in self.expt.save_after
        display = it in self.expt.show_after

        if save:
            if self.expt.save_particles:
                storage.dump(trainer.fantasy_particles, self.expt.pcd_particles_file(it))
            storage.dump(rbm, self.expt.rbm_file(it))
            if hasattr(trainer, 'avg_rbm'):
                storage.dump(trainer.avg_rbm, self.expt.avg_rbm_file(it))
            storage.dump(time.time() - self.t0, self.expt.time_file(it))

        if 'particles' in self.subset and (save or display):
            fig = rbm_vis.show_particles(rbm, trainer.fantasy_particles, self.expt.dataset, display=display,
                                         figtitle='PCD particles ({} updates)'.format(it))
            if display:
                pylab.gcf().canvas.draw()
            if save:
                misc.save_image(fig, self.expt.pcd_particles_figure_file(it))

        if 'gibbs_chains' in self.subset and (save or display):
            fig = diagnostics.show_chains(rbm, trainer.fantasy_particles, self.expt.dataset, display=display,
                                          figtitle='Gibbs chains (iteration {})'.format(it))
            if save:
                misc.save_image(fig, self.expt.gibbs_chains_figure_file(it))

        if 'objective' in self.subset:
            self.log_prob_tracker.update(rbm, trainer.fantasy_particles)

        if display:
            pylab.gcf().canvas.draw()

def p_h(data):
    P.gcf().set_size_inches(10, 5)
    x = np.arange(1, len(data[0]) + 1, 1)
    for y in data:
        P.plot(x, y)
    P.ylabel("Number of heads - Number of tails")
    P.xlabel("Number of coin tosses")

def plot_collinearity(motifs, best_Z):
    """Plot the cooccurrences of motifs.
    """
    import scipy.cluster.hierarchy as hier
    # from scipy.stats import pearsonr
    M = len(motifs)
    cooccurrences = numpy.ones((M, M))
    for m1 in xrange(M):
        for m2 in xrange(M):
            # both = sum(numpy.logical_and(m1seqs, m2seqs))
            # cooccurrences[m1,m2] = both/float(sum(m2seqs))
            cooccurrences[m1, m2] = \
                numpy.sqrt(sum(best_Z[m1] * best_Z[m2])) \
                / numpy.linalg.norm(best_Z[m2])
            # rho, _ = pearsonr(best_Z[m1], best_Z[m2])
            # cooccurrences[m1, m2] = rho
    Y = hier.centroid(cooccurrences)
    index = hier.fcluster(Y, -1) - 1
    cooccurrences = cooccurrences[index, :]
    cooccurrences = cooccurrences[:, index]
    pylab.pcolor(cooccurrences)
    pylab.colorbar()
    ax = pylab.gca()
    ax.set_xticks([])
    # ax.set_xticks(.5 + numpy.arange(M))
    # ax.set_xticklabels(motifs)
    ax.set_yticks(.5 + numpy.arange(M))
    ax.set_yticklabels(numpy.asarray(motifs)[index])
    ax.set_xlim((0, M))
    ax.set_ylim((0, M))
    for line in ax.yaxis.get_ticklines():
        line.set_markersize(0)
    pylab.gcf().subplots_adjust(left=.27, bottom=.02, top=.98, right=.99)

def analyze_year():
    fname = path + "/stretch/data_stretching_oneref_%s_%03d_5swin.npz"
    k = 0
    corrs = []
    datess = []
    stretchs = []
    while k <= 60:
        if not os.path.exists(fname % (station, k)):
            k += 1
            continue
        print "k", k
        npzfile = np.load(fname % (station, k))
        corr, stretch, dates = npzfile["corr"], npzfile["stretch"], npzfile["dates"]
        corr = corr[:, 1]
        stretch = stretch[:, 1]
        corrs.append(corr)
        datess.append(dates)
        stretchs.append(stretch)
        k += 1
    datess = np.hstack(datess)
    corrs = np.hstack(corrs)
    stretchs = np.hstack(stretchs)
    stretchs = smooth2(stretchs, 48 * 3)
    # stretchs, datess = smooth(stretchs, datess, 48)
    plt.plot(datess, stretchs)
    plt.gcf().autofmt_xdate()
    xfmt = mpl.dates.DateFormatter("%y-%m-%d %H:%M")
    plt.gca().xaxis.set_major_formatter(xfmt)
    plt.savefig(path + "/plots/year_%s_TIQ_5swin.pdf" % station)
    plt.show()

def plotDepth():
    sns.set_style("whitegrid", {"grid.color": "1", 'axes.linewidth': .5, "grid.linewidth": ".09"})
    sns.set_context("notebook", font_scale=1.4, rc={"lines.linewidth": 2.5})
    d = pd.read_pickle(utl.outpath + 'real/CD.F59.df').xs('D', level='READ', axis=1)
    (d.min(1) > 50).sum()

    (d > 50).sum().sum()

    z = pd.Series(np.ndarray.flatten(d.values))
    fontsize = 6
    mpl.rcParams.update({'font.size': fontsize})
    plt.figure(figsize=(6, 4), dpi=300);
    plt.subplot(2, 2, 1);
    z.value_counts().sort_index().plot()
    plt.xlim([0, 200]);
    plt.xlabel('Depth');
    plt.ylabel('Number of Measurments' + '\n (out of {:.1f}M)'.format(z.shape[0] / 1e6));
    plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 0))
    plt.title('Scaled PDF')
    pplt.annotate('(A)', xpad=0.85, ypad=0.45, fontsize=fontsize)
    plt.axvline(50, linestyle='--', linewidth=1, color='k')
    pplt.setSize(plt.gca(), fontsize)
    plt.subplot(2, 2, 2);

    z.value_counts().sort_index().cumsum().plot()
    plt.xlim([0, 200])
    plt.ylim([-3e5, 2.05 * 1e7])
    plt.xlabel('Depth');
    plt.title('Scaled CDF')
    pplt.annotate('(B)', xpad=0.85, ypad=0.45, fontsize=fontsize)
    plt.axvline(50, linestyle='--', linewidth=1, color='k')
    pplt.setSize(plt.gca(), fontsize)
    plt.subplot(2, 2, 3);
    d.min(1).value_counts().sort_index().plot()
    plt.xlim([0, 100]);
    plt.xlabel('Minimum Depth of each Variant');
    plt.ylabel('Number of Variants' + '\n (out of {:.1f}M)'.format(d.shape[0] / 1e6));
    plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 0))
    plt.rc('font', size=fontsize)
    pplt.annotate('(C)', xpad=0.85, ypad=0.45, fontsize=fontsize)
    plt.axvline(50, linestyle='--', linewidth=1, color='k')
    pplt.setSize(plt.gca(), fontsize)
    plt.subplot(2, 2, 4);
    d.min(1).value_counts().sort_index().cumsum().plot()
    plt.xlim([0, 60])
    plt.ylim([0.25 * -1e5, plt.ylim()[1]])
    plt.xlabel('Minimum Depth of each Variant');
    plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 0))
    pplt.annotate('(D)', xpad=0.85, ypad=0.45, fontsize=fontsize)
    plt.axvline(50, linestyle='--', linewidth=1, color='k')
    pplt.setSize(plt.gca(), fontsize)
    plt.gcf().subplots_adjust(bottom=0.15)
    plt.gcf().tight_layout(h_pad=0.1)
    fontsize = 6
    mpl.rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size': fontsize});
    mpl.rc('text', usetex=True)
    mpl.rcParams.update({'font.size': 1})

    pplt.savefig('depth', 300)
    plt.show()

def plotDepthHeterogenocity():
    dpi = 300
    sns.set_style("whitegrid", {"grid.color": "0.9", 'axes.linewidth': .5, "grid.linewidth": ".09"})
    _, ax = plt.subplots(2, 2, sharex=True, figsize=(6, 4), dpi=dpi)

    d = pd.read_pickle(utl.outpath + 'real/CD.F59.df').xs('D', level='READ', axis=1)
    std = d.std(1)
    loc = [std.idxmax(), (std == std.quantile(0.52)).replace({False: None}).dropna().index[0],
           (std == std.median()).replace({False: None}).dropna().index[-1],
           (std == std.quantile(0.8)).replace({False: None}).dropna().index[0]]
    ax = ax.reshape(-1)
    fontsize = 6
    for i, pos in enumerate(loc):
        eg = d.loc[pos]
        [eg[r].dropna().plot(marker='o', ax=ax[i], markersize=5) for r in range(3)];
        plt.xticks(d.columns.get_level_values('GEN').unique());
        plt.xlabel('');
        plt.ylabel('')
        print 'position={}:{}'.format(eg.name[0], eg.name[1]), get_axis_limits()

        if i in [0, 2]: ax[i].set_ylabel('Read Depth')
        if i > 1: ax[i].set_xlabel('Generation')
        if i == 0: ax[i].legend(['Replicate 1', 'Replicate 2', 'Replicate 3'], loc='upper center',
                                prop={'size': fontsize})
        yrang = pplt.get_axis_limits(upper=True, ax=ax[i])[1] - pplt.get_axis_limits(upper=False, ax=ax[i])[1]
        ax[i].set_ylim([min(0, ax[i].get_ylim()[0] - 0.05 * yrang), ax[i].get_ylim()[1] + 0.03 * yrang])
        ax[i].set_xlim([-2, 61]);
        ax[i].set_title('{}:{}'.format(eg.name[0], eg.name[1]))
        pplt.setSize(ax[i], fontsize)

    mpl.rc('font', **{'family': 'serif', 'serif': ['Computer Modern']});
    mpl.rc('text', usetex=True)
    plt.gcf().subplots_adjust(bottom=0.15)
    pplt.savefig('depthHetero', dpi)
    plt.show()

def plot_signal(t, u, fig_title='', file_name=''):
    """
    Plot a signal.

    Parameters
    ----------
    t : ndarray of floats
        Times (in s) at which the signal is defined.
    u : ndarray of floats
        Signal samples.
    fig_title : string
        Plot title.
    file_name : string
        File in which to save the plot.

    """

    p.clf()

    # Set the plot window title:
    p.gcf().canvas.set_window_title(fig_title)
    p.plot(t, u)
    p.xlabel('t (s)')
    p.ylabel('u(t)')
    p.title(fig_title)

    # Make the x axis exactly correspond to the plotted signal's time range:
    p.gca().set_xlim(min(t), max(t)) 
    p.draw_if_interactive()
    if file_name:
        p.savefig(file_name)

def demo_envelopes():

    for i, demo in enumerate(demos):
        try:
            prev = subplot(len(demos), 1, i+1, sharex=prev, sharey=prev)
        except NameError:
            prev = subplot(len(demos), 1, i+1)
        title(repr(demo).strip('{}').replace("'", ""))
        prev.yaxis.set_visible(False)
        prev.xaxis.set_visible(False)

        # Generate carrier of appropriate length for the trial
        #carrier = np.random.normal(size=demo['trial_dur']*fs)
        carrier = np.random.uniform(low=-1, high=1, size=demo['trial_dur']*fs)
        t, waveform, trial_trigs, set_trigs = generate_waveform(carrier, **demo)
        #fill_between(t, waveform, -waveform)
        plot(t, waveform, 'k')
        for trig in set_trigs:
            axvline(trig/fs, color='b', lw=5)
        #for trig in trial_trigs:
        #    axvline(trig/fs, color='b', lw=2.5)
        trig_times = np.true_divide(trial_trigs, fs)
        plot(trig_times, np.ones(len(trig_times)), 'ro')
    prev.xaxis.set_visible(True)
    gcf().subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.95,
                          hspace=0.5)
    axis(xmin=-1, ymin=-1.1, ymax=1.1)
    show()

def show_barlines(page):
    import pylab
    for i, barlines in enumerate(page.barlines):
        sd = page.staves.staff_dist[i]
        for j, barline_range in enumerate(barlines):
            barline_x = int(barline_range.mean())
            staff_y = page.staves.staff_y(i, barline_x)
            repeats = page.repeats[i][j]
            if repeats:
                # Draw thick bar
                pylab.fill_between([barline_x - sd/4,
                                    barline_x + sd/4],
                                   staff_y - sd*2,
                                   staff_y + sd*2,
                                   color='g')
                for letter, sign in (('L', -1), ('R', +1)):
                    if letter in repeats:
                        # Draw thin bar
                        bar_x = barline_x + sign * sd/2
                        pylab.plot([bar_x, bar_x],
                                   [staff_y - sd*2,
                                    staff_y + sd*2],
                                   color='g')
                        for y in (-1, +1):
                            circ = pylab.Circle((bar_x + sign*sd/2,
                                                 staff_y + y*sd/2),
                                                sd/4,
                                                color='g')
                            pylab.gcf().gca().add_artist(circ)
            else:
                pylab.plot([barline_x, barline_x],
                           [staff_y - sd*2,
                            staff_y + sd*2],
                           color='g')

def timeplot(*args,**kwargs):
    import matplotlib.pyplot as plt 
    import matplotlib.ticker 
    
    def HMSFormatter(value, loc): 
        h = value // 3600 
        m = (value - h * 3600) // 60 
        s = value % 60 
        
        return "%02d:%02d:%02d" % (h,m,s) 

    def HMSFormatter2(value, loc): 
        h = value // 3600 
        m = (value - h * 3600) // 60 
        s = value % 60 
        ms=value%1
        
        return "%02d:%02d.%03d" % (m,s,ms*1000) 


    
    t=args[0]
    if max(t)<10:  # use ms
        gca().xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(HMSFormatter2)) 
    else:
        gca().xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(HMSFormatter)) 
    
    
    
    plot(*args,**kwargs) 
    gcf().autofmt_xdate()

def testeps(d):
    M.clf()
    M.pcolor(d)
    M.axis('tight')
    M.colorbar()
    M.gcf().set_size_inches((7.5,6.))
    M.savefig('test.png',dpi=240)

def show_chains(rbm, state, dataset, num_particles=20, num_samples=20, show_every=10, display=True,
                figname='Gibbs chains', figtitle='Gibbs chains'):
    samples = gnp.zeros((num_particles, num_samples, state.v.shape[1]))
    state = state[:num_particles, :, :]

    for i in range(num_samples):
        samples[:, i, :] = rbm.vis_expectations(state.h)
        
        for j in range(show_every):
            state = rbm.step(state)

    npix = dataset.num_rows * dataset.num_cols
    rows = [vm.hjoin([samples[i, j, :npix].reshape((dataset.num_rows, dataset.num_cols)).as_numpy_array()
                      for j in range(num_samples)],
                     normalize=False)
            for i in range(num_particles)]
    grid = vm.vjoin(rows, normalize=False)

    if display:
        pylab.figure(figname)
        pylab.matshow(grid, cmap='gray', fignum=False)
        pylab.title(figtitle)
        pylab.gcf().canvas.draw()

    return grid

def hpSNR(metrics, charType):
    for k, v in metrics.iteritems():
        title = str(k)
        e = v[charType].strip().split(",")
        hp = []
        snr = []
        for i in e:
            k = i.split(":")
            key = k[0].strip()
            value = k[1].strip()
            hp.append(key)
            snr.append(float(value))
        pylab.clf()
        pylab.xlabel("HP")
        pylab.ylabel("SNR")
        pylab.title("%s_%s" % (charType, title))
        pylab.gcf().set_size_inches((8, 4))
        xaxis = [i for i in range(len(hp))]
        ymax = int(max(snr)) + 1
        xlen = int(len(hp)) + 1
        pylab.axis([0, xlen, 0, ymax])
        pylab.xticks(range(len(hp)))
        pylab.bar(xaxis, snr, width=0.5, linewidth=0, color="r", alpha=0.5, align="center")

        # pylab.legend()
        pylab.savefig(path.join(os.getcwd(), "%s_%s.png" % (charType, title)))
        pylab.clf()