/[escript]/trunk/doc/inversion/CostFunctions.tex
ViewVC logotype

Contents of /trunk/doc/inversion/CostFunctions.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4139 - (show annotations)
Tue Jan 15 00:13:28 2013 UTC (6 years, 7 months ago) by caltinay
File MIME type: application/x-tex
File size: 16467 byte(s)
Added L-BFGS reference and some minor changes.

1 \chapter{Cost Function}\label{chapter:ref:inversion cost function}
2 The general form of the cost function minimized in the inversion is given in the form (see also Chapter~\ref{chapter:ref:Drivers})
3 \begin{equation}\label{REF:EQU:DRIVE:10}
4 J(m) = J^{reg}(m) + \sum_{f} \mu^{data}_{f} \cdot J^{f}(p^f)
5 \end{equation}
6 where $m$ represents the level set function, $J^{reg}$ is the regularization term, see Chapter~\ref{Chp:ref:regularization},
7 and $J^{f}$ are a set of forward problems, see Chapter~\ref{Chp:ref:forward models} depending of
8 physical parameters $p^f$. The physical parameters $p^f$ are known functions
9 of the level set function $m$ which is the unknown to be calculated by the optimization process.
10 $\mu^{data}_{f}$ are trade-off factors. It is pointed out that the regularization term includes additional trade-off factors
11 The \class{InversionCostFunction} is class to define cost functions of an inversion. It is pointed out that
12 the \class{InversionCostFunction} class implements the \class{CostFunction} template class, see Chapter~\ref{chapter:ref:Minimization}.
13
14 In the simplest case there is a single forward model using a single physical parameter which is
15 derived form single-values level set function. The following script snippet shows the creation of the
16 \class{InversionCostFunction} for the case of a gravity inversion:
17 \begin{verbatim}
18 p=DensityMapping(...)
19 f=GravityModel(...)
20 J=InversionCostFunction(Regularization(...), \
21 mappings=p, \
22 forward_models=f)
23 \end{verbatim}
24 The argument \verb|...| refers to an appropriate argument list.
25
26 If two forward models are coming into play using two different physical parameters
27 the \member{mappings} and \member{forward_models} are defined as lists in the following form:
28 \begin{verbatim}
29 p_rho=DensityMapping(...)
30 p_k=SusceptibilityMapping(...)
31 f_mag=MagneticModel(...)
32 f_grav=GravityModel(...)
33
34 J=InversionCostFunction(Regularization(...), \
35 mappings=[p_rho, p_k], \
36 forward_models=[(f_mag, 1), (f_grav,0)])
37 \end{verbatim}
38 Here we define a joint inversion of gravity and magnetic data. \member{forward_models} is given as a list of
39 a tuple of a forward model and an index which referring to parameter in the \member{mappings} list to be used as an input.
40 The magnetic forward model \member{f_mag} is using the the second parameter (=\member{p_k}) in \member{mappings} list.
41 In this case the physical parameters are defined by a single-valued level set function. It is also possible
42 to link physical parameters to components of a level set function:
43 \begin{verbatim}
44 p_rho=DensityMapping(...)
45 p_k=SusceptibilityMapping(...)
46 f_mag=MagneticModel(...)
47 f_grav=GravityModel(...)
48
49 J=InversionCostFunction(Regularization(numLevelSets=2,...), \
50 mappings=[(p_rho,0), (p_k,1)], \
51 forward_models=[[(f_mag, 1), (f_grav,0)])
52 \end{verbatim}
53 The \member{mappings} argument is now a list of pairs where the first pair entry specifies the parameter mapping and
54 the second pair entry specifies the index of the component of the level set function to be used to evaluate the parameter.
55 In this case the level set function has two components, where the density mapping uses the first component of the level set function
56 while the susceptibility mapping uses the second component.
57
58 The \class{InversionCostFunction} API is defined as follows:
59
60 \begin{classdesc}{InversionCostFunction}{regularization, mappings, forward_models}
61 Constructor for the inversion cost function. \member{regularization} sets the regularization to be used, see Chapter~\ref{Chp:ref:regularization}.
62 \member{mappings} is a list of pairs where each pair comprises of a
63 physical parameter mapping (see Chapter~\ref{Chp:ref:mapping}) and an index which refers to the component of level set function
64 defined by the \member{regularization} to be used to calculate the corresponding physical parameter. If
65 the level set function has a single component the index can be omitted. If in addition there is a single physical parameter
66 the mapping can be given instead of a list. \member{forward_models} is a list of pairs where the
67 first pair component is a forward model ( see Chapter~\ref{Chp:ref:forward models}) and the second
68 pair component referring to the physical parameter in the \member{mappings} list providing the physical parameter for the model.
69 If a single physical parameter is present the index can be omitted. If in addition a single forward model is used this
70 forward model can be assigned to \member{forward_models} in replacement of a list.
71 \end{classdesc}
72
73
74 \begin{methoddesc}[InversionCostFunction]{getDomain}{}
75 """
76 returns the domain of the cost function
77 :rtype: 'Domain`
78 """
79 self.regularization.getDomain()
80 \end{methoddesc}
81
82 \begin{methoddesc}[InversionCostFunction]{getNumTradeOffFactors}{}
83 """
84 returns the number of trade-off factors being used including the
85 trade-off factors used in the regularization component.
86
87 :rtype: ``int``
88 """
89 return self.__num_tradeoff_factors
90 \end{methoddesc}
91
92 \begin{methoddesc}[InversionCostFunction]{getForwardModel}{idx=None}
93 """
94 returns the *idx*-th forward model.
95
96 :param idx: model index. If cost function contains one model only `idx`
97 can be omitted.
98 :type idx: ``int``
99 """
100 if idx==None: idx=0
101 f=self.forward_models[idx]
102 if isinstance(f, ForwardModel):
103 F=f
104 else:
105 F=f[0]
106 return F
107 \end{methoddesc}
108
109 \begin{methoddesc}[InversionCostFunction]{getRegularization}{}
110 """
111 returns the regularization
112 """
113 return self.regularization
114 \end{methoddesc}
115
116
117 \begin{methoddesc}[InversionCostFunction]{setTradeOffFactorsModels}{mu=None}
118 """
119 sets the trade-off factors for the forward model components.
120
121 :param mu: list of the trade-off factors. If not present ones are used.
122 :type mu: ``float`` in case of a single model or a ``list`` of ``float``
123 with the length of the number of models.
124 """
125 if mu==None:
126 self.mu_model=np.ones((self.numModels, ))
127 else:
128 if self.numModels > 1:
129 mu=np.asarray(mu)
130 if min(mu) > 0:
131 self.mu_model= mu
132 else:
133 raise ValueError("All value for trade-off factor mu must be positive.")
134 else:
135 mu=float(mu)
136 if mu > 0:
137 self.mu_model= [mu, ]
138 else:
139 raise ValueError("Trade-off factor must be positive.")
140 \end{methoddesc}
141
142 \begin{methoddesc}[InversionCostFunction]{setTradeOffFactorsRegularization}{mu=None, mu_c=None}
143 """
144 sets the trade of factors for the regularization component of the cost
145 function, see `Regularization` for details.
146
147 :param mu: trade-off factors for the level-set variation part
148 :param mu_c: trade-off factors for the cross gradient variation part
149 """
150 self.regularization.setTradeOffFactorsForVariation(mu)
151 self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
152 \end{methoddesc}
153
154 \begin{methoddesc}[InversionCostFunction]{setTradeOffFactors}{mu=None}
155 """
156 sets the trade-off factors for the forward model and regularization
157 terms.
158
159 :param mu: list of trade-off factors.
160 :type mu: ``list`` of ``float``
161 """
162 if mu is None:
163 mu=np.ones((self.__num_tradeoff_factors,))
164 self.setTradeOffFactorsModels(mu[:self.numModels])
165 self.regularization.setTradeOffFactors(mu[self.numModels:])
166 \end{methoddesc}
167
168 \begin{methoddesc}[InversionCostFunction]{createLevelSetFunction}{*props}
169 """
170 returns an instance of an object used to represent a level set function
171 initialized with zeros. Components can be overwritten by physical
172 properties 'props'. If present entries must correspond to the
173 `mappings` arguments in the constructor. Use `None` for properties for
174 which no value is given.
175 """
176 m=self.regularization.getPDE().createSolution()
177 if len(props) > 0:
178 for i in xrange(self.numMappings):
179 if props[i]:
180 mm=self.mappings[i]
181 if isinstance(mm, Mapping):
182 m=mm.getInverse(props[i])
183 elif len(mm) == 1:
184 m=mm[0].getInverse(props[i])
185 else:
186 m[mm[1]]=mm[0].getInverse(props[i])
187 return m
188 \end{methoddesc}
189
190 \begin{methoddesc}[InversionCostFunction]{getProperties}{m, return_list=False}
191 """
192 returns a list of the physical properties from a given level set
193 function *m* using the mappings of the cost function.
194
195 :param m: level set function
196 :type m: `Data`
197 :param return_list: if True a list is returned.
198 def _:type return_list: `bool`
199 :rtype: `list` of `Data`
200 """
201 props=[]
202 for i in xrange(self.numMappings):
203 mm=self.mappings[i]
204 if isinstance(mm, Mapping):
205 p=mm.getValue(m)
206 elif len(mm) == 1:
207 p=mm[0].getValue(m)
208 else:
209 p=mm[0].getValue(m[mm[1]])
210 props.append(p)
211 if self.numMappings > 1 or return_list:
212 return props
213 else:
214 return props[0]
215 \end{methoddesc}
216
217 \begin{methoddesc}[InversionCostFunction]{getDualProduct}{x, r}
218 """
219 Returns the dual product, see `Regularization.getDualProduct`
220
221 :type x: `Data`
222 :type r: `ArithmeticTuple`
223 :rtype: `float`
224 """
225 return self.regularization.getDualProduct(x, r)
226 \end{methoddesc}
227
228 \begin{methoddesc}[InversionCostFunction]{getArguments}{m}
229 """
230 returns pre-computed values that are shared in the calculation of
231 *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
232 mapped value of ``m``, the arguments from the forward model and the
233 arguments from the regularization.
234
235 :param m: current approximation of the level set function
236 :type m: `Data`
237 :return: tuple of of values of the parameters, pre-computed values for the forward model and
238 pre-computed values for the regularization
239 :rtype: `tuple`
240 """
241 args_reg=self.regularization.getArguments(m)
242
243 props=self.getProperties(m, return_list=True)
244 args_f=[]
245 for i in xrange(self.numModels):
246 f=self.forward_models[i]
247 if isinstance(f, ForwardModel):
248 aa=f.getArguments(props[0])
249 elif len(f) == 1:
250 aa=f[0].getArguments(props[0])
251 else:
252 idx = f[1]
253 f=f[0]
254 if isinstance(idx, int):
255 aa=f.getArguments(props[idx])
256 else:
257 pp=tuple( [ props[i] for i in idx] )
258 aa=f.getArguments(*pp)
259 args_f.append(aa)
260
261 return props, args_f, args_reg
262 \end{methoddesc}
263
264 \begin{methoddesc}[InversionCostFunction]{getValue}{m, *args}
265 """
266 Returns the value *J(m)* of the cost function at *m*.
267 If the pre-computed values are not supplied `getArguments()` is called.
268
269 :param m: current approximation of the level set function
270 :type m: `Data`
271 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
272 pre-computed values for the regularization
273 :rtype: `float`
274 """
275
276 if len(args)==0:
277 args=self.getArguments(m)
278
279 props=args[0]
280 args_f=args[1]
281 args_reg=args[2]
282
283 J = self.regularization.getValue(m, *args_reg)
284 print "J_reg = %e"%J
285
286 for i in xrange(self.numModels):
287
288 f=self.forward_models[i]
289 if isinstance(f, ForwardModel):
290 J_f = f.getValue(props[0],*args_f[i])
291 elif len(f) == 1:
292 J_f=f[0].getValue(props[0],*args_f[i])
293 else:
294 idx = f[1]
295 f=f[0]
296 if isinstance(idx, int):
297 J_f = f.getValue(props[idx],*args_f[i])
298 else:
299 args=tuple( [ props[j] for j in idx] + args_f[i])
300 J_f = f.getValue(*args)
301 print "J_f[%d] = %e"%(i, J_f)
302 print "mu_model[%d] = %e"%(i, self.mu_model[i])
303 J += self.mu_model[i] * J_f
304
305 return J
306 \end{methoddesc}
307
308 \begin{methoddesc}[InversionCostFunction]{getGradient}{m, *args}
309 """
310 returns the gradient of the cost function at *m*.
311 If the pre-computed values are not supplied `getArguments()` is called.
312
313 :param m: current approximation of the level set function
314 :type m: `Data`
315 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
316 pre-computed values for the regularization
317
318 :rtype: `ArithmeticTuple`
319 """
320 if len(args)==0:
321 args = self.getArguments(m)
322
323 props=args[0]
324 args_f=args[1]
325 args_reg=args[2]
326
327 g_J = self.regularization.getGradient(m, *args_reg)
328 p_diffs=[]
329 for i in xrange(self.numMappings):
330 mm=self.mappings[i]
331 if isinstance(mm, Mapping):
332 dpdm = mm.getDerivative(m)
333 elif len(mm) == 1:
334 dpdm = mm[0].getDerivative(m)
335 else:
336 dpdm = mm[0].getDerivative(m[mm[1]])
337 p_diffs.append(dpdm)
338
339 Y=g_J[0]
340 for i in xrange(self.numModels):
341 mu=self.mu_model[i]
342 f=self.forward_models[i]
343 if isinstance(f, ForwardModel):
344 Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
345 if self.numLevelSets == 1 :
346 Y +=Ys
347 else:
348 Y[0] +=Ys
349 elif len(f) == 1:
350 Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
351 if self.numLevelSets == 1 :
352 Y +=Ys
353 else:
354 Y[0] +=Ys
355 else:
356 idx = f[1]
357 f=f[0]
358 if isinstance(idx, int):
359 Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu
360 if self.numLevelSets == 1 :
361 if idx == 0:
362 Y+=Ys
363 else:
364 raise IndexError("Illegal mapping index.")
365 else:
366 Y[idx] += Ys
367 else:
368 args=tuple( [ props[j] for j in idx] + args_f[i])
369 Ys = f.getGradient(*args)
370 for ii in xrange(len(idx)):
371 Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]] * mu
372
373 return g_J
374 \end{methoddesc}
375
376
377 \begin{methoddesc}[InversionCostFunction]{getInverseHessianApproximation}{m, r, *args}
378 """
379 returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function
380 for a given gradient type *r* at a given location *m*: *H(m) p = r*
381
382 :param m: level set approximation where to calculate Hessian inverse
383 :type m: `Data`
384 :param r: a given gradient
385 :type r: `ArithmeticTuple`
386 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
387 pre-computed values for the regularization
388 :rtype: `Data`
389 :note: in the current implementation only the regularization term is
390 considered in the inverse Hessian approximation.
391
392 """
393 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
394 return m
395
396 \end{methoddesc}
397
398 \begin{methoddesc}[InversionCostFunction]{getNorm}{m}
399 """
400 returns the norm of ``m``
401
402 :param m: level set function
403 :type m: `Data`
404 :rtype: ``float``
405 """
406
407 \end{methoddesc}
408

  ViewVC Help
Powered by ViewVC 1.1.26