park.modelling.parameter

1 # This program is public domain 2 """ 3 Parameters and parameter sets. 4 5 Parameter defines an individual parameter, and ParameterSet groups them 6 into a hierarchy. 7 8 Individual models need to provide a parameter set with the correct 9 properties, either by using park.ParameterSet in their model definition, 10 or by providing a wrapper which can translate assignment to parameter.value 11 into the appropriate change in the wrapped model. See wrapper.py for 12 an example. 13 """ 14 15 __all__ = ['Parameter', 'ParameterSet'] 16 17 import numpy 18 import expression 19

20 -class Pnormal(object):

21 """ 22 Negative log likelihood function for a parameter from a Gaussian 23 distribution. 24 25 Given P(v;mu,sigma) = exp(-1/2 (mu-v)^2/sigma^2)/sqrt(2 pi sigma^2) 26 then -log(P) = 1/2 (mu-v)^2/sigma^2 + 1/2 log(2*pi*sigma^2) 27 28 Assuming that parameter P is selected from a normal distribution, 29 then P.likelihood = Pnormal(mu,sigma) 30 """

31 - def __init__(self, mean, std):

32 self.mean = mean 33 self.std = std 34 self.const = math.log(2*math.pi*std**2)/2

35 - def __call__(self, value):

36 return ((value-self.mean)/self.std)**2/2 + self.const

37 38 inf = numpy.inf

39 -class Parameter(object):

40 """ 41 A parameter is a box for communicating with the fitting service. 42 Parameters can have a number of properties, 43 44 Parameters have a number of properties: 45 46 name "string" 47 name of the parameter within the parameter set. 48 49 The name is read only. You can rename a parameter but only 50 in the context of the parameter set which contains it, using 51 parameterset.rename(par,name). This will change all expressions 52 containing the named parameter. 53 54 path 55 dotted name of the parameter within the set of models. The 56 dotted name is automatically generated by the parameter set 57 before expressions are parsed and evaluated. There are 58 some operations on parameter sets (such as renaming the 59 layer containing a parameter) which will force an adjustment 60 of all the underlying parameter names, as well as any 61 expressions in which they are referenced. 62 63 limits (low, high) 64 hard limits on the range of allowed parameter values, dictated 65 by physics rather than by knowledge of the particular system. 66 For example, thickness parameters would have limits (0,inf) 67 because negative thicknesses are unphysical. These limits 68 are enforced when setting range for the fit. 69 70 units "string" 71 units for the parameter. This should be a string, but 72 be parsable by whatever units package your application 73 supports. 74 75 tip "string" 76 parameter description, suitable for use in a tooltip 77 78 value double 79 current value of the parameter, either fixed, fitted or computed 80 81 range (low, high) 82 range of expected values for the parameter in the model 83 84 expression "string" 85 expression for the parameter in the model. This is a string 86 containing a formula for updating the parameter value based 87 on the values of other parameters in the system. The expression 88 is ignored if 'calculated' is False. 89 90 Note: the list of symbols available to the expression evaluator 91 defaults to the contents of the math module. The caller will be 92 able to override this within the fitting fitting class. 93 94 status 'fixed'|'computed'|'fitted' 95 the parameter type. Choose 'fixed' if the values is to 96 remain fixed throughout the fit, even if a range and an 97 expression have been specified. Choose 'computed' if the 98 value of the parameter is to be computed each time the 99 parameters are updated. Choose 'fitted' if an optimization 100 algorithm is supposed to search the parameter space. 101 102 likelihood 103 function to return the negative log likelihood of seeing a 104 particular parameter value. 2*likelihood(value) will be added 105 to the total cost function for the particular parameter set 106 during the fit. This will be on par with the probabilty 107 of seeing the particular theory function given the observed 108 datapoints when performing the fit (the residual term is 109 closely related to the log likelihood of the normal distribution). 110 111 Note: we are minimizing chi^2 = sum [ ((y-f(x;p))/dy)^2 ] rather 112 than -log P = sum [ ((y-f(x;p))/dy)^2/2 + log(2 pi dy^2) ], 113 where P is the probability of seeing f(x;p) given y,dy as the 114 mean and standard deviation of a normal distribution. Because 115 chi^2_p = - 2 * log P_p + constant, the minima of p are the same 116 for chi^2 and negative log likelihood. However, to weight the 117 likelihood properly when adding likelihood values to chisq, we 118 need the extra factor of 2 mentioned above. The usual statistical 119 implications of normalized chi^2 will of course be suspect, both 120 because the assumption of independence between the points in 121 chi^2 (which definitely do not hold for the new 'point' p_k), and 122 because of the additional 2 log(2 pi dp_k^2) constant, but given 123 the uncertainty in the estimate of the distribution parameters, 124 this is likely a minor point. 125 """ 126 # Protect parameter name from user modification 127 _name = "unknown"

128 - def _getname(self): return self._name

129 name = property(_getname,doc="parameter name") 130 131 # Do checking on the parameter range to make sure the model stays physical 132 _range = (-inf,inf)

133 - def _setrange(self, r):

134 if self.limits[0]<=r[0]<=r[1] <=self.limits[1]: 135 self._range = r 136 else: 137 raise ValueError("invalid range %s for %s"%(r,self.name))

138 - def _getrange(self): return self._range

139 range = property(_getrange,_setrange,doc="parameter range") 140 141 path = "" 142 value = 0. 143 limits = (-inf,inf) 144 expression = "" 145 status = 'fixed' # fixed, computed or fitted 146 likelihood = None 147 units = "" 148 tip = "Fitting parameter" 149 deriv = False

150 - def __init__(self, name="unknown", **kw):

151 self._name = name 152 for k,v in kw.items(): 153 if hasattr(self,k): setattr(self,k,v) 154 else: raise AttributeError("Unknown attribute %s"%k)

155 - def __str__(self): return self.path if self.path != '' else self.name

156 - def __repr__(self): return "Parameter('%s')"%self.name

157 - def summarize(self):

158 """ 159 Return parameter range string. 160 161 E.g., " Gold .....|.... 5.2043 in [2,7]" 162 """ 163 range = ['.']*10 164 lo,hi = p.range 165 portion = (p.value-lo)/(hi-lo) 166 if portion < 0: portion = 0. 167 elif portion >= 1: portion = 0.99999999 168 bar = math.floor(portion*len(range)) 169 range[bar] = '|' 170 range = "".join(range) 171 return "%25s %s %g in [%g,%g]" % (p.name,range,p.value,lo,hi)

172

173 - def isfitted(self): return self.status == 'fitted'

174 - def iscomputed(self): return self.status == 'computed'

175 - def isfixed(self): return self.status == 'fixed'

176 - def isrestrained(self): return self.likelihood is not None

177 - def isfeasible(self, value):

178 """Return true if the value is in the range""" 179 return self._range[0] <= value <= self._range[1]

180 - def setprefix(self, prefix):

181 """ 182 Set the full path to the parameter as used in expressions involving 183 the parameter name. 184 """ 185 self.path = prefix+self.name

186 - def get(self):

187 """ 188 Return the current value for a parameter. 189 """ 190 return self.value

191 - def set(self, value):

192 """ 193 Set a parameter to a value, a range or an expression. If it is a value, 194 the parameter will be fixed for the fit. If it is a range, the value 195 will be varying for the fit. If it is an expression, the parameter will 196 be calculated from the values of other parameters in the fit. 197 198 Raises ValueError if the value could not be interpreted. 199 """ 200 201 # Expression 202 if isinstance(value,basestring): 203 self.expression = value 204 self.status = 'computed' 205 return 206 207 # Fixed value 208 if numpy.isscalar(value): 209 self.value = value 210 self.status = 'fixed' 211 return 212 213 # Likelihood 214 if hasattr(value,'__call__'): 215 self.range = value.range 216 self.likelihood = value 217 self.status = 'fitted' 218 return 219 220 # Range 221 try: 222 lo,hi = numpy.asarray(value) 223 if not numpy.isscalar(lo) or not numpy.isscalar(hi): 224 raise Exception 225 self.range = (lo,hi) 226 self.status = 'fitted' 227 return 228 except: 229 pass 230 231 raise ValueError,\ 232 "parameter %s expects value, expression or range: %s"%(self.name,value)

233

234 -class ParameterSet(object):

235 """ 236 The set of parameters used to fit theory to data. 237 238 ParameterSet forms a hierarchy of parameters. The parameters themselves 239 are referred to by the path through the hierarchy, usually:: 240 241 fitname.component.parameter 242 243 Though more or fewer levels are permitted. Parameters are assumed to 244 have a unique label throughout the fit. This is required so that 245 expressions tying the results of one fit to another can uniquely 246 reference a parameter. 247 248 Attributes: 249 250 name 251 the name of the parameter set 252 path 253 the full dotted name of the parameter set 254 context 255 a dictionary providing additional context for evaluating parameters; 256 Note that this namespace is shared with other theory functions, so 257 populate it carefully. 258 259 Operations: 260 261 insert 262 add a new item to the set 263 append 264 add a new item to the end of the set 265 move 266 reorder the parameters in the set 267 """ 268 # Protect parameter set name from user modification

269 - def _getname(self): return self._name

270 - def _setname(self, name):

271 raise NotImplementedError("name is protected; use parameter.rename() to change it")

272 name = property(_getname,doc="parameter name") 273 path = ""

274 - def __init__(self, name="unknown", pars=[]):

275 #super(ParameterSet,self).__init__() 276 self._dict = {} 277 self._list = [] 278 for p in pars: 279 self._dict[p.name] = p 280 self._list.append(p) 281 self._name = name 282 self.context = {}

283

284 - def __len__(self):

285 return len(self._list)

286

287 - def _byname(self, name):

288 for i,p in enumerate(self._list): 289 if p.name == name: return i 290 raise KeyError(name+" not found in parameter set")

291

292 - def __getattr__(self, idx):

293 """Allow dotted-name parameter access""" 294 if idx in self._dict: 295 return self._dict[idx] 296 else: 297 raise AttributeError(idx)

298

299 - def __getitem__(self, idx):

300 """Allow indexing by name or by number""" 301 if isinstance(idx,int): 302 return self._list[idx] 303 else: 304 return self._dict[idx]

305

306 - def __delitem__(self, idx):

307 """Allow delete by name or by number""" 308 if isinstance(idx,int): 309 del self._dict[self._list[idx].name] 310 del self._list[idx] 311 else: 312 del self._list[self._byname(idx)] 313 del self._dict[idx]

314

315 - def __setitem__(self, idx, p):

316 """Allow replacement by number""" 317 if isinstance(idx,int): 318 # Replacing a specific position; make sure the corresponding 319 # item is removed from the dictionary 320 del self._dict[self._list[idx].name] 321 self._list[idx] = p 322 self._dict[p.name] = p 323 else: 324 raise KeyError("cannot insert parameters by name")

325

326 - def __iter__(self):

327 """ 328 Force iteration order 329 """ 330 for p in self._list: yield p

331

332 - def append(self, p):

333 """ 334 Add a parameter to the dictionary 335 """ 336 self._list.append(p) 337 self._dict[p.name] = p

338

339 - def insert(self, idx, p):

340 """ 341 Add a parameter to the dictionary 342 """ 343 self._list.insert(idx, p) 344 self._dict[p.name] = p

345

346 - def move(self, idx1, idx2):

347 """ 348 Move a parameter to a new position in the dictionary 349 """ 350 self._list.insert(idx2,self._list[idx1]) 351 if idx1 < idx2: 352 del self._list[idx1] 353 else: 354 del self._list[idx2+1]

355

356 -def fixed(self):

357 """ 358 Return the subset of the parameters which are fixed 359 """ 360 return [p for p in flatten(self) if p.isfixed()]

361

362 -def fitted(self):

363 """ 364 Return the subset of the paramters which are varying 365 """ 366 return [p for p in flatten(self) if p.isfitted()]

367

368 -def computed(self):

369 """ 370 Return the subset of the parameters which are calculated 371 """ 372 return [p for p in flatten(self) if p.iscomputed()]

373

374 -def restrained(self):

375 """ 376 Return the subset of the parameters which have a likelihood 377 function associated with them. 378 """ 379 return [p for p in flatten(self) if p.isrestrained()]

380

381 -def _lookup(self, parts):

382 """lookup helper""" 383 if len(parts) == 1: return self 384 for p in self: 385 if parts[1] == p.name: 386 if len(parts) == 2: 387 return p 388 elif isinstance(p, ParameterSet): 389 return _lookup(p, parts[1:]) 390 return None

391

392 -def find(self, name):

393 """Lookup parameter from dotted path""" 394 parts = name.split('.') 395 if parts[0] == self.name: 396 p = _lookup(self, name.split('.')) 397 if p: return p 398 raise KeyError("parameter %s not in parameter set"%name)

399

400 -def rename(self, target, name):

401 """ 402 Rename the parameter to something new. 403 Called from root of the parameter hierarchy, rename the particular 404 parameter object to something else. 405 406 This changes the internal name of the parameter, as well as all 407 expressions in which it occurs. If the parameter is actually 408 a parameter set, then it renames all parameters in the set. 409 """ 410 411 # Must run from root for self.setprefix and self.computed to work 412 if self.path != "": 413 raise ValueError,"rename must be called from root parameter set" 414 415 # Change the name of the node 416 target._name = name 417 418 # Determine which nodes (self and children) are affected 419 if isinstance(target,ParameterSet): 420 changeset = flatten(target) 421 else: 422 changeset = [target] 423 424 # Map the old names into the new names 425 old = [p.path for p in changeset] 426 setprefix(self) # Reset the path names of all parameters 427 new = [p.path for p in changeset] 428 mapping = dict(zip(old,new)) 429 430 # Perform the substitution into all of the expressions 431 exprs = computed(self) 432 for p in exprs: 433 p.expression = expression.substitute(p.expression, mapping)

434

435 -def flatten(self):

436 """ 437 Iterate over the elements in depth first order. 438 """ 439 L = [] 440 for p in self: 441 if isinstance(p, ParameterSet): 442 L += flatten(p) 443 else: 444 L.append(p) 445 return L 446 """ 447 # Iterators are cute but too hard to use since you can 448 # only use them in a [p for p in generator()] once. 449 for p in self: 450 if isinstance(p, ParameterSet): 451 for subp in p.flatten(): 452 yield subp 453 else: 454 yield p 455 """

456

457 -def gather_context(self):

458 """ 459 Gather all additional symbols that can be used in expressions. 460 461 For example, if reflectometry provides a volume fraction 462 function volfrac(rho1,rho2,frac) to compute densities, then 463 this function can be added as a context dictionary to the 464 reflectometry parameter set. Note that there is no guarantee 465 which function will be used if the same function exists in 466 two separate contexts. 467 """ 468 context = {} 469 context.update(self.context) 470 for p in self: 471 if isinstance(p, ParameterSet): 472 context.update(gather_context(p)) 473 return context

474

475 -def setprefix(self,prefix=None):

476 """ 477 Fill in the full path name for all the parameters in the tree. 478 479 Note: this function must be called from the root parameter set 480 to build proper path names. 481 482 This is required before converting parameter expressions into 483 function calls. 484 """ 485 if not isinstance(self, ParameterSet): 486 self.setprefix(prefix) 487 return 488 if prefix == None: 489 # We are called from root, so we don't have a path 490 self.path = "" 491 prefix = "" 492 else: 493 self.path = prefix+self.name 494 prefix = self.path+'.' 495 for p in self: 496 if isinstance(p, ParameterSet): 497 #print "setting prefix for",p,prefix 498 setprefix(p,prefix) 499 else: 500 p.setprefix(prefix)

501 502 503 504

505 -def test():

506 # Check single parameter 507 a = Parameter('a') 508 assert a.name == 'a' 509 a.value = 5 510 assert a.value == 5 511 # Check the setters 512 a.set(7) 513 assert a.value == 7 and a.status == 'fixed' and a.isfixed() 514 a.set([3,5]) 515 assert a.value == 7 and a.range[0] == 3 and a.range[1]==5 and a.isfitted() 516 a.set('3*M.b') 517 assert a.iscomputed() and a.expression == '3*M.b' 518 519 # Check limits 520 a.limits = (0,inf) 521 try: 522 a.range = (-1,1) 523 raise Exception,"Failed to check range in limits" 524 except ValueError: pass # Correct failure 525 526 # Check that we can't change name directly 527 try: 528 a.name = 'Q' 529 raise Exception,"Failed to protect name" 530 except AttributeError: pass # Correct failure 531 532 assert str(a) == 'a' # Before setpath, just print name 533 setprefix(a, 'M.') 534 assert a.path == 'M.a' 535 assert str(a) == 'M.a' # After setpath, print path 536 assert a.units == '' 537 a.units = 'kg' 538 assert a.units == 'kg' 539 assert repr(a) == "Parameter('a')" 540 541 # Check parameter set 542 b,c = Parameter('b'),Parameter('c') 543 M1 = ParameterSet('M1',[a,b,c]) 544 #assert M1[0] is a # No longer ordered 545 assert M1['a'] is a 546 a.set(5) # one value 547 b.set([3,5]) 548 c.set('3*M1.a') 549 assert computed(M1) == [c] 550 assert fitted(M1) == [b] 551 assert fixed(M1) == [a] 552 d = Parameter('d') 553 M1.insert(0,d) 554 assert set(fixed(M1)) == set([d,a]) 555 e = Parameter('e') 556 M1.append(e) 557 assert set(fixed(M1)) == set([d,a,e]) 558 559 a2,b2,c2 = [Parameter(s) for s in ('a','b','c')] 560 a2.set(15) 561 M2 = ParameterSet('M2',[a2,b2,c2]) 562 # Adjust parameter in set 563 b2.set([3,5]) 564 assert fitted(M2) == [b2] 565 566 # Hierarchical parameter sets 567 r = Parameter('r') 568 root = ParameterSet('root',[M1,r,M2]) 569 assert root.M1.a == a 570 assert find(root,'root.r') is r 571 assert fixed(root) == [d,a,e,r,a2,c2] 572 assert fitted(root) == [b,b2] 573 assert computed(root) == [c] 574 setprefix(root) 575 assert a2.path == "M2.a" 576 # Rename individual parameter 577 rename(root, a,'a1') 578 assert a.path == "M1.a1" 579 assert c.expression == "3*M1.a1" 580 # Rename parameter set 581 rename(root, M1, 'm1') 582 assert c.path == "m1.c" 583 assert c.expression == "3*m1.a1" 584 585 # Integration test: parameter and expression working together. 586 fn = expression.build_eval(flatten(root), gather_context(root)) 587 #import dis; dis.dis(fn) 588 fn() 589 assert c.value == 3*a.value 590 591 # Test context 592 M2.context['plus2'] = lambda x: x+2 593 c2.set('plus2(M2.a)') 594 assert computed(root) == [c,c2] 595 fn = expression.build_eval(flatten(root), gather_context(root)) 596 #print dis.dis(fn) 597 fn() 598 assert c2.value == a2.value+2 599 600 # Multilevel hierarchy 601 # Forming M3.a.x, M3.a.y, M3.b with M3.a.y = 2*M3.b+M3.a.x 602 x = Parameter('x'); x.set(15) 603 y = Parameter('y'); y.set('2*M3.b+M3.a.x') 604 b = Parameter('b'); b.set(10) 605 a = ParameterSet('a',[x,y]) 606 M3 = ParameterSet('M3',[a,b]) 607 root = ParameterSet('root',[M3]) 608 setprefix(root) 609 fn = expression.build_eval(flatten(root), gather_context(root)) 610 #import dis; dis.dis(fn) 611 fn() 612 #print "y-value:",y.value 613 assert y.value == 35

614 615 616 617 618 if __name__ == "__main__": test() 619

Source Code for Module park.modelling.parameter