## Thursday, June 10, 2010

### Vinay's Cheat Sheet

Vinay has been writing a Sherpa Cheat Sheet on AstroStat Slog. This is a good list of basic "how to" tutorials for Sherpa.

### 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.

versus

See the script below for more details

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

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.

Subscribe to:
Posts (Atom)