/[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 4154 - (show annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 17636 byte(s)
Round 1 of copyright fixes
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 """Collection of cost functions for inversion"""
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 class InversionCostFunction(MeteredCostFunction):
35 """
36 Class to define cost function *J(m)* for inversion with one or more forward model
37 based on a multi-valued level set function *m*:
38
39 *J(m) = J_reg(m) + sum_f mu_f * J_f(p)*
40
41 where *J_reg(m)* is the regularization and cross gradient component of the cost function applied
42 to a level set function *m*, *J_f(p)* are the data defect cost function involving a
43 physical forward model using the physical parameter(s) *p* and mu_f is trade-off factors for model f.
44
45
46 A forward model depends on a set of physical parameters *p* which are constructed from
47 components of the level set *m* function via mappings.
48
49 Example 1 (single forward model):
50 m=Mapping()
51 f=ForwardModel()
52 J=InversionCostFunction(Regularization(), m, f)
53
54 Example 2 (two forward models on a single valued level set)
55 m0=Mapping()
56 m1=Mapping()
57 f0=ForwardModel()
58 f1=ForwardModel()
59
60 J=InversionCostFunction(Regularization(), mappings=[m0, m1], forward_models=[(f0, 0), (f1,1)])
61
62 Example 2 (two forward models on 2-valued level set)
63 m0=Mapping()
64 m1=Mapping()
65 f0=ForwardModel()
66 f1=ForwardModel()
67
68 J=InversionCostFunction(Regularization(numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)])
69
70 :var provides_inverse_Hessian_approximation: if true the class provides an
71 approximative inverse of the Hessian operator.
72 """
73 provides_inverse_Hessian_approximation=True
74
75 def __init__(self, regularization, mappings, forward_models):
76 """
77 constructor for the cost function.
78 Stores the supplied object references and sets default weights.
79
80 :param regularization: the regularization part of the cost function
81 :type regularization: `Regularization`
82 :param mappings: the mappings to calculate physical parameters from the
83 regularization. This is a list of 2-tuples *(map, i)*
84 where the first component map defines a `Mapping` and
85 the second component *i* defines the index of the
86 component of level set function to be used to
87 calculate the mapping. An item in the list can be just
88 a `Mapping` object then the entire level set function
89 function is fed into the `Mapping` (typically used for
90 a single-component level set function.
91 :type mappings: ``list`` where each item is a ``tuple`` of `Mapping`
92 and ``int`` or just a `Mapping`
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 a
100 `mappings` list as a single entry.
101 :param forward_models: ``list`` where each item is ``tuple`` of
102 `ForwardModel` and ``list`` of ``int`` or is
103 `ForwardModel`.
104 """
105 super(InversionCostFunction, self).__init__()
106 self.regularization=regularization
107
108 if isinstance(mappings, Mapping):
109 self.mappings = [mappings ]
110 else:
111 self.mappings = mappings
112
113 if isinstance(forward_models, ForwardModel):
114 self.forward_models = [ forward_models ]
115 else:
116 self.forward_models=forward_models
117
118 self.numMappings=len(self.mappings)
119 self.numModels=len(self.forward_models)
120 self.numLevelSets = self.regularization.getNumLevelSets()
121 self.__num_tradeoff_factors = self.regularization.getNumTradeOffFactors() + self.numModels
122 self.setTradeOffFactorsModels()
123
124 def getDomain(self):
125 """
126 returns the domain of the cost function
127 :rtype: 'Domain`
128 """
129 self.regularization.getDomain()
130
131 def getNumTradeOffFactors(self):
132 """
133 returns the number of trade-off factors being used including the
134 trade-off factors used in the regularization component.
135
136 :rtype: ``int``
137 """
138 return self.__num_tradeoff_factors
139
140 def getForwardModel(self, idx=None):
141 """
142 returns the *idx*-th forward model.
143
144 :param idx: model index. If cost function contains one model only `idx`
145 can be omitted.
146 :type idx: ``int``
147 """
148 if idx==None: idx=0
149 f=self.forward_models[idx]
150 if isinstance(f, ForwardModel):
151 F=f
152 else:
153 F=f[0]
154 return F
155
156 def getRegularization(self):
157 """
158 returns the regularization
159
160 :rtype: `Regularization`
161 """
162 return self.regularization
163
164
165 def setTradeOffFactorsModels(self,mu=None):
166 """
167 sets the trade-off factors for the forward model components.
168
169 :param mu: list of the trade-off factors. If not present ones are used.
170 :type mu: ``float`` in case of a single model or a ``list`` of ``float``
171 with the length of the number of models.
172 """
173 if mu==None:
174 self.mu_model=np.ones((self.numModels, ))
175 else:
176 if self.numModels > 1:
177 mu=np.asarray(mu)
178 if min(mu) > 0:
179 self.mu_model= mu
180 else:
181 raise ValueError("All value for trade-off factor mu must be positive.")
182 else:
183 mu=float(mu)
184 if mu > 0:
185 self.mu_model= [mu, ]
186 else:
187 raise ValueError("Trade-off factor must be positive.")
188
189 def getTradeOffFactorsModels(self):
190 """
191 returns the trade-off factors for the forward models
192
193 :rtype: ``float`` or ``list`` of ``float``
194 """
195 if self.numModels>1:
196 return self.mu_model
197 else:
198 return self.mu_model[0]
199
200 def setTradeOffFactorsRegularization(self,mu=None, mu_c=None):
201 """
202 sets the trade of factors for the regularization component of the cost
203 function, see `Regularization` for details.
204
205 :param mu: trade-off factors for the level-set variation part
206 :param mu_c: trade-off factors for the cross gradient variation part
207 """
208 self.regularization.setTradeOffFactorsForVariation(mu)
209 self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
210
211 def setTradeOffFactors(self, mu=None):
212 """
213 sets the trade-off factors for the forward model and regularization
214 terms.
215
216 :param mu: list of trade-off factors.
217 :type mu: ``list`` of ``float``
218 """
219 if mu is None:
220 mu=np.ones((self.__num_tradeoff_factors,))
221 self.setTradeOffFactorsModels(mu[:self.numModels])
222 self.regularization.setTradeOffFactors(mu[self.numModels:])
223
224 def getTradeOffFactors(self, mu=None):
225 """
226 returns a list of the trade-off factors
227 :rtype: ``list`` of ``float``
228 """
229 mu1=self.getTradeOffFactorsModels(mu[:self.numModels])
230 mu2=self.regularization.getTradeOffFactors()
231 return [ m for m in mu1] + [ m for m in mu2]
232
233
234 def createLevelSetFunction(self, *props):
235 """
236 returns an instance of an object used to represent a level set function
237 initialized with zeros. Components can be overwritten by physical
238 properties 'props'. If present entries must correspond to the
239 `mappings` arguments in the constructor. Use `None` for properties for
240 which no value is given.
241 """
242 m=self.regularization.getPDE().createSolution()
243 if len(props) > 0:
244 for i in xrange(self.numMappings):
245 if props[i]:
246 mm=self.mappings[i]
247 if isinstance(mm, Mapping):
248 m=mm.getInverse(props[i])
249 elif len(mm) == 1:
250 m=mm[0].getInverse(props[i])
251 else:
252 m[mm[1]]=mm[0].getInverse(props[i])
253 return m
254
255 def getProperties(self, m, return_list=False):
256 """
257 returns a list of the physical properties from a given level set
258 function *m* using the mappings of the cost function.
259
260 :param m: level set function
261 :type m: `Data`
262 :param return_list: if True a list is returned.
263 :type return_list: `bool`
264 :rtype: `list` of `Data`
265 """
266 props=[]
267 for i in xrange(self.numMappings):
268 mm=self.mappings[i]
269 if isinstance(mm, Mapping):
270 p=mm.getValue(m)
271 elif len(mm) == 1:
272 p=mm[0].getValue(m)
273 else:
274 p=mm[0].getValue(m[mm[1]])
275 props.append(p)
276 if self.numMappings > 1 or return_list:
277 return props
278 else:
279 return props[0]
280
281 def _getDualProduct(self, x, r):
282 """
283 Returns the dual product, see `Regularization.getDualProduct`
284
285 :type x: `Data`
286 :type r: `ArithmeticTuple`
287 :rtype: `float`
288 """
289 return self.regularization.getDualProduct(x, r)
290
291 def _getArguments(self, m):
292 """
293 returns pre-computed values that are shared in the calculation of
294 *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
295 mapped value of ``m``, the arguments from the forward model and the
296 arguments from the regularization.
297
298 :param m: current approximation of the level set function
299 :type m: `Data`
300 :return: tuple of of values of the parameters, pre-computed values for the forward model and
301 pre-computed values for the regularization
302 :rtype: `tuple`
303 """
304 args_reg=self.regularization.getArguments(m)
305 # cache for physical parameters:
306 props=self.getProperties(m, return_list=True)
307 args_f=[]
308 for i in xrange(self.numModels):
309 f=self.forward_models[i]
310 if isinstance(f, ForwardModel):
311 aa=f.getArguments(props[0])
312 elif len(f) == 1:
313 aa=f[0].getArguments(props[0])
314 else:
315 idx = f[1]
316 f=f[0]
317 if isinstance(idx, int):
318 aa=f.getArguments(props[idx])
319 else:
320 pp=tuple( [ props[i] for i in idx] )
321 aa=f.getArguments(*pp)
322 args_f.append(aa)
323
324 return props, args_f, args_reg
325
326 def _getValue(self, m, *args):
327 """
328 Returns the value *J(m)* of the cost function at *m*.
329 If the pre-computed values are not supplied `getArguments()` is called.
330
331 :param m: current approximation of the level set function
332 :type m: `Data`
333 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
334 pre-computed values for the regularization
335 :rtype: `float`
336 """
337 # if there is more than one forward_model and/or regularization their
338 # contributions need to be added up. But this implementation allows
339 # only one of each...
340 if len(args)==0:
341 args=self.getArguments(m)
342
343 props=args[0]
344 args_f=args[1]
345 args_reg=args[2]
346
347 J = self.regularization.getValue(m, *args_reg)
348 print "J_reg = %e"%J
349
350 for i in xrange(self.numModels):
351
352 f=self.forward_models[i]
353 if isinstance(f, ForwardModel):
354 J_f = f.getValue(props[0],*args_f[i])
355 elif len(f) == 1:
356 J_f=f[0].getValue(props[0],*args_f[i])
357 else:
358 idx = f[1]
359 f=f[0]
360 if isinstance(idx, int):
361 J_f = f.getValue(props[idx],*args_f[i])
362 else:
363 args=tuple( [ props[j] for j in idx] + args_f[i])
364 J_f = f.getValue(*args)
365 print "J_f[%d] = %e"%(i, J_f)
366 print "mu_model[%d] = %e"%(i, self.mu_model[i])
367 J += self.mu_model[i] * J_f
368
369 return J
370
371 def _getGradient(self, m, *args):
372 """
373 returns the gradient of the cost function at *m*.
374 If the pre-computed values are not supplied `getArguments()` is called.
375
376 :param m: current approximation of the level set function
377 :type m: `Data`
378 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
379 pre-computed values for the regularization
380
381 :rtype: `ArithmeticTuple`
382 """
383 if len(args)==0:
384 args = self.getArguments(m)
385
386 props=args[0]
387 args_f=args[1]
388 args_reg=args[2]
389
390 g_J = self.regularization.getGradient(m, *args_reg)
391 p_diffs=[]
392 for i in xrange(self.numMappings):
393 mm=self.mappings[i]
394 if isinstance(mm, Mapping):
395 dpdm = mm.getDerivative(m)
396 elif len(mm) == 1:
397 dpdm = mm[0].getDerivative(m)
398 else:
399 dpdm = mm[0].getDerivative(m[mm[1]])
400 p_diffs.append(dpdm)
401
402 Y=g_J[0]
403 for i in xrange(self.numModels):
404 mu=self.mu_model[i]
405 f=self.forward_models[i]
406 if isinstance(f, ForwardModel):
407 Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
408 if self.numLevelSets == 1 :
409 Y +=Ys
410 else:
411 Y[0] +=Ys
412 elif len(f) == 1:
413 Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
414 if self.numLevelSets == 1 :
415 Y +=Ys
416 else:
417 Y[0] +=Ys
418 else:
419 idx = f[1]
420 f=f[0]
421 if isinstance(idx, int):
422 Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu
423 if self.numLevelSets == 1 :
424 if idx == 0:
425 Y+=Ys
426 else:
427 raise IndexError("Illegal mapping index.")
428 else:
429 Y[idx] += Ys
430 else:
431 args=tuple( [ props[j] for j in idx] + args_f[i])
432 Ys = f.getGradient(*args)
433 for ii in xrange(len(idx)):
434 Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]] * mu
435
436 return g_J
437
438
439 def _getInverseHessianApproximation(self, m, r, *args):
440 """
441 returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function
442 for a given gradient type *r* at a given location *m*: *H(m) p = r*
443
444 :param m: level set approximation where to calculate Hessian inverse
445 :type m: `Data`
446 :param r: a given gradient
447 :type r: `ArithmeticTuple`
448 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
449 pre-computed values for the regularization
450 :rtype: `Data`
451 :note: in the current implementation only the regularization term is
452 considered in the inverse Hessian approximation.
453
454 """
455 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
456 return m
457
458 def updateHessian(self):
459 """
460 notifies the class that the Hessian operator needs to be updated.
461 """
462 self.regularization.updateHessian()
463
464 def _getNorm(self, m):
465 """
466 returns the norm of ``m``
467
468 :param m: level set function
469 :type m: `Data`
470 :rtype: ``float``
471 """
472 return self.regularization.getNorm(m)

  ViewVC Help
Powered by ViewVC 1.1.26