def demo():
    import pylab

    # The module normalize is not part of the osrefl code base.
    from reflectometry.reduction import normalize

    from .examples import ng7 as dataset
    spec = dataset.spec()[0]
    water = WaterIntensity(D2O=20,probe=spec.probe)
    spec.apply(normalize())
    theory = water.model(spec.Qz,spec.detector.wavelength)

    pylab.subplot(211)
    pylab.title('Data normalized to water scattering (%g%% D2O)'%water.D2O)
    pylab.xlabel('Qz (inv Ang)')
    pylab.ylabel('Reflectivity')
    pylab.semilogy(spec.Qz,theory,'-',label='expected')
    scale = theory[0]/spec.R[0]
    pylab.errorbar(spec.Qz,scale*spec.R,scale*spec.dR,fmt='.',label='measured')

    spec.apply(water)
    pylab.subplot(212)
    #pylab.title('Intensity correction factor')
    pylab.xlabel('Slit 1 opening (mm)')
    pylab.ylabel('Incident intensity')
    pylab.yscale('log')
    pylab.errorbar(spec.slit1.x,spec.R,spec.dR,fmt='.',label='correction')

    pylab.show()

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 plot_sed(fluxes, backgrounds, errors, **kwargs):
    """
    Trivial SED plotting
    """
    pl.errorbar(band_waves.values(),fluxes-backgrounds,yerr=errors,marker='s', **kwargs)
    pl.xlabel('$\lambda$ (mm)')
    pl.ylabel('mJy/beam')

def _show_rates(rate, wo, wt, attenuator, tau_NP, tau_P):
    import pylab

    #pylab.figure()
    pylab.errorbar(rate, wt[0], yerr=wt[1], fmt='g.', label='attenuated')
    pylab.errorbar(rate, wo[0], yerr=wo[1], fmt='b.', label='unattenuated')

    pylab.xscale('log')
    pylab.yscale('log')
    pylab.xlabel('incident rate (counts/second)')
    pylab.ylabel('observed rate (counts/second)')
    pylab.legend(loc='best')
    pylab.grid(True)
    pylab.plot(rate, rate/attenuator, 'g-', label='target')
    pylab.plot(rate, rate, 'b-', label='target')

    Ipeak, Rpeak = peak_rate(tau_NP=tau_NP, tau_P=tau_P)
    if rate[0] <= Ipeak <= rate[-1]:
        pylab.axvline(x=Ipeak, ls='--', c='b')
        pylab.text(x=Ipeak, y=0.05, s=' %g'%Ipeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_xaxis_transform())
    if False:
        pylab.axhline(y=Rpeak, ls='--', c='b')
        pylab.text(y=Rpeak, x=0.05, s=' %g\n'%Rpeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_yaxis_transform())

	def plot(self, params, errors=None,label=''):
		params=[max(1e-100,p) for p in params]
		E=np.concatenate(([self._ERange[0]],self._splitE,[self._ERange[1]]))
		pl.plot(reduce(lambda a,b:a+b,[[e,e] for e in E]),[1e-10]+reduce(lambda a,b:a+b,[[p,p] for p in params])+[1e-10],label=label)
		if errors!=None:
			for i in range(len(E)-1):
				pl.errorbar([np.sqrt(E[i]*E[i+1])],[params[i]],yerr=[errors[i]],fmt='r')

def show_table(table_name,ls="none", fmt="o", legend=False, name="m", do_half=0):
	bt = fi.FITS(table_name)[1].read()
	rgpp = (np.unique(bt["rgp_lower"])+np.unique(bt["rgp_upper"]))/2
	nbins = rgpp.size

	plt.xscale("log")
	colours=["purple", "forestgreen", "steelblue", "pink", "darkred", "midnightblue", "gray", "sienna", "olive", "darkviolet"]
	pts = ["o", "D", "x", "^", ">", "<", "1", "s", "*", "+", "."]
	for i,r in enumerate(rgpp):
		sel = (bt["i"]==i)
		snr = 10** ((np.log10(bt["snr_lower"][sel]) + np.log10(bt["snr_upper"][sel]))/2)

		if do_half==1 and i>nbins/2:
			continue
		elif do_half==2 and i<nbins/2:
			continue
		if legend:
			plt.errorbar(snr, bt["%s"%name][i*snr.size:(i*snr.size)+snr.size], bt["err_%s"%name][i*snr.size:(i*snr.size)+snr.size], color=colours[i], ls=ls, fmt=pts[i], lw=2.5, label="$R_{gpp}/R_p = %1.2f-%1.2f$"%(np.unique(bt["rgp_lower"])[i],np.unique(bt["rgp_upper"])[i]))
		else:
			plt.errorbar(snr, bt["%s"%name][i*snr.size:(i*snr.size)+snr.size], bt["err_%s"%name][i*snr.size:(i*snr.size)+snr.size], color=colours[i], ls=ls, fmt=pts[i], lw=2.5)

	plt.xlim(10,300)
	plt.axhline(0, lw=2, color="k")
	
	plt.xlabel("Signal-to-Noise $SNR_w$")
	if name=="m":
		plt.ylim(-0.85,0.05)
		plt.ylabel("Multiplicative Bias $m \equiv (m_1 + m_2)/2$")
	elif name=="alpha":
		plt.ylabel(r"PSF Leakage $\alpha \equiv (\alpha _1 + \alpha _2)/2$")
		plt.ylim(-0.5,2)



	plt.legend(loc="lower right")

def plot_data(yRange=None):
    '''
    Plots and saves the cell measurement data.  Returns nothing.
    '''
    fig = plt.figure(figsize=(18,12))
    ax = plt.subplot(111)
    plt.errorbar(range(len(avgCells.index)), avgCells[column], yerr=stdCells[column], fmt='o')
    ax = plt.gca()
    ax.set(xticks=range(len(avgCells.index)), xticklabels=avgCells.index)
    xlims = ax.get_xlim()
    ax.set_xlim([lim-1 for lim in xlims])
    # adjust yRange if it was specified
    if yRange!=None:
        ax.set_ylim(yRange)
        fileName = column + ' exlcuding outliers'
    else:
        fileName = column
    plt.subplots_adjust(bottom=0.2, right=0.98, left=0.05)
    plt.title(column)
    plt.ylabel('mm')
    locs, labels = plt.xticks()
    plt.setp(labels, rotation=90)
    mng = plt.get_current_fig_manager()
    mng.window.state('zoomed')
    #plt.show()
    path1 = 'Y:/Test data/ACT02/vision inspection/plot_100_cells/'
    path2 = 'Y:/Nate/git/nuvosun-python-lib/vision system/plot_100_cells/'
    fig.savefig(path1 + fileName, bbox_inches = 'tight')
    fig.savefig(path2 + fileName, bbox_inches = 'tight')
    plt.close()

def p2dscatter(self, log=False, color=None, label=None, orientation='horizontal', **kwargs):
    """ use pylab.errorplot to visualize these scatter points
        
        Parameters:
          log        : if true create logartihmic plot

          (all other kwargs will be passed to pylab.errobar)
    """
    if len(self.x) == 0:
        return

    ax = p.gca()
    if color is None:
        color = next(ax._get_lines.color_cycle)
    
    kw = {"xerr" : self.xerr, "yerr" : self.yerr, "fmt" : "k", "capsize" : 0., "linestyle" : 'None', "color" : color}
    kw.update(kwargs)
    
    if orientation == 'vertical':
        x, y = self.y, self.x
        kw["xerr"], kw["yerr"] = kw["yerr"], kw["xerr"]
        axis_name = 'x'
    else:
        x, y = self.x, self.y
        axis_name = 'y'
    
    _set_logscale(ax, log, axis=axis_name)
    
    p.errorbar(x, y, **kw) 
        
    if not hasattr(ax, "_legend_proxy"):
        ax._legend_proxy = LegendProxy(ax)
    ax._legend_proxy.add_scatter(label=label, color=color)
    
    _h2label(self, orientation)

    def _plot_aggr_random(self, span, Nmax, marker='o', color='r', markersize=6):
        # those are the best submitter. Nothing to recompute, can be extracted
        # from the df itself.
        iauc = [self.df.ix[x].mean_auc for x in range(0, Nmax)]

        pylab.clf()
        pylab.plot([x for x in span], iauc, marker+color, markersize=markersize,
                   label="AUC (individual submissions)".format(self.mode))
        pylab.grid(True)
        #pylab.plot()
        pylab.xlabel("N", fontsize=20)
        pylab.ylabel("AUROC", fontsize=20)
        pylab.title("Aggregated AUROC (random case)", fontsize=20)

        pylab.errorbar(span, self.results.mean(axis=0), self.results.std(axis=0),
                       label="{} aggregation (over N submissions)".format(self.mode))
        pylab.legend(loc="lower left")

        self._random_results = {}
        self._random_results['x'] = span
        self._random_results['individual'] = iauc
        self._random_results['aggregation_mean'] = list(self.results.mean(axis=0))
        self._random_results['aggregation_std'] = list(self.results.std(axis=0))
        self._random_results['aggregation_all'] = [list(x) for x in self.results]

        xmax = pylab.xlim()[1]
        pylab.ylim([0.35, 0.86])
        pylab.xlim(0.5, xmax)

def plot_dmsq2(HOpions, OOpions, title=None, save=False, name=''):
    "Plot m^2_{vs} - m^2_{vv}/2."
    
    # Set up figure.
    fig = p.figure()
    p.rc('text', usetex=True)
    p.rc('font', size=16)
    p.rc('axes', linewidth=0.5)
    p.xlabel('$am_{s}$')
    p.ylabel('$m^2_{vs} - m^2_{vv}/2$')
    legend = ()

    xr = np.linspace(0.0,0.06)

    # First data set.
    hopions = [HOpions[1], HOpions[5]]
    r = OOpions[3]
    xs = [q.m1 for q in hopions]
    ys = [(q.msq - r.msq/2) for q in hopions]
    es = [nerror(q.sig_msq, r.sig_msq) for q in hopions]
    fit = line_fit2(zip(xs,ys,es))
    legend += p.errorbar(xs, ys, fmt='bo')[0],
    # Fit results
    p.errorbar(xr, fit.a+fit.b*xr, fmt='b-')
    print fit.a, fit.sig_a

    if save:
        p.savefig(name)
    else:
        p.show()

    def plot_results(self, results, xloc, color, ls, label):
        iter_counts = sorted(set([it for it, av in results.keys() if av == self.average]))
        sorted_results = [results[it, self.average] for it in iter_counts]

        avg = np.array([r.train_logprob() for r in sorted_results])
        if hasattr(r, 'train_logprob_interval'):
            lower = np.array([r.train_logprob_interval()[0] for r in sorted_results])
            upper = np.array([r.train_logprob_interval()[1] for r in sorted_results])

        if self.logscale:
            plot_cmd = pylab.semilogx
        else:
            plot_cmd = pylab.plot

        xloc = xloc[:len(avg)]

        lw = 2.

        if label not in self.labels:
            plot_cmd(xloc, avg, color=color, ls=ls, lw=lw, label=label)
        else:
            plot_cmd(xloc, avg, color=color, ls=ls, lw=lw)

        self.labels.add(label)

        pylab.xticks(fontsize='xx-large')
        pylab.yticks(fontsize='xx-large')

        try:
            pylab.errorbar(xloc, (lower+upper)/2., yerr=(upper-lower)/2., fmt='', ls='None', ecolor=color)
        except:
            pass

 def __primativePlotTGNs__(self,bare=bool(False)):
   """
   Is a macro of plotting commands that takes a list of TGNs that
   plots each of these individually as a collection of points.
   Creates a figure plotting the thread of list of TGNs using the
   centroid and an X,Y error bars.  Take a optional boolean to
   make the plot not include a title and legend
   """
   #Determine the index that corresponds to X and Y quantities
   xIndex=0
   yIndex=0
   xLabel="NULL"
   yLabel="NULL"
   (xIndex,xLabel,yIndex,yLabel)=self.__getIndexAndLabels__()
   plotValues=list()
   gpsTimesInList=list()
   for thisTGN in self.tgnList:
       label=str(thisTGN.getID())
       #Get the X,Y property
       (xC,xE)=thisTGN.getCentroidErrorViaIndex(xIndex)
       (yC,yE)=thisTGN.getCentroidErrorViaIndex(yIndex)
       plotValues.append([xC,yC,xE,yE,label])
       gpsTimesInList.append(thisTGN.getGPS())
   for x,y,ex,ey,txtLabel in plotValues:
       pylab.errorbar(x,y,xerr=ex,yerr=ey,label=txtLabel,marker='o')
   pylab.xlabel(str(xLabel))
   pylab.ylabel(str(yLabel))
   if not bare:
     pylab.title("TGNs: %i"%(min(gpsTimesInList)))
     pylab.legend()

    def plot_fitness(self, show=True, save=False):
        df = pd.DataFrame([res['t1opt']['results']['Best_score'].values
            for res in self.results.allRes])

        df = df.astype(float)

        pylab.clf()
        for res in self.results.allRes:
            pylab.plot(res['t1opt']['results']['Best_score'], '--', color='grey')
        pylab.grid()
        pylab.xlabel("Generation")
        pylab.ylabel("Score")
        #pylab.plot(df.mean().values, 'kx--', lw=3, label='Mean Score')

        y = df.mean().values
        x = range(0, len(y))
        yerr = df.std().values
        pylab.errorbar(x, y, yerr=yerr, xerr=None, fmt='-', label='Mean Score',
                color='k', lw=3)
        pylab.legend()

        if save is True:
            self._report.savefig("fitness.png")

        if show is False:
            pylab.close()

def nishiyama09(wavelength, AKs, makePlot=False):
    # Data pulled from Nishiyama et al. 2009, Table 1

    filters = ['V', 'J', 'H', 'Ks', '[3.6]', '[4.5]', '[5.8]', '[8.0]']
    wave =      np.array([0.551, 1.25, 1.63, 2.14, 3.545, 4.442, 5.675, 7.760])
    A_AKs =     np.array([16.13, 3.02, 1.73, 1.00, 0.500, 0.390, 0.360, 0.430])
    A_AKs_err = np.array([0.04,  0.04, 0.03, 0.00, 0.010, 0.010, 0.010, 0.010])

    # Interpolate over the curve
    spline_interp = interpolate.splrep(wave, A_AKs, k=3, s=0)

    A_AKs_at_wave = interpolate.splev(wavelength, spline_interp)
    A_at_wave = AKs * A_AKs_at_wave

    if makePlot:
        py.clf()
        py.errorbar(wave, A_AKs, yerr=A_AKs_err, fmt='bo', 
                    markerfacecolor='none', markeredgecolor='blue',
                    markeredgewidth=2)
        
        # Make an interpolated curve.
        wavePlot = np.arange(wave.min(), wave.max(), 0.1)
        extPlot = interpolate.splev(wavePlot, spline_interp)
        py.loglog(wavePlot, extPlot, 'k-')

        # Plot a marker for the computed value.
        py.plot(wavelength, A_AKs_at_wave, 'rs',
                markerfacecolor='none', markeredgecolor='red',
                markeredgewidth=2)
        py.xlabel('Wavelength (microns)')
        py.ylabel('Extinction (magnitudes)')
        py.title('Nishiyama et al. 2009')

    
    return A_at_wave

def scatter_stats(db, s1, s2, f1=None, f2=None, **kwargs):
    if f1 == None:
        f1 = lambda x: x  # constant function

    if f2 == None:
        f2 = f1

    x = []
    xerr = []

    y = []
    yerr = []

    for k in db:
        x_k = [f1(x_ki) for x_ki in db[k].__getattribute__(s1).gettrace()]
        y_k = [f2(y_ki) for y_ki in db[k].__getattribute__(s2).gettrace()]

        x.append(pl.mean(x_k))
        xerr.append(pl.std(x_k))

        y.append(pl.mean(y_k))
        yerr.append(pl.std(y_k))

        pl.text(x[-1], y[-1], " %s" % k, fontsize=8, alpha=0.4, zorder=-1)

    default_args = {"fmt": "o", "ms": 10}
    default_args.update(kwargs)
    pl.errorbar(x, y, xerr=xerr, yerr=yerr, **default_args)
    pl.xlabel(s1)
    pl.ylabel(s2)

def plot_one_ppc(model, t):
    """ plot data and posterior predictive check
    
    :Parameters:
      - `model` : data.ModelData
      - `t` : str, data type of 'i', 'r', 'f', 'p', 'rr', 'm', 'X', 'pf', 'csmr'
    
    """
    stats = model.vars[t]['p_pred'].stats()
    if stats == None:
        return

    pl.figure()
    pl.title(t)

    x = model.vars[t]['p_obs'].value.__array__()
    y = x - stats['quantiles'][50]
    yerr = [stats['quantiles'][50] - pl.atleast_2d(stats['95% HPD interval'])[:,0],
            pl.atleast_2d(stats['95% HPD interval'])[:,1] - stats['quantiles'][50]]
    pl.errorbar(x, y, yerr=yerr, fmt='ko', mec='w', capsize=0,
                label='Obs vs Residual (Obs - Pred)')

    pl.xlabel('Observation')
    pl.ylabel('Residual (observation-prediction)')

    pl.grid()
    l,r,b,t = pl.axis()
    pl.hlines([0], l, r)
    pl.axis([l, r, y.min()*1.1 - y.max()*.1, -y.min()*.1 + y.max()*1.1])

def compare_models(db, stoch="itn coverage", stat_func=None, plot_type="", **kwargs):
    if stat_func == None:
        stat_func = lambda x: x

    X = {}
    for k in sorted(db.keys()):
        c = k.split("_")[2]
        X[c] = []

    for k in sorted(db.keys()):
        c = k.split("_")[2]
        X[c].append([stat_func(x_ki) for x_ki in db[k].__getattribute__(stoch).gettrace()])

    x = pl.array([pl.mean(xc[0]) for xc in X.values()])
    xerr = pl.array([pl.std(xc[0]) for xc in X.values()])
    y = pl.array([pl.mean(xc[1]) for xc in X.values()])
    yerr = pl.array([pl.std(xc[1]) for xc in X.values()])

    if plot_type == "scatter":
        default_args = {"fmt": "o", "ms": 10}
        default_args.update(kwargs)
        for c in X.keys():
            pl.text(pl.mean(X[c][0]), pl.mean(X[c][1]), " %s" % c, fontsize=8, alpha=0.4, zorder=-1)
        pl.errorbar(x, y, xerr=xerr, yerr=yerr, **default_args)
        pl.xlabel("First Model")
        pl.ylabel("Second Model")
        pl.plot([0, 1], [0, 1], alpha=0.5, linestyle="--", color="k", linewidth=2)

    elif plot_type == "rel_diff":
        d1 = sorted(100 * (x - y) / x)
        d2 = sorted(100 * (xerr - yerr) / xerr)
        pl.subplot(2, 1, 1)
        pl.title("Percent Model 2 deviates from Model 1")

        pl.plot(d1, "o")
        pl.xlabel("Countries sorted by deviation in mean")
        pl.ylabel("deviation in mean (%)")

        pl.subplot(2, 1, 2)
        pl.plot(d2, "o")
        pl.xlabel("Countries sorted by deviation in std err")
        pl.ylabel("deviation in std err (%)")
    elif plot_type == "abs_diff":
        d1 = sorted(x - y)
        d2 = sorted(xerr - yerr)
        pl.subplot(2, 1, 1)
        pl.title("Percent Model 2 deviates from Model 1")

        pl.plot(d1, "o")
        pl.xlabel("Countries sorted by deviation in mean")
        pl.ylabel("deviation in mean")

        pl.subplot(2, 1, 2)
        pl.plot(d2, "o")
        pl.xlabel("Countries sorted by deviation in std err")
        pl.ylabel("deviation in std err")
    else:
        assert 0, "plot_type must be abs_diff, rel_diff, or scatter"

    return pl.array([x, y, xerr, yerr])

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 plotCorr(self,pars=None,SHOW=True,SAVE=True): 
   skipfits=True
   pb.clf()
   for irn in range(len(self.ranges)):
     for ist in range(len(self.params)/3):
       pb.errorbar(np.arange(len(self.Pcorr[irn,ist,:]))/self.Fs,self.Pcorr[irn,ist,:],yerr=self.PcorrERR[irn,ist,:],color=colours[ist])
     pb.ylim([0,1])
     pb.grid(True)
     pb.xlabel('Time (s)',fontsize=20)
     pb.ylabel('Probabilities',fontsize=20)
     pb.xticks(fontsize=16)
     pb.yticks(fontsize=16)
     if SAVE:
       fn=os.path.basename(self.filename)
       pb.savefig(self.figdir+'corr_'+'range_'+str(irn)+fn[:-4]+'.eps')
       f=open(self.datdir+fn[:-4]+'.dat','a')
       f.seek(0,2)
       f.write('#corr_range_'+str(irn)+'\n')
       #f.write('BinCentre(pA)\tBinMin(pA)\tBinMax(pA)\tCounts\n')
       #for i in range(len(self.n)):
       #  f.write(str(self.x[i])+'\t'+str(self.bins[i])+'\t'+str(self.bins[i+1])+'\t'+str(self.n[i])+'\n')
       f.close()
     if SHOW:pb.show()
     pb.clf()
   return

def blocks_per_trial(experiment, neutral_blocks = False, clf = True, fig_no = 1, last_n = 6):
    days = set([s.day for s in experiment.get_sessions('all', 'all')])  
    residual_trials = np.zeros(len(experiment.subject_IDs)) # Number of trials in last (uncompleted) block of session.
    mean_bpt, sd_bpt = ([],[]) # Lists to hold mean and standard deviation of blocks per trial for each day.
    for day in days:
        day_blocks_per_trial = []
        sessions = experiment.get_sessions('all', day)
        for session in sessions:
            assert hasattr(session,'blocks'), 'Session does not have block info.'
            ax = experiment.subject_IDs.index(session.subject_ID) # Index used for residual trials array.
            blocks_per_trial, residual_trials[ax] = _session_blocks_per_trial(session, residual_trials[ax], neutral_blocks)
            day_blocks_per_trial.append(blocks_per_trial)
        mean_bpt.append(ut.nanmean(day_blocks_per_trial))
        sd_bpt.append  (np.sqrt(ut.nanvar(np.array(day_blocks_per_trial))))
    days = np.array(list(days))-min(days) + 1
    p.figure(fig_no)
    if clf: p.clf()
    p.subplot(2,1,1)
    p.errorbar(days, mean_bpt, sd_bpt/np.sqrt(len(experiment.subject_IDs)))
    p.xlim(0.5, max(days) + 0.5)
    p.ylim(ymin = 0)
    p.ylabel('Blocks per trial')
    p.subplot(2,1,2)
    p.plot(days,1 / np.array(mean_bpt))
    p.xlabel('Day')
    p.ylabel('Trials per block')