/[escript]/branches/subworld2/downunder/py_src/inversioncostfunctions.py
ViewVC logotype

Contents of /branches/subworld2/downunder/py_src/inversioncostfunctions.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26