/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (19 months, 2 weeks ago) by jfenwick
File MIME type: text/x-python
File size: 21895 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2018 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Apache License, version 2.0
9 # http://www.apache.org/licenses/LICENSE-2.0
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-2018 by The University of Queensland
22 http://www.uq.edu.au
23 Primary Business: Queensland, Australia"""
24 __license__="""Licensed under the Apache License, version 2.0
25 http://www.apache.org/licenses/LICENSE-2.0"""
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 import Data, inner, interpolate
34 import numpy as np
35
36
37 class InversionCostFunction(MeteredCostFunction):
38 """
39 Class to define cost function *J(m)* for inversion with one or more
40 forward models based on a multi-valued level set function *m*:
41
42 *J(m) = J_reg(m) + sum_f mu_f * J_f(p)*
43
44 where *J_reg(m)* is the regularization and cross gradient component of the
45 cost function applied to a level set function *m*, *J_f(p)* are the data
46 defect cost functions involving a physical forward model using the
47 physical parameter(s) *p* and *mu_f* is the trade-off factor for model f.
48
49 A forward model depends on a set of physical parameters *p* which are
50 constructed from components of the level set function *m* via mappings.
51
52 Example 1 (single forward model):
53 m=Mapping()
54 f=ForwardModel()
55 J=InversionCostFunction(Regularization(), m, f)
56
57 Example 2 (two forward models on a single valued level set)
58 m0=Mapping()
59 m1=Mapping()
60 f0=ForwardModel()
61 f1=ForwardModel()
62
63 J=InversionCostFunction(Regularization(), mappings=[m0, m1], forward_models=[(f0, 0), (f1,1)])
64
65 Example 3 (two forward models on 2-valued level set)
66 m0=Mapping()
67 m1=Mapping()
68 f0=ForwardModel()
69 f1=ForwardModel()
70
71 J=InversionCostFunction(Regularization(self.numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)])
72
73 :note: If provides_inverse_Hessian_approximation is true, then the class
74 provides an approximative inverse of the Hessian operator.
75 """
76 provides_inverse_Hessian_approximation=True
77
78 def __init__(self, regularization, mappings, forward_models):
79 """
80 constructor for the cost function.
81 Stores the supplied object references and sets default weights.
82
83 :param regularization: the regularization part of the cost function
84 :type regularization: `Regularization`
85 :param mappings: the mappings to calculate physical parameters from the
86 regularization. This is a list of 2-tuples *(map, i)*
87 where the first component map defines a `Mapping` and
88 the second component *i* defines the index of the
89 component of level set function to be used to
90 calculate the mapping. Items in the list may also be
91 just `Mapping` objects in which case the entire level
92 set function is fed into the `Mapping` (typically used
93 for a single-component level set function.
94 :type mappings: `Mapping` or ``list``
95 :param forward_models: the forward models involved in the calculation
96 of the cost function. This is a list of 2-tuples
97 *(f, ii)* where the first component f defines a
98 `ForwardModel` and the second component *ii* a
99 list of indexes referring to the physical
100 parameters in the `mappings` list. The 2-tuple
101 can be replaced by a `ForwardModel` if the
102 `mappings` list has a single entry.
103 :param forward_models: `ForwardModel` or ``list``
104 """
105 super(InversionCostFunction, self).__init__()
106 self.regularization=regularization
107 self.numLevelSets = self.regularization.getNumLevelSets()
108
109 if isinstance(mappings, Mapping):
110 mappings = [ mappings ]
111
112 self.mappings=[]
113 for i in range(len(mappings)):
114 mm=mappings[i]
115 if isinstance(mm, Mapping):
116 m=mm
117 if self.numLevelSets>1:
118 idx=[ p for p in range(self.numLevelSets)]
119 else:
120 idx=None
121 elif len(mm) == 1:
122 m=mm[0]
123 if self.numLevelSets>1:
124 idx=[ p for p in range(self.numLevelSets)]
125 else:
126 idx=None
127 else:
128 m=mm[0]
129 if isinstance(mm[1], int):
130 idx=[mm[1]]
131 else:
132 idx=list(mm[1])
133 if self.numLevelSets>1:
134 for k in idx:
135 if k < 0 or k > self.numLevelSets-1:
136 raise ValueError("level set index %s is out of range."%(k,))
137
138 else:
139 if idx[0] != 0:
140 raise ValueError("Level set index %s is out of range."%(idx[0],))
141 else:
142 idx=None
143 self.mappings.append((m,idx))
144 self.numMappings=len(self.mappings)
145
146
147 if isinstance(forward_models, ForwardModel):
148 forward_models = [ forward_models ]
149 self.forward_models=[]
150 for i in range(len(forward_models)):
151 f=forward_models[i]
152 if isinstance(f, ForwardModel):
153 idx=[0]
154 fm=f
155 elif len(f) == 1:
156 idx=[0]
157 fm=f[0]
158 else:
159 if isinstance(f[1],int):
160 idx=[f[1]]
161 else:
162 idx=list(f[1])
163 for k in idx:
164 if k<0 or k> self.numMappings:
165 raise ValueError("mapping index %s in model %s is out of range."%(k,i))
166 fm=f[0]
167 self.forward_models.append((fm,idx))
168 self.numModels=len(self.forward_models)
169
170 trafo = self.regularization.getCoordinateTransformation()
171 for m in self.forward_models:
172 if not m[0].getCoordinateTransformation() == trafo:
173 raise ValueError("Coordinate transformation for regularization and model don't match.")
174
175 self.__num_tradeoff_factors = self.regularization.getNumTradeOffFactors() + self.numModels
176 self.setTradeOffFactorsModels()
177
178 def getDomain(self):
179 """
180 returns the domain of the cost function
181
182 :rtype: `Domain`
183 """
184 self.regularization.getDomain()
185
186 def getNumTradeOffFactors(self):
187 """
188 returns the number of trade-off factors being used including the
189 trade-off factors used in the regularization component.
190
191 :rtype: ``int``
192 """
193 return self.__num_tradeoff_factors
194
195 def getForwardModel(self, idx=None):
196 """
197 returns the *idx*-th forward model.
198
199 :param idx: model index. If cost function contains one model only `idx`
200 can be omitted.
201 :type idx: ``int``
202 """
203 if idx==None: idx=0
204 return self.forward_models[idx][0]
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, dtype=float)
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 mp, idx=self.mappings[i]
296 m2=mp.getInverse(props[i])
297 if idx:
298 if len(idx) == 1:
299 m[idx[0]]=m2
300 else:
301 for k in range(idx): m[idx[k]]=m2[k]
302 else:
303 if isinstance(m2, Data):
304 m=interpolate(m2, m.getFunctionSpace())
305 else:
306 m=Data(m2, m.getFunctionSpace())
307 return m
308
309 def getProperties(self, m, return_list=False):
310 """
311 returns a list of the physical properties from a given level set
312 function *m* using the mappings of the cost function.
313
314 :param m: level set function
315 :type m: `Data`
316 :param return_list: if ``True`` a list is returned.
317 :type return_list: ``bool``
318 :rtype: ``list`` of `Data`
319 """
320
321 props=[]
322 for i in range(self.numMappings):
323 mp, idx=self.mappings[i]
324 if idx:
325 if len(idx)==1:
326 p=mp.getValue(m[idx[0]])
327 else:
328 m2=Data(0.,(len(idx),),m.getFunctionSpace())
329 for k in range(len(idx)): m2[k]=m[idx[k]]
330 p=mp.getValue(m2)
331 else:
332 p=mp.getValue(m)
333 props.append(p)
334 if self.numMappings > 1 or return_list:
335 return props
336 else:
337 return props[0]
338
339 def _getDualProduct(self, x, r):
340 """
341 Returns the dual product, see `Regularization.getDualProduct`
342
343 :type x: `Data`
344 :type r: `ArithmeticTuple`
345 :rtype: ``float``
346 """
347 return self.regularization.getDualProduct(x, r)
348
349 def _getArguments(self, m):
350 """
351 returns pre-computed values that are shared in the calculation of
352 *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
353 mapped value of ``m``, the arguments from the forward model and the
354 arguments from the regularization.
355
356 :param m: current approximation of the level set function
357 :type m: `Data`
358 :return: tuple of of values of the parameters, pre-computed values
359 for the forward model and pre-computed values for the
360 regularization
361 :rtype: ``tuple``
362 """
363 args_reg=self.regularization.getArguments(m)
364 # cache for physical parameters:
365 props=self.getProperties(m, return_list=True)
366 args_f=[]
367 for i in range(self.numModels):
368 f, idx=self.forward_models[i]
369 pp=tuple( [ props[k] for k in idx] )
370 aa=f.getArguments(*pp)
371 args_f.append(aa)
372
373 return props, args_f, args_reg
374
375 def _getValue(self, m, *args):
376 """
377 Returns the value *J(m)* of the cost function at *m*.
378 If the pre-computed values are not supplied `getArguments()` is called.
379
380 :param m: current approximation of the level set function
381 :type m: `Data`
382 :param args: tuple of values of the parameters, pre-computed values
383 for the forward model and pre-computed values for the
384 regularization
385 :rtype: ``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 = self.regularization.getValue(m, *args_reg)
395 self.logger.debug("J_R (incl. trade-offs) = %e"%J)
396
397 for i in range(self.numModels):
398 f, idx=self.forward_models[i]
399 args=tuple( [ props[k] for k in idx] + list( args_f[i] ) )
400 J_f = f.getDefect(*args)
401 self.logger.debug("J_f[%d] = %e, mu_model[%d] = %e"%(i, J_f, i, self.mu_model[i]))
402 J += self.mu_model[i] * J_f
403
404 return J
405
406 def getComponentValues(self, m, *args):
407 return self._getComponentValues(m, *args)
408
409 def _getComponentValues(self, m, *args):
410 """
411 returns the values of the individual cost functions that make up *f(x)*
412 using the precalculated values for *x*.
413
414 :param x: a solution approximation
415 :type x: x-type
416 :rtype: ``list<<float>>``
417 """
418 if len(args)==0:
419 args=self.getArguments(m)
420
421 props=args[0]
422 args_f=args[1]
423 args_reg=args[2]
424
425 J_reg = self.regularization.getValue(m, *args_reg)
426 result = [J_reg]
427
428 for i in range(self.numModels):
429 f, idx=self.forward_models[i]
430 args=tuple( [ props[k] for k in idx] + list( args_f[i] ) )
431 J_f = f.getValue(*args)
432 self.logger.debug("J_f[%d] = %e, mu_model[%d] = %e"%(i, J_f, i, self.mu_model[i]))
433
434 result += [J_f] # self.mu_model[i] * ??
435
436 return result
437
438 def _getGradient(self, m, *args):
439 """
440 returns the gradient of the cost function at *m*.
441 If the pre-computed values are not supplied `getArguments()` is called.
442
443 :param m: current approximation of the level set function
444 :type m: `Data`
445 :param args: tuple of values of the parameters, pre-computed values
446 for the forward model and pre-computed values for the
447 regularization
448
449 :rtype: `ArithmeticTuple`
450
451 :note: returns (Y^,X) where Y^ is the gradient from regularization plus
452 gradients of fwd models. X is the gradient of the regularization
453 w.r.t. gradient of m.
454 """
455 if len(args)==0:
456 args = self.getArguments(m)
457
458 props=args[0]
459 args_f=args[1]
460 args_reg=args[2]
461
462 g_J = self.regularization.getGradient(m, *args_reg)
463 p_diffs=[]
464 for i in range(self.numMappings):
465 mm, idx=self.mappings[i]
466 if idx and self.numLevelSets > 1:
467 if len(idx)>1:
468 m2=Data(0,(len(idx),),m.getFunctionSpace())
469 for k in range(len(idx)): m2[k]=m[idx[k]]
470 dpdm = mm.getDerivative(m2)
471 else:
472 dpdm = mm.getDerivative(m[idx[0]])
473 else:
474 dpdm = mm.getDerivative(m)
475 p_diffs.append(dpdm)
476
477 Y=g_J[0] # Because g_J==(Y,X) Y_k=dKer/dm_k
478 for i in range(self.numModels):
479 mu=self.mu_model[i]
480 f, idx_f=self.forward_models[i]
481 args=tuple( [ props[k] for k in idx_f] + list( args_f[i] ) )
482 Ys = f.getGradient(*args) # this d Jf/d props
483 # in this case f depends on one parameter props only but this can
484 # still depend on several level set components
485 if Ys.getRank() == 0:
486 # run through all level sets k prop j is depending on:
487 idx_m=self.mappings[idx_f[0]][1]
488 # tmp[k] = dJ_f/d_prop * d prop/d m[idx_m[k]]
489 tmp=Ys * p_diffs[idx_f[0]] * mu
490 if idx_m:
491 if tmp.getRank()== 0:
492 for k in range(len(idx_m)):
493 Y[idx_m[k]]+=tmp # dJ_f /d m[idx_m[k]] = tmp
494 else:
495 for k in range(len(idx_m)):
496 Y[idx_m[k]]+=tmp[k] # dJ_f /d m[idx_m[k]] = tmp[k]
497 else:
498 Y+=tmp # dJ_f /d m[idx_m[k]] = tmp
499 else:
500 s=0
501 # run through all props j forward model f is depending on:
502 for j in range(len(idx_f)):
503 # run through all level sets k prop j is depending on:
504 idx_m=self.mappings[j][1]
505 if p_diffs[idx_f[j]].getRank() == 0 :
506 if idx_m: # this case is not needed (really?)
507 self.logger.error("something wrong A")
508 # tmp[k] = dJ_f/d_prop[j] * d prop[j]/d m[idx_m[k]]
509 tmp=Ys[s]*p_diffs[idx_f[j]] * mu
510 for k in range(len(idx_m)):
511 Y[idx_m[k]]+=tmp[k] # dJ_f /d m[idx_m[k]] = tmp[k]
512 else:
513 Y+=Ys[s]*p_diffs[idx_f[j]] * mu
514 s+=1
515 elif p_diffs[idx_f[j]].getRank() == 1 :
516 l=p_diffs[idx_f[j]].getShape()[0]
517 # tmp[k]=sum_j dJ_f/d_prop[j] * d prop[j]/d m[idx_m[k]]
518 tmp=inner(Ys[s:s+l], p_diffs[idx_f[j]]) * mu
519 if idx_m:
520 for k in range(len(idx_m)):
521 Y[idx_m[k]]+=tmp # dJ_f /d m[idx_m[k]] = tmp[k]
522 else:
523 Y+=tmp
524 s+=l
525 else: # rank 2 case
526 l=p_diffs[idx_f[j]].getShape()[0]
527 Yss=Ys[s:s+l]
528 if idx_m:
529 for k in range(len(idx_m)):
530 # dJ_f /d m[idx_m[k]] = tmp[k]
531 Y[idx_m[k]]+=inner(Yss, p_diffs[idx_f[j]][:,k])
532 else:
533 Y+=inner(Yss, p_diffs[idx_f[j]]) * mu
534 s+=l
535 return g_J
536
537 def _getInverseHessianApproximation(self, m, r, *args):
538 """
539 returns an approximative evaluation *p* of the inverse of the Hessian
540 operator of the cost function for a given gradient type *r* at a
541 given location *m*: *H(m) p = r*
542
543 :param m: level set approximation where to calculate Hessian inverse
544 :type m: `Data`
545 :param r: a given gradient
546 :type r: `ArithmeticTuple`
547 :param args: tuple of values of the parameters, pre-computed values
548 for the forward model and pre-computed values for the
549 regularization
550 :rtype: `Data`
551 :note: in the current implementation only the regularization term is
552 considered in the inverse Hessian approximation.
553 """
554 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
555 return m
556
557 def updateHessian(self):
558 """
559 notifies the class that the Hessian operator needs to be updated.
560 """
561 self.regularization.updateHessian()
562
563 def _getNorm(self, m):
564 """
565 returns the norm of `m`
566
567 :param m: level set function
568 :type m: `Data`
569 :rtype: ``float``
570 """
571 return self.regularization.getNorm(m)
572

  ViewVC Help
Powered by ViewVC 1.1.26