/[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 4561 - (show annotations)
Wed Dec 4 05:19:32 2013 UTC (5 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 20656 byte(s)
Removed spurious print and some minor changes in inversion example.

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
132 self.__Q={}
133 for i in range(self.numModels):
134 # find the index of the physical parameters used to evaluate a model
135 f=self.forward_models[i]
136 if isinstance(f, ForwardModel):
137 idx=[0]
138 elif len(f) == 1:
139 idx=[0]
140 else:
141 idx = f[1]
142 if isinstance(idx, int):
143 idx=[idx]
144 else:
145 pass
146
147 q=[]
148 for k in idx:
149 mm=self.mappings[k]
150 if isinstance(mm, Mapping):
151 q.append(0)
152 elif len(mm) == 1:
153 q.append(0)
154 else:
155 if isinstance(mm[1], int):
156 q.append(mm[1])
157 else:
158 q+=list(mm[1])
159 self.__Q[i]=q
160 # for each parameter find the component of the levelset function beeing used to define it:
161 self.__parameter_to_level_set={}
162 for k in xrange(self.numMappings):
163 mm=self.mappings[k]
164 if isinstance(mm, Mapping):
165 q=0
166 elif len(mm) == 1:
167 q=0
168 else:
169 if isinstance(mm[1], int):
170 q=mm[1]
171 self.__parameter_to_level_set[k]=q
172
173 def getDomain(self):
174 """
175 returns the domain of the cost function
176
177 :rtype: `Domain`
178 """
179 self.regularization.getDomain()
180
181 def getNumTradeOffFactors(self):
182 """
183 returns the number of trade-off factors being used including the
184 trade-off factors used in the regularization component.
185
186 :rtype: ``int``
187 """
188 return self.__num_tradeoff_factors
189
190 def getForwardModel(self, idx=None):
191 """
192 returns the *idx*-th forward model.
193
194 :param idx: model index. If cost function contains one model only `idx`
195 can be omitted.
196 :type idx: ``int``
197 """
198 if idx==None: idx=0
199 f=self.forward_models[idx]
200 if isinstance(f, ForwardModel):
201 F=f
202 else:
203 F=f[0]
204 return F
205
206 def getRegularization(self):
207 """
208 returns the regularization
209
210 :rtype: `Regularization`
211 """
212 return self.regularization
213
214 def setTradeOffFactorsModels(self, mu=None):
215 """
216 sets the trade-off factors for the forward model components.
217
218 :param mu: list of the trade-off factors. If not present ones are used.
219 :type mu: ``float`` in case of a single model or a ``list`` of
220 ``float`` with the length of the number of models.
221 """
222 if mu==None:
223 self.mu_model=np.ones((self.numModels, ))
224 else:
225 if self.numModels > 1:
226 mu=np.asarray(mu)
227 if min(mu) > 0:
228 self.mu_model= mu
229 else:
230 raise ValueError("All values for trade-off factor mu must be positive.")
231 else:
232 mu=float(mu)
233 if mu > 0:
234 self.mu_model= [mu, ]
235 else:
236 raise ValueError("Trade-off factor must be positive.")
237
238 def getTradeOffFactorsModels(self):
239 """
240 returns the trade-off factors for the forward models
241
242 :rtype: ``float`` or ``list`` of ``float``
243 """
244 if self.numModels>1:
245 return self.mu_model
246 else:
247 return self.mu_model[0]
248
249 def setTradeOffFactorsRegularization(self, mu=None, mu_c=None):
250 """
251 sets the trade-off factors for the regularization component of the
252 cost function, see `Regularization` for details.
253
254 :param mu: trade-off factors for the level-set variation part
255 :param mu_c: trade-off factors for the cross gradient variation part
256 """
257 self.regularization.setTradeOffFactorsForVariation(mu)
258 self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
259
260 def setTradeOffFactors(self, mu=None):
261 """
262 sets the trade-off factors for the forward model and regularization
263 terms.
264
265 :param mu: list of trade-off factors.
266 :type mu: ``list`` of ``float``
267 """
268 if mu is None:
269 mu=np.ones((self.__num_tradeoff_factors,))
270 self.setTradeOffFactorsModels(mu[:self.numModels])
271 self.regularization.setTradeOffFactors(mu[self.numModels:])
272
273 def getTradeOffFactors(self, mu=None):
274 """
275 returns a list of the trade-off factors.
276
277 :rtype: ``list`` of ``float``
278 """
279 mu1=self.getTradeOffFactorsModels(mu[:self.numModels])
280 mu2=self.regularization.getTradeOffFactors()
281 return [ m for m in mu1] + [ m for m in mu2]
282
283 def createLevelSetFunction(self, *props):
284 """
285 returns an instance of an object used to represent a level set function
286 initialized with zeros. Components can be overwritten by physical
287 properties `props`. If present entries must correspond to the
288 `mappings` arguments in the constructor. Use ``None`` for properties
289 for which no value is given.
290 """
291 m=self.regularization.getPDE().createSolution()
292 if len(props) > 0:
293 for i in range(self.numMappings):
294 if props[i]:
295 mm=self.mappings[i]
296 if isinstance(mm, Mapping):
297 m=mm.getInverse(props[i])
298 elif len(mm) == 1:
299 m=mm[0].getInverse(props[i])
300 else:
301 m[mm[1]]=mm[0].getInverse(props[i])
302 return m
303
304 def getProperties(self, m, return_list=False):
305 """
306 returns a list of the physical properties from a given level set
307 function *m* using the mappings of the cost function.
308
309 :param m: level set function
310 :type m: `Data`
311 :param return_list: if ``True`` a list is returned.
312 :type return_list: ``bool``
313 :rtype: ``list`` of `Data`
314 """
315 props=[]
316 for i in range(self.numMappings):
317 mm=self.mappings[i]
318 if isinstance(mm, Mapping):
319 p=mm.getValue(m)
320 elif len(mm) == 1:
321 p=mm[0].getValue(m)
322 else:
323 p=mm[0].getValue(m[mm[1]])
324 props.append(p)
325 if self.numMappings > 1 or return_list:
326 return props
327 else:
328 return props[0]
329
330 def _getDualProduct(self, x, r):
331 """
332 Returns the dual product, see `Regularization.getDualProduct`
333
334 :type x: `Data`
335 :type r: `ArithmeticTuple`
336 :rtype: ``float``
337 """
338 return self.regularization.getDualProduct(x, r)
339
340 def _getArguments(self, m):
341 """
342 returns pre-computed values that are shared in the calculation of
343 *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
344 mapped value of ``m``, the arguments from the forward model and the
345 arguments from the regularization.
346
347 :param m: current approximation of the level set function
348 :type m: `Data`
349 :return: tuple of of values of the parameters, pre-computed values
350 for the forward model and pre-computed values for the
351 regularization
352 :rtype: ``tuple``
353 """
354 args_reg=self.regularization.getArguments(m)
355 # cache for physical parameters:
356 props=self.getProperties(m, return_list=True)
357 args_f=[]
358 for i in range(self.numModels):
359 f=self.forward_models[i]
360 if isinstance(f, ForwardModel):
361 aa=f.getArguments(props[0])
362 elif len(f) == 1:
363 aa=f[0].getArguments(props[0])
364 else:
365 idx = f[1]
366 f=f[0]
367 if isinstance(idx, int):
368 aa=f.getArguments(props[idx])
369 else:
370 pp=tuple( [ props[i] for i in idx] )
371 aa=f.getArguments(*pp)
372 args_f.append(aa)
373
374 return props, args_f, args_reg
375
376 def _getValue(self, m, *args):
377 """
378 Returns the value *J(m)* of the cost function at *m*.
379 If the pre-computed values are not supplied `getArguments()` is called.
380
381 :param m: current approximation of the level set function
382 :type m: `Data`
383 :param args: tuple of values of the parameters, pre-computed values
384 for the forward model and pre-computed values for the
385 regularization
386 :rtype: ``float``
387 """
388 if len(args)==0:
389 args=self.getArguments(m)
390
391 props=args[0]
392 args_f=args[1]
393 args_reg=args[2]
394
395 J = self.regularization.getValue(m, *args_reg)
396
397 for i in range(self.numModels):
398 f=self.forward_models[i]
399 if isinstance(f, ForwardModel):
400 J_f = f.getDefect(props[0],*args_f[i])
401 elif len(f) == 1:
402 J_f=f[0].getDefect(props[0],*args_f[i])
403 else:
404 idx = f[1]
405 f=f[0]
406 if isinstance(idx, int):
407 J_f = f.getDefect(props[idx],*args_f[i])
408 else:
409 args=tuple( [ props[j] for j in idx] + args_f[i])
410 J_f = f.getDefect(*args)
411 self.logger.debug("J_f[%d] = %e, mu_model[%d] = %e"%(i, J_f, i, self.mu_model[i]))
412 J += self.mu_model[i] * J_f
413
414 return J
415
416 def getComponentValues(self, m, *args):
417 return self._getComponentValues(m, *args)
418
419 def _getComponentValues(self, m, *args):
420 """
421 returns the values of the individual cost functions that make up *f(x)*
422 using the precalculated values for *x*.
423
424 :param x: a solution approximation
425 :type x: x-type
426 :rtype: ``list<<float>>``
427 """
428 if len(args)==0:
429 args=self.getArguments(m)
430
431 props=args[0]
432 args_f=args[1]
433 args_reg=args[2]
434
435 J_reg = self.regularization.getValue(m, *args_reg)
436 result = [J_reg]
437
438 for i in range(self.numModels):
439 f=self.forward_models[i]
440 if isinstance(f, ForwardModel):
441 J_f = f.getValue(props[0],*args_f[i])
442 elif len(f) == 1:
443 J_f=f[0].getValue(props[0],*args_f[i])
444 else:
445 idx = f[1]
446 f=f[0]
447 if isinstance(idx, int):
448 J_f = f.getValue(props[idx],*args_f[i])
449 else:
450 args=tuple( [ props[j] for j in idx] + args_f[i])
451 J_f = f.getValue(*args)
452 self.logger.debug("J_f[%d] = %e, mu_model[%d] = %e"%(i, J_f, i, self.mu_model[i]))
453
454 result += [J_f] # self.mu_model[i] * ??
455
456 return result
457
458
459 def _getGradient(self, m, *args):
460 """
461 returns the gradient of the cost function at *m*.
462 If the pre-computed values are not supplied `getArguments()` is called.
463
464 :param m: current approximation of the level set function
465 :type m: `Data`
466 :param args: tuple of values of the parameters, pre-computed values
467 for the forward model and pre-computed values for the
468 regularization
469
470 :rtype: `ArithmeticTuple`
471
472 :note: returns (Y^,X) where Y^ is gradient from regularisation plus gradients of fwd models.
473 X is the gradient of the regularisation wrt gradient of m.
474 """
475 if len(args)==0:
476 args = self.getArguments(m)
477
478 props=args[0]
479 args_f=args[1]
480 args_reg=args[2]
481
482 g_J = self.regularization.getGradient(m, *args_reg)
483 p_diffs=[]
484 for i in range(self.numMappings):
485 mm=self.mappings[i]
486 if isinstance(mm, Mapping):
487 dpdm = mm.getDerivative(m)
488 elif len(mm) == 1:
489 dpdm = mm[0].getDerivative(m)
490 else:
491 dpdm = mm[0].getDerivative(m[mm[1]])
492 p_diffs.append(dpdm)
493
494 Y=g_J[0] # Beacause g_J==(Y,X) Y_k=dKer/dm_k
495 for i in range(self.numModels):
496 mu=self.mu_model[i]
497 f=self.forward_models[i]
498 if isinstance(f, ForwardModel):
499 Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
500 if self.numLevelSets == 1 :
501 Y +=Ys
502 else:
503 Y[self.__parameter_to_level_set[0]] +=Ys
504 elif len(f) == 1:
505 Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
506 if self.numLevelSets == 1 :
507 Y +=Ys
508 else:
509 Y[self.__parameter_to_level_set[0]] +=Ys
510 else:
511 idx = f[1]
512 f = f[0]
513 if isinstance(idx, int):
514 Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu
515 if self.numLevelSets == 1 :
516 Y+=Ys
517 else:
518 Y[self.__parameter_to_level_set[idx]] += Ys
519 else:
520 args = tuple( [ props[j] for j in idx] + args_f[i])
521 Ys = f.getGradient(*args)
522 for j in range(idx):
523 if self.numLevelSets == 1 :
524 Y+=Ys[j]* p_diffs[idx[j]] * mu
525 else:
526 Y[self.__parameter_to_level_set[j]]+=Ys[j]* p_diffs[idx[j]] * mu
527 return g_J
528
529 def _getInverseHessianApproximation(self, m, r, *args):
530 """
531 returns an approximative evaluation *p* of the inverse of the Hessian
532 operator of the cost function for a given gradient type *r* at a
533 given location *m*: *H(m) p = r*
534
535 :param m: level set approximation where to calculate Hessian inverse
536 :type m: `Data`
537 :param r: a given gradient
538 :type r: `ArithmeticTuple`
539 :param args: tuple of values of the parameters, pre-computed values
540 for the forward model and pre-computed values for the
541 regularization
542 :rtype: `Data`
543 :note: in the current implementation only the regularization term is
544 considered in the inverse Hessian approximation.
545 """
546 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
547 return m
548
549 def updateHessian(self):
550 """
551 notifies the class that the Hessian operator needs to be updated.
552 """
553 self.regularization.updateHessian()
554
555 def _getNorm(self, m):
556 """
557 returns the norm of `m`
558
559 :param m: level set function
560 :type m: `Data`
561 :rtype: ``float``
562 """
563 return self.regularization.getNorm(m)
564

  ViewVC Help
Powered by ViewVC 1.1.26