def plot_elecs_and_neurons(neuron_dict, ext_sim_dict, neural_sim_dict):
    pl.close('all')
    fig_all = pl.figure(figsize=[15,15])
    ax_all = fig_all.add_axes([0.1, 0.1, 0.8, 0.8], frameon=False)
    for elec in xrange(len(ext_sim_dict['elec_z'])):
        ax_all.plot(ext_sim_dict['elec_z'][elec], ext_sim_dict['elec_y'][elec], color='b',\
                marker='$E%i$'%elec, markersize=20 )    
    legends = []
    for i, neur in enumerate(neuron_dict):
        folder = os.path.join(neural_sim_dict['output_folder'], neuron_dict[neur]['name'])
        coor = np.load(os.path.join(folder,'coor.npy'))
        x,y,z = coor
        n_compartments = len(x)
        fig = pl.figure(figsize=[10, 10])
        ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], frameon=False)
        # Plot the electrodes
        for elec in xrange(len(ext_sim_dict['elec_z'])):
            ax.plot(ext_sim_dict['elec_z'][elec], ext_sim_dict['elec_y'][elec], color='b',\
                   marker='$%i$'%elec, markersize=20 )
        # Plot the neuron
        xmid, ymid, zmid = np.load(folder + '/coor.npy')
        xstart, ystart,zstart = np.load(folder + '/coor_start.npy')
        xend, yend, zend = np.load(folder + '/coor_end.npy')
        diam = np.load(folder + '/diam.npy')
        length = np.load(folder + '/length.npy')
        n_compartments = len(diam)
        for comp in xrange(n_compartments):
            if comp == 0:
                xcoords = pl.array([xmid[comp]])
                ycoords = pl.array([ymid[comp]])
                zcoords = pl.array([zmid[comp]])
                diams = pl.array([diam[comp]])    
            else:
                if zmid[comp] < 0.400 and zmid[comp] > -.400:  
                    xcoords = pl.r_[xcoords, pl.linspace(xstart[comp],
                                                         xend[comp], length[comp]*3*1000)]   
                    ycoords = pl.r_[ycoords, pl.linspace(ystart[comp],
                                                         yend[comp], length[comp]*3*1000)]   
                    zcoords = pl.r_[zcoords, pl.linspace(zstart[comp],
                                                         zend[comp], length[comp]*3*1000)]   
                    diams = pl.r_[diams, pl.linspace(diam[comp], diam[comp],
                                                length[comp]*3*1000)]
        argsort = pl.argsort(-xcoords)
        ax.scatter(zcoords[argsort], ycoords[argsort], s=20*(diams[argsort]*1000)**2,
                       c=xcoords[argsort], edgecolors='none', cmap='gray')
        ax_all.plot(zmid[0], ymid[0], marker='$%i$'%i, markersize=20, label='%i: %s' %(i, neur))
        #legends.append('%i: %s' %(i, neur))
        ax.axis(ext_sim_dict['plot_range'])
        ax.axis('equal')
        ax.axis(ext_sim_dict['plot_range'])
        ax.set_xlabel('z [mm]')
        ax.set_ylabel('y [mm]')
        fig.savefig(os.path.join(neural_sim_dict['output_folder'],\
                  'neuron_figs', '%s.png' % neur))
    ax_all.axis('equal')
    ax.axis(ext_sim_dict['plot_range'])
    ax_all.set_xlabel('z [mm]')
    ax_all.set_ylabel('y [mm]')
    ax_all.legend()
    fig_all.savefig(os.path.join(neural_sim_dict['output_folder'], 'fig.png'))

 def testTelescope(self):
     import matplotlib
     matplotlib.use('AGG')
     import matplotlib.mlab as ml
     import pylab as pl
     import time        
     w0 = 8.0
     k = 2*np.pi/3.0
     gb = GaussianBeam(w0, k)
     lens = ThinLens(150, 150)
     gb2 = lens*gb
     self.assertAlmostEqual(gb2._z0, gb._z0 + 2*150.0)
     lens2 = ThinLens(300, 600)
     gb3 = lens2*gb2
     self.assertAlmostEqual(gb3._z0, gb2._z0 + 2*300.0)
     self.assertAlmostEqual(gb._w0, gb3._w0/2.0)
     z = np.arange(0, 150)
     z2 = np.arange(150, 600)
     z3 = np.arange(600, 900)
     pl.plot(z, gb.w(z, k), z2, gb2.w(z2, k), z3, gb3.w(z3, k))
     pl.grid()
     pl.xlabel('z')
     pl.ylabel('w')
     pl.savefig('testTelescope1.png')
     time.sleep(0.1)
     pl.close('all')        

def plotEventTime(library, num, eventNames, sizes, times, events, filename = None):
  from pylab import close, legend, plot, savefig, show, title, xlabel, ylabel
  import numpy as np

  close()
  arches = sizes.keys()
  bs     = events[arches[0]].keys()[0]
  data   = []
  names  = []
  for event, color in zip(eventNames, ['b', 'g', 'r', 'y']):
    for arch, style in zip(arches, ['-', ':']):
      if event in events[arch][bs]:
        names.append(arch+'-'+str(bs)+' '+event)
        data.append(sizes[arch][bs])
        data.append(np.array(events[arch][bs][event])[:,0])
        data.append(color+style)
      else:
        print 'Could not find %s in %s-%d events' % (event, arch, bs)
  print data
  plot(*data)
  title('Performance on '+library+' Example '+str(num))
  xlabel('Number of Dof')
  ylabel('Time (s)')
  legend(names, 'upper left', shadow = True)
  if filename is None:
    show()
  else:
    savefig(filename)
  return

    def plot(self, filesuffix=('.png',)):

        pylab.figure()

        kPluses = 10**numpy.linspace(0, 3, 100)
        kMinuses = 10**numpy.linspace(6, 9, 100)
        figureOfMerits = numpy.zeros((len(kPluses), len(kMinuses), 4), 'd')
        for i, kPlus in enumerate(kPluses):
            for j, kMinus in enumerate(kMinuses):
                figureOfMerits[i, j, :] = self.figureOfMerit(self.generateData({'kPlus' : kPlus, 'kMinus' : kMinus}))
        data = self.generateData({'kPlus' : kPluses[0], 'kMinus' : kMinuses[0]})

        self._contourf(kMinuses, kPluses, figureOfMerits, data)

        pylab.xticks((10**6, 10**7, 10**8, 10**9), fontsize=self.fontsize)
        pylab.yticks((10**0, 10**1, 10**2, 10**3), fontsize=self.fontsize)
        pylab.xlabel(r'$k^-$ $\left(1\per\metre\right)$', fontsize=self.fontsize)
        pylab.ylabel(r'$k^+$ $\left(\power{\metre}{3}\per\mole\cdot\second\right)$', fontsize=self.fontsize)

        pylab.text(2 * 10**6, 7 * 10**2, r'I', fontsize=self.fontsize)
        pylab.text(3 * 10**7, 7 * 10**2, r'II', fontsize=self.fontsize)
        pylab.text(6 * 10**8, 7 * 10**2, r'III', fontsize=self.fontsize)
        pylab.text(6 * 10**8, 7 * 10**1, r'IV', fontsize=self.fontsize)

        for fP, kPlus, paxeslabel in ((1.143,3.51e+00, False), (0.975, 9.33e+00, False), (0.916, 3.51e+01, False), (0.89, 9.33e+01, False), (0.87, 3.51e+02, True)):
            for fM, kMinus, maxeslabel in ((1.4, 2.48e+06, False), (1.07, 7.05e+6, False), (0.96, 2.48e+07, False), (0.91, 7.05e+7, False), (0.88, 2.48e+08, True)):
                xpos = (numpy.log10(kMinus) - 6.) / 3. *  fM
                ypos = numpy.log10(kPlus) / 3. * fP        
                self.makeBackGroundPlot({'kPlus' : kPlus, 'kMinus' : kMinus}, xpos, ypos, axeslabel=paxeslabel and maxeslabel)

        for fs in filesuffix:
            pylab.savefig('kPlusVkMinus' + fs)
            pylab.close('all')

    def log_posterior(self,theta):

        model_g1 , model_g2, limit_mask , _ , _ = self.draw_model(theta)

        likelihood = self.log_likelihood(model_g1,model_g2,limit_mask)
        prior = self.log_prior(theta)
        if not np.isfinite(prior):
            posterior = -np.inf
        else:
            # use no info from prior for now
            posterior = likelihood 

        if logger.level == logging.DEBUG:
            n_progress = 10
        elif logger.level == logging.INFO:
            n_progress = 1000
        if self.n_model_evals % n_progress == 0:

            logger.info('%7d post=% 2.8e like=% 2.8e prior=% 2.4e M200=% 6.3e ' % (self.n_model_evals,posterior,likelihood,prior,theta[0]))

        if np.isnan(posterior):
            import pdb; pdb.set_trace()

        if self.save_all_models:

            self.plot_residual_g1g2(model_g1,model_g2,limit_mask)

            pl.suptitle('model post=% 10.8e M200=%5.2e' % (posterior,theta[0]) )
            filename_fig = 'models/res2.%04d.png' % self.n_model_evals
            pl.savefig(filename_fig)
            logger.debug('saved %s' % filename_fig)
            pl.close()


        return posterior

def plot_anat(brain):
    import os.path
    import pylab as pl
    from nibabel import load
    from nipy.labs import viz   
    import numpy as np

    img = load(brain)
    data = img.get_data()
    data[np.isnan(data)] = 0
    affine = img.get_affine() 
    viz.plot_anat(anat=data, anat_affine=affine, draw_cross=False, slicer='x')
    
    x_view = os.path.abspath('x_view.png')
    y_view = os.path.abspath('y_view.png')
    z_view = os.path.abspath('z_view.png')
    
    pl.savefig(x_view,bbox_inches='tight')
    
    viz.plot_anat(anat=data, anat_affine=affine, draw_cross=False, slicer='y')
    pl.savefig(y_view,bbox_inches='tight')
    
    viz.plot_anat(anat=data, anat_affine=affine, draw_cross=False, slicer='z')
    pl.savefig(z_view,bbox_inches='tight')
    
    images = [x_view, y_view, z_view]
    pl.close()
    return images

def plot_gc_distribution(pp, data):
    names = data.keys()

    # Plot the 2D histogram of coverage vs gc
    for name in names:
        x = [ i * 100 for i in data[name][GC_DISTRIBUTION_NAME]['gc_samples'] ]
        y = data[name][GC_DISTRIBUTION_NAME]['cov_samples']
        

        # Use the median to determine the range to show and round
        # to nearest 100 to avoid aliasing artefacts 
        m = np.median(y)
        y_limit = math.ceil( 2*m / 100) * 100
        hist,xedges,yedges = np.histogram2d(x,y, bins=[20, 50], range=[ [0, 100.0], [0, y_limit] ])

        # draw the plot
        extent = [xedges[0], xedges[-1], yedges[0], yedges[-1] ]
        pl.imshow(hist.T,extent=extent,interpolation='nearest',origin='lower', aspect='auto')

        pl.colorbar()
        pl.title(name + ' GC Bias')
        pl.xlabel("GC %")
        pl.ylabel("k-mer coverage")
        pl.savefig(pp, format='pdf')
        pl.close()

def plot_fragment_sizes(pp, data):

    # Trim outliers from the histograms
    names = data.keys()
    for name in names:
        h = data[name][FRAGMENT_SIZE_NAME]['sizes']
        sizes = {}
        for i in h:
            if i not in sizes:
                sizes[i] = 1
            else:
                sizes[i] += 1
        n = len(h)
        x = list()
        y = list()
        sum  = 0
        for i,j in sorted(sizes.items()):
            f = float(j) / n
            x.append(i)
            y.append(f)
            sum += f
        pl.plot(x, y)

    pl.xlim([0, 1000])
    pl.xlabel("Fragment Size (bp)")
    pl.ylabel("Proportion")
    pl.title("Estimated Fragment Size Histogram")
    pl.legend(names)
    pl.savefig(pp, format='pdf')
    pl.close()    

def test(cv,model,data,user,code,comp):
    test_power=np.array([float(r[2+comp])/max_power for r in data ])
    times=[datetime.datetime.strptime(r[0],'%Y-%m-%d %H:%M:%S UTC') for r in data]
    features=np.array([d[8:] for d in data],dtype=np.float)
    features[:,0]=features[:,0]/time_scale
    jobs=list(set([(r[1],r[2]) for r in data]))
    name_features=cv.transform([d[2] for d in data]).toarray()
    features=np.hstack((features,name_features))
    job_ids=[r[1] for r in data]
    prediction=model.predict(features)
    rmse=math.sqrt(np.average(((prediction-test_power)*max_power)**2))
    nrmse=math.sqrt(np.average(((prediction-test_power)/test_power)**2))
    corr=np.corrcoef(prediction,test_power)[0,1]
    r2=1-(sum((prediction-test_power)**2)/sum((test_power-np.average(test_power))**2))
    pl.figure(figsize=(6,7))
    pl.subplot(211)
    pl.plot(prediction*max_power,test_power*max_power,'+')
    if math.isnan(corr) or  math.isnan(r2) or math.isnan(rmse): 
        pl.title("RMSPE="+str(nrmse)+"RMSE="+str(rmse)+" Corr="+str(corr)+" R2="+str(r2))
    else:
        pl.title("RMSPE="+str(int(nrmse*1000)/1000.0)+" RMSE="+str(int(rmse*1000)/1000.0)+" Corr="+str(int(corr*1000)/1000.0)+" R2="+str(int(r2*1000)/1000.0))
    pl.xlabel('Predicted power')
    pl.ylabel('Real power')
    pl.plot([max(pl.xlim()[0],pl.ylim()[0]),min(pl.xlim()[1],pl.ylim()[1])],[max(pl.xlim()[0],pl.ylim()[0]),min(pl.xlim()[1],pl.ylim()[1])])
    pl.subplot(212)
    pl.plot(test_power*max_power)
    pl.plot(prediction*max_power)
    pl.ylabel('Power')
    pl.xlabel('Data point')
    #pl.legend(('Real power','Predicted power'))
    pl.subplots_adjust(hspace=0.35)
    pl.savefig('results'+str(month)+'global'+str(min_train)+'/'+user+code+'.pdf')
    pl.close()
    pkl.dump((nrmse,rmse,corr,r2,prediction*max_power,test_power*max_power,times,job_ids),file=gzip.open('results'+str(month)+'global'+str(min_train)+'/'+user+'test'+code+'.pkl.gz','w'))
    return prediction*max_power

def on_key_press(event):
    global old_t
    new_t=time.time()
    print old_t-new_t
    old_t=new_t
    if event.key == '+':
        a = axis()
        w = a[1] - a[0]
        axis([a[0] + w * .2, a[1] - w * .2, a[2], a[3]])
        draw()

    if event.key in ['-', '\'']:
        a = axis()
        w = a[1] - a[0]
        axis([a[0] - w / 3.0, a[1] + w / 3.0, a[2], a[3]])
        draw()

    if event.key == '.':
        a = axis()
        w = a[1] - a[0]
        axis([a[0] + w * .2, a[1] + w * .2, a[2], a[3]])
        draw()

    if event.key == ',':
        a = axis()
        w = a[1] - a[0]
        axis([a[0] - w * .2, a[1] - w * .2, a[2], a[3]])
        draw()

    if event.key == 'q':
        close()

def main():
    base_path = "/caps2/tsupinie/1kmf-control/"
    temp = goshen_1km_temporal(start=14400, end=14400)
    grid = goshen_1km_grid()
    n_ens_members = 40

    np.seterr(all='ignore')

    ens = loadEnsemble(base_path, [ 11 ], temp.getTimes(), ([ 'pt', 'p' ], computeDensity))
    ens = ens[0, 0]

    zs = decompressVariable(nio.open_file("%s/ena001.hdfgrdbas" % base_path, mode='r', format='hdf').variables['zp'])
    xs, ys = grid.getXY()
    xs = xs[np.newaxis, ...].repeat(zs.shape[0], axis=0)
    ys = ys[np.newaxis, ...].repeat(zs.shape[0], axis=0)

    eff_buoy = effectiveBuoyancy(ens, (zs, ys, xs), plane={'z':10})
    print eff_buoy

    pylab.figure()
    pylab.contourf(xs[0], ys[0], eff_buoy[0], cmap=matplotlib.cm.get_cmap('RdBu_r'))
    pylab.colorbar()

    grid.drawPolitical()

    pylab.suptitle("Effective Buoyancy")
    pylab.savefig("eff_buoy.png")
    pylab.close()
    return

def plot_values(X, Y, xlabel, ylabel, suffix, ptype='plot'):
    output_filename = constants.ATTRACTIVENESS_FOLDER_NAME + constants.DATASET + '_' + suffix

    X1 = [X[i] for i in range(len(X)) if X[i]>0 and Y[i]>0]
    Y1 = [Y[i] for i in range(len(X)) if X[i]>0 and Y[i]>0]
    X = X1
    Y = Y1
    
    pylab.close("all")
    
    pylab.figure(figsize=(8, 7))

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

    pylab.scatter(X, Y)
    
    #pylab.axis(vis.get_bounds(X, Y, False, False))

    #pylab.xscale('log')
    pylab.yscale('log')

    pylab.xlabel(xlabel)
    pylab.ylabel(ylabel)   
    #pylab.xlim(0.1,1)
    #pylab.ylim(ymin=0.01)
    #pylab.tight_layout()

    pylab.savefig(output_filename + '.pdf')

    def test_varying_inclination(self):
        #""" Test that the waveform is consistent for changes in inclination
        #"""
        sigmas = []
        incs = numpy.arange(0, 21, 1.0) * lal.PI / 10.0

        for inc in incs:
            # WARNING: This does not properly handle the case of SpinTaylor*
            # where the spin orientation is not relative to the inclination
            hp, hc = get_waveform(self.p, inclination=inc)
            s = sigma(hp, low_frequency_cutoff=self.p.f_lower)        
            sigmas.append(s)
         
        f = pylab.figure()
        pylab.axes([.1, .2, 0.8, 0.70])   
        pylab.plot(incs, sigmas)
        pylab.title("Vary %s inclination, $\\tilde{h}$+" % self.p.approximant)
        pylab.xlabel("Inclination (radians)")
        pylab.ylabel("sigma (flat PSD)")
        
        info = self.version_txt
        pylab.figtext(0.05, 0.05, info)
        
        if self.save_plots:
            pname = self.plot_dir + "/%s-vary-inclination.png" % self.p.approximant
            pylab.savefig(pname)

        if self.show_plots:
            pylab.show()
        else:
            pylab.close(f)

        self.assertAlmostEqual(sigmas[-1], sigmas[0], places=7)
        self.assertAlmostEqual(max(sigmas), sigmas[0], places=7)
        self.assertTrue(sigmas[0] > sigmas[5])

def CostVariancePlot(funct,args):
	pl=args[0]
	x=np.array([])
	y=np.array([])
	f=np.array([])
	z=np.array([])
	x=np.append(x,funct.rmsSet[:,0])
	y=np.append(y,funct.rmsSet[:,3])
	f=np.append(f,funct.rmsSet[:,5])
	z=np.append(z,funct.rmsSet[:,4])
	v=np.array([])
	v=np.append(v,[0])
	i=0
	while i<len(x):
		v=np.append(v,(z[i]/f[i])-(y[i]/f[i])*(y[i]/f[i]))
		i+=1
	v=np.delete(v,0)
	if centroidP(x,v):
		pl.set_yscale('log')
		pl.set_xscale('log') 
	else:
		pl.ticklabel_format(axis='both', style='sci', scilimits=(-2,5),pad=5,direction="bottom")
	pl.axis([0, np.amax(x)+(10*np.amax(x)/100), 0, np.amax(v)+(10*np.amax(v)/100)])
	pl.set_xlabel("read memory size",fontsize=8)
	pl.set_ylabel("cost",fontsize=8)
	pl.set_title("Variance Cost",fontsize=14)
	pl.grid(True)
	pl.tick_params(axis='x', labelsize=7)
	pl.tick_params(axis='y', labelsize=7)
	sc=pl.scatter(x,v,c=f,s=6,marker = 'o',lw=0.0,cmap=cmap,norm=norm)
	pylab.close()		

def plot_quality_scores(pp, data):
    names = data.keys()

    # Plot mean quality
    for name in names:

        mean_quality = data[name][QUALITY_SCORE_NAME]['mean_quality']
        indices = range(0, len(mean_quality))
        pl.plot(indices, mean_quality, linewidth=2)

    pl.ylim([0, 40])
    pl.xlabel("Base position")
    pl.ylabel("Mean Phred Score")
    pl.title("Mean quality score by position")
    pl.legend(names, loc="lower left")
    pl.savefig(pp, format='pdf')
    pl.close()

    # Plot >q30 fraction
    for name in names:

        q30_fraction = data[name][QUALITY_SCORE_NAME]['fraction_q30']
        indices = range(0, len(q30_fraction))
        pl.plot(indices, q30_fraction)

    pl.xlabel("Base position")
    pl.ylabel("Fraction at least Q30")
    pl.title("Fraction of bases at least Q30")
    pl.legend(names, loc="lower left")
    pl.savefig(pp, format='pdf')
    pl.close()    

def plot(x, z):             # plot input x versus convolution z
    pylab.close(1)
    pylab.figure(1)
    pylab.plot (x, 'b', label="sine")
    pylab.plot ([abs(i/max(z)) for i in z], 'r', label="convolution")
    pylab.title ("convolution of pure sine")
    pylab.legend()

def arcRVs(booSave = False, booShow = True, booFit = False):
    arcRVs = np.load('npy/arcRVs.npy')
    MJDs = np.load('npy/JDs.npy')
    
    colors = ['b','g','r','c']

    for epoch,MJD in enumerate(MJDs):
        for cam in range(4):
            y = arcRVs[:,epoch,cam]
            plt.plot(y,'.'+colors[cam])
            if booFit==True:
                x = np.arange(len(y))
                p = np.polyfit(x[-np.isnan(y)],y[-np.isnan(y)],1)
                plt.plot(x,x*p[0]+p[1])
                
        plt.title('MJD='+str(MJD))
        if booSave==True: 
            try:
                plotName = 'plots/arcRVs_'+str(epoch)
                print 'Attempting to save', plotName
                plt.savefig(plotName)

            except Exception,e: 
                print str(e)
                print 'FAILED'
        if booShow==True: plt.show()
        plt.close()        

 def save(self, name=None, format="png", dirc=None):
     """Saves Bloch sphere to file of type ``format`` in directory ``dirc``.
     
     Parameters
     ----------
     
     name : str 
         Name of saved image. Must include path and format as well.
         i.e. '/Users/Paul/Desktop/bloch.png'
         This overrides the 'format' and 'dirc' arguments.
     format : str 
         Format of output image.  
     dirc : str
         Directory for output images. Defaults to current working directory.
     
     Returns
     ------- 
     File containing plot of Bloch sphere.
     
     """
     self.make_sphere()
     if dirc:
         if not os.path.isdir(os.getcwd() + "/" + str(dirc)):
             os.makedirs(os.getcwd() + "/" + str(dirc))
     if name == None:
         if dirc:
             savefig(os.getcwd() + "/" + str(dirc) + "/bloch_" + str(self.savenum) + "." + format)
         else:
             savefig(os.getcwd() + "/bloch_" + str(self.savenum) + "." + format)
     else:
         savefig(name)
     self.savenum += 1
     if self.fig:
         close(self.fig)

def generate_glassbrain_image(image_pk):
    from neurovault.apps.statmaps.models import Image
    import neurovault
    import matplotlib as mpl
    mpl.rcParams['savefig.format'] = 'jpg'
    my_dpi = 50
    fig = plt.figure(figsize=(330.0/my_dpi, 130.0/my_dpi), dpi=my_dpi)
    
    img = Image.objects.get(pk=image_pk)    
    f = BytesIO()
    try:
        glass_brain = plot_glass_brain(img.file.path, figure=fig)
        glass_brain.savefig(f, dpi=my_dpi)
    except:
        # Glass brains that do not produce will be given dummy image
        this_path = os.path.abspath(os.path.dirname(__file__))
        f = open(os.path.abspath(os.path.join(this_path,
                                              "static","images","glass_brain_empty.jpg"))) 
        raise
    finally:
        plt.close('all')
        f.seek(0)
        content_file = ContentFile(f.read())
        img.thumbnail.save("glass_brain_%s.jpg" % img.pk, content_file)
        img.save()

    def plot(self, outputDirectory):
        """
        Plot both the raw kinetics data and the Arrhenius fit versus 
        temperature. The plot is saved to the file ``kinetics.pdf`` in the
        output directory. The plot is not generated if ``matplotlib`` is not
        installed.
        """
        # Skip this step if matplotlib is not installed
        try:
            import pylab
        except ImportError:
            return

        Tlist = 1000.0/numpy.arange(0.4, 3.35, 0.05)
        klist = numpy.zeros_like(Tlist)
        klist2 = numpy.zeros_like(Tlist)
        for i in range(Tlist.shape[0]):
            klist[i] = self.reaction.calculateTSTRateCoefficient(Tlist[i])
            klist2[i] = self.reaction.kinetics.getRateCoefficient(Tlist[i])

        order = len(self.reaction.reactants)
        klist *= 1e6 ** (order-1)
        klist2 *= 1e6 ** (order-1)

        pylab.semilogy(1000.0 / Tlist, klist, 'ok')
        pylab.semilogy(1000.0 / Tlist, klist2, '-k')
        pylab.xlabel('1000 / Temperature (1000/K)')
        pylab.ylabel('Rate coefficient ({0})'.format(self.kunits))
        pylab.savefig(os.path.join(outputDirectory, 'kinetics.pdf'))
        pylab.close()