/[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 4132 - (show annotations)
Fri Jan 11 05:46:49 2013 UTC (6 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 16985 byte(s)
Range of epydoc fixes.

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 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-2012 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 return self.regularization
161
162
163 def setTradeOffFactorsModels(self,mu=None):
164 """
165 sets the trade-off factors for the forward model components.
166
167 :param mu: list of the trade-off factors. If not present ones are used.
168 :type mu: ``float`` in case of a single model or a ``list`` of ``float``
169 with the length of the number of models.
170 """
171 if mu==None:
172 self.mu_model=np.ones((self.numModels, ))
173 else:
174 if self.numModels > 1:
175 mu=np.asarray(mu)
176 if min(mu) > 0:
177 self.mu_model= mu
178 else:
179 raise ValueError("All value for trade-off factor mu must be positive.")
180 else:
181 mu=float(mu)
182 if mu > 0:
183 self.mu_model= [mu, ]
184 else:
185 raise ValueError("Trade-off factor must be positive.")
186
187 def setTradeOffFactorsRegularization(self,mu=None, mu_c=None):
188 """
189 sets the trade of factors for the regularization component of the cost
190 function, see `Regularization` for details.
191
192 :param mu: trade-off factors for the level-set variation part
193 :param mu_c: trade-off factors for the cross gradient variation part
194 """
195 self.regularization.setTradeOffFactorsForVariation(mu)
196 self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
197
198 def setTradeOffFactors(self, mu=None):
199 """
200 sets the trade-off factors for the forward model and regularization
201 terms.
202
203 :param mu: list of trade-off factors.
204 :type mu: ``list`` of ``float``
205 """
206 if mu is None:
207 mu=np.ones((self.__num_tradeoff_factors,))
208 self.setTradeOffFactorsModels(mu[:self.numModels])
209 self.regularization.setTradeOffFactors(mu[self.numModels:])
210
211 def createLevelSetFunction(self, *props):
212 """
213 returns an instance of an object used to represent a level set function
214 initialized with zeros. Components can be overwritten by physical
215 properties 'props'. If present entries must correspond to the
216 `mappings` arguments in the constructor. Use `None` for properties for
217 which no value is given.
218 """
219 m=self.regularization.getPDE().createSolution()
220 if len(props) > 0:
221 for i in xrange(self.numMappings):
222 if props[i]:
223 mm=self.mappings[i]
224 if isinstance(mm, Mapping):
225 m=mm.getInverse(props[i])
226 elif len(mm) == 1:
227 m=mm[0].getInverse(props[i])
228 else:
229 m[mm[1]]=mm[0].getInverse(props[i])
230 return m
231
232 def getProperties(self, m, return_list=False):
233 """
234 returns a list of the physical properties from a given level set
235 function *m* using the mappings of the cost function.
236
237 :param m: level set function
238 :type m: `Data`
239 :param return_list: if True a list is returned.
240 :type return_list: `bool`
241 :rtype: `list` of `Data`
242 """
243 props=[]
244 for i in xrange(self.numMappings):
245 mm=self.mappings[i]
246 if isinstance(mm, Mapping):
247 p=mm.getValue(m)
248 elif len(mm) == 1:
249 p=mm[0].getValue(m)
250 else:
251 p=mm[0].getValue(m[mm[1]])
252 props.append(p)
253 if self.numMappings > 1 or return_list:
254 return props
255 else:
256 return props[0]
257
258 def _getDualProduct(self, x, r):
259 """
260 Returns the dual product, see `Regularization.getDualProduct`
261
262 :type x: `Data`
263 :type r: `ArithmeticTuple`
264 :rtype: `float`
265 """
266 return self.regularization.getDualProduct(x, r)
267
268 def _getArguments(self, m):
269 """
270 returns pre-computed values that are shared in the calculation of
271 *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
272 mapped value of ``m``, the arguments from the forward model and the
273 arguments from the regularization.
274
275 :param m: current approximation of the level set function
276 :type m: `Data`
277 :return: tuple of of values of the parameters, pre-computed values for the forward model and
278 pre-computed values for the regularization
279 :rtype: `tuple`
280 """
281 args_reg=self.regularization.getArguments(m)
282 # cache for physical parameters:
283 props=self.getProperties(m, return_list=True)
284 args_f=[]
285 for i in xrange(self.numModels):
286 f=self.forward_models[i]
287 if isinstance(f, ForwardModel):
288 aa=f.getArguments(props[0])
289 elif len(f) == 1:
290 aa=f[0].getArguments(props[0])
291 else:
292 idx = f[1]
293 f=f[0]
294 if isinstance(idx, int):
295 aa=f.getArguments(props[idx])
296 else:
297 pp=tuple( [ props[i] for i in idx] )
298 aa=f.getArguments(*pp)
299 args_f.append(aa)
300
301 return props, args_f, args_reg
302
303 def _getValue(self, m, *args):
304 """
305 Returns the value *J(m)* of the cost function at *m*.
306 If the pre-computed values are not supplied `getArguments()` is called.
307
308 :param m: current approximation of the level set function
309 :type m: `Data`
310 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
311 pre-computed values for the regularization
312 :rtype: `float`
313 """
314 # if there is more than one forward_model and/or regularization their
315 # contributions need to be added up. But this implementation allows
316 # only one of each...
317 if len(args)==0:
318 args=self.getArguments(m)
319
320 props=args[0]
321 args_f=args[1]
322 args_reg=args[2]
323
324 J = self.regularization.getValue(m, *args_reg)
325 print "J_reg = %e"%J
326
327 for i in xrange(self.numModels):
328
329 f=self.forward_models[i]
330 if isinstance(f, ForwardModel):
331 J_f = f.getValue(props[0],*args_f[i])
332 elif len(f) == 1:
333 J_f=f[0].getValue(props[0],*args_f[i])
334 else:
335 idx = f[1]
336 f=f[0]
337 if isinstance(idx, int):
338 J_f = f.getValue(props[idx],*args_f[i])
339 else:
340 args=tuple( [ props[j] for j in idx] + args_f[i])
341 J_f = f.getValue(*args)
342 print "J_f[%d] = %e"%(i, J_f)
343 print "mu_model[%d] = %e"%(i, self.mu_model[i])
344 J += self.mu_model[i] * J_f
345
346 return J
347
348 def _getGradient(self, m, *args):
349 """
350 returns the gradient of the cost function at *m*.
351 If the pre-computed values are not supplied `getArguments()` is called.
352
353 :param m: current approximation of the level set function
354 :type m: `Data`
355 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
356 pre-computed values for the regularization
357
358 :rtype: `ArithmeticTuple`
359 """
360 if len(args)==0:
361 args = self.getArguments(m)
362
363 props=args[0]
364 args_f=args[1]
365 args_reg=args[2]
366
367 g_J = self.regularization.getGradient(m, *args_reg)
368 p_diffs=[]
369 for i in xrange(self.numMappings):
370 mm=self.mappings[i]
371 if isinstance(mm, Mapping):
372 dpdm = mm.getDerivative(m)
373 elif len(mm) == 1:
374 dpdm = mm[0].getDerivative(m)
375 else:
376 dpdm = mm[0].getDerivative(m[mm[1]])
377 p_diffs.append(dpdm)
378
379 Y=g_J[0]
380 for i in xrange(self.numModels):
381 mu=self.mu_model[i]
382 f=self.forward_models[i]
383 if isinstance(f, ForwardModel):
384 Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
385 if self.numLevelSets == 1 :
386 Y +=Ys
387 else:
388 Y[0] +=Ys
389 elif len(f) == 1:
390 Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
391 if self.numLevelSets == 1 :
392 Y +=Ys
393 else:
394 Y[0] +=Ys
395 else:
396 idx = f[1]
397 f=f[0]
398 if isinstance(idx, int):
399 Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu
400 if self.numLevelSets == 1 :
401 if idx == 0:
402 Y+=Ys
403 else:
404 raise IndexError("Illegal mapping index.")
405 else:
406 Y[idx] += Ys
407 else:
408 args=tuple( [ props[j] for j in idx] + args_f[i])
409 Ys = f.getGradient(*args)
410 for ii in xrange(len(idx)):
411 Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]] * mu
412
413 return g_J
414
415
416 def _getInverseHessianApproximation(self, m, r, *args):
417 """
418 returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function
419 for a given gradient type *r* at a given location *m*: *H(m) p = r*
420
421 :param m: level set approximation where to calculate Hessian inverse
422 :type m: `Data`
423 :param r: a given gradient
424 :type r: `ArithmeticTuple`
425 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
426 pre-computed values for the regularization
427 :rtype: `Data`
428 :note: in the current implementation only the regularization term is
429 considered in the inverse Hessian approximation.
430
431 """
432 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
433 return m
434
435 def updateHessian(self):
436 """
437 notifies the class that the Hessian operator needs to be updated.
438 """
439 self.regularization.updateHessian()
440
441 def _getNorm(self, m):
442 """
443 returns the norm of ``m``
444
445 :param m: level set function
446 :type m: `Data`
447 :rtype: ``float``
448 """
449 return self.regularization.getNorm(m)

  ViewVC Help
Powered by ViewVC 1.1.26