source: sasview/park-1.2.1/park/model.py @ b11e127

# This program is public domain
"""
Define a park fitting model.

Usage
-----

The simplest sort of fitting model is something like the following::

    import numpy
    import park
    def G(x,mu,sigma):
        return numpy.exp(-0.5*(x-mu)**2/sigma**2)

    class Gauss(park.Model):
        parameters = ['center','width','scale']
        def eval(self,x):
            return self.scale * G(x,self.center,self.width)

It has a function which is evaluated at a series of x values and
a set of adjustable parameters controlling the shape of f(x).

You can check your module with something like the following::

    $ ipython -pylab

    from gauss import Gauss

    g = Gauss(center=5,width=1,scale=10)
    x = asarray([1,2,3,4,5])
    y = g(x)
    plot(x,y)

This should produce a plot of the Gaussian peak.

You will then want to try your model with some data.  Create a file
with some dummy data, such as gauss.dat::

    # x y
    4 2
    5 6
    6 7
    7 3

In order to match the model to data, you need to define a fitness
object.  This shows the difference between the model and the data,
which you can then plot or sum to create a weighted chisq value::

    f = park.Fitness(g, 'gauss.dat')
    plot(f.data.fit_x, f.residuals())

The data file can have up to four columns, either x,y or x,y,dy
or x,dx,y,dy where x,y are the measurement points and values,
dx is the instrument resolution in x and dy is the uncertainty
in the measurement value y.  You can pass any keyword arguments
to data.load that are accepted by numpy.loadtxt.  For example,
you can reorder columns or skip rows.  Additionally, you can modify
data attributes directly x,y,dx,dy and calc_x.  See help on park.Data1D
for details.

Once you have your model debugged you can use it in a fit::

    g = Gauss(center=[3,5],scale=7.2,width=[1,3])
    result = park.fit((g, 'gauss.dat'))
    result.plot()

In this example, center and width are allowed to vary but scale is fixed.

Existing models can be readily adapted to Park::

    class Gauss(object):
        "Existing model"
        center,width,scale = 0,1,1
        def __init__(self,**kw):
            for k,v in kw.items(): setattr(self,k,v)
        def eval(self,x):
            return self.scale *G(x,self.center,self.width)

    class GaussAdaptor(Gauss,Model):
        "PARK adaptor"
        parameters = ['center','width','scale']
        def __init__(self,*args,**kw):
            Model.__init__(self,*args)
            Gauss.__init__(self,**kw)

    g = GaussAdaptor(center=[3,5],scale=7.2,width=[1,3])
    result = park.fit((g, 'gauss.dat'))
    result.plot()

Models can become much more complex than the ones described above,
including multilevel models where fitting parameters can be added
and removed dynamically.

In many cases the park optimizer will need an adaptor for pre-existing
models.  The adaptor above relies on python properties to translate
model.par access into model._par.get() and model._par.set() where _par
is the internal name for par.  This technique works for simple static
models, but will not work well for sophisticated models which have,
for example, a dynamic parameter set where the model parameters cannot
be set as properties.  A solution to this problem is to subclass the
park.Parameter and override the value attribute as a property.

Once models are defined they can be used in a variety of contexts, such
as simultaneous fitting with constraints between the parameters.  With
some care in defining the model, computationally intensive fits can
be distributed across multiple processors.  We provide a simple user
interface for interacting with the model parameters and managing fits.
This can be extended with model specialized model editors which format
the parameters in a sensible way for the model, or allow direct manipulation
of the model structure.  The underlying fitting engine can also be
used directly from your own user interface.

"""

__all__ = ['Model']

import park
from park.parameter import Parameter, ParameterSet
from copy import copy, deepcopy


class ParameterProperty(object):
    """
    Attribute accessor for a named parameter.

    Objects of class ParameterProperty act similarly to normal property
    objects in that assignment to object.attr triggers the __set__
    method of ParameterProperty and queries of the value object.attr
    triggers the __get__ method.  These methods look up the actual
    parameter in the dictionary model.parameterset, delegating the
    the set/get methods of the underlying parameter object.

    For example::

        model.P1 = 5
        print model.P1

    is equivalent to::

        model.parameterset['P1'].set(5)
        print model.parameterset['P1'].get()

    To enable this behaviour, the model developer must assign a
    ParameterProperty object to a class attribute of the model.  It
    must be a class attribute.  Properties assigned as an instance
    attribute in the __init__ of the class will not be recognized
    as properties.

    An example model will look something like the following::

        class MyModel:
            # A property must be assigned as a class attribute!
            P1 = ParameterProperty('P1')
            P2 = ParameterProperty('P2')
            def __init__(self, P1=None, P2=None):
                parP1 = Parameter('P1')
                if P1 is not None: parP1.set(P1)
                parP2 = Parameter('P2')
                if P2 is not None: parP2.set(P2)
                self.parameterset = { 'P1': parP1, 'P2': parP2 }

    This work is done implicitly by MetaModel, and by subclassing
    the class Model, the model developer does not ever need to
    use ParameterProperty.
    """
    def __init__(self,name,**kw):
        self.name = name
    def __get__(self,instance,cls):
        return instance.parameterset[self.name].get()
    def __set__(self,instance,value):
        instance.parameterset[self.name].set(value)


class MetaModel(type):
    """
    Interpret parameters.

    This is a meta class, and the usual meta class rules apply.  That is,
    the Model base class should be defined like::

        class Model(object):
            __metaclass__ = MetaModel
            ...

    The MetaModel.__new__ method is called whenever a new Model
    class is created.  The name of the model, its superclasses and
    its attributes are passed to MetaModel.__new__ which creates the
    actual class.  It is not called for new instances of the model.

    MetaModel is designed to simplify the definition of parameters for
    the model.  When processing the class, MetaModel defines
    parameters which is a list of names of all the parameters in the
    model, and parameterset, which is a dictionary of the actual
    parameters and any parameter sets in the model.

    Parameters can be defined using a list of names in the parameter
    attribute, or by defining the individual parameters as separate
    attributes of class Parameter.

    For example, the following defines parameters P1, P2, and P3::

        class MyModel(Model):
            parameters = ['P1','P2']
            P2 = Parameter(range=[0,inf])
            P3 = Parameter()

    For each parameter, MetaModel will define a parameter accessor,
    and add the parameter definition to the parameter set. The accessor
    delegates query and assignment to the Parameter get/set methods. The
    attributes of the particular parameter instance can be
    adjusted using model.parameterset['name'].attribute = value.
    """
    def __new__(cls, name, bases, vars):
        #print 'calling model meta',vars

        # Find all the parameters for the model.  The listed parameters
        # are defined using::
        #    parameters = ['P1', 'P2', ...]
        # The remaining parameters are defined individually::
        #    P3 = Parameter()
        #    P4 = Parameter()
        # The undefined parameters are those that are listed but not defined.
        listed = vars.get('parameters',[])
        defined = [k for k,v in vars.items() if isinstance(v,Parameter)]
        undefined = [k for k in listed if k not in defined]
        #print listed,defined,undefined

        # Define a getter/setter for every parameter so that the user
        # can say model.name to access parameter name.
        #
        # Create a parameter object for every undefined parameter.
        # Check if the base class defines a default value, and use
        # that for the initial value.  We don't want to do this for
        # parameters explicitly defined since the user may have
        # given them a default value already.
        #
        # For predefined parameters we must set the name explicitly.
        # This saves the user from having to use, e.g.::
        #     P1 = Parameter('P1')
        pars = []
        for p in undefined:
            # Check if the attribute value is initialized in a base class
            for b in bases:
                if hasattr(b,p):
                    value = getattr(b,p)
                    break
            else:
                value = 0.
            #print "looking up value in base as",value
            pars.append(Parameter(name=p,value=value))
            vars[p] = ParameterProperty(p)
        for p in defined:
            # Make sure parameter name is correct.  Note that we are using
            # _name rather than .name to access the name attribute since
            # name is a read-only parameter.
            vars[p]._name = p
            pars.append(vars[p])
            vars[p] = ParameterProperty(p)

        # Construct the class attributes 'parameters' and 'parameterset'.
        # Listed parameters are given first, with unlisted parameters
        # following alphabetically.  For hierarchical parameter sets,
        # we also need to include the defined sets.
        unlisted = list(set(defined+undefined) - set(listed))
        unlisted.sort()
        parameters = listed + unlisted
        parsets = [k for k,v in vars.items() if isinstance(v,ParameterSet)]
        vars['parameters'] = parameters
        vars['parameterset'] = ParameterSet(pars=pars+parsets)
        #print 'pars',vars['parameters']
        #print 'parset',vars['parameterset']

        # Return a new specialized instance of the class with parameters
        # and parameterset made explicit.
        return type.__new__(cls, name, bases, vars)

class Model(object):
    """
    Model definition.

    The model manages attribute access to the fitting parameters and
    also manages the dataset.

    derivatives ['p1','p2',...]
        List of parameters for which the model can calculate
        derivatives.  The derivs
        The model function can compute the derivative with respect
        to this parameter.  The function model.derivs(x,[p1,p2,...])
        will return (f(x),df/dp1(x), ...).  The parameters and their
        order are determined by the fitting engine.

        Note: This is a property of the model, not the fit.
        Numerical derivatives will be used if the parameter is
        used in an expression or if no analytic derivative is
        available for the parameter.  Automatic differentiation
        on parameter expressions is possible, but beyond the scope
        of this project.

    eval(x)
        Evaluate the model at x.  This must be defined by the subclass.

    eval_deriv(x,pars=[])
        Evaluate the model and the derivatives at x.  This must be
        defined by the subclass.

    parameters
        The names of the model parameters.  If this is not provided, then
        the model will search the subclass for park.Parameter attributes
        and construct the list of names from that.  Any parameters in the
        list not already defined as park.Parameter attributes will be
        defined as parameters with a default of 0.

    parameterset
        The set of parameters defined by the model.  These are the
        parameters themselves, gathered into a park.ParameterSet.

    The minimum fittng model if you choose not to subclass park.Model
    requires parameterset and a residuals() method.
    """
    __metaclass__ = MetaModel
    derivatives = []
    def __init__(self,*args,**kw):
        """
        Initialize the model.

        Model('name',P1=value, P2=Value, ...)

        When overriding __init__ in the subclass be sure to call
        Model.__init__(self, *args, **kw).  This makes a private
        copy of the parameterset for the model and initializes
        any parameters set using keyword arguments.
        """
        #print 'calling model init on',id(self)
        # To avoid trashing the Model.parname = Parameter() template
        # when we create an instance and accessor for the parameter
        # we need to make sure that we are using our own copy of the
        # model dictionary stored in vars
        if len(args)>1: raise TypeError("expect name as model parameter")
        name = args[0] if len(args)>=1 else 'unknown'
        self.parameterset = deepcopy(self.parameterset)
        for k,v in kw.items():
            self.parameterset[k].set(v)
        self.parameterset._name = name
        #print 'done',id(self)

    def __call__(self, x, pars=[]):
        return self.eval(x)

    def eval(self, x):
        """
        Evaluate the model at x.

        This method needs to be specialized in the model to evaluate the
        model function.  Alternatively, the model can implement is own
        version of residuals which calculates the residuals directly
        instead of calling eval.
        """
        raise NotImplementedError('Model needs to implement eval')

    def eval_derivs(self, x, pars=[]):
        """
        Evaluate the model and derivatives wrt pars at x.

        pars is a list of the names of the parameters for which derivatives
        are desired.

        This method needs to be specialized in the model to evaluate the
        model function.  Alternatively, the model can implement is own
        version of residuals which calculates the residuals directly
        instead of calling eval.
        """
        raise NotImplementedError('Model needs to implement eval_derivs')

    def set(self, **kw):
        """
        Set the initial value for a set of parameters.

        E.g., model.set(width=3,center=5)
        """
        # This is a convenience funciton for the user.
        #
        for k,v in kw.items():
            self.parameterset[k].set(v)


def add_parameter(model, name, **kw):
    """
    Add a parameter to an existing park model.

    Note: this may not work if it is operating on a BaseModel.
    """
    setattr(model.__class__, name, park.model.ParameterProperty(name))
    model.parameterset.append(park.Parameter(name=name, **kw))


def test():
    # Gauss theory function
    import numpy
    eps,inf = numpy.finfo('d').eps,numpy.inf
    def G(x,mu,sigma):
        mu,sigma = numpy.asarray(mu),numpy.asarray(sigma)
        return numpy.exp(-0.5*(x-mu)**2/sigma**2)

    # === Minimal model ===
    # just list the fitting parameters and define the function.
    class Gauss(Model):
        parameters = ['center','width','scale']
        def eval(self,x):
            return self.scale * G(x,self.center,self.width)

    # Create a couple of models and make sure they don't conflict
    g1 = Gauss(center=5,width=1,scale=2)
    assert g1.center == 5
    g2 = Gauss(center=6,width=1)
    assert g1.center == 5
    assert g2.center == 6
    assert g1(5) == 2
    assert g2(6) == 0

    # === explore parameters ===
    # dedicated model using park parameters directly, and defining derivatives
    class Gauss(Model):
        center = Parameter(value=0,
                           tip='Peak center')
        width = Parameter(value=1,
                          limits=(eps,inf),
                          tip='Peak width (1-sigma)')
        scale = Parameter(value=1,
                          limits=(0,inf),
                          tip='Peak scale (>0)')
        def eval(self,x):
            return self.scale * G(x,self.center,self.width)

        # Derivatives
        derivatives = ['center','width','scale']
        def dscale(self, g, x):
            return g/self.scale
        def dcenter(self, g, x):
            return 2*(x-self.center)/self.width*g
        def dwidth(self, g, x):
            return 2*(x-self.center)/self.width**3*g
        dmap = dict(scale=dscale,center=dcenter,width=dwidth)
        def eval_derivs(self,x,pars=[]):
            """
            Calculate function value and derivatives wrt to the parameters.

            pars is a list of parameter names, possibly consisting of any
            parameter with deriv=True.
            """
            g = self.eval(x)
            dg = [self.dmap[p](self,g,x) for p in pars]
            return [g]+dg

    g1 = Gauss(center=5)
    g1.parameterset['center'].tip = 'This is the center'
    assert g1.center == 5
    g2 = Gauss(center=6)
    assert g1.center == 5
    assert g2.center == 6
    assert g1(5) == 1
    assert g2(6) == 1

    # ====== Test wrapper =======
    # delegate to existing model via inheritence
    class Gauss(object):
        """Pre-existing model"""
        center,width,scale = 0,1,1
        def __init__(self,**kw):
            #print "calling BaseGauss init on",id(self)
            for k,v in kw.items(): setattr(self,k,v)
            #print "done",id(self)
        def eval(self,x):
            return self.scale *G(x,self.center,self.width)

    class GaussAdaptor(Gauss,Model):
        """PARK wrapper"""
        parameters = ['center','width','scale']
        def __init__(self,*args,**kw):
            #print "calling Gauss init on",id(self)
            Model.__init__(self,*args)
            Gauss.__init__(self,**kw)
            #print "done",id(self)

    g1 = GaussAdaptor(center=5)
    g2 = GaussAdaptor(center=6)
    g3 = GaussAdaptor()
    assert g1.center == 5
    assert g2.center == 6
    assert g3.scale == 1
    assert g1(5) == 1
    assert g2(6) == 1

    g4 = GaussAdaptor('M3',center=6)
    assert g4.parameterset.name == 'M3'

    # dedicated multilevel model using park parameters directly
    class MultiGauss(Model):
        def add(self,name,model):
            pass

    # wrapped model using park parameters indirectly
    class BaseMultiGauss(object):
        def __init__(self):
            self.models = []
        def add(self,**kw):
            self.models.append(BaseGauss(**kw))
    class WrapMultiGauss(BaseMultiGauss,Model):
        def __init__(self):
            pass

if __name__ == "__main__": test()