/[escript]/trunk/downunder/py_src/inversioncostfunctions.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
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 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
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<<float>>``
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