def _brush(self,event,region,inverse=False):
        """
        This will loop through all the other subplots (without the brush region)
        and change the opacity of the points not associated with that region.
        
        when inverse is True, it will "unbrush" by resetting the opacity of the brushed points
        """
        opacity_fraction = self.opac
        # what variables are in the plot?
        plot_vars = [x[1] for x in self.axis_info if x[0] == event.inaxes][0]
        
        ## figure out the min max of this region
        minx, miny   = region.get_xy()
        maxx  = minx + region.get_width()
        maxy =  miny + region.get_height()
        
        ## now query the data to get all the sources that are inside this range
        if isinstance(self.data[plot_vars[0]][0],datetime.datetime):
            maxx = datetime.datetime.fromordinal(maxx)
            minx= datetime.datetime.fromordinal(minx)
        elif isinstance(self.data[plot_vars[0]][0],datetime.date):
            maxx = datetime.date.fromordinal(maxx)
            minx= datetime.date.fromordinal(minx)
            
        if isinstance(self.data[plot_vars[1]][0],datetime.datetime):
            maxy = datetime.datetime.fromordinal(maxx)
            miny= datetime.datetime.fromordinal(minx)
        elif isinstance(self.data[plot_vars[1]][0],datetime.date):
            maxy = datetime.date.fromordinal(maxy)
            miny= datetime.date.fromordinal(miny)
        
        inds = (self.data[plot_vars[0]]<= maxx) & (self.data[plot_vars[0]] > minx) & \
               (self.data[plot_vars[1]] <= maxy) & (self.data[plot_vars[1]] > miny)
        invinds = ~ inds  # get all indicies of those records not inside the region
        
        for a,pv in self.axis_info:
            # dont self brush!
            if a == event.inaxes:
                continue

            ## get the scatterplot color and alpha channel data
            self.t = a.collections[0]
            fc = self.t.get_facecolor() # this will be a 2d array
            '''Here we change the color and opacity of the points
            fc[index,0] = Red
            fc[index,1] = Green
            fc[index,2] = Blue
            fc[index,3] = Alpha
            
            default is  [ 0.4   ,  0.4   ,  1.    ,  1.0]
            '''
            if not inverse: 
                fc[invinds,2] /= 20. #reduce blue channel greatly
                fc[invinds,3] /= opacity_fraction 
            else:
                fc[invinds,2] *= 20.
                fc[invinds,3] *= opacity_fraction
            self.t.set_facecolor(fc)
            
        plt.draw()

def plot_spikes(time,voltage,APTimes,titlestr):
    """
    plot_spikes takes four arguments - the recording time array, the voltage
    array, the time of the detected action potentials, and the title of your
    plot.  The function creates a labeled plot showing the raw voltage signal
    and indicating the location of detected spikes with red tick marks (|)
    """
# Make a plot and markup
    plt.figure()
    plt.title(titlestr)
    plt.xlabel("Time (s)")
    plt.ylabel("Voltage (uV)") 

    plt.plot(time, voltage)
    
# Vertical positions for red marker
# The following attributes are configurable if required    
    vertical_markers_indent = 0.01 # 1% of Voltage scale height
    vertical_markers_height = 0.03 # 5% of Voltage scale height
    y_scale_height = 100 # Max of scale
    
    marker_ymin = 0.5 + ( max(voltage) / y_scale_height / 2 ) + vertical_markers_indent
    marker_ymax = marker_ymin + vertical_markers_height

# Drawing red markers for detected spikes
    for spike in APTimes:
        plt.axvline(spike, ymin=marker_ymin, ymax=marker_ymax, color='red')
    
    plt.draw()

    def test1():
        x = [0.5]*3
        xbounds = [(-5, 5) for y in x]


        GA = GenAlg(fitcalc1, x, xbounds, popMult=100, bitsPerGene=9, mutation=(1./9.), crossover=0.65, crossN=2, direction='min', maxGens=60, hammingDist=False)
        results = GA.run()
        print "*** DONE ***"
        #print results
        plt.ioff()
        #generate pareto frontier numerically
        x1_ = np.arange(-5., 0., 0.05)
        x2_ = np.arange(-5., 0., 0.05)
        x3_ = np.arange(-5., 0., 0.05)

        pfn = []
        for x1 in x1_:
            for x2 in x2_:
                for x3 in x3_:
                    pfn.append(fitcalc1([x1,x2,x3]))

        pfn.sort(key=lambda x:x[0])
        
        plt.figure()
        i = 0
        for x in results:
            plt.scatter(x[1][0], x[1][1], 20, c='r')

        plt.scatter([x[0] for x in pfn], [x[1] for x in pfn], 1.0, c='b', alpha=0.1)
        plt.xlim([-20,-1])
        plt.ylim([-12, 2])
        plt.draw()

def histogramac(histc):
    
    histc=histc.histogram(255)
    B.show()
    pylab.plot(histc)
    pylab.draw()
    pylab.pause(0.0001)

def ttplot(corf,sr,sl,norm,n, I1, I2,lag):
    global tplot_cum
    firstfile=int(input_info['n_first_image'])
    tplot=time.time()
    rchplot=int(ceil(log(n/nchannels)/log(2))+1)
    normplot=zeros((1,rcr),dtype=float32)
    
    for ir in xrange(rchplot):
       if ir==0:
           normplot[0,:nchannels]=1./arange(n-2,n-nchannels-2,-1)
       else:
           normplot[0,nchannels2*(ir+1):nchannels2*(ir+2)]=1./arange((n-1)/(2**ir)-nchannels2-1,(n-1)/(2**ir)-nchannels-1,-1)

    indt=int(nchannels+nchannels2*log(n/nchannels)/log(2))-2
    cc1=corf[0,:indt]/(sl[0,:indt]*sr[0,:indt])/normplot[0,:indt]
    cc2=corf[nq/2,:indt]/(sl[nq/2,:indt]*sr[nq/2,:indt])/normplot[0,:indt]
    t_axis=lag[0,:indt]
    t_axis2=tI_avg[0,:n]
    t_axis2b=tI_avg[0,:n]/dt+firstfile
    lm1.set_data(t_axis,cc1)
    lm2.set_data(t_axis2,I1)
    lm1b.set_data(t_axis,cc2)
    lm2b.set_data(t_axis2b,I2)
    ax1.set_xlim(min(t_axis),max(t_axis))
    ax1.set_ylim(min(cc1),max(cc1))
    ax1b.set_ylim(min(cc2),max(cc2))
    ax2.set_xlim(min(t_axis2),max(t_axis2))
    ax2b.set_xlim(min(t_axis2b),max(t_axis2b))
    ax2.set_ylim(min(I1),max(I1))
    ax2b.set_ylim(min(I2),max(I2))
    p.draw()
    tplot_cum+=time.time()-tplot

 def plot(self, outf=None, dosave=True, savedir="Plot/", show=True):
     if outf is None:
         outf = self.outf
         # print outf
     oo = mlab.csv2rec(outf, delimiter=" ")
     # print oo
     plt.errorbar(oo["time"] % self.period, oo["magnitude"], oo["error"], fmt="b.")
     plt.plot(oo["time"] % self.period, oo["model"], "ro")
     plt.title(
         "#%i P=%f d (chisq/dof = %f) r1+r2=%f"
         % (self.dotastro_id, self.period, self.outrez["chisq"], self.outrez.get("r1") + self.outrez.get("r2"))
     )
     ylim = plt.ylim()
     # print ylim
     if ylim[0] < ylim[1]:
         plt.ylim(ylim[1], ylim[0])
     plt.draw()
     if show:
         plt.show()
     if dosave:
         if not os.path.isdir(savedir):
             os.mkdir(savedir)
         plt.savefig("%splot%i.png" % (savedir, self.dotastro_id))  # ,self.period))
         print("Saved", "%splot%i.png" % (savedir, self.dotastro_id))  # ,self.period)
     plt.clf()

    def __init__(self, folder, **kwargs):  
        
        if not os.path.isdir(os.path.join(folder, 'plots')):
            os.mkdir(os.path.join(folder, 'plots'))
        plt.ioff()
        self.metrics_fig = plt.figure('Metrics')
        self.ax2 = self.metrics_fig.add_subplot(111)

        self.p1, = self.ax2.plot([], [], 'ro-', label='TEST: Pixel accuracy')
        self.p5, = self.ax2.plot([], [], 'rv-', label='TRAIN: Pixel accuracy')
        
        self.p2, = self.ax2.plot([], [], 'bo-', label='TEST: Mean-Per-Class accuracy')
        self.p6, = self.ax2.plot([], [], 'bv-', label='TRAIN:Mean-Per-Class accuracy')
        
        self.p3, = self.ax2.plot([], [], 'go-', label='TEST: Mean-Per-Class IU')
        self.p7, = self.ax2.plot([], [], 'gv-', label='TRAIN:Mean-Per-Class IU')
        
        self.p4, = self.ax2.plot([], [], 'ko-', label='TEST: Freq. weigh. mean IU')
        self.p8, = self.ax2.plot([], [], 'kv-', label='TRAIN:Freq. weigh. mean IU')
        


        plt.xlabel('iterations')
        self.handles2, self.labels2 = self.ax2.get_legend_handles_labels()
        self.lgd2 = self.ax2.legend(self.handles2, self.labels2, loc='upper center', bbox_to_anchor=(0.5,-0.2))
        self.ax2.grid(True)    
        plt.draw()

    def writeNudges(self, outfile='jitter.txt'):

        counters = np.arange(len(self.x))
        bjds = self.camera.counterToBJD(counters)
        time = bjds - np.min(bjds)
        plt.figure('jitter timeseries')
        gs = gridspec.GridSpec(2, 1, hspace=0.15)
        kw = dict(linewidth=2)
        ax = None

        for i, what in enumerate((self.x, self.y)):
            ax = plt.subplot(gs[i], sharex=ax, sharey=ax)
            ax.plot(time, what, **kw)
            ax.set_ylabel(['dRA (arcsec)', 'dDec (arcsec)'][i])
            if i == 0:
                ax.set_title('Jitter Timeseries from\n{}'.format(self.basename))

        plt.xlabel('Time from Observation Start (days)')
        plt.xlim(np.min(time), np.max(time))
        plt.draw()
        plt.savefig(outfile.replace('.txt', '.pdf'))

        data = [counters, bjds, self.x, self.y]
        names = ['imagenumber', 'bjd', 'arcsecnudge_ra', 'arcsecnudge_dec']

        t = astropy.table.Table(data=data, names=names)
        t.write(outfile.replace('.txt', '_amplifiedby{}.txt'.format(self.amplifyinterexposurejitter)),
                format='ascii.fixed_width', delimiter=' ')
        logger.info("save jitter nudge timeseries to {0}".format(outfile))

def matplotlib_set_plot(ax, plotter, outfile, default_camera=(14, -120),
                        hide_x=False, hide_y=False):
    ax.set_title(plotter.plot_title)

    tsize = 'medium'
    ax.set_xlabel(plotter.xaxis_label, fontsize=tsize)
    ax.set_ylabel(plotter.yaxis_label, fontsize=tsize)
    ax.set_zlabel(plotter.zaxis_label, fontsize=tsize)
    ax.ticklabel_format(axis='both', labelpad=150, useOffset=False)
    ax.set_xlim(*plotter.xaxis_range)
    ax.set_ylim(*plotter.yaxis_range)
    ax.set_zlim(*plotter.zaxis_range)
    ax.legend(fontsize='small')

    # getting a nice view over the whole mess in ppv
    ax.view_init(*default_camera)

    # hide axis-numbers:
    if hide_x:
        ax.get_xaxis().set_ticks([])
        ax.xaxis.set_visible(False)
        ax.get_xaxis().set_visible(False)
    if hide_y:
        ax.get_yaxis().set_ticks([])
        ax.yaxis.set_visible(False)
        ax.get_yaxis().set_visible(False)

    plt.draw()
    plt.savefig(outfile)
    plt.show()

def find_gates(mag1, mag2, param):
    col = mag1 - mag2

    lines = open(param, 'r').readlines()
    colmin, colmax = map(float, lines[4].split()[3:-1])
    mag1min, mag1max = map(float, lines[5].split()[:-1])
    #mag2min, mag2max = map(float, lines[5].split()[:-1])
    # click around
    fig, ax = plt.subplots()
    ax.plot(col, mag2, ',', color='k', alpha=0.2)
    ax.set_ylim(mag1max, mag1min)
    ax.set_xlim(colmin, colmax)

    ok = 1
    while ok == 1:
        print 'click '
        pts = np.asarray(plt.ginput(n=4, timeout=-1))
        exclude_gate = '1 {} 0 \n'.format(' '.join(['%.4f' % p for p in pts.flatten()]))
        pts = np.append(pts, pts[0]).reshape(5,2)
        ax.plot(pts[:,0], pts[:,1], color='r', lw=3, alpha=0.3)
        plt.draw()
        ok = move_on(0)
    lines[7] = exclude_gate
    # not so simple ... need them to be parallelograms.
    # PASS!

    # write new param file with exclude/include gate
    os.system('mv {0} {0}_bkup'.format(param))
    with open(param, 'w') as outp:
        [outp.write(l) for l in lines]
    print('wrote %s' % param)

def show_stat(net):
    plt.clf()

    f = plt.gcf()
    f.add_subplot('211')
    plt.title(net.checkpoint_name)
    plt.plot(net.stat['epoch'], net.stat['train']['error'], label='train')
    plt.plot(net.stat['epoch'], net.stat['val']['error'], label='val')
    plt.plot(net.stat['epoch'], net.stat['test']['error'], label='test')
    plt.legend(loc = 'lower left')
    plt.ylabel('error')
    plt.xlabel('epochs')
    plt.grid()

    f.add_subplot('212')
    plt.plot(net.stat['epoch'], net.stat['train']['cost'], label='train')
    plt.plot(net.stat['epoch'], net.stat['val']['cost'], label='val')
    plt.plot(net.stat['epoch'], net.stat['test']['cost'], label='test')
    plt.legend(loc = 'lower left')
    plt.ylabel('cost')
    plt.xlabel('epochs')
    plt.grid()

    plt.draw()
    plt.savefig(net.output_dir + 'stat.png')
    time.sleep(0.05)

def waterfall_plot(q,x,sampling=10,cmap=None,num_colors=100,outdir='./',outname='waterfall',format='eps',cbar_label='$|q| (a.u.)$'):
    plt.figure()
    plt.hold(True)
    colorVal = 'b'
    vmax = q[:,:].max()
    print vmax,len(q)
    for n in range(0,len(q),sampling):
        if cmap is not None:
            print q[n,:].max()
            colorVal = get_color(value=q[n,:].max(),cmap=cmap,vmax=vmax+.1,num_colors=num_colors)

        plt.plot(x,q[n,:]+n/10.0,label=str(n),color=colorVal,alpha=0.7)
    ax = plt.gca()
    for tic in ax.yaxis.get_major_ticks():
        tic.tick1On = tic.tick2On = False
        tic.label1On = tic.label2On = False

    if cmap is not None:
        scalar = get_smap(vmax=q[:,:].max()+.1,num_colors=sampling)
        cbar = plt.colorbar(scalar)

    plt.xlabel('$x\quad (a.u.)$')
    cbar.set_label(cbar_label)
    plt.draw()

    plt.savefig(os.path.join(outdir,outname+'.'+format),format=format,dpi=320,bbox_inches='tight')
    plt.close()
    return

def graphical_test(satisfactory=0):
    from matplotlib import cm, pylab
    def cons():
        return np.random.random(2)*4-2

    def foo(x,y,a,b):
        "banana function"
        tmp=a-x
        tmp*=tmp
        out=-x*x
        out+=y
        out*=out
        out*=b
        out+=tmp
        return out*(abs(np.cos((x-1)**2+(y-1)**2))+10.0/b)
    def f(params):
        return foo(params[0], params[1],1,100)

    optimizer=optimize(f, cons, verbose=False,its=1, hillWalks=0, satisfactory=satisfactory, finalWalk=0)
    
    bgx,bgy=np.mgrid[-2:2:1000j,-2:2:1000j]
    bg=foo(bgx,bgy, 1,100)
    for i in xrange(20):
        pylab.clf()
        pylab.imshow(bg, cmap=cm.RdBu,vmax=bg.mean()/10)
        for x in optimizer.pool: pylab.plot((x[2]+2)/4*1000,(x[1]+2)/4*1000, ('gx'))
        print optimizer.pool[0],optimizer.muterate
        pylab.gca().set_xbound(0,1000)
        pylab.gca().set_ybound(0,1000)
        pylab.draw()
        pylab.colorbar()
        optimizer.run()
        raw_input('enter to advance')
    return optimizer

def ttplot(corfp,srp,slp,n,I1,I2):
    global tplot_cum,dt,firstfile
    tplot=time.time()
    rchplot=int(ceil(log(n/chn)/log(2))+1)
    normplot=zeros((1,rcr),dtype=float32)
    
    for ir in xrange(rchplot):
       if ir==0:
           normplot[0,:chn]=1./arange(n-2,n-chn-2,-1)
       else:
           normplot[0,chn2*(ir+1.):chn2*(ir+2.)]=1./arange((n-1)/(2**ir)-chn2-1,(n-1)/(2**ir)-chn-1,-1)

    indt=int(chn+chn2*log(n/chn)/log(2))-2
    cc1=corfp[0,:indt]/(slp[0,:indt]*srp[0,:indt])/normplot[0,:indt]
    cc2=corfp[-1,:indt]/(slp[-1,:indt]*srp[-1,:indt])/normplot[0,:indt]
    t_axis=lag[0,:indt]
    t_axis2=tI_avg[0,:n]
    t_axis2b=tI_avg[0,:n]/dt+firstfile
    lm1.set_data(t_axis,cc1)
    lm2.set_data(t_axis2,I1)
    lm1b.set_data(t_axis,cc2)
    lm2b.set_data(t_axis2b,I2)
    ax1.set_xlim(min(t_axis),max(t_axis))
    ax1.set_ylim(min(cc1),max(cc1))
    ax1b.set_ylim(min(cc2),max(cc2))
    ax2.set_xlim(min(t_axis2),max(t_axis2))
    ax2b.set_xlim(min(t_axis2b),max(t_axis2b))
    ax2.set_ylim(min(I1),max(I1))
    ax2b.set_ylim(min(I2),max(I2))
    p.draw()
    tplot_cum+=time.time()-tplot
    return 

def plot(y, function):
    """ Show an animation of Poincare plot.

    --- arguments ---
    y: A list of initial values
    function: function which is argument of Runge-Kutta solver
    """
    h = dt
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.grid()
    time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
    plt.ion()

    for i in range(nmax + 1):
        for j in range(nstep):
            rk4 = RK.RK4(function)
            y = rk4.solve(y, j * h, h)
            # -pi <= theta <= pi
            while y[0] > pi:
                y[0] = y[0] - 2 * pi
            while y[0] < -pi:
                y[0] = y[0] + 2 * pi

        if ntransient <= i < nmax:          # <-- draw the poincare plots
            plt.scatter(y[0], y[1], s=2.0, marker='o', color='blue')
            time_text.set_text('n = %d' % i)
            plt.draw()

        if i == nmax:                       # <-- to stop the interactive mode
            plt.ioff()
            plt.scatter(y[0], y[1], s=2.0, marker='o', color='blue')
            time_text.set_text('n = %d' % i)
            plt.show()

def label_data(prefix, size=100, savename=None):
    from glob import glob
    from os.path import basename
    from PIL import Image
    from os.path import isfile
    if savename==None: savename=labelpath+'label_'+prefix+'.txt'
    # We want to avoid labeling an image twice, so keep track
    # of what we've labeled in previous labeling sessions.
    if isfile(savename):
        fileout = open(savename,'r')
        already_seen = [line.split(',')[0] for line in fileout]
        fileout.close()
    else: already_seen = []
    # Now reopen the file for appending.
    fileout = open(savename,'a')
    pl.ion()
    pl.figure(1,figsize=(9,9))
    files = glob(imgpath+prefix+'*.png')
    for file in np.random.choice(files, size=size, replace=False):
        if basename(file) in already_seen: continue
        pl.clf()
        pl.subplot(1,1,1)
        pl.imshow(np.array(Image.open(file)))
        pl.title(file)
        pl.axis('off')
        pl.draw()
        label = get_one_char()
        if label=='q': break
        fileout.write(basename(file)+','+label+'\n')
        print file,label
    fileout.close()
    return

def pick(event) :

    global data
    global X
    global Y
    global numpy_container
    if event.key == 'q' :
        if len(X) == 0 : return
        if not numpy_container :
            data = VectorDataSet(X)
        else :
            data = PyVectorDataSet(numpy.array(X))
        data.attachLabels(Labels(Y))
        X = []
        Y = []
        print 'done creating data.  close this window and use the decisionSurface function'
        pylab.disconnect(binding_id)
    if event.key =='1' or event.key == '2' :
        if event.inaxes is not None:
            print 'data coords', event.xdata, event.ydata
            X.append([event.xdata, event.ydata])
            Y.append(event.key)
            pylab.plot([event.xdata], [event.ydata], 
                       plotStr[int(event.key) - 1])
            pylab.draw()

    def plot(self):
        Ns = 1
        h = 0.1
        m = []
        for i in range(len(self.l)):
            for s in range(Ns):
                m.append(ml_to_xy((i, float(s)/Ns, 0), self.kappa, self.l, self.x0, self.theta0))
        m.append(ml_to_xy((len(self.l)-1, 1., 0), self.kappa, self.l, self.x0, self.theta0))

        plt.clf()
        plt.hold(True)
        x=[p[0] for p in m]
        y=[p[1] for p in m]
        
        bp, angles, tangents= breakpoints(self.kappa, self.l, self.x0, self.theta0)
        x=[p[0] for p in bp]
        y=[p[1] for p in bp]
        plt.plot(x[::10],y[::10],'k-')
#        xm = [p[0] for i, p in enumerate(bp) if self.markers[i]==1]
#        ym = [p[1] for i, p in enumerate(bp) if self.markers[i]==1]


#        plt.plot(xm,ym,'kx')

        plt.xlim((-100, 3000))
        plt.ylim((-3000, 100))
        plt.draw()

def runcand(doplot=False):
    
    import time
    l = os.listdir("BenLike/")
    m = open("features.csv","w")
    has_run = False
    for f in l:
        if f.find(".xml") != -1:
            fname = "BenLike/" + f
            print "working on", f
            x0,y,dy, name = _load_dotastro_data(fname)
            a = ebfeature(t=x0,m=y,merr=dy,name=name)
            a.gen_orbital_period(doplot=doplot)
            if doplot:
                plt.draw()
            if not has_run:
                ff = a.features.keys()
                ff.remove("run")
                ff.remove("p_pulse_initial")
                m.write("name," + ",".join(ff) + "\n")
                has_run = True
            
            m.write(os.path.basename(name) + "," + ",".join([str(a.features.get(s)) for s in ff]) + "\n")
            time.sleep(1)
    m.close()

def axes_animate_out(ax, iterations=10, max_animation_time=1):
    bottom, top = ax.get_ylim()

    alpha_objects = ax.findobj(has_alpha)
    alpha_divider = np.power(float(iterations-1), 2)

    def ease_out(cur_index, end_index, func):
        return 1. / (func(cur_index) / float(func(end_index)))

    start_time = time()
    for i in range(0, iterations):
        c = ease_out(iterations - i, iterations, roll_index_transformation)
        ax.set_ylim([bottom, c*top])

        for item in alpha_objects:
            item.set_alpha(alpha_index_transformation(iterations-i-1) / alpha_divider)

        if time() - start_time > max_animation_time:
            for item in alpha_objects:
                item.set_alpha(0)

            pl.draw()
            return

        pl.draw()