def test2_filtersb_test_real(self):
		"""Test if filter_sideband() on real data for inspection"""
		# NB: jpgs are required for correct [0,1] data range. PNG data range 
		# [0, 255] gives weird results?
		files = ['fringe_130622_154235Z_000082_img.jpg',
			'fringe_130622_154235Z_000139_img.jpg',
			'fringe_130622_154235Z_000220_img.jpg',
			'fringe_130622_154235Z_000330_img.jpg',
			'fringe_130622_154235Z_000412_img.jpg',
			'fringe_130622_154235Z_000505_img.jpg']
		fringes = [tim.file.read_file(pjoin(TESTDATAPATH,f)) for f in files]

		cfreq = fringe_cal(fringes, store_pow=False, do_embed=True)
		apt_mask = tim.im.mk_rad_mask(*fringes[0].shape) < 1

		self.fcache = {}
		phasepow_l = [filter_sideband(i, cfreq, 0.5, method='spectral', apt_mask=apt_mask, unwrap=True, wsize=-0.5, wfunc='cosine', do_embed=True, cache=self.fcache, ret_pow=True, get_complex=False, verb=0) for i in fringes]

		for im, phase in zip(fringes, phasepow_l):
			plt.figure(400, figsize=(4,4)); plt.clf()
			ax0 = plt.subplot2grid((2,2),(0, 0))
			ax0.set_title("fringes")
			ax0.imshow(im)
			ax1 = plt.subplot2grid((2,2),(0, 1))
			ax1.set_title("log side band power")
			ax1.imshow(np.log(phase[2]))
			ax2 = plt.subplot2grid((2,2),(1, 0))
			ax2.set_title("phase")
			ax2.imshow(phase[0])
			ax3 = plt.subplot2grid((2,2),(1, 1))
			ax3.set_title("amplitude")
			ax3.imshow(phase[1])
			if (raw_input("Continue [b=break]...") == 'b'): break
		plt.close()

def plot_simulation_results_initial_controls(time_points, y_sim):

    pl.rc("text", usetex = True)
    pl.rc("font", family="serif")

    pl.subplot2grid((4, 2), (0, 0))
    pl.plot(time_points, y_sim[:,0])
    pl.title("Simulation results for initial controls")
    pl.xlabel("t")
    pl.ylabel("X", rotation = 0, labelpad = 20)

    pl.subplot2grid((4, 2), (1, 0))
    pl.plot(time_points, y_sim[:,1])
    pl.xlabel("t")
    pl.ylabel("Y", rotation = 0, labelpad = 15)

    pl.subplot2grid((4, 2), (2, 0))
    pl.plot(time_points, y_sim[:,2])
    pl.xlabel("t")
    pl.ylabel(r"\phi", rotation = 0, labelpad = 15)

    pl.subplot2grid((4, 2), (3, 0))
    pl.plot(time_points, y_sim[:,3])
    pl.xlabel("t")
    pl.ylabel("v", rotation = 0, labelpad = 20)

    pl.subplot2grid((4, 2), (0, 1), rowspan = 4)
    pl.plot(y_sim[:,0], y_sim[:, 1])
    pl.title("Simulated race car path for initial controls")
    pl.xlabel("X")
    pl.ylabel("Y", rotation = 0, labelpad = 20)
    pl.show()

def plot_measurements(time_points, ydata):

    pl.rc("text", usetex = True)
    pl.rc("font", family="serif")

    pl.subplot2grid((4, 2), (0, 0))
    pl.plot(time_points, ydata[:,0])
    pl.title("Considered measurement data")
    pl.xlabel("t")
    pl.ylabel("X", rotation = 0, labelpad = 20)

    pl.subplot2grid((4, 2), (1, 0))
    pl.plot(time_points, ydata[:,1])
    pl.xlabel("t")
    pl.ylabel("Y", rotation = 0, labelpad = 15)

    pl.subplot2grid((4, 2), (2, 0))
    pl.plot(time_points, ydata[:,2])
    pl.xlabel("t")
    pl.ylabel(r"\phi", rotation = 0, labelpad = 15)

    pl.subplot2grid((4, 2), (3, 0))
    pl.plot(time_points, ydata[:,3])
    pl.xlabel("t")
    pl.ylabel("v", rotation = 0, labelpad = 20)

    pl.subplot2grid((4, 2), (0, 1), rowspan = 4)
    pl.plot(ydata[:,0], ydata[:, 1])
    pl.title("Considered racecar track (measured)")
    pl.xlabel("X")
    pl.ylabel("Y", rotation = 0, labelpad = 20)
    pl.show()

def weightedAffineInvCompWarpCost(targetIMG, templateIMG, templateWeight, curParameters, displayStuff, targetIMGMask = None):
	#function [ImageWarpedToTemplate, TX, TY, ErrorIMG, CostValue] = lk_weighted_run_affine_inv_comp_warpcost(targetIMG, templateIMG, TemplateW, CurParameters, displayStuff)
#	[ImageWarpedToTemplate, TX, TY] = lk_warp_image_affine(targetIMG, size(templateIMG), CurParameters);
	targetIMGToTemplate, TX, TY = warpImageAffine(targetIMG, templateIMG.shape, curParameters)
	if displayStuff == True:
		pylab.subplot(4, 3, 5);	CCSegUtils.showIMG(targetIMGToTemplate); pylab.title('template coordinates warped to image');
		pylab.subplot2grid((4, 3), (1, 2), rowspan = 2, colspan = 1); pylab.cla(); CCSegUtils.showIMG(targetIMG);
			
		pylab.plot(TX[:, 0], TY[:, 0], 'b-')
		pylab.plot(TX[0, :], TY[0, :], 'b-')
		pylab.plot(TX[:, -1], TY[:, -1], 'b-')
		pylab.plot(TX[-1, :], TY[-1, :], 'b-')
		pylab.title('Coordinates on target')
	#print "oiajfdoijadsf"

	errorIMG = targetIMGToTemplate - templateIMG

	
	LKCost = numpy.sum(errorIMG * templateWeight * errorIMG)
	
	# find out if any coordinates are not in the mask
	if targetIMGMask != None:
		T = CCSegUtils.interp2q(numpy.arange(1, targetIMG.shape[1] + 1), numpy.arange(1, targetIMG.shape[0] + 1), targetIMGMask, TX, TY, extrapval = 0)
		if numpy.any(T == 0):
			LKCost = numpy.inf
	
	return (targetIMGToTemplate, TX, TY, errorIMG, LKCost)

def plot_xin(viking,data,fit,output="plot.png",startrow=1,stoprow=-1,baseline=None):
    values=odict()
    #for fit functions
    
    values["vl1"] = odict()
    values["vl2"] = odict()
    
    for name in ["viking","data","fit"]:
        fname = getattr(args,name)
        try:
            values["vl1"][name]=read_data(fname, "vl1", args.startrow, args.stoprow)
            values["vl2"][name]=read_data(fname, "vl2", args.startrow, args.stoprow)
        except AttributeError as e:
            pass
    
    for v in ["vl1","vl2"]:
        values[v]["diff"] = dict(values[v]["viking"])
        values[v]["diff"]["y"] = values[v]["diff"]["y"] - values[v]["fit"]["y"]
    

    pylab.figure(figsize=(6,6))
    ylim=None
    
    plots = odict()
    top_plots = [  ["vl1",values["vl1"]["data"], "grey","-",1],
                    [None,values["vl1"]["fit"], "magenta","--",2],    
                    [None,values["vl1"]["viking"], "blue","--",2],
                  ["vl2",values["vl2"]["data"], "black","-",1],
                    [None,values["vl2"]["fit"], "green","--",2],    
                    [None,values["vl2"]["viking"], "red","--",2]]

    bottom_plots = [ ["vl1",values["vl1"]["diff"], "blue","-",2],
                     ["vl2",values["vl2"]["diff"], "red","-",2] ]

    ax = pylab.subplot2grid((9,1),(0, 0),rowspan=6)
    handles = odict()
    for n,p,c,l,w in top_plots:
        ylim,handle=segmented_plot(p["x"],p["y"],ylim=ylim, color=c, ls=l,lw=w, return_handle=True,alpha=1)
        if n:
            handles[n] = handle
    pylab.legend(handles,loc='upper left', framealpha=0.5)
    pylab.xlabel("Ls") ; pylab.xlim(0,360) ; pylab.xticks([0,60,120,180,240,300,360])
    pylab.ylabel("Pressure (Pa)") ; pylab.ylim(ylim)
    pylab.grid()
    
    
    ax = pylab.subplot2grid((9,1),(7, 0),rowspan=2)
    ylim=None
    handles = odict()
    for n,p,c,l,w in bottom_plots:
        ylim,handle=segmented_plot(p["x"],p["y"],ylim=ylim, color=c, ls=l,lw=w, return_handle=True,alpha=1)
        if n:
            handles[n] = handle
    pylab.legend(handles, loc="upper left", framealpha=0.5)
    pylab.xlabel("Ls") ; pylab.xlim(0,360) ; pylab.xticks([0,60,120,180,240,300,360])
    pylab.ylabel("P (Pa)") ; pylab.ylim(ylim)
    pylab.grid()
    
    pylab.savefig(filename=args.output, dpi=150)

def main( ):
    global loc_
    nGames = 20
    prevRows, prevCells = np.zeros( shape=(10,10) ), np.zeros( shape=(100,100))
    cellTransitionProbMean = [ ]
    rowTransitionProbMean = [ ]
    for i in range( nGames ):
        loc_ = 0
        while True:
            r1, c1 = int(loc_ / 10), loc_
            loc_ += roll_dice( )
            loc_ = check_for_snakes_and_ladders( loc_ )
            r2, c2 = int(loc_ / 10), loc_
            if loc_ < 100:
                cells_[c1,c2] += 1.0
                rows_[r1,r2] += 1.0
            if loc_ > 100:
                print( "Game %d is over" % i )
                loc_ = 0
                thisRows = renormalize( rows_ )
                thisCells = renormalize( cells_ )
                diffRow = thisRows - prevRows
                diffCell = thisCells - prevCells
                rowTransitionProbMean.append( diffRow.mean( ) )
                cellTransitionProbMean.append( diffCell.mean( ) )
                prevRows, prevCells = thisRows, thisCells
                # drawNow( img )
                break
                

    pCells = cells_ / cells_.sum( axis = 0 )
    pRows = rows_ / rows_.sum( axis = 0 )
    with np.errstate(divide='ignore', invalid='ignore'):
        pCells = np.true_divide(cells_, cells_.sum( axis = 0 ) )
        pRows = np.true_divide( rows_, rows_.sum( axis = 0 ) )
        pCells = np.nan_to_num( pCells )
        pRows = np.nan_to_num( pRows )
    pylab.figure( figsize=(10,6) )

    gridSize = (2, 2)
    ax1 = pylab.subplot2grid( gridSize, (0,0), colspan = 1 )
    ax2 = pylab.subplot2grid( gridSize, (0,1), colspan = 1 )
    ax3 = pylab.subplot2grid( gridSize, (1,0), colspan = 2 )
    ax1.set_title( 'Cell to cell transitions' )
    img1 = ax1.imshow( pCells, interpolation = 'none' )
    pylab.colorbar( img1, ax = ax1 ) #orientation = 'horizontal' )
    img2 = ax2.imshow( pRows, interpolation = 'none' )
    pylab.colorbar( img2, ax = ax2 ) # orientation = 'horizontal' )
    ax2.set_title( 'Row to row transitions' )

    ax3.plot( rowTransitionProbMean, '-o', label = 'Row probs' )
    ax3.plot( cellTransitionProbMean, '-*', label = 'Cell probs' )
    ax3.set_ylabel( 'diff of transition probs'  )
    ax3.set_xlabel( 'Number of games completed' )
    ax3.legend( )

    pylab.suptitle( 'Total games %d' % nGames )
    # pylab.tight_layout( )
    pylab.savefig( '%s.png' % sys.argv[0] )

    def _run_interface(self, runtime):
        import matplotlib

        matplotlib.use("Agg")
        import pylab as plt

        realignment_parameters = np.loadtxt(self.inputs.realignment_parameters)
        title = self.inputs.title

        F = plt.figure(figsize=(8.3, 11.7))
        F.text(0.5, 0.96, self.inputs.title, horizontalalignment="center")
        ax1 = plt.subplot2grid((2, 2), (0, 0), colspan=2)
        handles = ax1.plot(realignment_parameters[:, 0:3])
        ax1.legend(handles, ["x translation", "y translation", "z translation"], loc=0)
        ax1.set_xlabel("image #")
        ax1.set_ylabel("mm")
        ax1.set_xlim((0, realignment_parameters.shape[0] - 1))
        ax1.set_ylim(bottom=realignment_parameters[:, 0:3].min(), top=realignment_parameters[:, 0:3].max())

        ax2 = plt.subplot2grid((2, 2), (1, 0), colspan=2)
        handles = ax2.plot(realignment_parameters[:, 3:6] * 180.0 / np.pi)
        ax2.legend(handles, ["pitch", "roll", "yaw"], loc=0)
        ax2.set_xlabel("image #")
        ax2.set_ylabel("degrees")
        ax2.set_xlim((0, realignment_parameters.shape[0] - 1))
        ax2.set_ylim(
            bottom=(realignment_parameters[:, 3:6] * 180.0 / np.pi).min(),
            top=(realignment_parameters[:, 3:6] * 180.0 / np.pi).max(),
        )

        if isdefined(self.inputs.outlier_files):
            try:
                outliers = np.loadtxt(self.inputs.outlier_files)
            except IOError as e:
                if e.args[0] == "End-of-file reached before encountering data.":
                    pass
                else:
                    raise
            else:
                if outliers.size > 0:
                    ax1.vlines(outliers, ax1.get_ylim()[0], ax1.get_ylim()[1])
                    ax2.vlines(outliers, ax2.get_ylim()[0], ax2.get_ylim()[1])

        if title != "":
            filename = title.replace(" ", "_") + ".pdf"
        else:
            filename = "plot.pdf"

        F.savefig(filename, papertype="a4", dpi=self.inputs.dpi)
        plt.clf()
        plt.close()
        del F

        self._plot = filename

        runtime.returncode = 0
        return runtime

def SE_spectrum_analysis(sim_exp):
    #load matlab file where timeserie is stored
    #f = h5py.File('C:/Users/Squirel/Desktop/data_antoine/status_epilepticus/timeserie_SE_4j_rat8l_sr250Hz.mat', 'r')
    if sim_exp=='exp':
        f = h5py.File('C:/Users/Squirel/Desktop/data_antoine/status_epilepticus/EEG_rat8l.mat', 'r')
        data = f.get('EEG/data')
        data = pb.array(data) #transform to numpy array --> now each element of data is an array of one float
        sf = 250. # sampling freq (Hz)
        beg=120000; lgth = 20.; #stating time (s) and length (s) of the selected period
        spectrum_ymax = 0.2 #maximum value displayed in plot (to have all the same scale)
        ts_ylim = (-0.0020, 0.0010)
    else: #sim_exp=='sim'
        #radical='_40_HMR_40_ML_CpES1_0_CpES2_0_x0_35_noise3_0_noise2_30_noise1_60_gx1x2_20_gx2x1_20_gx1x1_20_gx2x2_20_r40_20s' #phase I
        #radical='_40_HMR_40_ML_CpES1_20_CpES2_20_x0_20_noise3_0_noise2_30_noise1_60_gx1x2_20_gx2x1_20_gx1x1_20_gx2x2_20_r40_20s' # phase II
        #radical='_40_HMR_40_ML_CpES1_80_CpES2_80_x0_20_noise3_0_noise2_30_noise1_60_gx1x2_20_gx2x1_20_gx1x1_20_gx2x2_20_r40_20s' #phase III
        radical='_40_HMR_40_ML_CpES1_40_CpES2_40_x0_35_noise3_0_noise2_30_noise1_60_gx1x2_20_gx2x1_20_gx1x1_20_gx2x2_20_r40_100s' #phase IV     !!! set beg=20 instead of 10s
        x1_plot = np.load('./traces_dec13/x1plot'+radical+'.npy')
        x2_plot = np.load('./traces_dec13/x2plot'+radical+'.npy')
        #downsampling to 1kHz
        x1_plot =  x1_plot[:,0:-1:10]
        x2_plot = x2_plot[:,0:-1:10]
        
        data = -(0.75*x1_plot.mean(axis=0) + 0.25*x2_plot.mean(axis=0)) #transform to numpy array --> now each element of data is an array of one float
        sf = 1000. # sampling freq (Hz)
        beg=20.; lgth = 80.; #stating time (s) and length (s) of the selected period
        spectrum_ymax = 300 #maximum value displayed in plot (to have all the same scale)
        ts_ylim = (-0.5, 1.75)
    
    minfreq=1.; maxfreq=50.; #Hz
    #pb.plot(range(lgth*sf), data[beg*sf:beg*sf+lgth*sf]); pb.show() #if splot needed 
    signal = data[beg*sf:beg*sf+lgth*sf]
    #signal = pb.sin(2*pb.pi*pb.arange(0,15,0.004))
    #pb.plot(signal); pb.show()
    spectrum = fft.fftn(signal)
    delta_f = 1.0 / ((lgth*sf)/2) * lgth
    #pb.plot(pb.arange(minfreq, maxfreq, delta_f),abs(spectrum[minfreq*delta_f:maxfreq*delta_f]), "."); pb.show()
    #pb.plot(pb.arange(minfreq, maxfreq, 1/lgth), abs(spectrum[minfreq*lgth:maxfreq*lgth]), "."); pb.show()
    
    spec = pb.arange(minfreq,maxfreq,1.0)
    for i in range(len(spec)):
        spec[i] = pb.mean(abs(spectrum[(minfreq+i)*lgth:(minfreq+i)*lgth + lgth]))
    #pb.plot(spec)
    #pb.show()
    
    #plots:
    fig = pb.figure()
    ax1 = pb.subplot2grid((2,2),(0,0),colspan=2)
    ax2 = pb.subplot2grid((2,2),(1,0))
    ax3 = pb.subplot2grid((2,2),(1,1))
    
    ax1.plot(signal)
    ax1.set_ylim(ts_ylim)
    ax2.plot(pb.arange(minfreq, maxfreq, 1/lgth), abs(spectrum[minfreq*lgth:maxfreq*lgth]), ".")
    ax2.set_ylim(0,spectrum_ymax)
    ax3.plot(spec)
    pb.show()

def plotPendulumWithThetaAndPhi(animate = True, \
    title = "PendulumNormalModes.png", t0 = p.t0, p0 = p.p0, nsteps = 120, massRat = rm):
    
    global rm
    old = rm
    rm = massRat
    
    fig = pylab.figure()
    if animate:
        ax1 = pylab.subplot2grid((1, 2), (0, 1), aspect = 1.0) # Pendulum
        ax2 = pylab.subplot2grid((1, 2), (0, 0)) # Theta/Phi plots
    else:
        ax1 = fig.gca(aspect = 1.0)
        fig2 = pylab.figure()
        ax2 = fig2.gca()
    
    pendSetup(ax1, titleHigh = True)
    ax1.set_axis_on()
    line1, = ax1.plot([], [], 'k', linestyle = '-', marker = 'o', \
        markeredgecolor = 'k', markerfacecolor = (0.5, 0.5, 1.0), markersize = 10)
    line1.set_markevery((1, 1))
    
    line2, = ax2.plot([], [], 'b-', label = r"$\theta$")
    line3, = ax2.plot([], [], 'r-', label = r"$\phi$")
    ax2.set_xlim((0, nsteps * p.dt))
    ax2.set_ylim((-30 * np.pi / 180.0, 30 * np.pi / 180))
    ax2.set_xlabel(r"simulation time $\tau$", fontsize = 14)
    ax2.set_ylabel("pendula angles", fontsize = 14)
    
    pendAnimator1 = PendulumAnimate(line1)
    angleAnimator = AngleAnimator(line2, line3)
    
    simulation = DopeSimulation([pendAnimator1, angleAnimator], \
        x0 = np.array([t0, p0, 0.0, 0.0]), nsteps = nsteps)
    
    if animate:
        anim = animation.FuncAnimation(fig, simulation, \
            blit = True, init_func = simulation.reset, \
            interval = 10, repeat = True)
        
        simulation.reset()
    
    for i in range(nsteps):
        simulation(i)
        
    pylab.legend(loc = 'lower right')
    pylab.tight_layout()
    
    if animate:
        pylab.show()
    else:
        fig.savefig("pendulum_" + title)
        fig2.savefig("angles_" + title)
    
    rm = old

def PlotMCMC(params,  EPIC, sampler, lc, GP=False, KepName=''):
    import corner
    if type(sampler)!=np.ndarray:
        samples = sampler.chain[:, -10000:, :].reshape((-1, len(params)))
    else:
        samples=sampler

    #Turning b into absolute...
    samples[:,1]=abs(samples[:,1])
    p.figure(1)

    #Clipping extreme values (top.bottom 0.1 percentiles)
    toclip=np.array([(np.percentile(samples[:,t],99.9)>samples[:,t]) // (samples[:,t]>np.percentile(samples[:,t],0.1)) for t in range(len(params))]).all(axis=0)
    samples=samples[toclip]

    #Earmarking the difference between GP and non
    labs = ["$T_{c}$","$b$","$v$","$Rp/Rs$","$u1$","$u2$","$\sigma_{white}$","$tau$","$a$"] if GP\
        else ["$T_{c}$","$b$","$v$","$Rp/Rs$","$F_0$","$u1$","$u2$"]

    #This plots the corner:
    fig = corner.corner(samples, label=labs, quantiles=[0.16, 0.5, 0.84], plot_datapoints=False)

    #Making sure the lightcuvre plot doesnt overstep the corner
    ndim=np.shape(samples)[1]
    rows=(ndim-1)/2
    cols=(ndim-1)/2

    #Printing Kepler name on plot
    p.subplot(ndim,ndim,ndim+3).axis('off')
    if KepName=='':
        if str(int(EPIC)).zfill(9)[0]=='2':
            KepName='EPIC'+str(EPIC)
        else:
            KepName='KIC'+str(EPIC)
    p.title(KepName, fontsize=22)

    #This plots the model on the same plot as the corner
    ax = p.subplot2grid((ndim,ndim), (0, ndim-cols), rowspan=rows-1, colspan=cols)
    modelfits=PlotModel(lc, samples, scale=1.0, GP=True)

    #plotting residuals beneath:
    ax = p.subplot2grid((ndim,ndim), (rows-1, ndim-cols), rowspan=1, colspan=cols)
    modelfits=PlotModel(lc,  samples, modelfits, residuals=True, GP=True, scale=1.0)

    #Saving as pdf. Will save up to 3 unique files.
    if os.path.exists(Namwd+'/Outputs/Corner_'+str(EPIC)+'_1.pdf'):
        if os.path.exists(Namwd+'/Outputs/Corner_'+str(EPIC)+'_2.pdf'):
            fname=Namwd+'/Outputs/Corner_'+str(EPIC)+'_3.pdf'
        else:
            fname=Namwd+'/Outputs/Corner_'+str(EPIC)+'_2.pdf'
    else:
        fname=Namwd+'/Outputs/Corner_'+str(EPIC)+'_1.pdf'
    p.savefig(fname,Transparent=True,dpi=300)

    return modelfits

 def __init__(self,overview):
     self.overview = overview
     self.fig = plt.figure(figsize=[14,12])
     self.prof_ax = plt.subplot2grid([5,9],[0,1],colspan=2)
     self.fold_ax = plt.subplot2grid([5,9],[1,1],colspan=2,rowspan=2,sharex=self.prof_ax)
     self.subs_ax = plt.subplot2grid([5,9],[1,0],rowspan=2,sharey=self.fold_ax)
     self.table_ax = plt.subplot2grid([5,9],[0,3],colspan=3,rowspan=3,frameon=False)
     self.dm_ax  = plt.subplot2grid([5,9],[0,6],colspan=2)
     self.acc_ax = plt.subplot2grid([5,9],[1,8],colspan=1,rowspan=2)
     self.dm_acc_ax = plt.subplot2grid([5,9],[1,6],colspan=2,rowspan=2,sharex=self.dm_ax,sharey=self.acc_ax)
     self.all_ax = plt.subplot2grid([6,9],[4,0],colspan=9,rowspan=3)
     self._plot_all_cands(self.all_ax)
     self.timers = {
         "read":Timer(),
         "prof":Timer(),
         "fold":Timer(),
         "stat":Timer(),
         "table":Timer(),
         "dm":Timer(),
         "acc":Timer(),
         "dmacc":Timer(),
         "write":Timer(),
         "clear":Timer()
         }
     self.header = self.overview._xml.find("header_parameters")
     self.fig.suptitle("Source name: %s"%self.header.find("source_name").text,fontsize=16)

def plot_rfs_and_theta(sim,neurons,connections):

    pre,post=neurons
    c=connections[0]

    weights=c.weights
    

    num_neurons=len(weights)
    fig=py.figure(figsize=(16,4*num_neurons))

    for i,w in enumerate(weights):
        try:  # check for channel
            neurons=pre.neuron_list
        except AttributeError:
            neurons=[pre]

        num_channels=len(neurons)

        count=0
        vmin,vmax=w.min(),w.max()
        for c,ch in enumerate(neurons):   
            try:
                rf_size=ch.rf_size
                if rf_size<0:
                    rf_size=py.sqrt(ch.N)
                    assert rf_size==int(rf_size)
                    rf_size=int(rf_size)

            except AttributeError:
                rf_size=py.sqrt(ch.N)
                assert rf_size==int(rf_size)
                rf_size=int(rf_size)


            py.subplot2grid((num_neurons,num_channels+1),(i, c),aspect='equal')
            subw=w[count:(count+rf_size*rf_size)]
            #py.pcolor(subw.reshape((rf_size,rf_size)),cmap=py.cm.gray)
            py.pcolormesh(subw.reshape((rf_size,rf_size)),cmap=py.cm.gray,
                vmin=vmin,vmax=vmax)
            py.xlim([0,rf_size]); 
            py.ylim([0,rf_size])
            py.axis('off')
            count+=rf_size*rf_size

    py.subplot2grid((num_neurons,num_channels+1),(0, num_channels))
    sim.monitors['theta'].plot()

    return fig

  def setup_main_screen(self):
    self.fig = pylab.figure(figsize=(8,8))
    pylab.subplots_adjust(top=.93, bottom=0.1, left=.15, right=.97, hspace=.01, wspace=.15)
    pylab.suptitle('Filter explore')
    self.ax = pylab.subplot2grid((2, 5), (0, 0), colspan=4) #Amplitude response
    self.ax2 = pylab.subplot2grid((2, 5), (1, 0), colspan=4) #Phase response

    meax = self.menuax = pylab.subplot2grid((2, 5), (0, 4))
    menu_names = ['butterworth', 'chebyshev I', 'chebyshev II', 'elliptic', 'bessel']
    self.filter_types = ['butter', 'cheby1', 'cheby2', 'ellip', 'bessel']
    self.menu_h = []
    for n, item in enumerate(menu_names):
      self.menu_h.append(meax.text(0.1, n, item, picker=5))
    self.menu_h[0].set_weight('bold')
    pylab.setp(meax, 'xticks', [], 'yticks', [], 'xlim',[0,1], 'ylim', [-.5, len(menu_names)-.5])

def pyplot(f, table, columns, resonance):
    import pylab

    global LINENUM
    first = LINENUM < 0
    LINENUM += 1
    if LINENUM >= len(colors):
        showplot(show=False)
        LINENUM = 0
    if LINENUM == 0: 
        if not first: FIGURES.append(pylab.gcf())
        pylab.figure(figsize=(18,3), tight_layout={'pad':0})
        pylab.subplot2grid((1,5),(0,0),colspan=4)
        #pylab.figure(figsize=(8,2.3), tight_layout={'pad':0})
        #pylab.subplot2grid((1,3),(0,0),colspan=2)

    #print "shape",table.shape
    name = os.path.splitext(os.path.basename(f))[0]
    if '-' in name: name = "-".join(name.split('-')[1:])
    label=name
    p = abundance(f)
    if p: label += " %.1f%%"%p
    #label += " res: %.2fA"%wavelength(resonance)
    color=colors[LINENUM]
    for i in (2,): #range(1,table.shape[0]):
        pylab.loglog(table[0,:].T*1e3, table[i,:], label=label, 
                     linestyle=lines[i-1], color=color)
        label='_nolegend_'


    # Table of relative total cross section
    if False: 
      if first:
        print " "*(8*table.shape[0])+"%7s %7s %7s %7s"%("0.5A","6A","15A","20A")
      TARGET=V2200
      y0 = [np.interp(TARGET,table[0,:],table[i,:]) 
          for i in range(1,table.shape[0])]
      b_c = np.sqrt(y0[0]/(0.01*4*np.pi))
      b_cL = np.sqrt(np.interp([L0p1,L0p2,L0p5,L6,L15,L20],table[0,:],table[1,:])/(0.01*4*np.pi))
      delta = (b_cL-b_c)/b_c*100
      #y0[0] = np.sqrt(y0[0]/(4*np.pi))
      print "%7s"%name," ".join("%7.3f"%vi for vi in y0)," ".join("%6.1f%%"%vi for vi in delta)
    else:
        if first:
            print "%7s %7s %15s"%("name","%","res. onset Ang/meV") 
        print "%7s %7s %7.2f %9.2f"%(
            name, ("%.3f"%p if p else "-"), 
            wavelength(resonance), resonance*1e3)

def singular(Ytra,Ytes,sdata):

  m=Ytra.shape[0]
  mt=Ytes.shape[0]
  
  xmean=np.mean(Ytra,axis=0)
  ytra0=Ytra-np.outer(np.ones(m),xmean)
  xmean=np.mean(Ytes,axis=0)
  ytes0=Ytes-np.outer(np.ones(mt),xmean)

  xnorm=np.sqrt(np.sum(ytra0**2,axis=0))
  Ctra=np.dot(ytra0.T,ytra0)
  xnorm=xnorm+(xnorm==0)
  Ctra=Ctra/np.outer(xnorm,xnorm)
  xnorm=np.sqrt(np.sum(ytes0**2,axis=0))
  Ctes=np.dot(ytes0.T,ytes0)
  xnorm=xnorm+(xnorm==0)
  Ctes=Ctes/np.outer(xnorm,xnorm)

  stra=np_lin.svd(Ctra)[1]
  stes=np_lin.svd(Ctes)[1]

  xlinewidth=3
  
  fig1=lab.figure(figsize=(12,6))
  ax1=lab.subplot2grid((10,2),(1,0), rowspan=9)
  ## ax=fig1.add_subplot(1,2,1)
  ax1.plot(stra,linewidth=xlinewidth, \
          linestyle='-',color='b')
  ax1.set_xlabel('Label indexes', fontsize=14)
  ax1.set_ylabel('Eigen values',fontsize=14)
  ax1.set_title('Training:',fontsize=16)
  ax1.grid(True)

  ## ax=fig1.add_subplot(1,2,2)
  ax2=lab.subplot2grid((10,2),(1,1), rowspan=9)
  ax2.plot(stes,linewidth=xlinewidth, \
          linestyle='-',color='b')
  ax2.set_xlabel('Label indexes', fontsize=14)
  ax2.set_ylabel('Eigen values',fontsize=14)
  ax2.set_title('Test:',fontsize=16)
  ax2.grid(True)

  fig1.suptitle('Eigen values of correlation matrixes: '+sdata,fontsize=18  )
    
  lab.show()

  return

def pca_spec(spec,W_red,mu,plot=False,plotname="spec.png"):
	spec = spec/(spec.dot(spec))-mu
	a = (W_red.T).dot(spec.T)
	model = W_red.dot(a)+mu
	spec+=mu
	if plot:
		from pylab import clf,plot,ylabel,subplot2grid,savefig
		waves=range(len(model))
		subplot2grid((3,3),(0,0),colspan=3,rowspan=2)
		ylabel('Scaled Flux')
		plot(waves,spec,'.',waves,model,'-r')
		subplot2grid((3,3),(2,0),colspan=3)
		plot(waves, (model-spec)/spec,".")
		ylabel('Residual (%)')
		savefig(plotname)
		clf()

	def plot_alt_az(self,grid,contour='',pfact=4):
		import pylab
		polar = pylab.subplot2grid(self.plotgrid,grid[:2],colspan=grid[2],rowspan=grid[3],projection='polar')
		polar.plot(np.radians(self.aa_az - 90),90 - self.aa_alt,'r.')
		polar.plot(np.radians(self.ar_az - 90),90 - self.ar_alt,'g.')
		polar.set_rmax(90)
		polar.set_xticklabels(['E','SE','S','SW','W','NW','N','NE'])
		if contour:
			X = np.radians(self.ar_az - 90)
			Y = 90 - self.ar_alt

			if contour == 'model':
				Z = self.diff_model_angular * 3600
				polar.set_title('Model differences')
			elif contour == 'real':
				Z = self.diff_angular * 3600
				polar.set_title('Real differences')

			xi = np.linspace(np.radians(-90),np.radians(271),num = 360 * pfact)
			yi = np.linspace(min(Y),max(Y),num = 90 * pfact)
			zi = pylab.griddata(X, Y, Z, xi, yi, interp='linear')
			ctf = polar.contourf(xi,yi,zi,cmap='hot')
			cbar = pylab.colorbar(ctf, orientation='horizontal', pad=0.05)
			cbar.set_ticks(range(0,int(max(Z)),60))
			cbar.ax.set_xticklabels(map("{0}'".format,range(0,int(max(Z) / 60))))
		else:
			polar.set_title('Alt-Az distribution')
                return polar

    def vecPlot2(self, row, col, r_m, mWx, mWy, data1, data2, dmr, frmSize, winSize):
        pl.clf()
        pl.ion()
        pWx = []
        l = dmr - row - col
        for f in range(l):
            pWx.append(mWx[row + f, col + f, :, 0])
        pWx = np.array(pWx)

        sqWx = np.sqrt(pWx * pWx)
        print sqWx
        r, c = sqWx.shape
        x = pl.arange(c + 1)
        y = pl.arange(r + 1)
        X, Y = pl.meshgrid(x, y)

        pl.subplot2grid((1, 2), (0, 0))
        pl.pcolor(X, Y, sqWx, vmin=0, vmax=1)
        pl.xlim(0, c)
        pl.ylim(0, r)
        pl.colorbar()
        pl.title("user_1 (t:" + str(row) + ")")
        pl.gray()

        pWy = []
        l = dmr - row - col
        for f in range(l):
            pWy.append(mWy[row + f, col + f, :, 0])
        pWy = np.array(pWy)

        sqWy = np.sqrt(pWy * pWy)
        print sqWy
        r, c = sqWy.shape
        x = pl.arange(c + 1)
        y = pl.arange(r + 1)
        X, Y = pl.meshgrid(x, y)

        pl.subplot2grid((1, 2), (0, 1))
        pl.pcolor(X, Y, sqWy, vmin=0, vmax=1)
        pl.xlim(0, c)
        pl.ylim(0, r)
        pl.colorbar()
        pl.title("user_2 (t:" + str(col) + ")")
        pl.gray()

        pl.tight_layout()
        pl.draw()

	def test1_filtersb_test_zern_sim(self):
		"""Test if filter_sideband() on simulated Zernike wavefront data"""
		### Make fake fringes
		cfreq0 = (18.81111508,  24.76214802)
		zndata = tim.zern.calc_zern_basis(10, min(self.sz)/2)
		rad_mask = tim.im.mk_rad_mask(*zndata['mask'].shape)
		apt_mask = zndata['mask']
		apt_mask9 = rad_mask<0.9

		np.random.seed(1337)
		zvecs = [np.random.random(10) for i in range(6)]
		zphases = [sum(zv*zmode for zv, zmode in zip(zvec, zndata['modes'])) for zvec in zvecs]
		zfringes = [sim_fringe(zph, cfreq0, noiseamp=0) for zph in zphases]

		cfreq = fringe_cal(zfringes, store_pow=False, do_embed=True, method='cog')
		# Should give approximately cfreq, difference should be less 
		# than 2% and pixel difference should be less than 0.3
		self.assertLess(np.abs(1-cfreq/np.r_[cfreq0]).mean(), 0.02)
		self.assertLess(np.abs(cfreq0 - cfreq).mean(), 0.4)

		self.fcache = {}
		for zphase, zfringe in zip(zphases, zfringes):
			phase, amp, ftpow = filter_sideband(zfringe, cfreq0, 0.5, method='spectral', apt_mask=apt_mask, unwrap=True, do_embed=True, wsize=-0.5, wfunc='cosine', ret_pow=True, cache=self.fcache, verb=0)

			dphase = (zphase*apt_mask-phase)
			dphase -= dphase[apt_mask].mean()

			plt.figure(400, figsize=(4,4)); plt.clf()
			ax0 = plt.subplot2grid((2,2),(0, 0))
			ax0.set_title("Input phase")
			im = ax0.imshow(zphase*apt_mask)
			ax1 = plt.subplot2grid((2,2),(0, 1))
			ax1.set_title("Input fringes")
			ax1.imshow(zfringe)
			plt.colorbar(im)
			ax2 = plt.subplot2grid((2,2),(1, 0))
			ax2.set_title("Rec. phase")
			im = ax2.imshow(phase)
			plt.colorbar(im)
			ax3 = plt.subplot2grid((2,2),(1, 1))
			ax3.set_title("In-rec. (inner 90%)")
			im = ax3.imshow(dphase*apt_mask9)
			plt.colorbar(im)

			if (raw_input("Continue [b=break]...") == 'b'): break
		plt.close()

def plot_initial_and_optimized_controls(time_points, \
    udata_init, udata_opt, umin, umax):

    pl.rc("text", usetex = True)
    pl.rc("font", family="serif")

    pl.subplot2grid((2, 1), (0, 0))
    pl.step(time_points[:-1], udata_init[:,0], label = "$\delta_{init}$")
    pl.step(time_points[:-1], udata_init[:,1], label = "$D_{init}$")

    pl.plot([time_points[0], time_points[-2]], [umin[0], umin[0]], \
        color = "b", linestyle = "dashed", label = "$\delta_{min}$")
    pl.plot([time_points[0], time_points[-2]], [umax[0], umax[0]], \
        color = "b", linestyle = "dotted", label = "$\delta_{max}$")

    pl.plot([time_points[0], time_points[-2]], [umin[1], umin[1]], \
        color = "g", linestyle = "dashed", label = "$D_{min}$")
    pl.plot([time_points[0], time_points[-2]], [umax[1], umax[1]], \
        color = "g", linestyle = "dotted", label = "$D_{max}$")

    pl.ylabel("$\delta,\,D$", rotation = 0)
    pl.ylim(-0.6, 1.1)
    pl.title("Initial and optimized controls")
    pl.legend(loc = "upper right")

    pl.subplot2grid((2, 1), (1, 0))
    pl.step(time_points[:-1], udata_opt[:,0], label = "$\delta_{opt,coll}$")
    pl.step(time_points[:-1], udata_opt[:,1], label = "$D_{opt,coll}$")

    pl.plot([time_points[0], time_points[-2]], [umin[0], umin[0]], \
        color = "b", linestyle = "dashed", label = "$\delta_{min}$")
    pl.plot([time_points[0], time_points[-2]], [umax[0], umax[0]], \
        color = "b", linestyle = "dotted", label = "$\delta_{max}$")

    pl.plot([time_points[0], time_points[-2]], [umin[1], umin[1]], \
        color = "g", linestyle = "dashed", label = "$D_{min}$")
    pl.plot([time_points[0], time_points[-2]], [umax[1], umax[1]], \
        color = "g", linestyle = "dotted", label = "$D_{max}$")

    pl.xlabel("$t$")
    pl.ylabel("$\delta,\,D$", rotation = 0)
    pl.ylim(-0.6, 1.1)
    pl.legend(loc = "upper right")

    pl.show()