source: sasview/park-1.2.1/park/parameter.py @ a189869

# This program is public domain
"""
Parameters and parameter sets.

Parameter defines an individual parameter, and ParameterSet groups them
into a hierarchy.

Individual models need to provide a parameter set with the correct
properties, either by using park.ParameterSet in their model definition,
or by providing a wrapper which can translate assignment to parameter.value
into the appropriate change in the wrapped model.  See wrapper.py for
an example.
"""

__all__ = ['Parameter', 'ParameterSet']

import numpy
import expression        

class Pnormal(object):
    """
    Negative log likelihood function for a parameter from a Gaussian
    distribution.
    
    Given P(v;mu,sigma) = exp(-1/2 (mu-v)^2/sigma^2)/sqrt(2 pi sigma^2)
    then -log(P) = 1/2 (mu-v)^2/sigma^2 + 1/2 log(2*pi*sigma^2)
    
    Assuming that parameter P is selected from a normal distribution,
    then P.likelihood = Pnormal(mu,sigma)
    """
    def __init__(self, mean, std):
        self.mean = mean
        self.std = std
        self.const = math.log(2*math.pi*std**2)/2
    def __call__(self, value):
        return ((value-self.mean)/self.std)**2/2 + self.const
    
inf = numpy.inf
class Parameter(object):
    """
    A parameter is a box for communicating with the fitting service.
    Parameters can have a number of properties, 
    
    Parameters have a number of properties:
    
    name "string"
        name of the parameter within the parameter set.  
        
        The name is read only.  You can rename a parameter but only 
        in the context of the parameter set which contains it, using
        parameterset.rename(par,name).  This will change all expressions
        containing the named parameter.
    
    path
        dotted name of the parameter within the set of models.  The
        dotted name is automatically generated by the parameter set
        before expressions are parsed and evaluated.  There are
        some operations on parameter sets (such as renaming the
        layer containing a parameter) which will force an adjustment
        of all the underlying parameter names, as well as any
        expressions in which they are referenced.
    
    limits (low, high)
        hard limits on the range of allowed parameter values, dictated
        by physics rather than by knowledge of the particular system.
        For example, thickness parameters would have limits (0,inf)
        because negative thicknesses are unphysical.  These limits
        are enforced when setting range for the fit.

    units "string"
        units for the parameter.  This should be a string, but
        be parsable by whatever units package your application
        supports.

    tip "string"
        parameter description, suitable for use in a tooltip
    
    value double
        current value of the parameter, either fixed, fitted or computed
    
    range (low, high)
        range of expected values for the parameter in the model
    
    expression "string"
        expression for the parameter in the model.  This is a string
        containing a formula for updating the parameter value based
        on the values of other parameters in the system.  The expression
        is ignored if 'calculated' is False.

        Note: the list of symbols available to the expression evaluator
        defaults to the contents of the math module.  The caller will be
        able to override this within the fitting fitting class.
    
    status 'fixed'|'computed'|'fitted'
        the parameter type.  Choose 'fixed' if the values is to
        remain fixed throughout the fit, even if a range and an
        expression have been specified.  Choose 'computed' if the
        value of the parameter is to be computed each time the
        parameters are updated.  Choose 'fitted' if an optimization
        algorithm is supposed to search the parameter space.

    likelihood
        function to return the negative log likelihood of seeing a
        particular parameter value.  2*likelihood(value) will be added
        to the total cost function for the particular parameter set
        during the fit.  This will be on par with the probabilty
        of seeing the particular theory function given the observed
        datapoints when performing the fit (the residual term is
        closely related to the log likelihood of the normal distribution).
        
        Note: we are minimizing chi^2 = sum [ ((y-f(x;p))/dy)^2 ] rather 
        than  -log P = sum [ ((y-f(x;p))/dy)^2/2 + log(2 pi dy^2) ],
        where P is the probability of seeing f(x;p) given y,dy as the
        mean and standard deviation of a normal distribution.  Because
        chi^2_p = - 2 * log P_p + constant, the minima of p are the same
        for chi^2 and negative log likelihood.  However, to weight the
        likelihood properly when adding likelihood values to chisq, we
        need the extra factor of 2 mentioned above.  The usual statistical
        implications of normalized chi^2 will of course be suspect, both
        because the assumption of independence between the points in
        chi^2 (which definitely do not hold for the new 'point' p_k), and
        because of the additional 2 log(2 pi dp_k^2) constant, but given
        the uncertainty in the estimate of the distribution parameters,
        this is likely a minor point.
    """
    # Protect parameter name from user modification
    _name = "unknown"
    def _getname(self): return self._name
    name = property(_getname,doc="parameter name")
    
    # Do checking on the parameter range to make sure the model stays physical
    _range = (-inf,inf)
    def _setrange(self, r):
        if self.limits[0]<=r[0]<=r[1] <=self.limits[1]:
            self._range = r
        else:
            raise ValueError("invalid range %s for %s"%(r,self.name))
    def _getrange(self): return self._range
    range = property(_getrange,_setrange,doc="parameter range")

    path = ""
    value = 0.
    limits = (-inf,inf)
    expression = ""
    status = 'fixed' # fixed, computed or fitted
    likelihood = None
    units = ""
    tip = "Fitting parameter"
    deriv = False
    def __init__(self, name="unknown", **kw):
        self._name = name
        for k,v in kw.items():
            if hasattr(self,k): setattr(self,k,v)
            else: raise AttributeError("Unknown attribute %s"%k)
    def __str__(self): return self.path if self.path != '' else self.name
    def __repr__(self): return "Parameter('%s')"%self.name
    def summarize(self):
        """
        Return parameter range string.
        
        E.g.,  "       Gold .....|.... 5.2043 in [2,7]"
        """
        range = ['.']*10
        lo,hi = p.range
        portion = (p.value-lo)/(hi-lo)
        if portion < 0: portion = 0.
        elif portion >= 1: portion = 0.99999999
        bar = math.floor(portion*len(range))
        range[bar] = '|'
        range = "".join(range)
        return "%25s %s %g in [%g,%g]"  % (p.name,range,p.value,lo,hi)

    def isfitted(self): return self.status == 'fitted'
    def iscomputed(self): return self.status == 'computed'
    def isfixed(self): return self.status == 'fixed'
    def isrestrained(self): return self.likelihood is not None
    def isfeasible(self, value):
        """Return true if the value is in the range"""
        return self._range[0] <= value <= self._range[1]
    def setprefix(self, prefix):
        """
        Set the full path to the parameter as used in expressions involving
        the parameter name.
        """
        self.path = prefix+self.name
    def get(self):
        """
        Return the current value for a parameter.
        """
        return self.value
    def set(self, value):
        """
        Set a parameter to a value, a range or an expression.  If it is a value,
        the parameter will be fixed for the fit.  If it is a range, the value
        will be varying for the fit.  If it is an expression, the parameter will
        be calculated from the values of other parameters in the fit.

        Raises ValueError if the value could not be interpreted.
        """
        
        # Expression
        if isinstance(value,basestring):
            self.expression = value
            self.status = 'computed'
            return
        
        # Fixed value
        if numpy.isscalar(value):
            self.value = value
            self.status = 'fixed'
            return
        
        # Likelihood
        if hasattr(value,'__call__'):
            self.range = value.range
            self.likelihood = value
            self.status = 'fitted'
            return
        
        # Range
        try:
            lo,hi = numpy.asarray(value)
            if not numpy.isscalar(lo) or not numpy.isscalar(hi):
                raise Exception
            self.range = (lo,hi)
            self.status = 'fitted'
            return
        except:
            pass
        
        raise ValueError,\
            "parameter %s expects value, expression or range: %s"%(self.name,value)

class ParameterSet(list):
    """
    The set of parameters used to fit theory to data.
    
    ParameterSet forms a hierarchy of parameters.  The parameters themselves
    are referred to by the path through the hierarchy, usually::
    
         fitname.component.parameter
         
    Though more or fewer levels are permitted.  Parameters are assumed to
    have a unique label throughout the fit.  This is required so that
    expressions tying the results of one fit to another can uniquely
    reference a parameter.
    
    Attributes:
    
    name
        the name of the parameter set
    path
        the full dotted name of the parameter set
    context
        a dictionary providing additional context for evaluating parameters;
        Note that this namespace is shared with other theory functions, so
        populate it carefully.
    """
    # Protect parameter set name from user modification
    def _getname(self): return self._name
    def _setname(self, name): 
        raise NotImplementedError("parameter.name is protected; use fit.rename_parameter() to change it")
    name = property(_getname,doc="parameter name")
    path = ""
    def __init__(self, name="unknown", pars=[]):
        super(ParameterSet,self).__init__(pars)
        self._name = name
        self.context = {}

    def _byname(self, parts):
        """byname recursion function"""
        if len(parts) == 1: return self
        for p in self:
            if parts[1] == p.name:
                if len(pars) == 2:
                    return p
                elif isinstance(p, ParameterSet):
                    return p._byname(parts[1:])
        return None

    def byname(self, name):
        """Lookup parameter from dotted path"""
        parts = name.split('.')
        if parts[0] == self.name:
            p =  _byname(self, name.split('.'))
            if p: return p
        raise KeyError("parameter %s not in parameter set"%name)
    
    def __getitem__(self, idx):
        """Allow indexing by name or by number"""
        if isinstance(idx,basestring):
            for p in self:
                if p.name == idx: 
                    return p
            raise KeyError("%s is not in the parameter set"%idx)
        else:
            return super(ParameterSet,self).__getitem__(idx)

    def flatten(self):
        """
        Iterate over the elements in depth first order.
        """
        import park
        L = []
        for p in self:
            # Yuck! I really only want to try flattening if p is a
            # ParameterSet but it seems we have parameter.ParameterSet
            # and park.parameter.ParameterSet as separate entities,
            # depending on how things were included.  The solution is
            # probably to force absolute include paths always.
            try:
                L += p.flatten()
            except:
                L.append(p)
        return L
        """
        # Iterators are cute but too hard to use since you can
        # only use them in a [p for p in generator()] once.
        for p in self:
            if isinstance(p, ParameterSet):
                for subp in p.flatten(): 
                    yield subp
            else: 
                yield p
        """

    def _fixed(self):
        """
        Return the subset of the parameters which are fixed
        """
        return [p for p in self.flatten() if p.isfixed()]
    fixed = property(_fixed,doc=_fixed.__doc__)

    def _fitted(self):
        """
        Return the subset of the paramters which are varying
        """
        return [p for p in self.flatten() if p.isfitted()]
    fitted = property(_fitted,doc=_fitted.__doc__)

    def _computed(self):
        """
        Return the subset of the parameters which are calculated
        """
        return [p for p in self.flatten() if p.iscomputed()]
    computed = property(_computed,doc=_computed.__doc__)

    def _restrained(self):
        """
        Return the subset of the parameters which have a likelihood
        function associated with them.
        """
        return [p for p in self.flatten() if p.isrestrained()]
    restrained = property(_restrained,doc=_restrained.__doc__)

    def setprefix(self,prefix=None):
        """
        Fill in the full path name for all the parameters in the tree.
        
        Note: this function must be called from the root parameter set
        to build proper path names.
        
        This is required before converting parameter expressions into
        function calls.
        """
        if prefix == None:
            # We are called from root, so we don't have a path
            self.path = ""
            prefix = ""
        else:
            self.path = prefix+self.name
            prefix = self.path+'.'
        for p in self:
            #print "setting prefix for",p,prefix
            p.setprefix(prefix)

    def rename(self, par, name):
        """
        Rename the parameter to something new.
        Called from root of the parameter hierarchy, rename the particular
        parameter object to something else.

        This changes the internal name of the parameter, as well as all
        expressions in which it occurs.  If the parameter is actually
        a parameter set, then it renames all parameters in the set.
        """
        
        # Must run from root for self.setprefix and self.computed to work
        if self.path != "":  
            raise RuntimeError,"rename must be called from root parameter set"

        # Change the name of the node
        par._name = name
        
        # Determine which nodes (self and children) are affected
        if isinstance(par,ParameterSet):
            changeset = par.flatten()
        else:
            changeset = [par]
        
        #  Map the old names into the new names
        old = [p.path for p in changeset]
        self.setprefix() # Reset the path names of all parameters
        new = [p.path for p in changeset]
        mapping = dict(zip(old,new))

        # Perform the substitution into all of the expressions
        exprs = self.computed
        for p in exprs:
            p.expression = expression.substitute(p.expression, mapping)

    def gather_context(self):
        """
        Gather all additional symbols that can be used in expressions.
        
        For example, if reflectometry provides a volume fraction
        function volfrac(rho1,rho2,frac) to compute densities, then
        this function can be added as a context dictionary to the
        reflectometry parameter set.  Note that there is no guarantee
        which function will be used if the same function exists in
        two separate contexts.
        """
        context = {} # Create a new dictionary
        context.update(self.context)
        for p in self:
            if hasattr(p,'gather_context'):
                context.update(p.gather_context())
        return context


def test():
    # Check single parameter
    a = Parameter('a')
    assert a.name == 'a'
    a.value = 5
    assert a.value == 5
    # Check the setters
    a.set(7)
    assert a.value == 7 and a.status == 'fixed' and a.isfixed()
    a.set([3,5])
    assert a.value == 7 and a.range[0] == 3 and a.range[1]==5 and a.isfitted()
    a.set('3*M.b')
    assert a.iscomputed() and a.expression == '3*M.b'

    # Check limits
    a.limits = (0,inf)
    try: 
        a.range = (-1,1)
        raise Exception,"Failed to check range in limits"
    except ValueError: pass # Correct failure
    
    # Check that we can't change name directly
    try:
        a.name = 'Q'
        raise Exception,"Failed to protect name"
    except AttributeError: pass # Correct failure

    assert str(a) == 'a'  # Before setpath, just print name
    a.setprefix('M.')
    assert a.path == 'M.a'
    assert str(a) == 'M.a'  # After setpath, print path
    assert a.units == ''
    a.units = 'kg'
    assert a.units == 'kg'
    assert repr(a) == "Parameter('a')"

    # Check parameter set
    b,c = Parameter('b'),Parameter('c')
    M1 = ParameterSet('M1',[a,b,c])
    assert M1[0] is a
    assert M1['a'] is a
    a.set(5) # one value
    b.set([3,5])
    c.set('3*M1.a')
    assert M1.computed == [c]
    assert M1.fitted == [b]
    assert M1.fixed == [a]
    d = Parameter('d')
    M1.insert(0,d)
    assert M1.fixed == [d,a]
    e = Parameter('e')
    M1.append(e)
    assert M1.fixed == [d,a,e]

    a2,b2,c2 = [Parameter(s) for s in ('a','b','c')]
    a2.set(15)
    M2 = ParameterSet('M2',[a2,b2,c2])
    # Adjust parameter in set
    b2.set([3,5])
    assert M2.fitted == [b2]
    
    # Hierarchical parameter sets
    r = Parameter('r')
    root = ParameterSet('root',[M1,r,M2])
    assert root.fixed == [d,a,e,r,a2,c2]
    assert root.fitted == [b,b2]
    assert root.computed == [c]
    root.setprefix()
    assert a2.path == "M2.a"
    # Rename individual parameter
    root.rename(a,'a1')
    assert a.path == "M1.a1"
    assert c.expression == "3*M1.a1"
    # Rename parameter set
    root.rename(M1,'m1')
    assert c.path == "m1.c"
    assert c.expression == "3*m1.a1"

    # Integration test: parameter and expression working together.
    fn = expression.build_eval(root.flatten(), root.gather_context())
    #import dis; dis.dis(fn)
    fn()
    assert c.value == 3*a.value

    # Test context
    M2.context['plus2'] = lambda x: x+2
    c2.set('plus2(M2.a)')
    assert set(root.computed) == set([c,c2])
    fn = expression.build_eval(root.flatten(), root.gather_context())
    #print dis.dis(fn)
    fn()
    assert c2.value == a2.value+2

    # Multilevel hierarchy
    # Forming M3.a.x, M3.a.y, M3.b with M3.a.y = 2*M3.b+M3.a.x
    x = Parameter('x'); x.set(15)
    y = Parameter('y'); y.set('2*M3.b+M3.a.x')
    b = Parameter('b'); b.set(10)
    a = ParameterSet('a',[x,y])
    M3 = ParameterSet('M3',[a,b])
    root = ParameterSet('root',[M3])
    root.setprefix()
    fn = expression.build_eval(root.flatten(), root.gather_context())
    #import dis; dis.dis(fn)
    fn()
    #print "y-value:",y.value
    assert y.value == 35
    



if __name__ == "__main__": test()