def plot_tracks(src, fakewcs, spa=None, **kwargs):
    # NOTE -- MAGIC 61 = monthly; this is ASSUMEd below.
    tt = np.linspace(2010., 2015., 61)
    t0 = TAITime(None, mjd=TAITime.mjd2k + 365.25*10)
    #rd0 = src.getPositionAtTime(t0)
    #print 'rd0:', rd0
    xx,yy = [],[]
    rr,dd = [],[]
    for t in tt:
        #print 'Time', t
        rd = src.getPositionAtTime(t0 + (t - 2010.)*365.25*24.*3600.)
        ra,dec = rd.ra, rd.dec
        rr.append(ra)
        dd.append(dec)
        ok,x,y = fakewcs.radec2pixelxy(ra,dec)
        xx.append(x - 1.)
        yy.append(y - 1.)

    if spa is None:
        spa = [None,None,None]
    for rows,cols,sub in spa:
        if sub is not None:
            plt.subplot(rows,cols,sub)
        ax = plt.axis()
        plt.plot(xx, yy, 'k-', **kwargs)
        plt.axis(ax)

    return rr,dd,tt

 def plot_vm(self):
     """Plot Vm for presynaptic compartment and soma - along with
     the same in NEURON simulation if possible."""
     pylab.subplot(211)
     pylab.title('Soma Vm')
     pylab.plot(self.tseries*1e3, self.somaVmTab.vec * 1e3,
                label='Vm (mV) - moose')
     pylab.plot(self.tseries*1e3, self.injectionTab.vec * 1e9,
                label='Stimulus (nA)')
     try:
         nrn_data = np.loadtxt('../nrn/data/%s_soma_Vm.dat' % \
                                   (self.celltype))
         nrn_indices = np.nonzero(nrn_data[:, 0] <= self.tseries[-1]*1e3)[0]                        
         pylab.plot(nrn_data[nrn_indices,0], nrn_data[nrn_indices,1], 
                    label='Vm (mV) - neuron')
     except IOError:
         print 'No neuron data found.'
     pylab.legend()
     pylab.subplot(212)
     pylab.title('Presynaptic Vm')
     pylab.plot(self.tseries*1e3, self.presynVmTab.vec * 1e3,
                label='Vm (mV) - moose')
     pylab.plot(self.tseries*1e3, self.injectionTab.vec * 1e9,
                label='Stimulus (nA)')
     try:
         nrn_data = np.loadtxt('../nrn/data/%s_presynaptic_Vm.dat' % \
                                   (self.celltype))
         nrn_indices = np.nonzero(nrn_data[:, 0] <= self.tseries[-1]*1e3)[0]
         pylab.plot(nrn_data[nrn_indices,0], nrn_data[nrn_indices,1], 
                    label='Vm (mV) - neuron')
     except IOError:
         print 'No neuron data found.'
     pylab.legend()
     pylab.show()

    def Xtest2(self):
        """
        Test from Kate Marvel
        As the following code snippet demonstrates, regridding a
        cdms2.tvariable.TransientVariable instance using regridTool='regrid2' 
        results in a new array that is masked everywhere.  regridTool='esmf' 
        and regridTool='libcf' both work as expected.

        This passes.
        """
        import cdms2 as cdms
        import numpy as np

        filename = cdat_info.get_sampledata_path() + '/clt.nc'
        a=cdms.open(filename)
        data=a('clt')[0,...]

        print data.mask #verify this data is not masked

        GRID= data.getGrid() # input = output grid, passes

        test_data=data.regrid(GRID,regridTool='regrid2')

        # check that the mask does not extend everywhere...
        self.assertNotEqual(test_data.mask.sum(), test_data.size)
        
        if PLOT:
            pylab.subplot(2, 1, 1)
            pylab.pcolor(data[...])
            pylab.title('data')
            pylab.subplot(2, 1, 2)
            pylab.pcolor(test_data[...])
            pylab.title('test_data (interpolated data)')
            pylab.show()

def trace(data, name, format='png', datarange=(None, None), suffix='', path='./', rows=1, columns=1, 
    num=1, last=True, fontmap = None, verbose=1):
    """
    Generates trace plot from an array of data.

    :Arguments:
        data: array or list
            Usually a trace from an MCMC sample.

        name: string
            The name of the trace.
            
        datarange: tuple or list
            Preferred y-range of trace (defaults to (None,None)).

        format (optional): string
            Graphic output format (defaults to png).

        suffix (optional): string
            Filename suffix.

        path (optional): string
            Specifies location for saving plots (defaults to local directory).
            
        fontmap (optional): dict
            Font map for plot.

    """

    if fontmap is None: fontmap = {1:10, 2:8, 3:6, 4:5, 5:4}

    # Stand-alone plot or subplot?
    standalone = rows==1 and columns==1 and num==1

    if standalone:
        if verbose>0:
            print_('Plotting', name)
        figure()

    subplot(rows, columns, num)
    pyplot(data.tolist())
    ylim(datarange)

    # Plot options
    title('\n\n   %s trace'%name, x=0., y=1., ha='left', va='top', fontsize='small')

    # Smaller tick labels
    tlabels = gca().get_xticklabels()
    setp(tlabels, 'fontsize', fontmap[rows/2])

    tlabels = gca().get_yticklabels()
    setp(tlabels, 'fontsize', fontmap[rows/2])

    if standalone:
        if not os.path.exists(path):
            os.mkdir(path)
        if not path.endswith('/'):
            path += '/'
        # Save to file
        savefig("%s%s%s.%s" % (path, name, suffix, format))

def plot_Barycenter(dataset_name, feat, unfeat, repo):

    if dataset_name==MNIST:
        _, _, test=get_data(dataset_name, repo, labels=True)
        xtest1,_,_, labels,_=test
    else:
        _, _, test=get_data(dataset_name, repo, labels=False)
        xtest1,_,_ =test
        labels=np.zeros((len(xtest1),))
    # get labels
    def bary_wdl2(index): return _bary_wdl2(index, xtest1, feat, unfeat)
    
    n=xtest1.shape[-1]
    
    num_class = (int)(max(labels)+1)
    barys=[bary_wdl2(np.where(labels==i)) for i in range(num_class)]
    pl.figure(1, (num_class, 1))
    for i in range(num_class):
        pl.subplot(1,10,1+i)
        pl.imshow(barys[i][0,0,:,:],cmap='Blues',interpolation='nearest')
        pl.xticks(())
        pl.yticks(())
        if i==0:
            pl.ylabel('DWE Bary.')
        if num_class >1:
            pl.title('{}'.format(i))
    pl.tight_layout(pad=0,h_pad=-2,w_pad=-2) 
    pl.savefig("imgs/{}_dwe_bary.pdf".format(dataset_name))

            def graphError(predictedPercentChanges, actualPercentChanges, title):
                # considering error and only considering it as error when the signs are different

                def computeSignedError(pred, actual):

                    if (pred > 0 and actual > 0) or (pred < 0 and actual < 0):
                        return 0

                    else:
                        error = abs(pred - actual)
                        # print 'pred: {0}, actual: {1}, error: {2}'.format(pred, actual, error)
                        return error

                signedError = map(
                    lambda pred, actual: computeSignedError(pred, actual), predictedPercentChanges, actualPercentChanges
                )
                pl.figure(2)
                pl.title(title + " Error")
                pl.subplot(211)
                pl.plot(signedError)
                pl.xlabel("Time step")
                pl.ylabel("Error (0 if signs are same and normal error if signs are different)")

                pl.figure(3)
                pl.title(title + " Actual vs Predictions")
                pl.subplot(211)
                pl.plot(
                    range(len(predictedPercentChanges)),
                    predictedPercentChanges,
                    "ro",
                    range(len(actualPercentChanges)),
                    actualPercentChanges,
                    "bs",
                )

def T2_cpmg_process(folder_to_process,plot='y'):
    """Given a folder of images will process cpmg data and return
    fitted T2 values and associated uncertainties"""
    data=img_roi_signal([folder_to_process],['EchoTime'])
    rois=data[0][0]
    TEs=data[2][0]
    mean_signal_mat=data[3]
    serr_signal_mat=data[4]
        
    if plot=='y':
        plt.figure()
    spin_echo_fits=[]
    for jj in np.arange(len(rois)-2):
        mean_sig=mean_signal_mat[0,jj,:]
        #serr_sig=serr_signal_mat[0,jj,:]
        mean_noise=np.mean(mean_signal_mat[0,-2,:])
        try:
            spin_echo_fit = SE_fit_new( np.array(TEs[0:]), mean_sig[0:], mean_noise, 'n' )
            if plot=='y':
                TE_full=np.arange(0,400,1)
                plt.subplot(4,4,jj+1)
                plt.plot(np.array(TEs[0:]), mean_sig[0:],'o')
                plt.plot(TE_full,spin_echo_fit(TE_full))
            
            spin_echo_fits.append(spin_echo_fit)
        except RuntimeError: 
            print 'RuntimeError'
            spin_echo=fitting.model('M0*exp(-x/T2)+a',{'M0':0,'T2':0,'a':0}) 
            spin_echo_fits.append(spin_echo)
    return spin_echo_fits   

def main():
    star = StarBinary(90.0, 0.5, nside=64, limb_law=-1, limb_coeff=[0.8])
    Flux0 = star.flux(0.0)
    phases = np.linspace(-0.5, 0.5, 100)

    flux1 = np.zeros(len(phases))
    for i in range(len(phases)):
        flux1[i] = star.flux(phases[i]) / Flux0

    # for tt in np.arange(0,360,80):
    # star.makeSpot(0,-45,10.,0.8)
    # star.makeSpot(45,+00,10.,0.8)
    # star.makeSpot(00, -90, 20., 0.)
    # for theta in range(0,360,45):
    #	star.makeSpot(theta,65,10.,0.8)

    # star.makeSpot(00,-90,20.,0.)
    # star.makeSpot(0,+45,10.,0.8)
    # star.makeSpot(270,-10.,10.,0.8)
    # star.makeSpot(180,-45.,10.,0.8)
    #	star.makeSpot(tt,65.,10.,0.5)

    flux2 = np.zeros(len(phases))
    for i in range(len(phases)):
        flux2[i] = star.flux(phases[i]) / Flux0

    H.mollview(star.I, sub=211, rot=(-90, 90))

    #ff = np.loadtxt('/tmp/cl.dat', unpack=True)

    py.subplot(212)
    py.plot(phases, flux1, '-')
    # py.plot(phases,flux2,'-')
    #py.plot(ff[0], ff[1], '.')
    py.show()

def sim_results(obs, modes, stars, model, data):


    synth = model.generate_data(modes)
    synth_stats = model.summary_stats(synth)
    obs_stats = model.summary_stats(obs)


    f = plt.figure(figsize=(15,3))
    plt.suptitle('Obs Cand.:{}; Sim Cand.:{}'.format(obs.size, synth.size))
    plt.rc('legend', fontsize='xx-small', frameon=False)
    plt.subplot(121)
    bins = opt_bin(obs_stats[0],synth_stats[0])
    plt.hist(obs_stats[0], bins=bins, histtype='step', label='Data', lw=2)
    plt.hist(synth_stats[0], bins=bins, histtype='step', label='Simulation', lw=2)
    plt.xlabel(r'$\xi$')
    plt.legend()

    plt.subplot(122)
    bins = opt_bin(obs_stats[1],synth_stats[1])
    plt.hist(obs_stats[1], bins=np.arange(bins.min()-0.5, bins.max()+1.5,
                                          1),
             histtype='step', label='Data', log=True, lw=2)
    plt.hist(synth_stats[1], bins=np.arange(bins.min()-0.5, bins.max()+1.5,
                                            1),
             histtype='step', label='Simulation', log=True, lw=2)
    plt.xlabel(r'$N_p$')
    plt.legend()

def graphSimpleResults(resultsDir):
    x=[]
    y=[]
    files=[open("%s/objectiveFunctionReport.txt" % resultsDir),
       open("%s/fitnessReport.txt" % resultsDir)]
    for f in files:
        x.append([])
        y.append([])
        i=len(x)-1
        for line in f:
            line=line.split(',')
            if line[0] != "gen":
                x[i].append(int(line[0]))
                y[i].append(float(line[1]))

    ylen=len(y[0])
    pl.subplot(2,1,1)
    pl.plot(x[0],y[0],'bo')
    pl.ylabel('Maximum x')
    pl.title('Maximizing x**2 with SGA')
    pl.annotate("{0:,}".format(y[0][0]),xy=(x[0][0],y[0][0]),  xycoords='data',
             xytext=(50, 30), textcoords='offset points',
             arrowprops=dict(arrowstyle="->") )
    pl.annotate("{0:,}".format(y[0][ylen-1]),xy=(x[0][ylen-1],y[0][ylen-1]),  xycoords='data',
             xytext=(-30, -30), textcoords='offset points',
             arrowprops=dict(arrowstyle="->") )

    pl.subplot(2,1,2)
    pl.plot(x[1],y[1],'go')
    pl.xlabel('Generation')
    pl.ylabel('Fitness')
    pl.savefig("%s/simple_result.png" % resultsDir)

def graphTSPResults(resultsDir,numberOfTours):
    x=[]
    y=[]
    files=[open("%s/objectiveFunctionReport.txt" % resultsDir),
       open("%s/fitnessReport.txt" % resultsDir)]
    for f in files:
        x.append([])
        y.append([])
        i=len(x)-1
        for line in f:
            line=line.split(',')
            if line[0] != "gen":
                x[i].append(int(line[0]))
                y[i].append(float(line[1] if i==1 else line[2]))
            
    ylen=len(y[0])
    pl.subplot(2,1,1)
    pl.plot(x[0],y[0],'bo')
    pl.ylabel('Minimum Distance')
    pl.title("TSP with a %s City Tour" % numberOfTours)
    pl.annotate("{0:,}".format(y[0][0]),xy=(x[0][0],y[0][0]),  xycoords='data',
             xytext=(30, -30), textcoords='offset points',
             arrowprops=dict(arrowstyle="->") )
    pl.annotate("{0:,}".format(y[0][ylen-1]),xy=(x[0][ylen-1],y[0][ylen-1]),  xycoords='data',
             xytext=(-30,30), textcoords='offset points',
             arrowprops=dict(arrowstyle="->") )

    pl.subplot(2,1,2)
    pl.plot(x[1],y[1],'go')
    pl.xlabel('Generation')
    pl.ylabel('Fitness')
    pl.savefig("%s/tsp_result.png" % resultsDir)
    pl.clf()

def main():
    src_cv_img_1 = cv.LoadImage("data/gorskaya/images/1/g126.jpg")
    src_cv_img_gray_1 = cv.LoadImage("data/gorskaya/images/1/g126.jpg", 
        cv.CV_LOAD_IMAGE_GRAYSCALE)

    src_cv_img_2 = cv.LoadImage("data/gorskaya/images/1/g127.jpg")
    src_cv_img_gray_2 = cv.LoadImage("data/gorskaya/images/1/g127.jpg", 
        cv.CV_LOAD_IMAGE_GRAYSCALE)

    (keypoints_1, descriptors_1) = \
        cv.ExtractSURF(src_cv_img_gray_1, None, cv.CreateMemStorage(), 
            (0, 30000, 3, 1))
    (keypoints_2, descriptors_2) = \
        cv.ExtractSURF(src_cv_img_gray_2, None, cv.CreateMemStorage(), 
            (0, 30000, 3, 1))

    print("Found {0} and {1} keypoints".format(
        len(keypoints_1), len(keypoints_2)))

    src_arr_1 = array(src_cv_img_1[:, :])[:, :, ::-1]
    src_arr_2 = array(src_cv_img_2[:, :])[:, :, ::-1]

    pylab.rc('image', interpolation='nearest')

    pylab.subplot(121)
    pylab.imshow(src_arr_1)
    pylab.plot(*zip(*[k[0] for k in keypoints_1]), 
        marker='.', color='r', ls='')

    pylab.subplot(122)
    pylab.imshow(src_arr_2)
    pylab.plot(*zip(*[k[0] for k in keypoints_2]), 
        marker='.', color='r', ls='')

    pylab.show()

def display_coeff(data=None):
    betaAll,betaErrAll, R2adjAll = measure_stamp_coeff(data = data, zernike_max_order=20)
    ind = np.arange(len(betaAll[0]))
    momname = ('M20','M22.Real','M22.imag','M31.real','M31.imag','M33.real','M33.imag')
    fmtarr = ['bo-','ro-','go-','co-','mo-','yo-','ko-']
    pl.figure(figsize=(17,13))
    for i in range(7):
        pl.subplot(7,1,i+1)
        pl.errorbar(ind,betaAll[i],yerr = betaErrAll[i],fmt=fmtarr[i])
        pl.grid()
        pl.xlim(-1,21)
        if i ==0:
            pl.ylim(-10,65)
        elif i ==1:
            pl.ylim(-5,6)
        elif i ==2:
            pl.ylim(-5,6)
        elif i == 3:
            pl.ylim(-0.1,0.1)
        elif i == 4:
            pl.ylim(-0.1,0.1)
        elif i ==5:
            pl.ylim(-100,100)
        elif i == 6:
            pl.ylim(-100,100)
        pl.xticks(ind,('','','','','','','','','','','','','','','','','','','',''))
        pl.ylabel(momname[i])
    pl.xticks(ind,('0','1','2','3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19'))
    pl.xlabel('Zernike Coefficients')
    return '--- done ! ----'

def traceplot(traces, thin, burn):
    '''
    Plot parameter estimates for different levels of the model
    into the same plots. Black lines are individual observers
    and red lines are mean estimates.
    '''
    variables = ['Slope1', 'Slope2', 'Offset', 'Split']
    for i, var in enumerate(variables):
        plt.subplot(2, 2, i + 1)
        vals = get_values(traces, var, thin, burn)
        dim = (vals.min() - vals.std(), vals.max() + vals.std())
        x = plt.linspace(*dim, num=1000)
        for v in vals.T:
            a = gaussian_kde(v)
            y = a.evaluate(x)
            y = y / y.max()
            plt.plot(x, y, 'k', alpha=.5)
        try:
            vals = get_values(traces, 'Mean_' + var, thin, burn)
            a = gaussian_kde(vals)
            y = a.evaluate(x)
            y = y / y.max()
            plt.plot(x, y, 'r', alpha=.75)
        except KeyError:
            pass
        plt.ylim([0, 1.1])
        plt.yticks([0])
        sns.despine(offset=5, trim=True)
        plt.title(var)

def display_orbit(input, amplitude=0.000001):
    import matplotlib.pyplot
    import pylab

    # Compare the power spectral density functions of the system with and
    #  without the input sequence.
    x = [0.4, 0.6]
    X = []
    for i in range(warmups):
        x = network(x)
        x[0] = x[0] + ( amplitude * input[0][i % len(input[0])] )
        x[1] = x[1] + ( amplitude * input[1][i % len(input[1])] )
    for i in range(measure):
        X.append(x[0])
        x = network(x)
        x[0] = x[0] + ( amplitude * input[0][i % len(input[0])] )
        x[1] = x[1] + ( amplitude * input[1][i % len(input[1])] )
    pylab.subplot(2,1,1)
    matplotlib.pyplot.psd(X,1024,32)

    x = [0.4, 0.6]
    X = []
    for i in range(warmups):
        x = network(x)
        # x[0] = x[0] + ( amplitude * input[i % len(input)] )
    for i in range(measure):
        X.append(x[0])
        x = network(x)
        # x[0] = x[0] + ( amplitude * input[i % len(input)] )
    pylab.subplot(2,1,2)
    matplotlib.pyplot.psd(X,1024,32)

    pylab.show()

def lookatresults(data, modes, theta=None, vert=False, labels=None):


    P = data[-1][0]
    n = P.shape[0]

    if labels == None:
        labels = [""] * n
    else:
        pass 

    if vert == True:
        subplots = range(n*100+11,n*100+n+11,1)
        figsize = (6, 3*n)
    elif vert == 'four':
        subplots = [221, 222, 223, 224]
        figsize = (10, 10)
    else:
        subplots = range(100+n*10+1,100+n*10+1+n,1)
        figsize = (5*n, 3)

    f = stats.gaussian_kde(data[-1][0])
    int_guess = np.mean(data[-1][0], axis=1)
    modes = minimize(neg, int_guess, args=(f)).x

    thetas = []
    P = data[-1][0]
    labelpad = 20

    for i in xrange(n):
        x = P[i]
        t = r'$\theta_{3:}$ {1:.2f} +{2:.2f}/-{0:.2f}'.format(
            modes[i]-stats.scoreatpercentile(x, 16),
            modes[i],
            stats.scoreatpercentile(x, 84)-modes[i], i+1)

        thetas.append(t)

    if P.shape[1] > 10:
        bins = np.sqrt(P.shape[1])
    else:
        bins=10
    fig = plt.figure(figsize=figsize)
    
    for i in xrange(n):
        print subplots[i]
        plt.subplot(int(subplots[i]))
        #plt.title(thetas[0])
        ker = stats.gaussian_kde(P[i])
        h = plt.hist(P[i], bins=bins, normed=True, alpha=1)
        x = np.linspace(h[1][0],h[1][-1],1000)
        plt.plot(x,ker(x))
        plt.xlabel(labels[i], labelpad=labelpad, fontsize=24)
        if theta != None:
            plt.axvline(theta[0])

    for t in thetas: 
        print t

    return fig

def display_head(set_x, set_y, n = 5):
    '''
    show some figures based on gray image matrixs
    
    @type set_x: TensorSharedVariable, 
    @param set_x: gray level value matrix of the
    
    @type set_y: TensorVariable, 
    @param set_y: label of the figures    
    
    @type n: int, 
    @param n: numbers of figure to be display, less than 10, default 5
    '''
    import pylab
    
    if n > 10: n = 10
    img_x = set_x.get_value()[0:n].reshape(n, 28, 28)
    img_y = set_y.eval()[0:n]
    
    for i in range(n): 
        pylab.subplot(1, n, i+1); 
        pylab.axis('off'); 
        pylab.title(' %d' % img_y[i])
        pylab.gray()
        pylab.imshow(img_x[i])

def main( runTime ):
    try:
        moose.delete('/acc92')
        print("Deleted old model")
    except Exception as  e:
        print("Could not clean. model not loaded yet")

    moose.loadModel('acc92_caBuff.g',loadpath,'gsl')  
    ca = moose.element(loadpath+'/kinetics/Ca')
    pr = moose.element(loadpath+'/kinetics/protein')
    clockdt = moose.Clock('/clock').dts 
    moose.setClock(8, 0.1)#simdt
    moose.setClock(18, 0.1)#plotdt
    print clockdt
    print " \t \t simdt ", moose.Clock('/clock').dts[8],"plotdt ",moose.Clock('/clock').dts[18]
    ori =  ca.concInit
    tablepath = loadpath+'/kinetics/Ca'
    tableele = moose.element(tablepath)
    table = moose.Table2(tablepath+'.con')
    x = moose.connect(table, 'requestOut', tablepath, 'getConc')
    tablepath1 = loadpath+'/kinetics/protein'
    tableele1 = moose.element(tablepath1)
    table1 = moose.Table2(tablepath1+'.con')
    x1 = moose.connect(table1, 'requestOut', tablepath1, 'getConc')

    ca.concInit = ori
    print("[INFO] Running for 4000 with Ca.conc %s " % ca.conc)
    moose.start(4000)

    ca.concInit = 5e-03
    print("[INFO] Running for 20 with Ca.conc %s " % ca.conc)
    moose.start(20)

    ca.concInit = ori
    moose.start( runTime ) #here give the interval time 

    ca.concInit = 5e-03
    print("[INFO] Running for 20 with Ca.conc %s " % ca.conc)
    moose.start(20)

    ca.concInit = ori
    print("[INFO] Running for 2000 with Ca.conc %s " % ca.conc)
    moose.start(2000)

    pylab.figure()
    pylab.subplot(2, 1, 1)
    t = numpy.linspace(0.0, moose.element("/clock").runTime, len(table.vector)) # sec
    pylab.plot( t, table.vector, label="Ca Conc (interval- 8000s)" )
    pylab.legend()
    pylab.subplot(2, 1, 2)
    t1 = numpy.linspace(0.0, moose.element("/clock").runTime, len(table1.vector)) # sec
    pylab.plot( t1, table1.vector, label="Protein Conc (interval- 8000s)" )
    pylab.legend()
    pylab.savefig( os.path.join( dataDir, '%s_%s.png' % (table1.name, runTime) ) )

    print('[INFO] Saving data to csv files in %s' % dataDir)
    tabPath1 = os.path.join( dataDir, '%s_%s.csv' % (table.name, runTime))
    numpy.savetxt(tabPath1, numpy.matrix([t, table.vector]).T, newline='\n')
    tabPath2 = os.path.join( dataDir, '%s_%s.csv' % (table1.name, runTime) )
    numpy.savetxt(tabPath2, numpy.matrix([t1, table1.vector]).T, newline='\n')

def plot_signal(x,y,title,labelx,labley,position):
    pylab.subplot (9, 1, position)
    pylab.plot(x,y)
    pylab.title(title)
    pylab.xlabel(labelx)
    pylab.ylabel(labley)
    pylab.grid(True)

    def draw(self):

        print self.edgeno

        pos = 0
        dy = 8
        edgeno = self.edgeno
        edge = self.edges[edgeno]
        edgeprev = self.edges[edgeno-1]
        p = np.round(edge["top"](1024))
        top = min(p+2*dy, 2048)
        bot = min(p-2*dy, 2048)
        self.cutout = self.flat[1][bot:top,:].copy()

        pl.figure(1)
        pl.clf()
        start = 0
        dy = 512
        for i in xrange(2048/dy):
            pl.subplot(2048/dy,1,i+1)
            pl.xlim(start, start+dy)

            if i == 0: pl.title("edge %i] %s|%s" % (edgeno,
                edgeprev["Target_Name"], edge["Target_Name"]))


            pl.subplots_adjust(left=.07,right=.99,bottom=.05,top=.95)

            pl.imshow(self.flat[1][bot:top,start:start+dy], extent=(start,
                start+dy, bot, top), cmap='Greys', vmin=2000, vmax=6000)

            pix = np.arange(start, start+dy)
            pl.plot(pix, edge["top"](pix), 'r', linewidth=1)
            pl.plot(pix, edgeprev["bottom"](pix), 'r', linewidth=1)
            pl.plot(edge["xposs_top"], edge["yposs_top"], 'o')
            pl.plot(edgeprev["xposs_bot"], edgeprev["yposs_bot"], 'o')


            hpp = edge["hpps"]
            pl.axvline(hpp[0],ymax=.5, color='blue', linewidth=5)
            pl.axvline(hpp[1],ymax=.5, color='red', linewidth=5)

            hpp = edgeprev["hpps"]
            pl.axvline(hpp[0],ymin=.5,color='blue', linewidth=5)
            pl.axvline(hpp[1],ymin=.5,color='red', linewidth=5)


            if False:
                L = top-bot
                Lx = len(edge["xposs"])
                for i in xrange(Lx):
                    xp = edge["xposs"][i]
                    frac1 = (edge["top"](xp)-bot-1)/L
                    pl.axvline(xp,ymin=frac1)

                for xp in edgeprev["xposs"]: 
                    frac2 = (edgeprev["bottom"](xp)-bot)/L
                    pl.axvline(xp,ymax=frac2)

            start += dy