Index: park-1.2.1/park/model.py
===================================================================
--- park-1.2.1/park/model.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/model.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,511 @@
+# 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()
