# 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 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 = 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) 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(parts) == 2: return p elif isinstance(p, ParameterSet): return p._byname(parts[1:]) else: raise return None 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) 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()