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

  ViewVC Help
Powered by ViewVC 1.1.26