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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4138 - (hide annotations)
Mon Jan 14 23:57:00 2013 UTC (6 years, 7 months ago) by gross
File MIME type: application/x-tex
File size: 16467 byte(s)
some rearrangement of part II
1 gross 4138 \chapter{Cost Function}\label{chapter:ref:inversion cost function}
2 gross 4133 The general form of the cost function minimized in the inversion is given in the form (see also Chapter~\ref{Chp:ref:introduction})
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 gross 4131
14 gross 4133 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 gross 4131
26 gross 4133 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 gross 4119 \end{classdesc}
72 gross 4133
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}

  ViewVC Help
Powered by ViewVC 1.1.26