def draw_hist_cdf(data,fig=None,nbins=None,subpl=None,nohist=False,**figargs):
    '''input data is a list of dicts with keys: data,histcolor,plotline
    '''

    bins = None

    if fig: pylab.figure(fig,**figargs)
    if subpl: pylab.subplot(subpl)
    
    results = []
    for d in data:
        if bins is not None:
            n,bins,patches=pylab.hist(d['data'],bins=bins,normed=1,fc=histcolors[d['plotline'][0]])
        elif nbins is not None:
            n,bins,patches=pylab.hist(d['data'],bins=nbins,normed=1,fc=histcolors[d['plotline'][0]])
        else:
            n,bins,patches=pylab.hist(d['data'],normed=1,fc=histcolors[d['plotline'][0]])
        results.append(n)

    if nohist:
        pylab.cla()
    else:
        pylab.twinx()

    for i,d in enumerate(data):
        y = pylab.plot(bins,numpy.cumsum(results[i])/sum(results[i]), d['plotline'], linewidth=2)

    def plot_dH_ds(self, mean=0):
        # self.mu = mean
        x = np.linspace(1E-2,5,1000)
        def h(s):
            self.sigma = np.sqrt(s)
            self.sigma2 = s
            self.compute_Z()
            return self.H()
        def dh(s):
            self.sigma = np.sqrt(s)
            self.sigma2 = s
            self.compute_Z()
            return self.dH_dvar()
        def dmean_dvar(s):
            self.sigma = np.sqrt(s)
            self.sigma2 = s
            self.compute_Z()
            return self.dmean_dvar()
        H = np.vectorize(h)
        dH = np.vectorize(dh)
        # MV = np.vectorize(dmean_dvar)
        for m in range(1):
#             self.mu = m
            self.compute_Z()
            pb.figure("mu:{:.4f} var:{:.4f} side:{:s}".format(self.mu, self.sigma, self.side))
            pb.clf()
            pb.plot(x,dH(x), label='dH', lw=1.5)
            #pb.plot(x,MV(x), label='dMV', lw=1.5)
            pb.plot(x[:-1], np.diff(H(x))/np.diff(x), label='numdH', lw=1.5)
            pb.legend()
            pb.twinx()
            pb.plot(x,H(x), label='H')

def show_plot(title, true1, true2, data1, data2, x, plot):

    if (not plot): return

    pylab.clf()

    pylab.title(title)

    l1 = "-"
    l2 = "--"
    colt = "r"  # true data
    cold = "b"  # data
    colr = "g"  # residual

    pylab.plot(x, true1, label='true-1', color=colt, linestyle=l1)
    pylab.plot(x, true2, label='true-2', color=colt, linestyle=l2)

    pylab.plot(x, data1, label='data-1', color=cold, linestyle=l1)
    pylab.plot(x, data2, label='data-2', color=cold, linestyle=l2)
    pylab.legend(loc='upper left')

    pylab.twinx()
    pylab.plot(x, (data1-true1)/true1, label='data-1 % err', color=colr, 
               linestyle=l1)
    pylab.plot(x, (data2-true2)/true2, label='data-2 % err', color=colr, 
               linestyle=l2)
    pylab.legend(loc='upper right')

    pylab.show()

    return

def plot_comparison(results_dir):
	'''INTERFACE: Plot the results of comparing two accounts for common ads.

	Args:
		results_dir: Directory path containing the comparison results file and
		also where the comparison plot should be saved.
	'''
	fd = open(results_dir + "/results.txt", "r")
	bases = []
	counts = []
	others = []
	misses = []
	founds = []

	for line in fd.readlines():
		if "BaseTrial" in line:
			continue

		base, count, other, missed = line.strip().split()
		base = int(base) + 1
		count = int(count)
		other = int(other) + 1
		missed = int(missed)
		found = 1 - (missed/float(count))

		if len(bases) > 0 and bases[-1] == base:
			others[-1] = other
			misses[-1] = missed
			founds[-1] = found
		else:
			bases.append(base)
			counts.append(int(count))
			others.append(other)
			misses.append(missed)
			founds.append(found)

	fd.close()

	pylab.ylim([0, 1])
	pylab.yticks(adLib.float_range(0, 1, 0.1))
	pylab.xlabel("Base Trial")
	pylab.ylabel("Fraction of Ads Found")
	pylab.plot(bases, founds, "b.", label="Found " + \
										str(round(sum(founds)/len(founds), 3)))
	pylab.legend(loc="upper left", prop={'size':10})
	pylab.title("Common Ads in Identical Accounts")

	pylab.twinx()
	pylab.xlim([0, 100])
	pylab.xticks(range(0, 100, 10))
	pylab.ylim([0, 15])
	pylab.yticks(range(0, 15, 1))
	pylab.ylabel("Number of Trials to Find Base Ads")
	pylab.plot(bases, others, "r.", label="In Trials " + \
								str(round(sum(others)/float(len(others)), 3)))
	pylab.legend(loc="lower right", prop={'size':10})

	pylab.savefig(results_dir + "/results.png")
	pylab.clf()
	print results_dir, sum(counts)/float(len(counts)), sum(misses)/float(len(misses))

def plot_error(x,x_hat):
    import pylab
    
    pylab.figure()
    pylab.plot(x,x)
    pylab.plot(x,x_hat)
    pylab.twinx()
    pylab.plot(x,x_hat-x)
    pylab.title('RMSE: %f'%rms(x-x_hat))

def test():

    starttime = time.time()

    # set up theta, phi, r, quantities
    rlo = 6.2e10
    rhi = 6.96e10
    thi = 40.0*numpy.pi/180.
    tlo = -40.0*numpy.pi/180.
    phi = 40.0*numpy.pi/180.
    plo = -40.0*numpy.pi/180.

    dr = (rhi - rlo)/100
    dth = (thi - tlo)/100
    dphi = (phi - plo)/100

    radius = numpy.arange(rlo, rhi+dr, dr)
    theta  = numpy.arange(tlo, thi+dth, dth)
    phi    = numpy.arange(plo, phi+dphi, dphi)

    # data has dimensions of (nth, nphi, nr, nq=1)
    data = func(radius, theta, phi)

    # avgdata has dimensions of (nr, nq)
    avgdata = shell_avg(data, theta, phi, method='simp')

    # only consider first quantity
    avgdata = avgdata[:,0]

    true = exact(radius, phi, theta)

    l2norm = numpy.sqrt(dr*numpy.sum((true-avgdata)**2))

    print "\nL2 Norm: ",l2norm
    print

    endtime = time.time()
    print "elapsed time: ", endtime-starttime
    print

    pylab.clf()

    pylab.plot(radius, true, label='true', color='r', linestyle='--')
    pylab.plot(radius, avgdata, label='avg', color='b', linestyle='-')
    pylab.legend(loc='lower left')

    pylab.twinx()
    pylab.plot(radius, (avgdata - true)/true, label='% err', 
               color='g', linestyle='--')
    pylab.legend(loc='lower right')

    pylab.show()

    return

def plot_quality(btquality):
    """Make diagnostic plot of output of computeBTQuality()."""
    import pylab as pl

    btq = btquality
    pl.figure()
    pl.plot(btq.index, (btq.N - btq.N.ix[0]) / float(btq.N.ix[0]), "bo")
    pl.ylabel("$(N - N_0) / N_0$", color="b")
    pl.xlabel("Frame")
    pl.twinx()
    pl.plot(btq.index, -(btq.Nconserved - btq.N.ix[0]) / float(btq.N.ix[0]), "r^")
    pl.ylabel("Fraction dropped", color="r")

    def optimize_lbfgsb(self, hessian_terms=10, plotfn=None):

        XX = []
        OO = []
        def objective(x, tractor, stepsizes, lnp0):
            res = lnp0 - tractor(x * stepsizes)
            print 'LBFGSB objective:', res
            if plotfn:
                XX.append(x.copy())
                OO.append(res)
            return res

        from scipy.optimize import fmin_l_bfgs_b

        stepsizes = np.array(self.getStepSizes())
        p0 = np.array(self.getParams())
        lnp0 = self.getLogProb()

        print 'Active parameters:', len(p0)

        print 'Calling L-BFGS-B ...'
        X = fmin_l_bfgs_b(objective, p0 / stepsizes, fprime=None,
                          args=(self, stepsizes, lnp0),
                          approx_grad=True, bounds=None, m=hessian_terms,
                          epsilon=1e-8, iprint=0)
        p1,lnp1,d = X
        print d
        print 'lnp0:', lnp0
        self.setParams(p1 * stepsizes)
        print 'lnp1:', self.getLogProb()

        if plotfn:
            import pylab as plt
            plt.clf()
            XX = np.array(XX)
            OO = np.array(OO)
            print 'XX shape', XX.shape
            (N,D) = XX.shape
            for i in range(D):
                OO[np.abs(OO) < 1e-8] = 1e-8
                neg = (OO < 0)
                plt.semilogy(XX[neg,i], -OO[neg], 'bx', ms=12, mew=2)
                pos = np.logical_not(neg)
                plt.semilogy(XX[pos,i], OO[pos], 'rx', ms=12, mew=2)
                I = np.argsort(XX[:,i])
                plt.plot(XX[I,i], np.abs(OO[I]), 'k-', alpha=0.5)
                plt.ylabel('Objective value')
                plt.xlabel('Parameter')
            plt.twinx()
            for i in range(D):
                plt.plot(XX[:,i], np.arange(N), 'r-')
                plt.ylabel('L-BFGS-B iteration number')
            plt.savefig(plotfn)

 def fit(self, X, n_epochs = 100, plot=True):
     p = self.p
     q = self.q 
     w_ma = np.random.randn(p)
     w_ar = np.random.randn(q)
     k = max(p,q)
     n = X.shape[0]
     p_offset = k - p
     q_offset = k - q 
     Y = np.random.randn(n)
     ma_changes = [] 
     ar_changes = [] 
     errors = [] 
     learning_rate = self.learning_rate 
     for i in xrange(n_epochs):
         if self.verbose: print "Epoch ", i 
         old_ma = w_ma.copy()
         old_ar = w_ar.copy()
         for j in np.random.permutation(n - k):    
             curr_idx = j+k
             x_prev = X[j+p_offset : curr_idx]
             y_prev = Y[j+q_offset : curr_idx]
             pred = np.dot(x_prev, w_ma) + np.dot(y_prev, w_ar)
             Y[curr_idx] = pred 
             err = X[curr_idx] - pred
             w_ma += err * x_prev * learning_rate 
             w_ar += err * y_prev * learning_rate 
         ma_change = np.linalg.norm(old_ma - w_ma)
         ma_changes.append(ma_change)
         ar_change = np.linalg.norm(old_ar - w_ar)
         ar_changes.append(ar_change)
         mean_abs_error = np.mean(np.abs(Y[k:] - X[k:])) 
         errors.append(mean_abs_error)
         if self.verbose:
             print "MA weight change", ma_change
             print "AR weight change", ar_change 
             print "Mean abs. error:", mean_abs_error
         
         if ar_change < self.tol and ma_change < self.tol: break 
     if plot:
         import pylab
         pylab.plot(ma_changes)
         pylab.plot(ar_changes)
         pylab.legend(['MA', 'AR'])
         pylab.twinx()
         pylab.plot(errors, 'r')
     self.w_ma = w_ma
     self.w_ar = w_ar 

def plot_two_values(X, Y1, Y2, xlabel, y1label, y2label, suffix):
    output_filename = constants.CHARTS_FOLDER_NAME + constants.DATASET + '_' + suffix

    pylab.figure(figsize=(15, 7))

    pylab.rcParams.update({'font.size': 20})

    pylab.xlabel(xlabel)

    ax1 = pylab.gca()
    ax1.plot(X, Y1, 'b')
    ax1.plot(X, [0.0 for x in X], 'b--')
    ax1.set_ylabel(y1label, color='b')   
    for tl in ax1.get_yticklabels():
        tl.set_color('b')
    

    ax2 = pylab.twinx()
    ax2.plot(X, Y2, 'r')
    ax2.plot(X, [0.0 for x in X], 'r--')
    ax2.set_ylabel(y2label, color='r')   
    for tl in ax2.get_yticklabels():
        tl.set_color('r')

    pylab.savefig(output_filename + '.pdf')

    def plot(self, derivative=0, xmin="auto", xmax="auto", steps=500, smooth=0, simple='auto', clear=True, yaxis='left'):
        if simple=='auto': simple = self.simple

        # get the min and max
        if xmin=="auto": xmin = self.xmin
        if xmax=="auto": xmax = self.xmax

        # get and clear the figure and axes
        f = _pylab.gcf()
        if clear and yaxis=='left': f.clf()

        # setup the right-hand axis
        if yaxis=='right':  a = _pylab.twinx()
        else:               a = _pylab.gca()

        # define a new simple function to plot, then plot it
        def f(x): return self.evaluate(x, derivative, smooth, simple)
        _pylab_help.plot_function(f, xmin, xmax, steps, clear, axes=a)

        # label it
        th = "th"
        if derivative == 1: th = "st"
        if derivative == 2: th = "nd"
        if derivative == 3: th = "rd"
        if derivative: self.ylabel = str(derivative)+th+" derivative of "+self.ylabel+" spline"
        a.set_xlabel(self.xlabel)
        a.set_ylabel(self.ylabel)
        a.figure.canvas.Refresh()

 def postProcess(self):
     coreName = self.study.cacheDir + os.sep + "hz-" + self.name
     f = open(coreName)
     vals = {}
     for l in f:
         toks = l.split("\t")
         vals[toks[0], toks[1]] = toks[2:]
     f.close()
     for pop in self.pops:
         pylab.clf()
         pylab.title(pop)
         labels = []
         ehzs = []
         ohzs = []
         ns = []
         for indiv in self.study.pops.getIndivs(pop):
             labels.append("\n".join(list(indiv)))
             o, e, n, f = vals[indiv]
             ehz = 1 - float(e) / int(n)
             ohz = 1 - float(o) / int(n)
             ehzs.append(ehz)
             ohzs.append(ohz)
             ns.append(int(n))
         pylab.xticks(rotation=90, fontsize="x-small")
         pylab.legend()
         pylab.plot(ehzs, "+", label="ExpHe")
         pylab.plot(ohzs, "+", label="ObsHe")
         a2 = pylab.twinx()
         a2.plot(ns, ".")
         pylab.xticks(list(range(len(labels))), labels)
         xmin, xmax = pylab.xlim()
         pylab.xlim(xmin - 1, xmax + 1)
         pylab.savefig(coreName + "-" + pop + ".png")

def plot_props(xlab, props, magbins, delta, flags_to_use,plot_info):
    magbinctr = (magbins[1:]+magbins[0:-1])/2

    ax1 = pl.gca()
    pl.xlabel(xlab)

    ax2 = pl.twinx()

    for flag in flags_to_use:
        ax2.plot(magbinctr, props[flag], marker = 'o', ms = plot_info[flag]['ms'], ls = '-', color = plot_info[flag]['color'], label = plot_info[flag]['label'])


        ax2.set_ylabel('fraction of galaxies', fontsize=8)
        ax1.set_ylabel('Total galaxies', fontsize=8)

        ax1.yaxis.tick_right()
        ax1.yaxis.set_label_position("right")
        ax2.yaxis.tick_left()
        ax2.yaxis.set_label_position("left")

        ax2.set_ylim(0,1.0)

    ax1.bar(magbins[:-1], props['total'], width = delta, color = plot_info['total']['color'], log = False, zorder = -100) 

        
    #for tick in ax1.yaxis.get_major_ticks():
    #    tick.set_label('%2.1e' %float(tick.get_label())) 
            
        #ticklabs = ax1.get_yticklabels()
        #print ticklabs
        #print ['%2.1e' %float(x.get_text()) for x in ticklabs]
        #ax1.set_yticklabels( ['%2.1e' %float(x) for x in ticklabs] )
    return ax1, ax2

def plottao():
  f = pl.figure(2)
  pl.clf()
  pl.pcolor(Rar,Rmr,tao_e)
  #pl.contourf(Rar,Rmr,tao_e,10)
  cb = pl.colorbar()
  pl.contour(Rar,Rmr,tao_e,10,colors='r')
  pl.ylim(Rms[0],Rms[1])
  cb.set_label('tao_e at x = 0 [ms]')
  pl.xlabel('Ra [Ohm*cm]')
  pl.ylabel('Rm [Ohm*cm^2]')
 
  pl.twinx()
  #Show real tao equivalent values right
  pl.ylim(Rms[0]*1e-3,Rms[1]*1e-3)  
  pl.ylabel('tao [ms]')

def plot_props(xlab, props, magbins, delta, flags_to_use,plot_info,do_ci=False):
    magbinctr = (magbins[1:]+magbins[0:-1])/2

    ax1 = pl.gca()
    pl.xlabel(xlab)

    ax2 = pl.twinx()

    max_yticksl = 6
    max_yticksr = 5
    max_xticks = 5
    ylocl = plt.MaxNLocator(max_yticksl,prune='lower')
    ylocr = plt.MaxNLocator(max_yticksr,prune='lower')
    xloc = plt.MaxNLocator(max_xticks)

    for flag in flags_to_use:
        if do_ci:
            ax2.errorbar(magbinctr, np.array(props[flag]),
                yerr=np.array(props[str(flag)+'_err']).T,
                ecolor=plot_info[flag]['color'], elinewidth=1, 
                capsize=3, marker = plot_info[flag]['marker'],
                ms = plot_info[flag]['ms'], ls = plot_info[flag]['ls'], 
                color = plot_info[flag]['color'], 
                label = plot_info[flag]['label'])
        else:
            ax2.plot(magbinctr, np.array(props[flag]), 
                 marker = plot_info[flag]['marker'],
                 ms = plot_info[flag]['ms'], ls = plot_info[flag]['ls'], 
                 color = plot_info[flag]['color'], 
                 label = plot_info[flag]['label'])



        ax2.set_ylabel('fraction of galaxies', fontsize=8)
        ax1.set_ylabel('Total galaxies', fontsize=8)

        ax1.yaxis.tick_right()
        ax1.yaxis.set_label_position("right")
        ax2.yaxis.tick_left()
        ax2.yaxis.set_label_position("left")


    ax1.bar(magbins[:-1], props['total'], width = delta, color = plot_info['total']['color'], log = False, zorder = -100) 

    ax1.xaxis.set_major_locator(xloc)
    ax2.xaxis.set_major_locator(xloc)
    ax1.yaxis.set_major_locator(ylocr)
    ax2.yaxis.set_major_locator(ylocl)
    ax2.set_ylim(0,1.0)

    #for tick in ax1.yaxis.get_major_ticks():
    #    tick.set_label('%2.1e' %float(tick.get_label())) 
            
        #ticklabs = ax1.get_yticklabels()
        #print ticklabs
        #print ['%2.1e' %float(x.get_text()) for x in ticklabs]
        #ax1.set_yticklabels( ['%2.1e' %float(x) for x in ticklabs] )
    ax1.yaxis.set_tick_params(labelsize=6)
    ax2.yaxis.set_tick_params(labelsize=6)
    return ax1, ax2

def handler(signum, frame):
    try:
        fd.close()  # might not be open yet
    except:
        pass

    a=inspect.stack()
    stacklevel=0
    for k in range(len(a)):
        if (string.find(a[k][1], 'profileplot.py') > 0):
            stacklevel=k
            break
    myf=sys._getframe(stacklevel).f_globals
#    myf=frame.f_globals
    t=myf['t']
    y11=myf['y11']
    y22=myf['y22']
    numfile=myf['numfile']
    pl.plot(t,y11,lw=2)
    pl.plot(t,y22,lw=2)
    if max(y11)>=max(y22):
        pl.axis([0.9*min(t),1.1*max(t),0.9*min(y11),1.1*max(y11)])
    else:
        pl.axis([0.9*min(t),1.1*max(t),0.9*min(y22),1.1*max(y22)])
    pl.xlabel('time (sec)')
    pl.ylabel('memory footprint') #note virtual vs. resident
    font=FontProperties(size='small')
    pl.legend(('virtual','resident'),loc=[0.7,0.85], prop=font)
    ax2 = pl.twinx()
    pl.ylabel('No of open File Descriptors')
    ax2.yaxis.tick_right()
    pl.plot(t,numfile, 'r-.',lw=2)
    pl.legend(['No. of Open FDs'],loc=[0.7,0.8], prop=font)
    pl.title('memory usage of casapy for '+testname)

    #s="test-plot.ps" change to PNG for web browser compatibility
    #s="test-plot.png"

    #set up special images directory?
    if os.path.basename(myf['RESULT_DIR']) == myf['testname']:
        png_path=myf['RESULT_DIR']
    else:
        png_path=myf['RESULT_DIR']+time.strftime('/%Y_%m_%d')

    png_filename=png_path +'/'+myf['testname']+'_profile.png'
    if not os.path.isdir(png_path):
        os.mkdir(png_path)

    pl.savefig(png_filename);
    ht=htmlPub(myf['webpage'], 'Memory profile of '+myf['testname'])
    body1=['<pre>Memory profile of run of test %s at %s </pre>'%(myf['testname'],time.strftime('%Y/%m/%d/%H:%M:%S'))]
    body2=['']
    ht.doBlk(body1, body2,myf['testname']+'_profile.png', 'Profile')
    ht.doFooter()

    sys.stdout.flush()

    sys.exit()
    return

def plot_derivatives(x,y,subpl,xmax, sm=20):
    x,y = x[:-sm/2],Util.smooth(y,sm)[:-sm/2]
    x0,y0 = derivative(x,y,False)
    x0,y0 = x0[:-sm/2],Util.smooth(y0,sm)[:-sm/2]
    x1,y1 = derivative(x0,y0,False)
    x1,y1 = x1[:-sm/2],Util.smooth(y1,sm)[:-sm/2]

    pylab.subplot(*subpl)
    #pylab.cla()
    pylab.plot(x,y,'b',label='fn')
    pylab.twinx()
    pylab.plot(x0,y0,'g',label='d0')
    pylab.yticks([])
    pylab.twinx()
    pylab.plot(x1,y1,'r',label='d1')
    pylab.yticks([])
    pylab.plot([x1[0],x1[-1]],[0,0],'k:')
    pylab.xlim(0,xmax)

 def _plotResults(self, path):
     '''
     Power and Incoming
     '''        
     fig = pylab.figure()
     fields = getColumns(MONITOR_NAME)
     data = file2data(fields, MONITOR_NAME)
     pylab.grid(True)
     ax1 = pylab.subplot(111)        
     l1 = pylab.plot(data[4], color='red')
     ax2 = pylab.twinx()        
     l2 = pylab.plot(data[1], color='blue')
    
     ax2.legend((l1,l2), ('Incoming', 'Power'), 'best')
     ax1.set_ylabel("Incoming (req/s)")
     ax2.set_ylabel("Power (W)")
     ax1.set_xlabel("Time (s)")
     ax2.set_yticks(range(0, 600, 100))
     pylab.ylim(-300, 600)
     
     fig.savefig(path + 'power.eps')
     
     '''
     Sat and Incoming
     '''
     fig = pylab.figure()
     fields = getColumns(MONITOR_NAME)
     data = file2data(fields, MONITOR_NAME)
     pylab.grid(True)
     ax1 = pylab.subplot(111)        
     l1 = pylab.plot(data[4], color='blue')
     ax2 = pylab.twinx()        
     l2 = pylab.plot(data[10], color='red')
     qos_target = self._perfRegulator.controller._setPoint        
     pylab.axhline(qos_target, color='orange')        
     
     ax2.legend((l1,l2), ('Incoming', 'Saturation'), 'best')
     ax1.set_ylabel("Incoming (req/s)")
     ax2.set_ylabel("Saturation")
     ax1.set_xlabel("Time (s)")
     ax2.set_yticks([x/10.0 for x in range(0, 20, 2)])
     pylab.ylim(-2.5, 2)
     
     fig.savefig(path + 'sat.eps')

def flow_plot(ax1,expt_points,data_points, regime):
	plt.figure(j)
	ax1.plot(expt_points,data_points,'k-')
	#ax1.set_xlabel('z (microm)')
	ax1.set_ylabel(regime)
	ax2 = twinx()
	ax2.plot(expt_points,radius_points, 'r-')
	ax2.axvline(x=0, ymin=0.01, ymax=0.99, color='b',linestyle=':')
	ax2.axvline(x=(d/1e-6), ymin=0.01, ymax=0.99, color='b',linestyle=':')
	plt.show

def plot_results(results, ai, ylab, xlab="Time (h)"):
    nplot = len(results)

    for i, (result_mean, result_std, epsilon, lab, y2lab, y2max, y2min) in \
    enumerate(results):

        subplot(nplot, 1, i + 1)

        if i == nplot-1:
            xlabel(xlab)

        title(lab)

        x = arange(0.0, result_mean.shape[1], 1.0)
        y = result_mean[ai, :]
        e = result_std[ai, :]
        y2 = epsilon[ai, :]

#        plot(x, result_mean[ai, :],
#             color=clr[ai % nc],
#             linestyle=ls[ai % ns],
#             label=lab)

        errorbar(x, y, yerr=e, fmt='ko', linestyle="None",
                 label=ylab,
                 capsize=0, markersize=3)#, linewidth=0.2)
        ylabel(ylab)

        l = legend(loc="upper right")
        l.get_frame().set_linewidth(0.5)

        # Exploration rate plot.
        twinx()
        plot(x, y2, color="black", label=y2lab)
        ylabel(y2lab)
        if y2max is not None:
            ylim(ymax=y2max)
        if y2min is not None:
            ylim(ymin=y2min)

        l = legend(loc="lower right")
        l.get_frame().set_linewidth(0.5)