Thursday, June 10, 2010

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.

No comments:

Post a Comment