# Contents of /trunk/downunder/py_src/inversioncostfunctions.py

Revision 4416 - (show annotations)
Fri May 17 02:32:33 2013 UTC (6 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 19212 byte(s)
```Renaming a thing and adding some doco.
```
 1 2 ############################################################################## 3 # 4 # Copyright (c) 2003-2013 by University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development since 2012 by School of Earth Sciences 13 # 14 ############################################################################## 15 16 """Cost functions for inversions with one or more forward models""" 17 18 __copyright__="""Copyright (c) 2003-2013 by University of Queensland 19 http://www.uq.edu.au 20 Primary Business: Queensland, Australia""" 21 __license__="""Licensed under the Open Software License version 3.0 22 23 __url__= 24 25 __all__ = [ 'InversionCostFunction'] 26 27 from .costfunctions import MeteredCostFunction 28 from .mappings import Mapping 29 from .forwardmodels import ForwardModel 30 from esys.escript.pdetools import ArithmeticTuple 31 from esys.escript import Data 32 import numpy as np 33 34 35 class InversionCostFunction(MeteredCostFunction): 36 """ 37 Class to define cost function *J(m)* for inversion with one or more 38 forward models based on a multi-valued level set function *m*: 39 40 *J(m) = J_reg(m) + sum_f mu_f * J_f(p)* 41 42 where *J_reg(m)* is the regularization and cross gradient component of the 43 cost function applied to a level set function *m*, *J_f(p)* are the data 44 defect cost functions involving a physical forward model using the 45 physical parameter(s) *p* and *mu_f* is the trade-off factor for model f. 46 47 A forward model depends on a set of physical parameters *p* which are 48 constructed from components of the level set function *m* via mappings. 49 50 Example 1 (single forward model): 51 m=Mapping() 52 f=ForwardModel() 53 J=InversionCostFunction(Regularization(), m, f) 54 55 Example 2 (two forward models on a single valued level set) 56 m0=Mapping() 57 m1=Mapping() 58 f0=ForwardModel() 59 f1=ForwardModel() 60 61 J=InversionCostFunction(Regularization(), mappings=[m0, m1], forward_models=[(f0, 0), (f1,1)]) 62 63 Example 3 (two forward models on 2-valued level set) 64 m0=Mapping() 65 m1=Mapping() 66 f0=ForwardModel() 67 f1=ForwardModel() 68 69 J=InversionCostFunction(Regularization(numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)]) 70 71 :cvar provides_inverse_Hessian_approximation: if true the class provides an 72 approximative inverse of the Hessian operator. 73 """ 74 provides_inverse_Hessian_approximation=True 75 76 def __init__(self, regularization, mappings, forward_models): 77 """ 78 constructor for the cost function. 79 Stores the supplied object references and sets default weights. 80 81 :param regularization: the regularization part of the cost function 82 :type regularization: `Regularization` 83 :param mappings: the mappings to calculate physical parameters from the 84 regularization. This is a list of 2-tuples *(map, i)* 85 where the first component map defines a `Mapping` and 86 the second component *i* defines the index of the 87 component of level set function to be used to 88 calculate the mapping. Items in the list may also be 89 just `Mapping` objects in which case the entire level 90 set function is fed into the `Mapping` (typically used 91 for a single-component level set function. 92 :type mappings: `Mapping` or ``list`` 93 :param forward_models: the forward models involved in the calculation 94 of the cost function. This is a list of 2-tuples 95 *(f, ii)* where the first component f defines a 96 `ForwardModel` and the second component *ii* a 97 list of indexes referring to the physical 98 parameters in the `mappings` list. The 2-tuple 99 can be replaced by a `ForwardModel` if the 100 `mappings` list has a single entry. 101 :param forward_models: `ForwardModel` or ``list`` 102 """ 103 super(InversionCostFunction, self).__init__() 104 self.regularization=regularization 105 106 if isinstance(mappings, Mapping): 107 self.mappings = [mappings ] 108 else: 109 self.mappings = mappings 110 111 if isinstance(forward_models, ForwardModel): 112 self.forward_models = [ forward_models ] 113 else: 114 self.forward_models=forward_models 115 116 trafo = regularization.getCoordinateTransformation() 117 for m in self.forward_models: 118 if isinstance(m, ForwardModel): 119 F=m 120 else: 121 F=m[0] 122 if not F.getCoordinateTransformation() == trafo: 123 raise ValueError("Coordinate transformation for regularization and model don't match.") 124 125 self.numMappings=len(self.mappings) 126 self.numModels=len(self.forward_models) 127 self.numLevelSets = self.regularization.getNumLevelSets() 128 self.__num_tradeoff_factors = self.regularization.getNumTradeOffFactors() + self.numModels 129 self.setTradeOffFactorsModels() 130 131 def getDomain(self): 132 """ 133 returns the domain of the cost function 134 135 :rtype: `Domain` 136 """ 137 self.regularization.getDomain() 138 139 def getNumTradeOffFactors(self): 140 """ 141 returns the number of trade-off factors being used including the 142 trade-off factors used in the regularization component. 143 144 :rtype: ``int`` 145 """ 146 return self.__num_tradeoff_factors 147 148 def getForwardModel(self, idx=None): 149 """ 150 returns the *idx*-th forward model. 151 152 :param idx: model index. If cost function contains one model only `idx` 153 can be omitted. 154 :type idx: ``int`` 155 """ 156 if idx==None: idx=0 157 f=self.forward_models[idx] 158 if isinstance(f, ForwardModel): 159 F=f 160 else: 161 F=f[0] 162 return F 163 164 def getRegularization(self): 165 """ 166 returns the regularization 167 168 :rtype: `Regularization` 169 """ 170 return self.regularization 171 172 def setTradeOffFactorsModels(self, mu=None): 173 """ 174 sets the trade-off factors for the forward model components. 175 176 :param mu: list of the trade-off factors. If not present ones are used. 177 :type mu: ``float`` in case of a single model or a ``list`` of 178 ``float`` with the length of the number of models. 179 """ 180 if mu==None: 181 self.mu_model=np.ones((self.numModels, )) 182 else: 183 if self.numModels > 1: 184 mu=np.asarray(mu) 185 if min(mu) > 0: 186 self.mu_model= mu 187 else: 188 raise ValueError("All values for trade-off factor mu must be positive.") 189 else: 190 mu=float(mu) 191 if mu > 0: 192 self.mu_model= [mu, ] 193 else: 194 raise ValueError("Trade-off factor must be positive.") 195 196 def getTradeOffFactorsModels(self): 197 """ 198 returns the trade-off factors for the forward models 199 200 :rtype: ``float`` or ``list`` of ``float`` 201 """ 202 if self.numModels>1: 203 return self.mu_model 204 else: 205 return self.mu_model[0] 206 207 def setTradeOffFactorsRegularization(self, mu=None, mu_c=None): 208 """ 209 sets the trade-off factors for the regularization component of the 210 cost function, see `Regularization` for details. 211 212 :param mu: trade-off factors for the level-set variation part 213 :param mu_c: trade-off factors for the cross gradient variation part 214 """ 215 self.regularization.setTradeOffFactorsForVariation(mu) 216 self.regularization.setTradeOffFactorsForCrossGradient(mu_c) 217 218 def setTradeOffFactors(self, mu=None): 219 """ 220 sets the trade-off factors for the forward model and regularization 221 terms. 222 223 :param mu: list of trade-off factors. 224 :type mu: ``list`` of ``float`` 225 """ 226 if mu is None: 227 mu=np.ones((self.__num_tradeoff_factors,)) 228 self.setTradeOffFactorsModels(mu[:self.numModels]) 229 self.regularization.setTradeOffFactors(mu[self.numModels:]) 230 231 def getTradeOffFactors(self, mu=None): 232 """ 233 returns a list of the trade-off factors. 234 235 :rtype: ``list`` of ``float`` 236 """ 237 mu1=self.getTradeOffFactorsModels(mu[:self.numModels]) 238 mu2=self.regularization.getTradeOffFactors() 239 return [ m for m in mu1] + [ m for m in mu2] 240 241 def createLevelSetFunction(self, *props): 242 """ 243 returns an instance of an object used to represent a level set function 244 initialized with zeros. Components can be overwritten by physical 245 properties `props`. If present entries must correspond to the 246 `mappings` arguments in the constructor. Use ``None`` for properties 247 for which no value is given. 248 """ 249 m=self.regularization.getPDE().createSolution() 250 if len(props) > 0: 251 for i in range(self.numMappings): 252 if props[i]: 253 mm=self.mappings[i] 254 if isinstance(mm, Mapping): 255 m=mm.getInverse(props[i]) 256 elif len(mm) == 1: 257 m=mm[0].getInverse(props[i]) 258 else: 259 m[mm[1]]=mm[0].getInverse(props[i]) 260 return m 261 262 def getProperties(self, m, return_list=False): 263 """ 264 returns a list of the physical properties from a given level set 265 function *m* using the mappings of the cost function. 266 267 :param m: level set function 268 :type m: `Data` 269 :param return_list: if ``True`` a list is returned. 270 :type return_list: ``bool`` 271 :rtype: ``list`` of `Data` 272 """ 273 props=[] 274 for i in range(self.numMappings): 275 mm=self.mappings[i] 276 if isinstance(mm, Mapping): 277 p=mm.getValue(m) 278 elif len(mm) == 1: 279 p=mm[0].getValue(m) 280 else: 281 p=mm[0].getValue(m[mm[1]]) 282 props.append(p) 283 if self.numMappings > 1 or return_list: 284 return props 285 else: 286 return props[0] 287 288 def _getDualProduct(self, x, r): 289 """ 290 Returns the dual product, see `Regularization.getDualProduct` 291 292 :type x: `Data` 293 :type r: `ArithmeticTuple` 294 :rtype: ``float`` 295 """ 296 return self.regularization.getDualProduct(x, r) 297 298 def _getArguments(self, m): 299 """ 300 returns pre-computed values that are shared in the calculation of 301 *J(m)* and *grad J(m)*. In this implementation returns a tuple with the 302 mapped value of ``m``, the arguments from the forward model and the 303 arguments from the regularization. 304 305 :param m: current approximation of the level set function 306 :type m: `Data` 307 :return: tuple of of values of the parameters, pre-computed values 308 for the forward model and pre-computed values for the 309 regularization 310 :rtype: ``tuple`` 311 """ 312 args_reg=self.regularization.getArguments(m) 313 # cache for physical parameters: 314 props=self.getProperties(m, return_list=True) 315 args_f=[] 316 for i in range(self.numModels): 317 f=self.forward_models[i] 318 if isinstance(f, ForwardModel): 319 aa=f.getArguments(props[0]) 320 elif len(f) == 1: 321 aa=f[0].getArguments(props[0]) 322 else: 323 idx = f[1] 324 f=f[0] 325 if isinstance(idx, int): 326 aa=f.getArguments(props[idx]) 327 else: 328 pp=tuple( [ props[i] for i in idx] ) 329 aa=f.getArguments(*pp) 330 args_f.append(aa) 331 332 return props, args_f, args_reg 333 334 def _getValue(self, m, *args): 335 """ 336 Returns the value *J(m)* of the cost function at *m*. 337 If the pre-computed values are not supplied `getArguments()` is called. 338 339 :param m: current approximation of the level set function 340 :type m: `Data` 341 :param args: tuple of values of the parameters, pre-computed values 342 for the forward model and pre-computed values for the 343 regularization 344 :rtype: ``float`` 345 """ 346 if len(args)==0: 347 args=self.getArguments(m) 348 349 props=args[0] 350 args_f=args[1] 351 args_reg=args[2] 352 353 J = self.regularization.getValue(m, *args_reg) 354 355 for i in range(self.numModels): 356 f=self.forward_models[i] 357 if isinstance(f, ForwardModel): 358 J_f = f.getDefect(props[0],*args_f[i]) 359 elif len(f) == 1: 360 J_f=f[0].getDefect(props[0],*args_f[i]) 361 else: 362 idx = f[1] 363 f=f[0] 364 if isinstance(idx, int): 365 J_f = f.getDefect(props[idx],*args_f[i]) 366 else: 367 args=tuple( [ props[j] for j in idx] + args_f[i]) 368 J_f = f.getDefect(*args) 369 self.logger.debug("J_f[%d] = %e"%(i, J_f)) 370 self.logger.debug("mu_model[%d] = %e"%(i, self.mu_model[i])) 371 J += self.mu_model[i] * J_f 372 373 return J 374 375 def getComponentValues(self, m, *args): 376 return self._getComponentValues(m, *args) 377 378 def _getComponentValues(self, m, *args): 379 """ 380 returns the values of the individual cost functions that make up *f(x)* 381 using the precalculated values for *x*. 382 383 :param x: a solution approximation 384 :type x: x-type 385 :rtype: ``list<>`` 386 """ 387 if len(args)==0: 388 args=self.getArguments(m) 389 390 props=args[0] 391 args_f=args[1] 392 args_reg=args[2] 393 394 J_reg = self.regularization.getValue(m, *args_reg) 395 result = [J_reg] 396 397 for i in range(self.numModels): 398 f=self.forward_models[i] 399 if isinstance(f, ForwardModel): 400 J_f = f.getValue(props[0],*args_f[i]) 401 elif len(f) == 1: 402 J_f=f[0].getValue(props[0],*args_f[i]) 403 else: 404 idx = f[1] 405 f=f[0] 406 if isinstance(idx, int): 407 J_f = f.getValue(props[idx],*args_f[i]) 408 else: 409 args=tuple( [ props[j] for j in idx] + args_f[i]) 410 J_f = f.getValue(*args) 411 self.logger.debug("J_f[%d] = %e"%(i, J_f)) 412 self.logger.debug("mu_model[%d] = %e"%(i, self.mu_model[i])) 413 414 result += [J_f] # self.mu_model[i] * ?? 415 416 return result 417 418 def _getGradient(self, m, *args): 419 """ 420 returns the gradient of the cost function at *m*. 421 If the pre-computed values are not supplied `getArguments()` is called. 422 423 :param m: current approximation of the level set function 424 :type m: `Data` 425 :param args: tuple of values of the parameters, pre-computed values 426 for the forward model and pre-computed values for the 427 regularization 428 429 :rtype: `ArithmeticTuple` 430 431 :note: returns (Y^,X) where Y^ is gradient from regularisation plus gradients of fwd models. 432 X is the gradient of the regularisation wrt gradient of m. 433 """ 434 if len(args)==0: 435 args = self.getArguments(m) 436 437 props=args[0] 438 args_f=args[1] 439 args_reg=args[2] 440 441 g_J = self.regularization.getGradient(m, *args_reg) 442 p_diffs=[] 443 for i in range(self.numMappings): 444 mm=self.mappings[i] 445 if isinstance(mm, Mapping): 446 dpdm = mm.getDerivative(m) 447 elif len(mm) == 1: 448 dpdm = mm[0].getDerivative(m) 449 else: 450 dpdm = mm[0].getDerivative(m[mm[1]]) 451 p_diffs.append(dpdm) 452 453 Y=g_J[0] # Beacause g_J==(Y,X) Y_k=dKer/dm_k 454 for i in range(self.numModels): 455 mu=self.mu_model[i] 456 f=self.forward_models[i] 457 if isinstance(f, ForwardModel): 458 Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu 459 if self.numLevelSets == 1 : 460 Y +=Ys 461 else: 462 Y[0] +=Ys 463 elif len(f) == 1: 464 Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0] * mu 465 if self.numLevelSets == 1 : 466 Y +=Ys 467 else: 468 Y[0] +=Ys 469 else: 470 idx = f[1] 471 f = f[0] 472 if isinstance(idx, int): 473 Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu 474 if self.numLevelSets == 1 : 475 if idx == 0: 476 Y+=Ys 477 else: 478 raise IndexError("Illegal mapping index.") 479 else: 480 Y[idx] += Ys 481 else: 482 args = tuple( [ props[j] for j in idx] + args_f[i]) 483 Ys = f.getGradient(*args) 484 for ii in range(len(idx)): 485 Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]] * mu 486 487 return g_J 488 489 def _getInverseHessianApproximation(self, m, r, *args): 490 """ 491 returns an approximative evaluation *p* of the inverse of the Hessian 492 operator of the cost function for a given gradient type *r* at a 493 given location *m*: *H(m) p = r* 494 495 :param m: level set approximation where to calculate Hessian inverse 496 :type m: `Data` 497 :param r: a given gradient 498 :type r: `ArithmeticTuple` 499 :param args: tuple of values of the parameters, pre-computed values 500 for the forward model and pre-computed values for the 501 regularization 502 :rtype: `Data` 503 :note: in the current implementation only the regularization term is 504 considered in the inverse Hessian approximation. 505 """ 506 m=self.regularization.getInverseHessianApproximation(m, r, *args[2]) 507 return m 508 509 def updateHessian(self): 510 """ 511 notifies the class that the Hessian operator needs to be updated. 512 """ 513 self.regularization.updateHessian() 514 515 def _getNorm(self, m): 516 """ 517 returns the norm of `m` 518 519 :param m: level set function 520 :type m: `Data` 521 :rtype: ``float`` 522 """ 523 return self.regularization.getNorm(m) 524

 ViewVC Help Powered by ViewVC 1.1.26