def plot_bars(pos_count, title='', max_pathway_length=8, legend_loc='upper right'):
    n_labels = len(pos_count)
    ind = np.arange(max_pathway_length)
    width = 0.2
    
    fig = pylab.figure()
    pylab.hold(True)
    ax = fig.add_subplot(111)

    colors = {'No known regulation':'grey', 'Activated':'green', 'Inhibited':'red', 'Mixed regulation':'blue'}
    plot_order = ['Inhibited', 'Mixed regulation', 'Activated', 'No known regulation']    
    
    i = 0
    for label in plot_order:
        curr_vals = pos_count[label][1:max_pathway_length+1]
        if (sum(curr_vals) < 20):
            n_labels -= 1
            continue
        ax.bar(ind + i * width, tuple([j * 1.0 /sum(curr_vals) for j in curr_vals]), width, color=colors[label], label=('%s (%d)' % (label, sum(curr_vals))))
        i += 1
    
    ax.set_ylabel('Fraction of reactions per type')
    ax.set_xlabel('Position in pathway')
    
    ax.set_xticks(ind+ width * n_labels/2)
    ax.set_xticklabels( ind + 1 )
    
    legendfont = matplotlib.font_manager.FontProperties(size=11)
    pylab.legend(loc=legend_loc, prop=legendfont)
    pylab.title(title)

    pylab.hold(False)
    
    return fig

def plot_coord_mapping(mapper,sheet,style='b-'):
    """
    Plot a coordinate mapping for a sheet.

    Given a CoordinateMapperFn (as for a CFProjection) and a sheet
    of the projection, plot a grid showing where the sheet's units
    are mapped.
    """

    from pylab import plot,hold,ishold

    xs = sheet.sheet_rows()
    ys = sheet.sheet_cols()

    hold_on = ishold()
    if not hold_on:
        plot()
    hold(True)

    for y in ys:
        pts = [mapper(x,y) for x in xs]
        plot([u for u,v in pts],
             [v for u,v in pts],
             style)

    for x in xs:
        pts = [mapper(x,y) for y in ys]
        plot([u for u,v in pts],
             [v for u,v in pts],
             style)

    hold(hold_on)

def _test_graph():
    i = 10000
    x = np.linspace(0,3.7*pi,i)
    y = (0.3*np.sin(x) + np.sin(1.3 * x) + 0.9 * np.sin(4.2 * x) + 0.06 *
    np.random.randn(i))
    y *= -1
    x = range(i)

    _max, _min = peakdetect(y,x,750, 0.30)
    xm = [p[0] for p in _max]
    ym = [p[1] for p in _max]
    xn = [p[0] for p in _min]
    yn = [p[1] for p in _min]

    plot = pylab.plot(x,y)
    pylab.hold(True)
    pylab.plot(xm, ym, 'r+')
    pylab.plot(xn, yn, 'g+')

    _max, _min = peak_det_bad.peakdetect(y, 0.7, x)
    xm = [p[0] for p in _max]
    ym = [p[1] for p in _max]
    xn = [p[0] for p in _min]
    yn = [p[1] for p in _min]
    pylab.plot(xm, ym, 'y*')
    pylab.plot(xn, yn, 'k*')
    pylab.show()

def estimateHarness( delta = 1e-1,
                     alpha = .0,
                     num_samples=10,
                     Tf_sample = Tf ):
    #Load all the simulated trajectories
    file_name = os.path.join(RESULTS_DIR,
                              'OU_Xs.a=%.3f_N=%d.npy'%(alpha,
                                                       num_samples));
    trajectoryBank = load(file_name)
    
    #Select an arbitrary trajectory: (here the 2nd)
    figure(); hold(True);
    n_thin = int(delta / dt); print n_thin
    N_sample = int(Tf_sample / dt) 
#    for idx in xrange(1,10):
#    for idx in [2]: #xrange(3,4):
    for idx in xrange(1,num_samples+1):
        ts, Xs = trajectoryBank[:N_sample,0], trajectoryBank[:N_sample,idx]
    
        #Select sampling rate:    
        #Generate sampled data, by sub-sampling the fine trajectory:    
        ts_thin = ts[::n_thin];
        Xs_thin = Xs[::n_thin];
        
        #Obtain estimator
#        est_params = estimateParams(Xs_thin, delta)
#        print 'est original: %.4f,%.4f, %.4f'%(est_params[0],est_params[1],est_params[2])
        est_params = estimateParamsBeta(Xs_thin, delta, alpha)
        print 'est reduced: %.4f,%.4f, %.4f'%(est_params[0],est_params[1],est_params[2])
        plot(ts_thin, Xs_thin);
         
    print 'true param values:', [mu, beta, sigma]

def plot_arm_speed(axis, startTime=-1):
    rootName = 'siemensSensors'
    f = netcdf.netcdf_file(rootName+'Data.nc', 'r')
    data1 = f.variables[rootName+'.data.'+'carouselSpeedSetpoint'].data[startSample:]
    data2 = f.variables[rootName+'.data.'+'carouselSpeedSmoothed'].data[startSample:]
    ts_trigger = f.variables[rootName+'.data.ts_trigger'].data[startSample:]*1.0e-9

    # Load the actual arm speed from the arm gyro
    rootName = 'armboneLisaSensors'
    fiile = netcdf.netcdf_file(rootName+'Data.nc', 'r')

    rawdata4 = fiile.variables['armboneLisaSensors.GyroState.gr'].data[startSample:]
    ts_trigger4 = fiile.variables['armboneLisaSensors.GyroState.ts_trigger'].data[startSample:]*1.0e-9
    #fullscale = 2000 # deg/sec
    #data4 = -1.0 * rawdata4 / (2**15) * fullscale * pi/180 - 0.0202 # Rad/s
    data4 = rawdata4

    if startTime == -1:
        startTime = ts_trigger[0]

    times = ts_trigger-startTime
    times4 = ts_trigger4-startTime

    pylab.hold(True)

    plot(times, data2, '.-', label='On Motor Side of Belt')
    plot(times4, data4,'.-',  label='From Gyro on Arm')
    plot(times, data1, '.-', label='Setpoint (Echoed)')
    ylabel('Arm rotation speed [Rad/s]')
    xlabel('Time [s]')
    #legend(['Setpoint (Echoed)', 'Setpoint (Sent)', 'On Motor Side of Belt', 'From Gyro on Arm'])
    title('Plot of Signals Related to Arm Speed')
    return startTime

def DRIVplot(folder,keys):
  T = 281
  APiterator = [5,10]
  AP = Analysis.AnalyseFile()
  P = Analysis.AnalyseFile()
  if folder[0]['IVtemp'] == T:
    scale = 1e6
    plt.hold(True)
    plt.title('NLIV in P and AP at ' + str(T) + 'K')
    plt.xlabel('Current ($\mu$A)')
    plt.ylabel('V$_{NL}$ ($\mu$V)')
    for f in folder:
      if f['iterator'] in APiterator:
        AP.add_column(f.Voltage,str(f['iterator']))
      else:
        P.add_column(f.Voltage,str(f['iterator']))        
    AP.apply(func,0,replace=False,header='Mean NLVoltage')
    P.apply(func,0,replace=False,header='Mean NLVoltage')    
    
    I = numpy.arange(-295e-6,295e-6,1e-6)
    
    ap = interpolate.interp1d(f.column('Current'),AP.column('Mean NLV'))    
    p = interpolate.interp1d(f.column('Current'),P.column('Mean NLV')) 
    
    print P
    plt.title(' ',verticalalignment='bottom')
    plt.xlabel('Current ($\mu$A)')
    #plt.ylabel('V$_{NL}$/|I| (V/A)')
    plt.ylabel('$\Delta$V$_{NL}$/|I| (mV/A)') 
    plt.plot(f.column('Current')*scale,1e3*(P.column('Mean NLV')-AP.column('Mean NLV'))/abs(f.column('Current')),label =''+str(T)+ ' K')
    #plt.plot(f.column('Current')*scale,1e3*(P.column('Mean NLV'))/abs(f.column('Current')),label ='P at '+str(T)+ ' K')
    #plt.plot(f.column('Current')*scale,1e3*(AP.column('Mean NLV'))/abs(f.column('Current')),label ='AP at '+str(T)+ ' K')        
    plt.legend(loc='upper left')
  else:
    return 1  

def NormDeltaRvT(folder,keys):
  if folder[0]['IVtemp']<250 and folder[0]['IVtemp']>5:
    APiterator = [5,10]
    AP = Analysis.AnalyseFile()
    P = Analysis.AnalyseFile()
    tsum = 0.0
    for f in folder:
      if f['iterator'] in APiterator:
        AP.add_column(f.column('Voltage'),str(f['iterator']))
      else:
        P.add_column(f.column('Voltage'),str(f['iterator']))
      tsum = tsum + f['Sample Temp']
      
    AP.apply(func,0,replace=False,header='Mean NLV')
    AP.add_column(f.Current,column_header = 'Current')
    P.apply(func,0,replace=False,header='Mean NLV')
    P.add_column(f.Current,column_header = 'Current')
    
    APfit= AP.curve_fit(quad,'Current','Mean NLV',bounds=lambda x,y:x,result=True,header='Fit',asrow=True)
    Pfit = P.curve_fit(quad,'Current','Mean NLV',bounds=lambda x,y:x,result=True,header='Fit',asrow=True)
    
    DeltaR = Pfit[2] - APfit[2]
    ErrDeltaR = numpy.sqrt((Pfit[3]**2)+(APfit[3]**2))
    Spinsig.append(DeltaR/Res_Cu(tsum/10))
    Spinsig_error.append(ErrDeltaR)
    
    Temp.append(tsum/10)
    
    plt.hold(True)
    plt.title('$\Delta$R$_s$ vs T from linear coef of\nNLIV fit for '+f['Sample ID'],verticalalignment='bottom')
    plt.xlabel('Temperture (K)')
    plt.ylabel(r'$\Delta$R$_s$/$\rho$')
    plt.errorbar(f['IVtemp'],1e3*DeltaR,1e3*ErrDeltaR,ecolor='k',marker='o',mfc='r', mec='k')
    #plt.plot(f['IVtemp'],ErrDeltaR,'ok')
    return Temp, Spinsig

def _smooth_demo():
    from numpy import linspace, sin, ones
    from pylab import subplot, plot, hold, axis, legend, title, show, randn

    t = linspace(-4, 4, 100)
    x = sin(t)
    xn = x + randn(len(t)) * 0.1
    y = smooth(x)

    ws = 31

    subplot(211)
    plot(ones(ws))

    windows = ["flat", "hanning", "hamming", "bartlett", "blackman"]

    hold(True)
    for w in windows[1:]:
        eval("plot(" + w + "(ws) )")

    axis([0, 30, 0, 1.1])

    legend(windows)
    title("The smoothing windows")
    subplot(212)
    plot(x)
    plot(xn)
    for w in windows:
        plot(smooth(xn, 10, w))
    l = ["original signal", "signal with noise"]
    l.extend(windows)

    legend(l)
    title("Smoothing a noisy signal")
    show()

 def plotResults(self, titlestr="", ylimits=[0.5,1.05], plotfunc = pl.semilogx, ylimitsB=[0,101],
                  legend_loc=3, show=True ):
     pl.figure(num=None, figsize=(15,5))
     xvals = range(1, (1+len(self.removed)) )
     #Two subplots. One the left is the test accuracy vs. iteration
     pl.subplot(1,2,1)
     plotfunc(xvals, self.test_acc_list, "b", label="Test Accuracy")    
     pl.hold(True)
     plotfunc(xvals, self.getRollingAvgTestAcc(window_size=10), "r", label="Test Acc (rolling avg)")
     plotfunc(xvals, self.getRollingAvgTrainAcc(window_size=10), "g--", label="Train Acc (rolling avg)")
     pl.ylim(ylimits)
     if titlestr == "":
         pl.title("Iterative Feature Removal")
     else:
         pl.title(titlestr)
     pl.ylabel("Test Accuracy")
     pl.xlabel("Iteration")
     pl.legend(loc=legend_loc) #3=lower left
     pl.hold(False)
     
     #second subplot. On the right is the number of features removed per iteration
     pl.subplot(1,2,2)
     Ns = [ len(lst) for lst in self.removed ]
     pl.semilogx(xvals, Ns, "bo", label="#Features per Iteration")
     pl.xlabel("Iteration")
     pl.ylabel("Number of Features Selected")
     pl.title("Number of Features Removed per Iteration")
     pl.ylim(ylimitsB)
     
     pl.subplots_adjust(left=0.05, bottom=0.15, right=0.95, top=0.90, wspace=0.20, hspace=0.20)
     if show: pl.show()

def plotRes_varyingTrees( data_dict, dataset_name, max_correct=3000 , show=True):
    '''
    Plots the results of a varyingNumTrees() experiment, using a dictionary
    structure to hold the data. See the loadRes_varyingTrees() comments on the
    dictionary layout.
    '''
    xvals = data_dict['NumTrees']
    
    #prox forest trials
    pf_avg = data_dict['PF'].mean(axis=0)
    pf_std = data_dict['PF'].std(axis=0)
    pf_95_conf = 1.96 * pf_std / math.sqrt(data_dict['PF'].shape[0])

    #kdt forest trials
    kdt_avg = data_dict['KDT'].mean(axis=0)
    kdt_std = data_dict['KDT'].std(axis=0)
    kdt_95_conf = 1.96 * kdt_std / math.sqrt(data_dict['KDT'].shape[0])
    
    #plot average results of each, bounded by lower and upper bounds of 95% conf intv
    pl.hold(True)
    pl.errorbar(xvals, pf_avg/max_correct, yerr=pf_95_conf/max_correct, fmt='-r', 
                label="PF")
    pl.errorbar(xvals, kdt_avg/max_correct, yerr=kdt_95_conf/max_correct, fmt='-.b',
                label="KDT")
    pl.ylim([0,1.05])
    pl.title(dataset_name)
    pl.xlabel("Number of Trees in Forest")
    pl.ylabel("Percent Correct")
    pl.legend(loc='lower right')
    if show: pl.show()

def test_radial_profiles():
    arr = random_periodic_upsample(128, 16, seed=0)
    mask = np.zeros(arr.shape, dtype=np.bool_)
    arr_x = vcalc.cderivative(arr, 'X_DIR')
    arr_y = vcalc.cderivative(arr, 'Y_DIR')
    arr_div = np.sqrt(arr_x**2 + arr_y**2)
    surf = _cp.TopoSurface(arr)
    rprofs = radial_profiles(surf, threshold=25, expand_regions=1, other_arr=arr_div, mask=mask)
    arr[mask] = 2 * arr.max()
    pl.imshow(arr, interpolation='nearest')
    pl.figure()
    pl.imshow(arr_div)
    pl.figure()
    pl.hold(True)
    linreg_xy = ([], [])
    for minmax, (rprof, region) in rprofs.items():
        # minmax_flux = arr_div[minmax]
        pts, fluxes, avg_fluxes, avg_fluxes_errs, avg_dists, avg_dists_errs = \
                zip(*rprof)
        linreg_xy[0].extend(fluxes)
        linreg_xy[1].extend(avg_fluxes)
        # fluxes = np.abs(np.array(fluxes) - minmax_flux)
        # avg_fluxes = np.abs(np.array(avg_fluxes) - minmax_flux)
        # pl.plot(avg_dists, avg_fluxes, 'd-')
        pl.plot(avg_dists, avg_fluxes, 'd-')
    pl.grid()
    slope, intercept, rval, pval, stderr = stats.linregress(*linreg_xy)
    print
    print "slope: %f" % slope
    print "intercept: %f" % intercept
    print "rval: %f" % rval
    print "pval: %f" % pval
    print "stderr: %f" % stderr
    import pdb; pdb.set_trace()

def Wave2DShow(ufield, ds, vel=None, vmin=None, vmax=None):
    r"""
    Show a 2D pressure field at some instant of time.
    As background is shown velocity field.
    Same dimension as ufield.

    * ufield    : 2d pressure field at an instant of time
    * ds        : space discretization
    * vel       : 2d background velocity field
    * vmin/vmax : vmin/vmax of imshow
    """
    #max index time and max index space
    maxt = np.shape(snapshots)[0]
    maxk = np.shape(snapshots)[1]
    maxi = np.shape(snapshots)[2]    

    print "vmin : ", vmin, "vmax : ", vmax
    # space axis starting at 0 in x and z (using y coz' plotting)
    # extents of the picture,
    xmin, xmax = 0, ds*maxi
    ymin, ymax = 0, ds*maxk
    extent= xmin, xmax, ymax, ymin
    py.hold(True)
    if not vel == None:
        py.imshow(vel, interpolation='bilinear', cmap=cm.jet, extent=extent,  origin='upper', aspect='auto')

    py.imshow(ufield, interpolation='bilinear', cmap=cm.Greys_r, alpha=0.8, extent=extent, origin='upper', aspect='auto', vmin=vmin, vmax=vmax)
    py.hold(False)
    # optional cmap=cm.jet, apect='auto' adjust aspect to the previous plot
    py.show()

def plot_the_overview(samples, i, j,  output_image_file):

    pylab.hold(True)
    pylab.scatter(samples[:,i], samples[:,j])
    pylab.draw()
    pylab.savefig(output_image_file, dpi=150)
    pylab.close()

def pinwheel_overlay(pinwheels, contours=None, style='wo',linewidth=1,mmap=None):
   """
   Plots the pinwheel locations and optionally the real and imaginary
   pinwheel contours. Designed to be overlayed over an OR map.
   """
   fig = plt.figure(frameon=False)
   fig.patch.set_alpha(0.0)
   ax = plt.subplot(111, aspect='equal', frameon=True)
   ax.patch.set_alpha(0.0)
   plt.hold(True)
    
   plt.imshow(mmap,cmap='hsv',extent=(0, 1.0, 0, 1.0))
   (recontours, imcontours) = contours if contours else ([],[])
   for recontour in recontours:
      plt.plot(recontour[:,0], recontour[:,1],'k',linewidth=linewidth)
   for imcontour in imcontours:
      plt.plot(imcontour[:,0], imcontour[:,1],'w', linewidth=linewidth)

   Xs, Ys = zip(*pinwheels)
   plt.plot(np.array(Xs), np.array(Ys), style)

   plt.xlim((0.0,1.0));         plt.ylim((0.0,1.0))
   ax.xaxis.set_ticks([]);      ax.yaxis.set_ticks([])
   ax.xaxis.set_ticklabels([]); ax.yaxis.set_ticklabels([])
   return fig

	def plot_groups_at_time_point(self,t, feat1, feat2):


		markers = ['ro', 'go', 'bo', 'yo', 'ko', 'mo', 'co']

		cp_list = self.cell_tracker.list_of_cell_profiles_per_timestamp[t].list_of_cell_profiles



		fig = pylab.figure( facecolor='white')

		counter = -1

		for group_name in self.groups.keys():

			counter +=1
			gr = self.groups[group_name][t]

			feat1_vals = []
			feat2_vals = []
			
			for idx in gr:
				feat1_vals.append(cp_list[idx].dict_of_features[feat1])
				feat2_vals.append(cp_list[idx].dict_of_features[feat2])

			pylab.plot(feat1_vals, feat2_vals, markers[counter], label = group_name)

			pylab.hold(True)

		fig.canvas.set_window_title("Time point %s" % t)

		pylab.legend(loc="best")
		pylab.xlabel(feat1)
		pylab.ylabel(feat2)
		pylab.grid()

def PlotTorques(robot, traj0, traj1, dt=0.001, taumin=[], taumax=[],
                figstart=0):
    from pylab import figure, clf, hold, gca, plot, axis, title, xlabel, ylabel
    colorcycle = ['r', 'g', 'b', 'm', 'c', 'y', 'k']
    colorcycle = colorcycle[0:traj0.dimension]
    Tmax = max(traj0.duration, traj1.duration)
    tvect0, tauvect0 = ComputeTorques(traj0, robot, dt)
    tvect1, tauvect1 = ComputeTorques(traj1, robot, dt)
    figure(figstart)
    clf()
    hold('on')
    ax = gca()
    ax.set_color_cycle(colorcycle)
    plot(tvect0, tauvect0, '--', linewidth=2)
    ax.set_color_cycle(colorcycle)
    plot(tvect1, tauvect1, linewidth=2)
    ax.set_color_cycle(colorcycle)
    for a in taumax:
        plot([0, Tmax], [a, a], '-.')
    ax.set_color_cycle(colorcycle)
    for a in taumin:
        plot([0, Tmax], [a, a], '-.')
    if len(taumax) > 0:
        axis([0, Tmax, 1.2 * min(taumin), 1.2 * max(taumax)])
    title('Joint torques', fontsize=20)
    xlabel('Time (s)', fontsize=18)
    ylabel('Joint torques (Nm)', fontsize=18)

def anim_slices(dir="bckreact_b5_n01_E1_M14_t4d3",out=[2,3,4,5,6,7,8,9],end=[1.,1.,1.],N=None,var='d',log=True,fig=1,over=False, rmax=10.0):
#def anim_slices(dir="backreaction_b5_3d3yr",out=[3,6],end=[1.,1.,1.],N=None,var='d',log=True,fig=1,over=False, rmax=15.0, mp4="dens_maps.mp4"):
    # reading time intervals
    tseq = read_times(dir=dir+'/')

    pylab.close(fig)
    figure = matplotlib.pyplot.figure(num=fig)
    #a = matplotlib.pyplot.gca()
    a = None
    
    ims = []
    for i in out:
        print 'output #',i, '  age:',tseq[i-1]
        connect(dir=dir,out=i,var=[var])
        load_slices()
        [map,img, jet] = show_slice(var=var,log=False,rmax=rmax,ax=a,time=tseq[i-1])
        ann = pylab.annotate(tseq[i-1]+' yrs', xy=(.8, .9),  xycoords='axes fraction', horizontalalignment='center', verticalalignment='center', color="white") 
	ims.append((img,ann,))
	a = jet.ax
	#ims.append((img,ann,jet,))
	# artificially expand the longer steps
        #if eval(tseq[i-1])>(5.e2-1.0): 
            #ims.append((img,ann))
        #    if eval(tseq[i-1])>1.5e3-1.0:  ims.append((img,ann))
	fname = '_tmp%03d.png'%i
	print 'Saving frame', fname
	figure.savefig(basedir+dir+'/'+fname)
	pylab.hold(False)

    matplotlib.pyplot.xlabel('r, pc')
    matplotlib.pyplot.ylabel('r, pc')

    im_ani = anim.ArtistAnimation(figure, ims, interval=200, repeat_delay=500, blit=True)

def demo():
    from pylab import hold, linspace, subplot, plot, legend, show
    hold(True)
    #y = [9,6,1,3,8,4,2]
    #y = [9,11,13,3,-2,0,2]
    y = [9, 11, 2, 3, 8, 0]
    #y = [9,9,1,3,8,2,2]
    x = linspace(0, 1, len(y))
    t = linspace(x[0], x[-1], 400)
    subplot(211)
    plot(t, bspline(y, t, clamp=False), '-.y',
         label="unclamped bspline")  # bspline
    # bspline
    plot(t, bspline(y, t, clamp=True), '-y', label="clamped bspline")
    plot(sorted(x), y, ':oy', label="control points")
    legend()
    #left, right = _derivs(t, bspline(y, t, clamp=False))
    #print(left, (y[1] - y[0]) / (x[1] - x[0]))

    subplot(212)
    xt, yt = pbs(x, y, t, clamp=False)
    plot(xt, yt, '-.b', label="unclamped pbs")  # pbs
    xt, yt = pbs(x, y, t, clamp=True)
    plot(xt, yt, '-b', label="clamped pbs")  # pbs
    #xt,yt = pbs(x,y,t,clamp=True, parametric=True)
    # plot(xt,yt,'-g') # pbs
    plot(sorted(x), y, ':ob', label="control points")
    legend()
    show()

def BetavsT_AP_P(folder,keys):
  AP = Analysis.AnalyseFile()
  P = Analysis.AnalyseFile()
  I = 300e-6

  for f in folder:
    a=Analysis.AnalyseFile(f)
    fit= a.curve_fit(quad,'Current','Voltage',bounds=lambda x,y:x,result=True,header='Fit',asrow=True)
   
    if f['iterator'] == 7:
      AP.add_column(fit,str(f['iterator']))
      AP['Sample Temp'] = f['Sample Temp']
      Spc = ((-0.01411*f['Sample Temp'])-0.11185)*1e-6  
      T_d = -fit[0]*(I*I)/Spc
      AP['DeltaTemp'] = T_d
    elif f['iterator'] == 6:
      P.add_column(fit,str(f['iterator']))
      P['Sample Temp'] = f['Sample Temp']
      Spc = ((-0.01411*f['Sample Temp'])-0.11185)*1e-6  
      T_d = -fit[0]*(I*I)/Spc
      P['DeltaTemp'] = T_d
  
  
    
  plt.hold(True)
  plt.title(r'$\beta$ coef of NLIV vs Temp')
  plt.xlabel('Temperture (K)')
  plt.ylabel(r'$\beta$ (V/A$^2$)')
  plt.plot(f['IVtemp'],P['DeltaTemp']-AP['DeltaTemp'],'ok')

def plot_quad_mesh(pts, quads):
    """Plot a quadrilateral mesh."""

    def plot_quad(pts):
        """Plot one quad."""
        plot(np.r_[pts[:,0], pts[0,0]],
             np.r_[pts[:,1], pts[0,1]],
             lw=1.5, color="black")

    hold(True)
    for k,q in enumerate(quads):
        plot_quad(pts[q])

    if len(quads) < 400:
        for k,q in enumerate(quads):
            center = pts[q].sum(axis=0) / 4.0
            text(center[0], center[1], "%d" % k, color="black",
                 fontsize=10,
                 horizontalalignment='center', verticalalignment='center')

    if len(pts) < 200:
        for k,p in enumerate(pts):
            text(p[0], p[1], "%d" % k, color="black",
                 fontsize=10,
                 bbox=dict(boxstyle = "round", fc = "white"),
                 horizontalalignment='center', verticalalignment='center')