User-defined Sherpa model types using Python Classes

Here is a quick example of a native Sherpa model (logparlaw1d) defined as a Python class. This class can then be registered as Sherpa model type using the add_model() function (found in Sherpa 4.2.2). The class attribute, name, can then be used as a model type within a standard Sherpa model expression.


To python users, instantiation of this class is obscured within the Sherpa model expression syntax. The add_model() function simply adds this class to Sherpa's function factory for creating model instances. The variable, mylogparlaw, in the local namespace is actually the Sherpa function factory. This variable's __getattr__() method has been overloaded to create instances of MyLogParLaw using the name of the attribute requested. Think of the following declarations as comparable.

lp1 = mylogparlaw.lp1

versus

lp1 = MyLogParLaw('lp1')

See the script below for more details

#!/usr/bin/env python

from sherpa.astro.ui import *
from sherpa.models import ArithmeticModel, Parameter
import numpy as np


class MyLogParaLaw(ArithmeticModel):
   """


               /  x \ - p[1] - p[2] * log10(x/p[0])
   f(x) = p[3] |----|
               \p[0]/

   """
   def __init__(self, name='mylogparlaw'):
       # p[0]
       self.ref = Parameter(name, 'ref', 1, alwaysfrozen=True)

       # p[1]
       self.c1 = Parameter(name, 'c1', 1)

       # p[2]
       self.c2 = Parameter(name, 'c2', 1)

       # p[3]
       self.ampl = Parameter(name, 'ampl', 1, 0)

       ArithmeticModel.__init__(self, name, (self.ref,self.c1,
                                             self.c2,self.ampl))


   def calc(self, p, xlo, xhi=None, *args, **kwargs):
       frac = (xlo/p[0])
       exponent = -p[1] - p[2] * np.log10(frac)
       return p[3] * np.power(frac, exponent)

if __name__ == '__main__':

   # register the class MyLogParLaw as a Sherpa model type
   add_model(MyLogParLaw)

   # Set the Sherpa model expression as the sum of MyLogParLaw
   # and Const1D
   set_model(mylogparlaw.lp1+const1d.c1)

   # Initialize the parameter values
   c1.c0 = 10
   lp1.ampl = 75
   lp1.c1 = 0.5
   lp1.c2 = 1.75
   lp1.ref = 2.1
   print lp1

   # Define a data space grid to be from [0.1, 10] in intervals
   # of 0.01
   dataspace1d(0.1,10,0.01, dstype=Data1D)

   # Evaluate the model over set of initialized parameters, add
   # poisson noise, and set as data space values.
   fake()

   # Use Chi2 with data variance as fit statistic
   set_stat('chi2datavar')

   # Set the statistical error of data space as 10% of faked
   # data values
   set_staterror(0.1, fractional=True)

   # Plot faked data space values with statistical errorbars
   plot_data()

The plot of faked data should look something like the image above.