def plotConn():
    # Create plot
    figh = figure(figsize=(8,6))
    figh.subplots_adjust(left=0.02) # Less space on left
    figh.subplots_adjust(right=0.98) # Less space on right
    figh.subplots_adjust(top=0.96) # Less space on bottom
    figh.subplots_adjust(bottom=0.02) # Less space on bottom
    figh.subplots_adjust(wspace=0) # More space between
    figh.subplots_adjust(hspace=0) # More space between
    h = axes()
    totalconns = zeros(shape(f.connprobs))
    for c1 in range(size(f.connprobs,0)):
        for c2 in range(size(f.connprobs,1)):
            for w in range(f.nreceptors):
                totalconns[c1,c2] += f.connprobs[c1,c2]*f.connweights[c1,c2,w]*(-1 if w>=2 else 1)
    imshow(totalconns,interpolation='nearest',cmap=bicolormap(gap=0))

    # Plot grid lines
    hold(True)
    for pop in range(f.npops):
        plot(array([0,f.npops])-0.5,array([pop,pop])-0.5,'-',c=(0.7,0.7,0.7))
        plot(array([pop,pop])-0.5,array([0,f.npops])-0.5,'-',c=(0.7,0.7,0.7))

    # Make pretty
    h.set_xticks(range(f.npops))
    h.set_yticks(range(f.npops))
    h.set_xticklabels(f.popnames)
    h.set_yticklabels(f.popnames)
    h.xaxis.set_ticks_position('top')
    xlim(-0.5,f.npops-0.5)
    ylim(f.npops-0.5,-0.5)
    clim(-abs(totalconns).max(),abs(totalconns).max())
    colorbar()

    def plot_field(self,x,y,u=None,v=None,F=None,contour=False,outdir=None,plot='quiver',figname='_field',format='eps'):
        outdir = self.set_dir(outdir)
        p = 64
        if F is None: F=self.calc_F(u,v)

        plt.close('all')
        plt.figure()
        #fig, axes = plt.subplots(nrows=1)
        if contour:
            plt.hold(True)
            plt.contourf(x,y,F)
        if plot=='quiver':
            plt.quiver(x[::p],y[::p],u[::p],v[::p],scale=0.1)
        
        if plot=='pcolor':
            plt.pcolormesh(x[::4],y[::4],F[::4],cmap=plt.cm.Pastel1)
            plt.colorbar()

        if plot=='stream':
            speed = F[::16]
            plt.streamplot(x[::16], y[::16], u[::16], v[::16], density=(1,1),color='k')

        plt.xlabel('$x$ (a.u.)')
        plt.ylabel('$y$ (a.u.)')
        
        plt.savefig(os.path.join(outdir,figname+'.'+format),format=format,dpi=320,bbox_inches='tight')
        plt.close()

def startAnim(x, m, th, Tsim, inter=1, Tstart=0, h=0.002):
#    fig, axM = subplo1ts()
#    axX = axM.twinx()
    fig = figure()
    axM = subplot(211)
    axX = subplot(212, sharex=axM)
    anim = MyAnim(x, m, th, Tsim, inter, Tstart/h, h)

    anim.line1, = axM.plot([], [], 'b')
    anim.line2, = axX.plot([], [], 'g')
    setp(axM.get_xticklabels(), visible=False)
    axM.set_xlim([-pi, pi])
    axX.set_xlim([-pi, pi])
#    axM.set_ylim([0., 3])
#    axX.set_ylim([0., 1.])
    axM.set_ylim([amin(m), amax(m)])
    axX.set_ylim([amin(x), amax(x)])

    axM.set_ylabel(r"$m$")
    axX.set_ylabel(r"$x$")
    axX.set_xlabel(r"$\theta$")
    hold(False)
    for tl in axM.get_yticklabels():
        tl.set_color('b')
    for tl in axX.get_yticklabels():
        tl.set_color('g')
    anim.axM = axM
    anim.axX = axX
    fig.canvas.mpl_connect('button_press_event', anim.onClick)
    a = FuncAnimation(fig, anim, frames=anim.dataGen, init_func=anim.init,
                 interval=0, blit=True, repeat=True)
    show()
    return a

def plot_waveforms(time,voltage,APTimes,titlestr):
    """
    plot_waveforms 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 waveforms for each
    detected action potential
    """

    plt.figure()

    ## Your Code Here
    for n in range(0,APTimes.size):
        ind = time[(time>(APTimes[n] - 0.003)) & (time<(APTimes[n] + 0.003))]
        sp = np.zeros(ind.size)
    
        for i in range(0,ind.size):
            sp[i] = plt.find(time == ind[i])
        sp=sp.astype(np.int64)
        x = np.linspace(-3.0e-3, 3.0e-3, sp.size)
        plt.plot(x,voltage[sp])
        plt.hold(True)
        
    plt.xlabel('Time (s)')
    plt.ylabel('Voltages (s)')
    plt.title(titlestr)    
    plt.show()

def saveCMDataFig(CMdata, trialDuration, trialReps, saveDir, plotName, timeStr):
    plt.figure(1, figsize=(14, 11), dpi=80)
    plt.clf()
    
    numFreq = len(CMdata.freqArray)
    numAmp = len(CMdata.ampArray)
    nGraphs = numFreq + 1
    numRows = int(np.ceil(nGraphs ** 0.5))
    numCols = int(np.ceil(nGraphs / numRows))
    numSD = 2.5
    for n in range(0, numFreq):
        CMresp = CMdata.CMResp[n, :]
        noise_mean = CMdata.noise_mean[n, :]
        noise_std = CMdata.noise_std[n, :]
        noise = noise_mean + (numSD * noise_std)
        
        plt.subplot(numRows, numCols, n+1)
        plt.plot(CMdata.ampArray, 1e6*CMresp,'-bo', label='Response (uV)')
        plt.hold('on')
        plt.plot(CMdata.ampArray,  1e6*(noise_mean + 3*noise_std), '-ro', label='Noise')   
        plt.plot(CMdata.ampArray, noise_mean, ':r', label='Noise mean')
        # plt.plot(CMdata.ampArray, noise, '-rs', label='+ %g SD' % numSD)
        
        plt.xlabel('Amplitude (dB)', fontsize=10)
        if n == (numFreq-1):
            plt.legend(loc='upper left', fontsize=10)
        #plt.pcolormesh(t, f, Sxx)
        plt.ylabel('dB SPL', fontsize=10)
        plt.title('%0.2f kHz' % (CMdata.freqArray[n]/1e3), x=0.15, fontsize=10)
        
    plt.show()
    fname = os.path.join(saveDir, plotName)
    plt.savefig(fname)

def ensemble_demo():
    utc = Calendar()
    t_start = utc.time(YMDhms(2011, 9, 1))
    t_fc_ens_start = utc.time(YMDhms(2015, 7, 26))
    disp_start = utc.time(YMDhms(2015, 7, 20))
    dt = deltahours(1)
    n_obs = int(round((t_fc_ens_start - t_start)/dt))
    n_fc_ens = 30
    n_disp = int(round(t_fc_ens_start - disp_start)/dt) + n_fc_ens + 24*7

    obs_time_axis = Timeaxis(t_start, dt, n_obs + 1)
    fc_ens_time_axis = Timeaxis(t_fc_ens_start, dt, n_fc_ens)
    display_time_axis = Timeaxis(disp_start, dt, n_disp)

    q_obs_m3s_ts = observed_tistel_discharge(obs_time_axis.total_period())
    ptgsk = create_tistel_simulator(PTGSKOptModel, tistel.geo_ts_repository(tistel.grid_spec.epsg()))
    initial_state = burn_in_state(ptgsk, t_start, utc.time(YMDhms(2012, 9, 1)), q_obs_m3s_ts)

    ptgsk.run(obs_time_axis, initial_state)
    current_state = adjust_simulator_state(ptgsk, t_fc_ens_start, q_obs_m3s_ts)
    q_obs_m3s_ts = observed_tistel_discharge(display_time_axis.total_period())
    ens_repos = tistel.arome_ensemble_repository(tistel.grid_spec)
    ptgsk_fc_ens = create_tistel_simulator(PTGSKModel, ens_repos)
    sims = ptgsk_fc_ens.create_ensembles(fc_ens_time_axis, t_fc_ens_start, current_state)
    for sim in sims:
        sim.simulate()
    plt.hold(1)
    percentiles = [10, 25, 50, 75, 90]
    plot_percentiles(sims, percentiles, obs=q_obs_m3s_ts)
    #plt.interactive(1)
    plt.show()

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 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 plot_waveforms(time,voltage,APTimes,titlestr):
    """
    plot_waveforms 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 waveforms for each
    detected action potential
    """
   
    plt.figure()
   
    ## Your Code Here 
    indices = []
    
    for x in range(len(APTimes)):
        for i in range(len(time)):
            if(time[i]==APTimes[x]):
                indices.append(i)
            

    ##print indices
    Xval = np.linspace(-.003,.003,200)
    print len(Xval)
    for x in range(len(APTimes)):
        plt.plot(Xval, voltage[indices[x]-100:indices[x]+100])
        plt.title(titlestr)
        plt.xlabel('Time (s)')
        plt.ylabel('Voltage (uV)')
        plt.hold(True)

    
    
    plt.show()

def plot_svc(X, y, mysvc, bounds=None, grid=50):
    if bounds is None:
        xmin = np.min(X[:, 0], 0)
        xmax = np.max(X[:, 0], 0)
        ymin = np.min(X[:, 1], 0)
        ymax = np.max(X[:, 1], 0)
    else:
        xmin, ymin = bounds[0], bounds[0]
        xmax, ymax = bounds[1], bounds[1]
    aspect_ratio = (xmax - xmin) / (ymax - ymin)
    xgrid, ygrid = np.meshgrid(np.linspace(xmin, xmax, grid),
                              np.linspace(ymin, ymax, grid))
    plt.gca(aspect=aspect_ratio)
    plt.xlim(xmin, xmax)
    plt.ylim(ymin, ymax)
    plt.xticks([])
    plt.yticks([])
    plt.hold(True)
    plt.plot(X[y == 1, 0], X[y == 1, 1], 'bo')
    plt.plot(X[y == -1, 0], X[y == -1, 1], 'ro')
    
    box_xy = np.append(xgrid.reshape(xgrid.size, 1), ygrid.reshape(ygrid.size, 1), 1)
    if mysvc is not None:
        scores = mysvc.decision_function(box_xy)
    else:
        print 'You must have a valid SVC object.'
        return None;
    
    CS=plt.contourf(xgrid, ygrid, scores.reshape(xgrid.shape), alpha=0.5, cmap='jet_r')
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[0], colors='k', linestyles='solid', linewidths=1.5)
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[-1,1], colors='k', linestyles='dashed', linewidths=1)
    plt.plot(mysvc.support_vectors_[:,0], mysvc.support_vectors_[:,1], 'ko', markerfacecolor='none', markersize=10)
    CB = plt.colorbar(CS)

def plot_fits(direction_rates,fit_curve,title):
    """
    This function takes the x-values and the y-values  in units of spikes/s 
    (found in the two columns of direction_rates and fit_curve) and plots the 
    actual values with circles, and the curves as lines in both linear and 
    polar plots.
    """
    curve_xs = np.arange(direction_rates[0,0], direction_rates[-1,0])
    fit_ys2 = normal_fit(curve_xs,fit_curve[0],fit_curve[1],fit_curve[2])
    
    
    plt.subplot(2,2,3)
    plt.plot(direction_rates[:,0],direction_rates[:,1],'o',hold=True)
    plt.plot(curve_xs,fit_ys2,'-')
    plt.xlabel('Direction of Motions (Degrees)')
    plt.ylabel('Firing Rates (Spikes/sec)')
    plt.title(title)
    plt.axis([0, 360, 0, 40])
    plt.xticks(direction_rates[:,0])
    
    fit_ys = normal_fit(direction_rates[:,0],fit_curve[0],fit_curve[1],fit_curve[2])
    plt.subplot(2,2,4,polar=True)
    spkiecount = np.append(direction_rates[:,1],direction_rates[0,1])
    plt.polar(np.arange(0,361,45)*np.pi/180,spkiecount,'o',label='Firing Rate (spike/s)')
    plt.hold(True)
    spkiecount_y = np.append(fit_ys,fit_ys[0])
    plt.plot(np.arange(0,361,45)*np.pi/180,spkiecount_y,'-')    
    plt.legend(loc=8)
    plt.title(title)
    
    fit_ys2 = np.transpose(np.vstack((curve_xs,fit_ys2)))    
    
    return(fit_ys2)

def mark_cross(center, **kwargs):
    """Mark a cross. Correct for matplotlib imshow funny coordinate system.
    """
    N = 20
    plt.hold(1)
    plt.axhline(y=center[1]-0.5, **kwargs)
    plt.axvline(x=center[0]-0.5, **kwargs)

def saveABRDataFig(ABRdata, ABRparams, saveDir, plotName, timeStr):
    plt.figure(1, figsize=(14, 11), dpi=80) 
    plt.clf()
    
    numFreq = len(ABRdata.freqArray)
    numAmp = len(ABRdata.ampArray)
    numSD = 5
    # clr_lbls = ['b', 'g', 'r', 'c', 'm',  'y', 'k']
    nGraphs = numFreq + 1
    numRows = int(np.ceil(nGraphs ** 0.5))
    numCols = int(np.ceil(nGraphs / numRows))
    for n in range(0, numFreq):
        plt.subplot(numRows, numCols, n+1)
        ABRresp = 1e6*ABRdata.ABRResp[n, :]
        noise = 1e6*numSD*ABRdata.noise_std[n, :]
        plt.plot(ABRdata.ampArray, ABRresp, '-bo', label='Signal')
        plt.hold('on')
        plt.plot(ABRdata.ampArray, noise, '-r', label='Noise (%0.1f SD)' % numSD)
        if n == (0):
            plt.ylabel('Resp PP (uV)', fontsize=10)
        if n == (numFreq-1):
            plt.legend(loc='upper left', fontsize=10)
            plt.xlabel('Amplitude (dB)', fontsize=10)
        #plt.pcolormesh(t, f, Sxx)
        
        plt.title('%0.2f kHz' % (ABRdata.freqArray[n]/1e3), x=0.15, fontsize=10)

    plt.show()
    fname = os.path.join(saveDir, plotName)
    plt.savefig(fname)
    
    plt.figure(2, figsize=(14, 11), dpi=80)
    plt.clf()
    
    # clr_lbls = ['b', 'g', 'r', 'c', 'm',  'y', 'k']
    t = ABRdata.t
    for n in range(0, numFreq):
        plt.subplot(numRows, numCols, n+1)
        offset = 0
        for a in range(numAmp):
            tr = 1e6*ABRdata.tracings[n, a, :]
            if a > 0:
                offset = offset + np.abs(np.min(tr))

            plt.plot(t*1e3, tr+offset, '-b', label='Signal')
            plt.hold('on')
            offset = offset + np.max(tr)
            
        if n == (0):
            plt.ylabel('Resp (uV)', fontsize=10)
        if n == (numFreq-1):
            plt.xlabel('Time (ms)', fontsize=10)
        #plt.pcolormesh(t, f, Sxx)
        
        plt.title('%0.2f kHz' % (ABRdata.freqArray[n]/1e3), x=0.15, fontsize=10)
    
    plt.show()
    fname = os.path.join(saveDir, plotName + ' tracings')
    plt.savefig(fname)

def gasket(pa, pb, pc, level, col):
    if level == 0:
        plt.fill([pa[0], pb[0], pc[0]], [pa[1], pb[1], pc[1]], col) 
        plt.hold(True)
    else:
        gasket(pa, (pa + pb) / 2., (pa + pc) / 2., level - 1, col) 
        gasket(pb, (pb + pa) / 2., (pb + pc) / 2., level - 1, col) 
        gasket(pc, (pc + pa) / 2., (pc + pb) / 2., level - 1, col)

def Sierpinski(a,b,c,k,iteration): 
    if iteration==0:
        plt.fill([a[0], b[0], c[0]], [a[1], b[1], c[1]],'b') 
        plt.hold(True)
    else:
        Sierpinski(a,bis(a,b,k),bis(a,c,k),k,iteration-1)
        Sierpinski(b,bis(a,b,k),bis(b,c,k),k,iteration-1)
        Sierpinski(c,bis(a,c,k),bis(b,c,k),k,iteration-1)
        plt.hold(True)

 def display_multiple_egg_carton(u0,s0,s1,ds,C,s):
     X = solve_multiple(C,u0,s0,s1,ds,s)
     plt.hold('on')
     # plot the parametrized data on to the egg carton
     u,v = X[:,0], X[:,2]
     x = u
     y = v
     z = np.sin(u)*np.cos(v)
     ax.plot(x,y,z,'--r')
     return plt

def plotPolygon(lat,lon,poly,regnum,plotFile):
    fig = plt.figure(figsize=(6,6))
    px = [pol[0] for pol in poly]
    py = [pol[1] for pol in poly]
    plt.plot(px,py,'b')
    plt.hold(True)
    plt.plot(lon,lat,'rx')
    plt.title('%.4f,%.4f in Region %i' % (lat,lon,regnum))
    plt.savefig(plotFile)
    plt.close()    

 def display_multiple_catenoid(u0,s0,s1,ds,C,s):
     X = solve_multiple(C,u0,s0,s1,ds,s)
     plt.hold('on')
     # plot the parametrized data on to the catenoid
     u,v = X[:,0], X[:,2]
     x = np.cos(u) - v*np.sin(u)
     y = np.sin(u) + v*np.cos(u)
     z = v
     ax.plot(x,y,z,'--r')
     return plt

def sierpinski(a, b, c, iterations):
    x=(random.random(),random.random(),random.random())
    if iterations  == 0:
        plt.fill([a[0], b[0], c[0]], [a[1], b[1], c[1]], color=x,alpha=0.9) 
        plt.hold(True)
    else:
        sierpinski(a, (a + b) / 2., (a + c) / 2., iterations  - 1) 
        sierpinski(b, (b + a) / 2., (b + c) / 2., iterations  - 1) 
        sierpinski(c, (c + a) / 2., (c + b) / 2., iterations  - 1)
        plt.fill([(a[0] + b[0]) / 2., (a[0] + c[0]) / 2., (b[0] + c[0]) / 2.], [(a[1] + b[1]) / 2., (a[1] + c[1]) / 2., (b[1] + c[1]) / 2.], color=x,alpha=0.9)

 def display_multiple_torus(u0,s0,s1,ds,C,s):
     X = solve_multiple(C,u0,s0,s1,ds,s)
     plt.hold('on')
     # plot the parametrized data on to the sphere
     u,v = X[:,0], X[:,2]
     x = (2 + 1*np.cos(v))*np.cos(u)
     y = (2 + 1*np.cos(v))*np.sin(u)
     z = np.sin(v)
     ax.plot(x,y,z,'--r')
     return plt