def plot(self, marker='o', color='red', m=0, M=6000):
        """Plots the position / height of the peaks in the well

        .. plot::
            :include-source:

            from fragment_analyser import Line, fa_data
            l = Line(fa_data("alternate/peaktable.csv"))
            well = l.wells[0]
            well.plot()

        """
        import pylab
        if len(self.df) == 0:
            print("Nothing to plot (no peaks)")
            return
        x = self.df['Size (bp)'].astype(float).values
        y = self.df['RFU'].astype(float).values

        pylab.stem(x, y, marker=marker, color=color)
        pylab.semilogx()
        pylab.xlim([1, M])
        pylab.ylim([0, max(y)*1.2])
        pylab.grid(True)
        pylab.xlabel("size (bp)")
        pylab.ylabel("RFU")
        return x, y

 def plot_weights(self):
     import pylab
     x, w = zip(*self.distribution)
     pylab.stem(x, 100*numpy.array(w))
     pylab.title('Weight distribution')
     pylab.xlabel(self.P.name)
     pylab.ylabel('Percentage')

def display_fit(xs, target, method, color=None, acronym=None, xmax=None, linewidth=2):
    if not hasattr(method, "__iter__"):
        method = (method,)
        if not color is None:
            color = (color,)
        if not acronym is None:
            acronym = (acronym,)
    if not xs.shape[0] == 1:
        raise ValueError("univariate data expected")
    if xmax == None:
        xmax = int(np.max(np.abs(xs.squeeze()))) + 1
    x = np.linspace(-xmax, xmax, 2 * xmax / 0.01)
    x = np.reshape(x, (1, x.size))
    fits = [m.fit(x) for m in method]
    if color is None:
        for fit in fits:
            plt.plot(x.squeeze(), fit, linewidth=linewidth)
    else:
        for fit, col in zip(fits, color):
            plt.plot(x.squeeze(), fit, col, linewidth=linewidth)
    if not acronym is None:
        plt.legend(acronym)
    target_xs = np.exp(target(xs.squeeze()))
    target_x = np.exp(target(x.squeeze()))
    plt.stem(xs.squeeze(), target_xs, linefmt="k-", markerfmt="ko", basefmt="k-")
    plt.plot(x.squeeze(), target_x, "k")
    x0, x1, y0, y1 = plt.axis()
    plt.plot((x0, x1), (0, 0), "k")
    plt.show()

def matlab_fourier_anal(data, sample_rate = 1.0):
 
  n = data.shape[0] 
  t = np.arange(n)

  P.figure()
  P.plot(t/sample_rate,data)
  P.xlabel('Time Units')
  P.ylabel('Data Magnitude')

  Y = fftpack.fft(data);
  freqs = np.arange(n/2.0) * sample_rate / n  
  power = np.abs(Y[:n/2.0])

  P.figure()
  markerline, stemlines, baseline = P.stem(freqs, power, '--')
  P.setp(markerline, 'markerfacecolor', 'b')
  # P.setp(baseline, 'color','r', 'linewidth', 2)
  P.xlabel('Cycles/ Unit Time')
  P.ylabel('Power')
  P.title('Periodogram')

  period = 1. / freqs
  
  k = np.arange(50)
  f = k / float(n)
  power2 = power[k];

  P.figure()
  markerline, stemlines, baseline = P.stem(f, power2, '--')
  P.setp(markerline, 'markerfacecolor', 'b')
  P.setp(baseline, 'color','r', 'linewidth', 2)
  P.xlabel('Unit Time / Cycle')
  P.ylabel('Power')
  P.title('Periodogram: First 50 periods')

 def plot_reflection(self):
     from pylab import stem, title, xlabel, ylabel
     if self.reflection is not None:
         stem(list(range(0, len(self.reflection))), abs(self.reflection))
         title('Reflection coefficient evolution')
         xlabel('Order')
         ylabel('Reflection Coefficient absolute values')
     else:
         logging.warning("Reflection coefficients not available or not yet computed.")

def decode_symbols(r, i, corr_index, r0, i0, Nbits, fs=48000, baud=300,
        plot=False):
    Ns = fs/baud
    r = r[corr_index:]
    i = i[corr_index:]

    r0 = r0[corr_index:]
    i0 = i0[corr_index:]

   #print len(r)
   #print len(r0)
   #print len(i0)

    r = r/np.max(r)*2.2 #change this 2 depending on the input amplitude
    i = i/np.max(i)*2.2 #change this 2 depending on the input amplitude

    if plot:
        fig = plt.figure(figsize = (16,4))
        plt.plot(2*r0)
        plt.plot(r)
        plt.title('Real part, raw input and normalized')

        fig = plt.figure(figsize = (16,4))
        plt.plot(2*i0)
        plt.plot(i)
        plt.title('Imaginary part, raw input and normalized')

    ####Decode
    idx = np.r_[Ns/2:len(r):Ns]

    #r = np.around(r)
    #i = np.around(i)

    r_dec = np.around(r[idx])
    i_dec = np.around(i[idx])

    if plot:
        fig = plt.figure(figsize = (16,4))
        plt.plot(2*r0)
        plt.plot(r)
        plt.plot(np.around(r))
        plt.stem(idx, r_dec)
        plt.title('Real part, decoded by sampling values as indicated')

        fig = plt.figure(figsize = (16,4))
        plt.plot(2*i0)
        plt.plot(i)
        plt.plot(np.around(i))
        plt.stem(idx, i_dec)
        plt.title('Imaginary part, decoded by sampling values as indicated')
        fig = plt.figure(figsize = (8,8))
        plt.scatter(r0*2, i0*2, c='r')
        plt.scatter(r_dec, i_dec, c='g')

        plt.title('Constellation of input message vs decoded symbols')

    return r_dec + 1j*i_dec

def print_spectrum(spectrum, filename):
    x = spectrum.mz
    y = spectrum.intensity
    pylab.stem(x, y, markerfmt='b,')
    pylab.xlabel('m/z')
    pylab.ylabel('Transformed intensity')
    pylab.grid()
    pylab.xlim(0.95*min(x), 1.05*max(x))
    pylab.savefig(filename)

 def plot_reflection(self):
     from pylab import stem, title, xlabel, ylabel
     if self.reflection is not None:
         stem(list(range(0, len(self.reflection))), abs(self.reflection))
         title('Reflection coefficient evolution')
         xlabel('Order')
         ylabel('Reflection Coefficient absolute values')
     else:
         import warnings
         warnings.warn("""Reflection coefficients not available with
             the current method.""")

    def plot_spectra(self):
        """ Plots the original spectrum and the sum of cisoids.
        """

        import pylab as plt
        plt.plot(self.psd[0], self.psd[1])
        freqs = [el[0] for el in self.soc]
        ampls = [el[1] for el in self.soc]
        plt.stem(freqs, ampls)
        plt.xlabel(r'$f$ in Hz')
        plt.ylabel(r'$S(f)$')
        plt.show()

 def nzcols_hist(P,file=None):
     if file==None: pylab.ion()
     else: pylab.ioff()
     nzcols = +P.nzcols
     Nmax = max(nzcols)
     y = matrix(0,(Nmax,1))
     for n in nzcols:
         y[n-1] += 1
     y = sparse(y)
     f = pylab.figure(); f.clf()
     pylab.stem(y.I+1,y.V)#; pylab.xlim(-1,Nmax+2)
     if file: pylab.savefig(P._pname + "_nzcols_hist.png")

def fourier_analysis(sig, time_step = 1.0, top_freqs = 5, zoomed_num = 15):

  time_vec = np.arange(0, len(sig), time_step)
  
  sample_freq = fftpack.fftfreq(sig.size, d=time_step)
  
  sig_fft = fftpack.fft(sig)
  
  # Only take first half of sampling frequency 
  pidxs = np.where(sample_freq > 0)
  freqs, power = sample_freq[pidxs], np.abs(sig_fft)[pidxs]
  
  # plot with zoomed in sub plot
  P.figure()
  P.plot(freqs, power)
  P.ylabel('Power')
  P.xlabel('Frequency [Hz]')
  axes = P.axes([0.3, 0.3, 0.5, 0.5])
  P.title('Peak frequency')
  P.stem(freqs[:zoomed_num], power[:zoomed_num])
  #P.setp(axes, yticks=[])
  #P.savefig('source/fftpack-frequency.png')
  
  # Find top x frequencies to use in reconstruction 
  full_sort_idx = power.argsort()
  
  find_idx = full_sort_idx >= top_freqs 
  
  
  sorted_power = power[sort_idx][::-1] # Sort decending
  component_freqs = freqs[sort_idx[:top_freqs]]
  
  # copy fft
  reduced_sig_fft = sig_fft.copy()
  
  # set all values not in component freqs to 0 
  L = np.array(reduced_sig_fft, dtype = bool)
  L[sort_idx[:top_freqs]] = False  
  reduced_sig_fft[L] = 0
  
  # Reconstruct signal
  reconstruct_sig = fftpack.ifft(reduced_sig_fft)
  
  # Plot original and reconstructed signal 
  P.figure()
  P.plot(time_vec, sig)
  P.plot(time_vec, reconstruct_sig, linewidth=2)
  P.ylabel('Amplitude')
  P.xlabel('Time [s]')
  
  return sig_fft, reduced_sig_fft, reconstruct_sig, freqs, power, component_freqs

def demo(text=None):
    from nltk.corpus import brown
    import pylab
    tt=TextTilingTokenizer(demo_mode=True)
    if text is None: text=brown.raw()[:10000]
    s,ss,d,b=tt.tokenize(text)
    pylab.xlabel("Sentence Gap index")
    pylab.ylabel("Gap Scores")
    pylab.plot(range(len(s)), s, label="Gap Scores")
    pylab.plot(range(len(ss)), ss, label="Smoothed Gap scores")
    pylab.plot(range(len(d)), d, label="Depth scores")
    pylab.stem(range(len(b)),b)
    pylab.legend()
    pylab.show()

def plot_inflections(x,y):
    """
    Plot inflection points in the spline curve.
    """
    m = (y[2:]-y[:-2])/(x[2:]-x[:-2])
    b = y[2:] - m*x[2:]
    delta = y[1:-1] - (m*x[1:-1] + b)
    t = linspace(x[0],x[-1],400)
    import pylab
    ax1=pylab.subplot(211)
    pylab.plot(t,monospline(x,y,t),'-b',x,y,'ob')
    pylab.subplot(212, sharex=ax1)
    delta_x = x[1:-1]
    pylab.stem(delta_x,delta)
    pylab.plot(delta_x[delta<0],delta[delta<0],'og')
    pylab.axis([x[0],x[-1],min(min(delta),0),max(max(delta),0)])

def impz(b,a=1):
  impulse = np.repeat(0.,50); impulse[0] =1.
  x = np.arange(0,50)
  response = signal.lfilter(b,a,impulse)
  pl.subplot(211)
  pl.stem(x, response)
  pl.ylabel('Amplitude') 
  pl.xlabel(r'n (samples)')
  pl.title(r'Impulse response')
  pl.subplot(212)
  step = np.cumsum(response)
  pl.stem(x, step)
  pl.ylabel('Amplitude') 
  pl.xlabel(r'n (samples)')
  pl.title(r'Step response')
  pl.subplots_adjust(hspace=0.5)

def hist_mtx(mtx, tstr=''):
    """
    Given a piano-roll matrix, 128 MIDI piches x beats, plot the pitch class histogram
    """
    i_min, i_max = np.where(mtx.mean(1))[0][[0,-1]]
    P.figure(figsize=(14.5,8))    
    P.stem(np.arange(i_max+1-i_min),mtx[i_min:i_max+1,:].sum(1))
    ttl = 'Note Frequency'
    if tstr: ttl+=': '+tstr
    P.title(ttl,fontsize=16)
    t=P.xticks(np.arange(0,i_max+1-i_min,3),pc_labels[i_min:i_max+1:3],fontsize=14)
    P.xlabel('Pitch Class', fontsize=14)
    P.ylabel('Frequency', fontsize=14)
    ax = P.axis()
    P.axis(xmin=-0.5)
    P.grid()

def impz(b,a=1):
    l = len(b)
    impulse = numpy.repeat(0.,l); impulse[0] =1.
    x = numpy.arange(0,l)
    response = scipy.signal.lfilter(b,a,impulse)
    pylab.subplot(211)
    pylab.stem(x, response)
    pylab.ylabel('Amplitude')
    pylab.xlabel(r'n (samples)')
    pylab.title(r'Impulse response')
    pylab.subplot(212)
    step = numpy.cumsum(response)
    pylab.stem(x, step)
    pylab.ylabel('Amplitude')
    pylab.xlabel(r'n (samples)')
    pylab.title(r'Step response')
    pylab.subplots_adjust(hspace=0.5)

    def plot_spectra(self):
        """ Plots the original Spectrum and the sum of sinusoids.
        """

        import matplotlib as mp
        mp.rc('font', family='serif', size=22)
        import pylab as plt
        freqs = [f[0] for f in self.sos]
        negfreqs = [-f[0] for f in self.sos]
        # convert the amplitudes of the sinusoids to a psd
        ampls = [(a[1] / 2) ** 2 for a in self.sos]
        plt.stem(freqs, ampls)
        plt.stem(negfreqs, ampls)
        plt.xlim(1.5 * negfreqs[-1], 1.5 * freqs[-1])
        plt.ylim(0, 1.5 * max(ampls))
        plt.xlabel(r'$f$ in Hz')
        plt.ylabel(r'$\| S(f) \|$')
        plt.show()

def impz(b,a=1):
    '''
    #Plot step and impulse response
    from http://mpastell.com/2010/01/18/fir-with-scipy/
    '''
    l = len(b)
    impulse = pylab.repeat(0.,l); impulse[0] =1.
    x = numpy.arange(0,l)
    response = signal.lfilter(b,a,impulse)
    pylab.subplot(211)
    pylab.stem(x, response)
    pylab.ylabel('Amplitude')
    pylab.xlabel(r'n (samples)')
    pylab.title(r'Impulse response')
    pylab.subplot(212)
    step = numpy.cumsum(response)
    pylab.stem(x, step)
    pylab.ylabel('Amplitude')
    pylab.xlabel(r'n (samples)')
    pylab.title(r'Step response')
    pylab.subplots_adjust(hspace=0.5)

def test_volterra():
    f=array([-1200.,-1000.,1000.,1200.])
    v=array([2.,2.,2.,2.])
    w=Waveform(f,v)
    v1,v2,v31,v32=lna_volterra(w)
    pylab.subplot(4,1,1)
    pylab.stem(v1.get_x()[0],v1.get_y())
    pylab.subplot(4,1,2)
    pylab.stem(v2.get_x()[0],v2.get_y())
    pylab.subplot(4,1,3)
    pylab.stem(v31.get_x()[0],v31.get_y())
    pylab.subplot(4,1,4)
    pylab.stem(v32.get_x()[0],v32.get_y())
    pylab.show()

def _add_stems_to_plot(interval, stem_bed, samples, plot):
    stems = _get_stems_by_callers(stem_bed.tabix_intervals(interval))
    callers = sorted(stems.keys())
    caller_colormap = _get_caller_colormap(callers)
    caller_heights = _get_caller_heights(callers, plot)
    for caller in callers:
        stem_color = caller_colormap[caller]
        caller_stems = stems[caller]
        stem_heights = list(repeat(caller_heights[caller], len(caller_stems)))
        markerline, _, baseline = stem(caller_stems, stem_heights, "-.", label=caller)
        setp(markerline, "markerfacecolor", stem_color)
        setp(baseline, "color", "r", "linewidth", 0)
        plt.legend()