Extending existing models (with a side-order of XSPEC)

I was recently asked by a colleage whether Sherpa allowed you to subtract a model. This is "easy" - since you can create model expressions using +, -, *, and / - but it turned out that what was really wanted fitted into my current documentation obsession, writing user models in Sherpa. It also lets me play Sherpa's support for XSPEC models.

The aim of this notebook is to show you how you can write models in Sherpa that use existing Sherpa models. In this particular case I will create a model that is the result of subtracting two XSPEC vvnei models, providing an interface where the Hydrogen abundance is no-longer a parameter (I make no claim that this is a physically-meaningful model), but the approach is valid for all Sherpa models. The notebook ends up with a brief discussion about optimising model, but this is meant more as a primer on the subject than an in-depth discussion.

I am going to assume that you have read at least Writing your own user model but please feel free to look through all the notebooks I've written in the series.

Author and disclaimer

This was written by Douglas Burke on June 16 2015. This notebook, and others that may be of interest, can be found on GitHub at https://github.com/DougBurke/sherpa-standalone-notebooks.

The information in this document is placed into the Publc Domain. It is not an official product of the Chandra X-ray Center, and I make no guarantee that it is not without bugs or embarassing typos. Please contact me via the GitHub repository or on Twitter - at @doug_burke - if you have any questions.

Setting up for XSPEC

If you are using the Sherpa provided with CIAO then you can skip this section, as the XSPEC extension is provided for you.

If you are using the standalone version of Sherpa, then you need to have the XSPEC libraries and then tell Sherpa where they are, before building Sherpa. The current (as of mid-June, 2015) documentation on this process is a bit lacking, but I have a pull request that should improve this. For today's notebook I am going to assume that you have built a version of XSPEC using HEASoft version 6.16, that is, XSPEC version 12.8.2 (with or without patches). The important parts of the installation are:

  • $HEADAS/lib - for building
  • $HEADAS/../spectral - for using the models

The xspec_config section of Sherpa's setup.cfg file must be edited so that it says (the HEADAS environment variable should be expanded, and the lib sub-directory checked to see if the version numbers on the library names need adjusting):

with-xspec=True
xspec_lib_dirs=$HEADAS/lib
xspec_libraries=XSFunctions XSModel XSUtil XS wcs-4.20
cfitsio_libraries=cfitsio_3.37
ccfits_libraries=CCfits_2.4
cfitsio_lib_dirs=$HEADAS/lib
ccfits_lib_dirs=$HEADAS/lib
gfortran_lib_dirs=$HEADAS/lib

After this, run

python setup.py build
python setup.py test
python setup.py install

This assumes that my patch has not been accepted yet, so that there is no wcslib_libraries configuration option. If this has been added (or something similar) then please follow the instructions provided with Sherpa rather than these ones!

Loading the XSPEC models

The XSPEC models are automatically loaded - if available - by the sherpa.astro.ui module. Before I load any of these, I am going to turn off the logging information that confuses the ipython notebooks; for use in scripts or within ipython on the command line this is not needed:

In [1]:
import logging
logging.getLogger('sherpa').propagate = 0

With this, we can now load Sherpa, and I also load in the XSPEC module directly, so that I can access several low-level routines:

In [2]:
import numpy as np
from matplotlib import pyplot as plt

from sherpa.astro import ui
from sherpa.astro import xspec
WARNING: imaging routines will not be available, 
failed to import sherpa.image.ds9_backend due to 
'RuntimeErr: DS9Win unusable: Could not find ds9 on your PATH'
In [3]:
%matplotlib inline

To check that things are working, we can find out the version of the XSPEC library in use (for this example I did not apply any of the patches that have been made to 12.8.2):

In [4]:
xspec.get_xsversion()
Out[4]:
'12.8.2'

If this command fails with the message

ImportError: XSPEC initialization failed; check HEADAS environment variable

then check that this variable is set correctly. It should be set to the parent directory of the lib directory used for the xspec_lib_dirs setting used to build Sherpa. For example

import os
os.environ['HEADAS'] = '/usr/local/headas-6.16/x86_64-unknown-linux-gnu-libc2.21-0/'

and then - assuming this is the correct setting ;-) - then the xspec.get_xsversion() command should work.

In [5]:
xspec.get_xschatter()
Out[5]:
0

The XSPEC chatter level - which can be read with xspec.get_xschatter and set with xspec.set_xschatter can be used to change the amount of output printed to the screen when X-Spec models are evaluated. The default value used by Sherpa - of 0 - means that very-little information is displayed. The default value used by XSPEC itself is 10 - and values around 25 can be used for debugging (e.g. to check that data files are being read correctly; although the amount of information provided depends on the model and the version of XSPEC in use).

The model

The model I am going to develop is a version of the vvnei model - which is a non-equilibrium ionisation model with variable abundances - which has no Hydrogen. The standard way for XSPEC to model variable abundances is to make them all relative to Hydrogen, but in this case we (or at least, my colleague), does not want any Hydrogren! One solution for this is to evaluate the model twice, once with a Hydrogen abundance of 1, and the required elemental abundances, and then a Hydrogen-only version. The second version is subtracted from the first (hence the original question I was asked) to return the required model.

This can be done directly in Sherpa, using two different instances of the xsvvnei model, but we can also write our own user model to wrap up this behavior, and provide a slightly-simpler interface to the user.

Before we write the user model, let's try and see the individual components. The following creates two instances of the vvnei model (Sherpa prepends xs to the XSPEC model names), called mdl_with_elems and mdl_h_only. I then set several of the elemental abundances in the mdl_with_elems component - using ui.set_par - and then clear out all the non-Hydrogen abundances in the mdl_h_only component (I happen to know that the check I wrote is sufficient for this model type to select only the abundances I wanted to change). This code shows two different ways of accessing model parameters (by name, with ui.set_par, or by iterating directly from the model component), but I don't want to talk too much about this today:

In [6]:
ui.xsvvnei.mdl_with_elems
ui.xsvvnei.mdl_h_only

# Set some abundances
abunds = { 'c': 2.3, 'fe': 1.8, 'o': 3.4, 'si': 2.1, 'ar': 1.2 }
for (k,v) in abunds.iteritems():
    ui.set_par('mdl_with_elems.' + k, v)
    
# Clear out the abundances for the Hydrogen-only version
for p in mdl_h_only.pars:
    if len(p.name) < 3 and p.name not in ['kT', 'H']:
        p.val = 0

This particular model has a lot of parameters!

In [7]:
print(mdl_with_elems)
xsvvnei.mdl_with_elems
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   mdl_with_elems.kT thawed            1       0.0808         79.9        keV
   mdl_with_elems.H frozen            1            0         1000           
   mdl_with_elems.He frozen            1            0         1000           
   mdl_with_elems.Li frozen            1            0         1000           
   mdl_with_elems.Be frozen            1            0         1000           
   mdl_with_elems.B frozen            1            0         1000           
   mdl_with_elems.C frozen          2.3            0         1000           
   mdl_with_elems.N frozen            1            0         1000           
   mdl_with_elems.O frozen          3.4            0         1000           
   mdl_with_elems.F frozen            1            0         1000           
   mdl_with_elems.Ne frozen            1            0         1000           
   mdl_with_elems.Na frozen            1            0         1000           
   mdl_with_elems.Mg frozen            1            0         1000           
   mdl_with_elems.Al frozen            1            0         1000           
   mdl_with_elems.Si frozen          2.1            0         1000           
   mdl_with_elems.P frozen            1            0         1000           
   mdl_with_elems.S frozen            1            0         1000           
   mdl_with_elems.Cl frozen            1            0         1000           
   mdl_with_elems.Ar frozen          1.2            0         1000           
   mdl_with_elems.K frozen            1            0         1000           
   mdl_with_elems.Ca frozen            1            0         1000           
   mdl_with_elems.Sc frozen            1            0         1000           
   mdl_with_elems.Ti frozen            1            0         1000           
   mdl_with_elems.V frozen            1            0         1000           
   mdl_with_elems.Cr frozen            1            0         1000           
   mdl_with_elems.Mn frozen            1            0         1000           
   mdl_with_elems.Fe frozen          1.8            0         1000           
   mdl_with_elems.Co frozen            1            0         1000           
   mdl_with_elems.Ni frozen            1            0         1000           
   mdl_with_elems.Cu frozen            1            0         1000           
   mdl_with_elems.Zn frozen            1            0         1000           
   mdl_with_elems.Tau thawed        1e+11        1e+08        5e+13     s/cm^3
   mdl_with_elems.redshift frozen            0       -0.999           10           
   mdl_with_elems.norm thawed            1            0        1e+24           

We can see that the 'Hydrogen only' one has been cleared out:

In [8]:
print(mdl_h_only)
xsvvnei.mdl_h_only
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   mdl_h_only.kT thawed            1       0.0808         79.9        keV
   mdl_h_only.H frozen            1            0         1000           
   mdl_h_only.He frozen            0            0         1000           
   mdl_h_only.Li frozen            0            0         1000           
   mdl_h_only.Be frozen            0            0         1000           
   mdl_h_only.B frozen            0            0         1000           
   mdl_h_only.C frozen            0            0         1000           
   mdl_h_only.N frozen            0            0         1000           
   mdl_h_only.O frozen            0            0         1000           
   mdl_h_only.F frozen            0            0         1000           
   mdl_h_only.Ne frozen            0            0         1000           
   mdl_h_only.Na frozen            0            0         1000           
   mdl_h_only.Mg frozen            0            0         1000           
   mdl_h_only.Al frozen            0            0         1000           
   mdl_h_only.Si frozen            0            0         1000           
   mdl_h_only.P frozen            0            0         1000           
   mdl_h_only.S frozen            0            0         1000           
   mdl_h_only.Cl frozen            0            0         1000           
   mdl_h_only.Ar frozen            0            0         1000           
   mdl_h_only.K frozen            0            0         1000           
   mdl_h_only.Ca frozen            0            0         1000           
   mdl_h_only.Sc frozen            0            0         1000           
   mdl_h_only.Ti frozen            0            0         1000           
   mdl_h_only.V frozen            0            0         1000           
   mdl_h_only.Cr frozen            0            0         1000           
   mdl_h_only.Mn frozen            0            0         1000           
   mdl_h_only.Fe frozen            0            0         1000           
   mdl_h_only.Co frozen            0            0         1000           
   mdl_h_only.Ni frozen            0            0         1000           
   mdl_h_only.Cu frozen            0            0         1000           
   mdl_h_only.Zn frozen            0            0         1000           
   mdl_h_only.Tau thawed        1e+11        1e+08        5e+13     s/cm^3
   mdl_h_only.redshift frozen            0       -0.999           10           
   mdl_h_only.norm thawed            1            0        1e+24           

The models can be evaluated on the grid 0.1 to 10 keV, with bin spacing 0.01 keV (XSPEC models expect the input grid to be in keV). The line emission from the elements is visible here:

In [9]:
egrid = np.arange(0.1, 10, 0.01)
plt.plot(egrid, mdl_with_elems(egrid))
plt.plot(egrid, mdl_h_only(egrid))
plt.yscale('log')
plt.xscale('log')

Subtracting these two will give the elemental emission for a situation where there is no Hydrogen. I also convert the plot values into those more commonly displayed: the XSPEC models return the flux in units of photon/cm^2/s for each bin, so we need to divide by the bin width to get the flux density. Since I gave the models a single argument, it was taken to be the lower value of each bin edge, so for n bins the return value is n-1 values, with the last one being 0 (since the input and output array lengths match, we add a 0 to the end of the model values); this is why the egrid and y values are suffixed with [:-1] when plotting, to remove this last element.

In [10]:
de = egrid[1:] - egrid[:-1]
y = mdl_with_elems(egrid) - mdl_h_only(egrid) 
plt.plot(egrid[:-1], y[:-1] / de)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Energy (keV)')
plt.ylabel('photon/cm$^2$/s/keV')
Out[10]:
<matplotlib.text.Text at 0x7f5d84ff7210>

Aside

Here I evaluate one of the models on the grid 0.1 to 0.2, 0.2 to 0.3, 0.3 to 0.4 keV. Note that a 0 is added onto the end to ensure that the return array has the same length as the input array:

In [11]:
mdl_with_elems([0.1, 0.2, 0.3, 0.4])
Out[11]:
array([ 0.43303607,  0.54599511,  0.15643566,  0.        ])

If I explicitly give arrays for the lower and upper edges then we do not get this extra 0 element:

In [12]:
mdl_with_elems([0.1,0.2,0.3], [0.2,0.3,0.4])
Out[12]:
array([ 0.43303607,  0.54599511,  0.15643566])

End aside

Using the existing models

Since we can write model expressions using subtraction, then the "Hydrogen free" variant of xsvvnei could be handled explicitly, by using the model combination

(ui.xsvvnei.mdl_elems - ui.xsvvnei.mdl_honly)

in the source expression. It then requires

  • linking the kT, Tau, redshift, and norm parameters: for example

     ui.link(mdl_honly.kt, mdl_elems.kt)
     ui.link(mdl_honly.tau, mdl_elems.tau)
     ui.link(mdl_honly.redshift, mdl_elems.redshift)
     ui.link(mdl_honly.norm, mdl_elems.norm)
  • setting the mdl_honly elemental abundances to 0:

     for p in mdl_honly.pars:
         if p.name.lower() not in ['kt', 'h', 'tau', redshift', 'norm']:
             p.val = 0
             p.frozen = True # this is actually the default, but just check
  • fixing the Hydrogen abundances

     mdl_honly.h = 1
     mdl_elems.h = 1
     ui.freeze(mdl_honly.h, mdl_elems.h)
  • thaw the elemental abundances in mdl_elems that you want to fit.

This can be done, but if you started using this model a lot, you'd quickly find it annoying (and easy to miss a step). What is needed is a way to create a user model - as described in previous notebooks - but instead of writing the function, use the existing Sherpa code to evaluate the models.

Writing a simple user model

So, can we simplify the above by writing a simple user model? First of all, let's look at evaluating the model for a given set of parameters. The evaluations above have required an instance of the model - which I created by calling

ui.xsvvnei.mdl_with_elems

but could have also used

ui.create_model_component('xsvvnei', 'mdl_with_elems')

I could use this approach in a user model, but it would involve coming up with unique model names which is messy and inelegant. What I really want to be able to do is evaluate the model without creating a component. This can be done by using parts of Sherpa directly. In this case, the models are defined in the sherpa.astro.xspec module, and we can create an object directly:

In [13]:
mdl = xspec.XSvvnei()

The model can be given a name by sending in a string to the constructor, but this is only used in screen output which we do not care about, so we let the default value be used. The only real difference to calling this, rather than ui.create_model_component (or writing ui.xsvvnei.xxx), is that this model is not added to the Sherpa session, so it does not conflict with anything (e.g. ui.list_model_components() will not include it). This makes it ideal for our purposes, since it is essentially an anonymous model.

An anonymous model function is a start, but I also need an easy way to call the model with different parameters: one set with the elemental abundances, and the second time with Hydrogen only. I could set one set of parameters, evaluate the model, and then swap in the second set and re-evaluate. However, the Sherpa model class defines the calc method, which has the signature:

In [14]:
help(mdl.calc)
Help on method calc in module sherpa.models.model:

calc(cls, pars, xlo, *args, **kwargs) method of sherpa.astro.xspec.XSvvnei instance

The pars argument is the array of parameter values (in the order displayed from print(mdl) or returned by mdl.pars) and the remaining arguments are the independent axes on which to evaluate the model. This is therefore the second part to writing the user model.

In [15]:
def xsvvnei_no_h(pars, *args):
    "pars is kT, list of 29 element abundances (excluding H), tau, redshift, norm"

    # It's a bit excessive to create this each iteration
    mdl = xspec.XSvvnei()
    
    pall = [pars[0], 1.0]
    pall.extend(pars[1:])
    
    phonly = [pars[0], 1.0]
    phonly.extend([0.0] * 29)
    phonly.extend(pars[-3:])

    # Evaluate the "Hydrogen-free" model
    return mdl.calc(pall, *args) - mdl.calc(phonly, *args)

As a check, does this give the same answer as evaluating the separate models?

In [16]:
pars = [p.val for p in mdl_with_elems.pars if p.name != "H"]
ynew = xsvvnei_no_h(pars, egrid)
y = mdl_with_elems(egrid) - mdl_h_only(egrid) 
np.all(ynew == y)
Out[16]:
True

This can be added to Sherpa using ui.load_user_model to create a model instance (i.e. component). Since I want to copy over the parameter settings from the xsvvnei model (at least, apart from for H), I have written a small function that both loads the model - with ui.load_user_model - and sets the parameter values - with ui.add_user_pars:

In [17]:
def make_xsvvnei_noh_model(name):
    """Create an instance of the 'xsvvnei' model which has
    no Hydrogen."""
    
    # Create the model instance
    ui.load_user_model(xsvvnei_no_h, name)
    
    # Get a copy of the default parameter settings for the vvnei
    # model (excluding Hydrogen).
    pars = [p for p in xspec.XSvvnei().pars if p.name != 'H']
    
    parnames = [p.name for p in pars]
    parvals = [p.val for p in pars]
    parmins = [p.min for p in pars]
    parmaxs = [p.max for p in pars]
    parunits = [p.units for p in pars]
    parf = [p.frozen for p in pars]
    
    # Set up the parameters for this instance
    ui.add_user_pars(name, parnames, parvals=parvals, parmins=parvals,
                     parmaxs=parmaxs, parunits=parunits, parfrozen=parf)
    
    # Return the model
    return ui.get_model_component(name)
In [18]:
make_xsvvnei_noh_model('noh')
Out[18]:
<UserModel model instance 'usermodel.noh'>
In [19]:
print(noh)
usermodel.noh
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   noh.kT       thawed            1            1         79.9        keV
   noh.He       frozen            1            1         1000           
   noh.Li       frozen            1            1         1000           
   noh.Be       frozen            1            1         1000           
   noh.B        frozen            1            1         1000           
   noh.C        frozen            1            1         1000           
   noh.N        frozen            1            1         1000           
   noh.O        frozen            1            1         1000           
   noh.F        frozen            1            1         1000           
   noh.Ne       frozen            1            1         1000           
   noh.Na       frozen            1            1         1000           
   noh.Mg       frozen            1            1         1000           
   noh.Al       frozen            1            1         1000           
   noh.Si       frozen            1            1         1000           
   noh.P        frozen            1            1         1000           
   noh.S        frozen            1            1         1000           
   noh.Cl       frozen            1            1         1000           
   noh.Ar       frozen            1            1         1000           
   noh.K        frozen            1            1         1000           
   noh.Ca       frozen            1            1         1000           
   noh.Sc       frozen            1            1         1000           
   noh.Ti       frozen            1            1         1000           
   noh.V        frozen            1            1         1000           
   noh.Cr       frozen            1            1         1000           
   noh.Mn       frozen            1            1         1000           
   noh.Fe       frozen            1            1         1000           
   noh.Co       frozen            1            1         1000           
   noh.Ni       frozen            1            1         1000           
   noh.Cu       frozen            1            1         1000           
   noh.Zn       frozen            1            1         1000           
   noh.Tau      thawed        1e+11        1e+11        5e+13     s/cm^3
   noh.redshift frozen            0            0           10           
   noh.norm     thawed            1            1        1e+24           

A quick check shows how the Oxygen abundance changes the results for the default temperature:

In [20]:
plt.plot(egrid, noh(egrid))
plt.xscale('log')
plt.yscale('log')

noh.o = 5
plt.plot(egrid, noh(egrid))
Out[20]:
[<matplotlib.lines.Line2D at 0x7f5d84dc6e90>]

Using a class-based approach

As with the 'User models in Sherpa' notebook, we can also create a class for this model:

In [21]:
from sherpa.models import model

class XSvvneiNoH(model.ArithmeticModel):
    """A version of XSvvnei which has no Hydrogen. The
    abundances of the other elements are given relative to
    solar abundance of Hydrogen.
    """
    
    def __init__(self, name='xxvvneinoh'):
        mdl = xspec.XSvvnei()
        
        pars = []
        for par in mdl.pars:
            if par.name == 'H':
                continue
                
            newpar = model.Parameter(name, par.name, par.val, min=par.min, max=par.max,
                                     hard_min=par.hard_min, hard_max=par.hard_max,
                                     units=par.units, frozen=par.frozen,
                                     alwaysfrozen=par.alwaysfrozen, hidden=par.hidden)
            setattr(self, par.name, newpar)
            pars.append(newpar)
            
        model.ArithmeticModel.__init__(self, name, tuple(pars))
        
    @model.modelCacher1d
    def calc(self, pars, x, *args, **kwargs):
        return xsvvnei_no_h(pars, x, *args)

# Tell Sherpa about the model
ui.add_model(XSvvneiNoH)

This class can then be used to create model components: in the following I show the lines from Carbon, Oxygen, and Argon lines.

In [22]:
ui.create_model_component('xsvvneinoh', 'mdl')

# Clear out all the elemental abundances
for par in mdl.pars:
    if len(par.name) < 3 and par.name != 'kT':
        par.val = 0

# Colors are taken from ColorBrewer2, or rather,
# this StackOverflow answer: http://stats.stackexchange.com/a/118183
        
# Ensure there's some continuum
mdl.he = 1

# Carbon
mdl.c = 1
plt.plot(egrid, mdl(egrid), c='#e41a1c')

# Oxygen
mdl.c = 0
mdl.o = 1
plt.plot(egrid, mdl(egrid), c='#377eb8')

# Argon
mdl.o = 0
mdl.ar = 1
plt.plot(egrid, mdl(egrid), c='#4daf4a')

# Continuum
mdl.ar = 0
plt.plot(egrid, mdl(egrid), c='#984ea3')

plt.legend(['C', 'O', 'Ar', 'Continuum'])
plt.xscale('log')
plt.yscale('log')

Can we make this faster?

If I were to be using this for actual analysis, I would probably re-write xsvvnei_no_h so that it is sent the model instance to use, rather than re-create it each time the model is evaluated, since the latter feels wasteful. However, it is really important to actually measure the time these functions take to run to make sure changes you make are worth it, so I'll do that too.

In [23]:
# I repeat the original version as a reference.
def xsvvnei_no_h_1(pars, *args):
    "pars is kT, list of 29 element abundances (excluding H), tau, redshift, norm"

    # It's a bit excessive to create this each iteration
    mdl = xspec.XSvvnei()
    
    pall = [pars[0], 1.0]
    pall.extend(pars[1:])
    
    phonly = [pars[0], 1.0]
    phonly.extend([0.0] * 29)
    phonly.extend(pars[-3:])

    # Evaluate the "Hydrogen-free" model
    return mdl.calc(pall, *args) - mdl.calc(phonly, *args)

# This version requires a model instance to be sent in
def xsvvnei_no_h_2(mdl, pars, *args):
    "pars is kT, list of 29 element abundances (excluding H), tau, redshift, norm"

    pall = [pars[0], 1.0]
    pall.extend(pars[1:])
    
    phonly = [pars[0], 1.0]
    phonly.extend([0.0] * 29)
    phonly.extend(pars[-3:])

    # Evaluate the "Hydrogen-free" model
    return mdl.calc(pall, *args) - mdl.calc(phonly, *args)    

The IPython interpreter provides the %timeit command, which is a basic - but useful, as I'll show - way to time functions. In this case, I want to evaluate the two functions, which means I need a parameter array, a model (for xsvvnei_no_h_2), and I'll use egrid for the grid.

In [24]:
pars = [p.val for p in mdl.pars]
mdl_2 = xspec.XSvvnei()

The %timeit command takes the function to evaluate, and runs it multiple times, and reports its best estimate for the run time (it's an adaptive scheme, in that it tailors the number of loops to how long the function takes, and reports the time taken in "sensible" units).

In [25]:
%timeit xsvvnei_no_h_1(pars, egrid)
10 loops, best of 3: 20.2 ms per loop
In [26]:
%timeit xsvvnei_no_h_2(mdl_2, pars, egrid)
100 loops, best of 3: 11.9 ms per loop

As we can see in this example, the second function is roughly two times faster than the first one (which surprised me). Now, there's plenty of things to be wary of when using this, since

  • I have not actualy checked that the two functions return the same values (which you really must do before trying to work out which is faster);

  • it relies on the function always doing the same amount of work each iteration (it might be that the function loads in data files for particular parameter values, or sets up memory or loads in data the first time it is run);

  • I have only used a single set of parameters and one energy grid and the difference in run time could well vary depending on these values;

  • when used in a fit, there is plenty of other code being run, so there comes a point at which reducing the run time of a model component is no longer worth it (perhaps there's another model being evaluated which is significantly slower).

So, this is a start, but it's not the end of the story. There are other commands and modules in Python that can help you optimise the run time (or memory usage), of Python code (e.g. Timing and Profiling, SciPy notes on Optimizing code, Jake VanderPlas' example of optimizing Python code, SnakeViz for visualizing the cProfile output, and plenty of other resources out there on the internet).

Just one more thing

The code above uses the generic Sherpa model interface to evaluate a model; that is, it uses the calc method of the Model class to evaluate the model, and can be used for all Sherpa models. I happen to know that, for XSPEC models, this method is defined as

class XSModel(ArithmeticModel):

    @modelCacher1d
    def calc(self, *args, **kwargs):
        return self._calc(*args, **kwargs)

and that the XSvvnei class defines _calc to be:

class XSvvnei(XSAdditiveModel):

    _calc =  _xspec.C_vvnei
...

which means that we can call _calc instead of calc, saving even-more time. Let's see how much this really helps:

In [27]:
from sherpa.astro.xspec import _xspec

def xsvvnei_no_h_3(mdl, pars, *args):
    "pars is kT, list of 29 element abundances (excluding H), tau, redshift, norm"

    pall = [pars[0], 1.0]
    pall.extend(pars[1:])
    
    phonly = [pars[0], 1.0]
    phonly.extend([0.0] * 29)
    phonly.extend(pars[-3:])

    # Evaluate the "Hydrogen-free" model
    return mdl._calc(pall, *args) - mdl._calc(phonly, *args)

def xsvvnei_no_h_4(pars, *args):
    "pars is kT, list of 29 element abundances (excluding H), tau, redshift, norm"

    pall = [pars[0], 1.0]
    pall.extend(pars[1:])
    
    phonly = [pars[0], 1.0]
    phonly.extend([0.0] * 29)
    phonly.extend(pars[-3:])

    # Evaluate the "Hydrogen-free" model
    return _xspec.C_vvnei(pall, *args) - _xspec.C_vvnei(phonly, *args)    
In [28]:
%timeit xsvvnei_no_h_1(pars, egrid)
%timeit xsvvnei_no_h_2(mdl_2, pars, egrid)
%timeit xsvvnei_no_h_3(mdl_2, pars, egrid)
%timeit xsvvnei_no_h_4(pars, egrid)
10 loops, best of 3: 20.3 ms per loop
100 loops, best of 3: 11.9 ms per loop
100 loops, best of 3: 11.8 ms per loop
100 loops, best of 3: 11.9 ms per loop

So, in this case, the extra "optimisations" do not really save much time. Given these results, I would be reluctant to use them, since they rely on internals of the code that could be changed at any time.