# 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 math
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
[docs]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
[docs] def summarize(self):
"""
Return parameter range string.
E.g., " Gold .....|.... 5.2043 in [2,7]"
"""
range = ['.']*10
lo,hi = self.range
portion = (self.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]" % (self.name,range,self.value,lo,hi)
[docs] def isfitted(self): return self.status == 'fitted'
[docs] def iscomputed(self): return self.status == 'computed'
[docs] def isfixed(self): return self.status == 'fixed'
[docs] def isrestrained(self): return self.likelihood is not None
[docs] def isfeasible(self, value):
"""Return true if the value is in the range"""
return self._range[0] <= value <= self._range[1]
[docs] def setprefix(self, prefix):
"""
Set the full path to the parameter as used in expressions involving
the parameter name.
"""
self.path = prefix+self.name
[docs] def get(self):
"""
Return the current value for a parameter.
"""
return self.value
[docs] 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)
[docs]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(parts) == 2:
return p
elif isinstance(p, ParameterSet):
return p._byname(parts[1:])
else:
raise
return None
[docs] def byname(self, name):
"""Lookup parameter from dotted path"""
parts = name.split('.')
if parts[0] == self.name:
p = self._byname(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)
[docs] 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__)
[docs] 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)
[docs] 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)
[docs] 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()