park.modelling.model

1 # This program is public domain 2 """ 3 Define a park fitting model. 4 5 Usage 6 ===== 7 8 The simplest sort of fitting model is something like the following:: 9 10 import numpy 11 import park 12 def G(x,mu,sigma): 13 return numpy.exp(-0.5*(x-mu)**2/sigma**2) 14 15 class Gauss(park.Model): 16 parameters = ['center','width','scale'] 17 def eval(self,x): 18 return self.scale * G(x,self.center,self.width) 19 20 It has a function which is evaluated at a series of x values and 21 a set of adjustable parameters controlling the shape of f(x). 22 23 You can check your module with something like the following:: 24 25 $ ipython -pylab 26 27 from gauss import Gauss 28 29 g = Gauss(center=5,width=1,scale=10) 30 x = asarray([1,2,3,4,5]) 31 y = g(x) 32 plot(x,y) 33 34 This should produce a plot of the Gaussian peak. 35 36 You will then want to try your model with some data. Create a file 37 with some dummy data, such as gauss.dat:: 38 39 # x y 40 4 2 41 5 6 42 6 7 43 7 3 44 45 In order to match the model to data, you need to define a fitness 46 object. This shows the difference between the model and the data, 47 which you can then plot or sum to create a weighted chisq value:: 48 49 f = park.Fitness(g, 'gauss.dat') 50 plot(f.data.fit_x, f.residuals()) 51 52 The data file can have up to four columns, either x,y or x,y,dy 53 or x,dx,y,dy where x,y are the measurement points and values, 54 dx is the instrument resolution in x and dy is the uncertainty 55 in the measurement value y. You can pass any keyword arguments 56 to data.load that are accepted by numpy.loadtxt. For example, 57 you can reorder columns or skip rows. Additionally, you can modify 58 data attributes directly x,y,dx,dy and calc_x. See help on park.Data1D 59 for details. 60 61 Once you have your model debugged you can use it in a fit:: 62 63 g = Gauss(center=[3,5],scale=7.2,width=[1,3]) 64 result = park.fit((g, 'gauss.dat')) 65 result.plot() 66 67 In this example, center and width are allowed to vary but scale is fixed. 68 69 Existing models can be readily adapted to Park:: 70 71 class Gauss(object): 72 "Existing model" 73 center,width,scale = 0,1,1 74 def __init__(self,**kw): 75 for k,v in kw.items(): setattr(self,k,v) 76 def eval(self,x): 77 return self.scale *G(x,self.center,self.width) 78 79 class GaussAdaptor(Gauss,Model): 80 "PARK adaptor" 81 parameters = ['center','width','scale'] 82 def __init__(self,*args,**kw): 83 Model.__init__(self,*args) 84 Gauss.__init__(self,**kw) 85 86 g = GaussAdaptor(center=[3,5],scale=7.2,width=[1,3]) 87 result = park.fit((g, 'gauss.dat')) 88 result.plot() 89 90 Models can become much more complex than the ones described above, 91 including multilevel models where fitting parameters can be added 92 and removed dynamically. 93 94 95 You can access model parameters directly:: 96 97 M1 = Gauss('M1',center=[3,5],scale=7.2,width=[1,3]) 98 print M1.center # Access the value for parameter 'center' 99 M1.center = 7 # Set the value for parameter 'center' 100 M1.center = [5,9] # Set the fit range for parameter 'center' 101 M1.center = "2*M1.scale" # Set parameter 'center' to an expression 102 103 Note that a strange side-effect of the convenience functions used to 104 access model parameters is that assignment statements do not act like 105 you expect. The following sets the initial value of the parameter to 106 7 then sets the range to [5,9]:: 107 108 M1.center = 7 # Set the value for parameter 'center' 109 M1.center = [5,9] # Set the fit range for parameter 'center' 110 111 If these statements were reversed, the initial value would still be 7 and 112 the range [5,9] but the parameter would be set to 'fixed' rather than 'fitted'. 113 True global optimizers will not care what the initial value is, so from the 114 perspective of the fit this should not change the results. Most optimizers 115 will however use the initial value as part of the search. 116 117 Since we may choose to improve this interface in future, programmatic access 118 to models should go through parameter sets:: 119 120 M1.parameterset['center'].value 121 M1.parameterset['center'].range 122 ... 123 124 Once models are defined they can be used in a variety of contexts, such 125 as simultaneous fitting with constraints between the parameters. With 126 some care in defining the model, computationally intensive fits can 127 be distributed across multiple processors. We provide a simple user 128 interface for interacting with the model parameters and managing fits. 129 This can be extended with model specialized model editors which format 130 the parameters in a sensible way for the model, or allow direct manipulation 131 of the model structure. The underlying fitting engine can also be 132 used directly from your own user interface. 133 134 For those situations where the properties are dynamic or programmatically 135 generated the meta class machinery which translates parameters into 136 properties is likely not going to be useful. In this case you can 137 inherit directly from park.model.BaseModel, or simply define the 138 required API by hand. 139 """ 140 141 __all__ = ['Model'] 142 143 import inspect 144 from copy import copy, deepcopy 145 146 import park 147 from park.modelling.parameter import Parameter, ParameterSet 148 149

150 -class BaseModel(object):

151 """ 152 Model definition. 153 154 The model manages attribute access to the fitting parameters and 155 also manages the dataset. 156 157 derivatives ['p1','p2',...] 158 List of parameters for which the model can calculate 159 derivatives. The eval_derivs() method will be called 160 to calculate the derivatives with respect to the parameters. 161 Parameters with no analytic derivatives available will be 162 approximated by numerical differentiation. 163 164 eval(x) 165 Evaluate the model at x. 166 167 eval_deriv(x,pars=['p1','p2',...]) 168 Evaluate the model and the derivatives at x. 169 170 parameterset 171 The set of parameters defined by the model. 172 173 The minimum fitting model if you choose not to subclass park.BaseModel 174 requires parameterset and an eval() method. 175 176 Note that BaseModel is a simple class which 177 """ 178 derivatives = [] 179 parameters = [] 180 parameterset = None

181 - def __init__(self,*args,**kw):

182 """ 183 Initialize the model. 184 185 Model('name',P1=value, P2=Value, ...) 186 187 When overriding __init__ in the subclass be sure to call 188 Model.__init__(self, *args, **kw). This makes a private 189 copy of the parameterset for the model and initializes 190 any parameters set using keyword arguments. 191 """ 192 #print 'calling model init on',id(self) 193 # To avoid trashing the Model.parname = Parameter() template 194 # when we create an instance and accessor for the parameter 195 # we need to make sure that we are using our own copy of the 196 # model dictionary stored in vars 197 if len(args)>1: raise TypeError("expect name as model parameter") 198 name = args[0] if len(args)>=1 else 'unknown' 199 # We are copying the parameters from the model class to the 200 # model instance. If we don't do this, then parameters from different 201 # instances of the same model will collide. 202 self.parameterset = deepcopy(self.parameterset) 203 for k,v in kw.items(): 204 self.parameterset[k].set(v) 205 self.parameterset._name = name

206 #print 'done',id(self) 207

208 - def eval(self, *args):

209 """ 210 Evaluate the model at x. 211 212 This method needs to be specialized in the model to evaluate the 213 model function. 214 """ 215 raise NotImplementedError('Model needs to implement eval')

216

217 - def eval_derivs(self, x, pars=[]):

218 """ 219 Evaluate the model and derivatives wrt pars at x. 220 221 pars is a list of the names of the parameters for which derivatives 222 are desired. The fitting engine may be working with a subset of the 223 parameters stored in the model, and the parameters may be requested 224 in a different order from what they were stored. 225 226 Returns (f(x), df/dp1(x), df/dp2(x), ...). 227 228 Note that the analytic derivatives need not be defined for all 229 parameters; numerical derivatives can be calculated when needed. 230 If when analytic derivatives are available, they may not be 231 used because the parameter in question is used as part of an 232 expression for other parameters. Automatic differentiation 233 on parameter expressions is possible, but beyond the scope 234 of this project. 235 236 This method needs to be specialized in the model to evaluate the 237 model function. 238 """ 239 raise NotImplementedError('Model needs to implement eval_derivs')

240 241 # User convenience functions

242 - def set(self, **kw):

243 """ 244 Set the initial value for a set of parameters. 245 246 E.g., model.set(width=3,center=5) 247 """ 248 for k,v in kw.items(): 249 self.parameterset[k].set(v)

250

251 - def __getitem__(self, p):

252 """ 253 Return the value for parameter p (by name or number). 254 """ 255 return self.parameterset[p].get()

256

257 - def __setitem__(self, p, value):

258 """ 259 Set the value for parameter p (by name or number). 260 """ 261 return self.parameterset[p].set(value)

262

263 - def __call__(self, *args):

264 """ 265 Evaluate the model at x. 266 """ 267 return self.eval(*args)

268 269 270 # The following classes make it easy to define a model from scratch 271 # by specifying the parameters at the class level and evaluating them 272 # as attributes. For example: 273 # class Scale(Model): 274 # a = Parameter() 275 # def eval(self, x): return self.a*x 276 # This involves several pieces: 277 # ParameterProperty defines the accessor 278 # MetaModel converts class parameters to accessors and parameter sets 279 # Model is the public interface

280 -class ParameterProperty(object):

281 """ 282 Attribute accessor for a named parameter. 283 284 Objects of class ParameterProperty act similarly to normal property 285 objects in that assignment to object.attr triggers the __set__ 286 method of ParameterProperty and queries of the value object.attr 287 triggers the __get__ method. These methods look up the actual 288 parameter in the dictionary model.parameterset, delegating the 289 the set/get methods of the underlying parameter object. 290 291 For example:: 292 293 model.P1 = 5 294 print model.P1 295 296 is equivalent to:: 297 298 model.parameterset['P1'].set(5) 299 print model.parameterset['P1'].get() 300 301 To enable this behaviour, the model developer must assign a 302 ParameterProperty object to a class attribute of the model. It 303 must be a class attribute. Properties assigned as an instance 304 attribute in the __init__ of the class will not be recognized 305 as properties. 306 307 An example model will look something like the following:: 308 309 class MyModel: 310 # A property must be assigned as a class attribute! 311 P1 = ParameterProperty('P1') 312 P2 = ParameterProperty('P2') 313 def __init__(self, P1=None, P2=None): 314 parP1 = Parameter('P1') 315 if P1 is not None: parP1.set(P1) 316 parP2 = Parameter('P2') 317 if P2 is not None: parP2.set(P2) 318 self.parameterset = { 'P1': parP1, 'P2': parP2 } 319 320 This work is done implicitly by MetaModel, and by subclassing 321 the class Model, the model developer does not ever need to 322 use ParameterProperty. 323 """

324 - def __init__(self,name,**kw):

325 self.name = name

326 - def __get__(self,instance,cls):

327 return instance.parameterset[self.name].get()

328 - def __set__(self,instance,value):

329 instance.parameterset[self.name].set(value)

330 331

332 -class MetaModel(type):

333 """ 334 Place class level parameter definitions into the parameter set. 335 336 This is a meta class, and the usual meta class rules apply. That is, 337 the Model base class should be defined like:: 338 339 class Model(object): 340 __metaclass__ = MetaModel 341 ... 342 343 The MetaModel.__new__ method is called whenever a new Model 344 class is created. The name of the model, its superclasses and 345 its attributes are passed to MetaModel.__new__ which creates the 346 actual class. It is not called for new instances of the model. 347 348 MetaModel is designed to simplify the definition of parameters for 349 the model. When processing the class, MetaModel defines 350 parameters which is a list of names of all the parameters in the 351 model, and parameterset, which is a dictionary of the actual 352 parameters and any parameter sets in the model. 353 354 Parameters can be defined using a list of names in the parameter 355 attribute, or by defining the individual parameters as separate 356 attributes of class Parameter. 357 358 For example, the following defines parameters P1, P2, and P3:: 359 360 class MyModel(Model): 361 parameters = ['P1','P2'] 362 P2 = Parameter(range=[0,inf]) 363 P3 = Parameter() 364 365 For each parameter, MetaModel will define a parameter accessor, 366 and add the parameter definition to the parameter set. The accessor 367 delegates query and assignment to the Parameter get/set methods. The 368 attributes of the particular parameter instance can be 369 adjusted using model.parameterset['name'].attribute = value. 370 """

371 - def __new__(cls, name, bases, vars):

372 #print 'calling model meta',vars 373 374 # Find all the parameters for the model. The listed parameters 375 # are defined using:: 376 # parameters = ['P1', 'P2', ...] 377 # The remaining parameters are defined individually:: 378 # P3 = Parameter() 379 # P4 = Parameter() 380 # The undefined parameters are those that are listed but not defined. 381 listed = vars.get('parameters',[]) 382 defined = [k for k,v in vars.items() if isinstance(v,Parameter)] 383 undefined = [k for k in listed if k not in defined] 384 #print listed,defined,undefined 385 386 # Define a getter/setter for every parameter so that the user 387 # can say model.name to access parameter name. 388 # 389 # Create a parameter object for every undefined parameter. 390 # Check if the base class defines a default value, and use 391 # that for the initial value. We don't want to do this for 392 # parameters explicitly defined since the user may have 393 # given them a default value already. 394 # 395 # For predefined parameters we must set the name explicitly. 396 # This saves the user from having to use, e.g.:: 397 # P1 = Parameter('P1') 398 pars = [] 399 for p in undefined: 400 # Check if the attribute value is initialized in a base class 401 for b in bases: 402 if hasattr(b,p): 403 value = getattr(b,p) 404 break 405 else: 406 value = 0. 407 #print "looking up value in base as",value 408 pars.append(Parameter(name=p,value=value)) 409 vars[p] = ParameterProperty(p) 410 for p in defined: 411 # Make sure parameter name is correct. Note that we are using 412 # _name rather than .name to access the name attribute since 413 # name is a read-only parameter. 414 vars[p]._name = p 415 pars.append(vars[p]) 416 vars[p] = ParameterProperty(p) 417 418 # Construct the class attributes 'parameters' and 'parameterset'. 419 # Listed parameters are given first, with unlisted parameters 420 # following alphabetically. For hierarchical parameter sets, 421 # we also need to include the defined sets. 422 unlisted = list(set(defined+undefined) - set(listed)) 423 unlisted.sort() 424 parameters = listed + unlisted 425 parsets = [k for k,v in vars.items() if isinstance(v,ParameterSet)] 426 vars['parameters'] = parameters 427 vars['parameterset'] = ParameterSet(pars=pars+parsets) 428 #print 'pars',vars['parameters'] 429 #print 'parset',vars['parameterset'] 430 431 # Return a new specialized instance of the class with parameters 432 # and parameterset made explicit. 433 return type.__new__(cls, name, bases, vars)

434

435 -class Model(BaseModel):

436 """ 437 Model definition. 438 439 The model manages attribute access to the fitting parameters and 440 also manages the dataset. 441 442 derivatives ['p1','p2',...] 443 List of parameters for which the model can calculate 444 derivatives. The eval_derivs() method will be called 445 to calculate the derivatives with respect to the parameters. 446 Parameters with no analytic derivatives available will be 447 approximated by numerical differentiation. 448 449 eval(x) 450 Evaluate the model at x. This must be defined by the subclass. 451 452 eval_deriv(x,pars=[]) 453 Evaluate the model and the derivatives at x. This must be 454 defined by the subclass. 455 456 parameters 457 The names of the model parameters. If this is not provided, then 458 the model will search the subclass for park.Parameter attributes 459 and construct the list of names from that. Any parameters in the 460 list not already defined as park.Parameter attributes will be 461 defined as parameters with a default of 0. 462 463 parameterset 464 The set of parameters defined by the model. These are the 465 parameters themselves, gathered into a park.ParameterSet. 466 467 The minimum fittng model if you choose not to subclass park.Model 468 requires parameterset and a residuals() method. 469 """ 470 __metaclass__ = MetaModel

471

472 -def add_parameter(model, name, **kw):

473 """ 474 Add a parameter to an existing park model. 475 476 Note: this may not work if it is operating on a BaseModel. 477 """ 478 setattr(model.__class__, name, park.model.ParameterProperty(name)) 479 model.parameterset.append(park.Parameter(name=name, **kw))

480

481 -def _extract_parameters(fn):

482 """ 483 Find the keyword arguments associated with a function and return them 484 as parameters. 485 """ 486 names,arg,kw,values = inspect.getargspec(fn) 487 if len(values)<1: 488 raise ValueError("Need keyword args to model function") 489 parameters = names[-len(values):] 490 eval_args = names[:len(names)-len(values)] 491 return [park.Parameter(k,value=v) for k,v in zip(parameters,values)]

492

493 -class FunctionModel(BaseModel):

494 """ 495 Define a model based on a function. The fitting parameters should 496 be defined as keyword arguments to the function. 497 """

498 - def __init__(self, fn, *args, **kw):

499 self.fn = fn 500 self.parameterset = park.ParameterSet('',_extract_parameters(fn)) 501 self.set(**kw)

502

503 - def eval(self, *args):

504 kw = dict([(p.name,p.value) for p in self.parameterset]) 505 return self.fn(*args, **kw)

506 507

508 -def test():

509 # Gauss theory function 510 import numpy 511 eps,inf = numpy.finfo('d').eps,numpy.inf 512 def G(x,mu,sigma): 513 mu,sigma = numpy.asarray(mu),numpy.asarray(sigma) 514 return numpy.exp(-0.5*(x-mu)**2/sigma**2)

515 516 # === Minimal model === 517 # just list the fitting parameters and define the function. 518 class Gauss(Model): 519 parameters = ['center','width','scale'] 520 def eval(self,x): 521 return self.scale * G(x,self.center,self.width) 522 523 # Create a couple of models and make sure they don't conflict 524 g1 = Gauss(center=5,width=1,scale=2) 525 assert g1.center == 5 526 g2 = Gauss(center=6,width=1) 527 assert g1.center == 5 528 assert g2.center == 6 529 assert g1(5) == 2 530 assert g2(6) == 0 531 532 # === explore parameters === 533 # dedicated model using park parameters directly, and defining derivatives 534 class Gauss(Model): 535 center = Parameter(value=0, 536 tip='Peak center') 537 width = Parameter(value=1, 538 limits=(eps,inf), 539 tip='Peak width (1-sigma)') 540 scale = Parameter(value=1, 541 limits=(0,inf), 542 tip='Peak scale (>0)') 543 def eval(self,x): 544 return self.scale * G(x,self.center,self.width) 545 546 # Derivatives 547 derivatives = ['center','width','scale'] 548 def dscale(self, g, x): 549 return g/self.scale 550 def dcenter(self, g, x): 551 return 2*(x-self.center)/self.width*g 552 def dwidth(self, g, x): 553 return 2*(x-self.center)/self.width**3*g 554 dmap = dict(scale=dscale,center=dcenter,width=dwidth) 555 def eval_derivs(self,x,pars=[]): 556 """ 557 Calculate function value and derivatives wrt to the parameters. 558 559 pars is a list of parameter names, possibly consisting of any 560 parameter with deriv=True. 561 """ 562 g = self.eval(x) 563 dg = [self.dmap[p](self,g,x) for p in pars] 564 return [g]+dg 565 566 g1 = Gauss(center=5) 567 g1.parameterset['center'].tip = 'This is the center' 568 assert g1.center == 5 569 g2 = Gauss(center=6) 570 assert g1.center == 5 571 assert g2.center == 6 572 assert g1(5) == 1 573 assert g2(6) == 1 574 575 # ====== Test wrapper ======= 576 # delegate to existing model via inheritence 577 class Gauss(object): 578 """Pre-existing model""" 579 center,width,scale = 0,1,1 580 def __init__(self,**kw): 581 #print "calling BaseGauss init on",id(self) 582 for k,v in kw.items(): setattr(self,k,v) 583 #print "done",id(self) 584 def eval(self,x): 585 return self.scale *G(x,self.center,self.width) 586 587 class GaussAdaptor(Gauss,Model): 588 """PARK wrapper""" 589 parameters = ['center','width','scale'] 590 def __init__(self,*args,**kw): 591 #print "calling Gauss init on",id(self) 592 Model.__init__(self,*args) 593 Gauss.__init__(self,**kw) 594 #print "done",id(self) 595 596 g1 = GaussAdaptor(center=5) 597 g2 = GaussAdaptor(center=6) 598 g3 = GaussAdaptor() 599 assert g1.center == 5 600 assert g2.center == 6 601 assert g3.scale == 1 602 assert g1(5) == 1 603 assert g2(6) == 1 604 605 g4 = GaussAdaptor('M3',center=6) 606 assert g4.parameterset.name == 'M3' 607 608 # define a model from a keyword function 609 def line(x, a=1, b=2): return a*x+b 610 linemodel = FunctionModel(line,a=5) 611 assert linemodel['a'] == 5 612 assert linemodel(3) == 5*3+2 613 614 # dedicated multilevel model using park parameters directly 615 class MultiGauss(Model): 616 def add(self,name,model): 617 pass 618 619 # wrapped model using park parameters indirectly 620 class BaseMultiGauss(object): 621 def __init__(self): 622 self.models = [] 623 def add(self,**kw): 624 self.models.append(BaseGauss(**kw)) 625 class WrapMultiGauss(BaseMultiGauss,Model): 626 def __init__(self): 627 pass 628 629 if __name__ == "__main__": test() 630

Source Code for Module park.modelling.model