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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4132 by gross, Fri Jan 11 03:54:16 2013 UTC revision 4133 by gross, Fri Jan 11 06:41:10 2013 UTC
# Line 14  template for inversion drivers. By defau Line 14  template for inversion drivers. By defau
14  \end{classdesc*}  \end{classdesc*}
15    
16  \begin{methoddesc}[InversionDriver]{getCostFunction}{}  \begin{methoddesc}[InversionDriver]{getCostFunction}{}
17  returns the cost function of the inversion. This will be an instance of the \class{InversionCostFunction} class, see section~\ref{XXX}.  returns the cost function of the inversion. This will be an instance of the \class{InversionCostFunction} class, see section~\ref{chapter:ref:inversion cost function}.
18  Use this method to access or alter attribute or methods of the underlying cost function.  Use this method to access or alter attribute or methods of the underlying cost function.
19  \end{methoddesc}  \end{methoddesc}
20    
# Line 95  For examples of usage please see Chapter Line 95  For examples of usage please see Chapter
95  Driver class to perform an inversion of magnetic anomaly data. This class  Driver class to perform an inversion of magnetic anomaly data. This class
96  is a sub-class of the \class{InversionDriver} class. The class uses the standard  is a sub-class of the \class{InversionDriver} class. The class uses the standard
97  \class{Regularization} class for a single level set function, see Chapter~\ref{Chp:ref:regularization},  \class{Regularization} class for a single level set function, see Chapter~\ref{Chp:ref:regularization},
98  \class{SusceptibilityMapping} mapping, see Section~\label{Chp:ref:mapping susceptibility}, and the linear  \class{SusceptibilityMapping} mapping, see Section~\ref{Chp:ref:mapping susceptibility}, and the linear
99  magnetic forward model \class{MagneticModel}, see Section~\ref{sec:forward magnetic}.  magnetic forward model \class{MagneticModel}, see Section~\ref{sec:forward magnetic}.
100  \end{classdesc}  \end{classdesc}
101    
# Line 180  set initial guesses for density and susc Line 180  set initial guesses for density and susc
180    
181    
182    
183  \section{Inversion Cost Function}  \section{Cost Function for Inversion}\label{chapter:ref:inversion cost function}
184    The general form of the cost function minimized in the inversion is given in the form (see also Chapter~\ref{Chp:ref:introduction})
185    \begin{equation}\label{REF:EQU:DRIVE:10}
186    J(m) = J^{reg}(m) + \sum_{f} \mu^{data}_{f} \cdot J^{f}(p^f)
187    \end{equation}
188    where $m$ represents the level set function, $J^{reg}$ is the regularization term, see Chapter~\ref{Chp:ref:regularization},
189    and $J^{f}$ are a set of forward problems, see Chapter~\ref{Chp:ref:forward models} depending of
190    physical parameters $p^f$.  The physical parameters $p^f$ are known functions
191    of the  level set function $m$ which is the unknown to be calculated by the optimization process.
192    $\mu^{data}_{f}$ are trade-off factors. It is pointed out that the regularization term includes additional trade-off factors
193    The \class{InversionCostFunction} is class to define cost functions of an inversion. It is pointed out that
194    the \class{InversionCostFunction} class implements the \class{CostFunction} template class, see Chapter~\ref{chapter:ref:Minimization}.
195    
196    In the simplest case there is a single forward model using a single physical parameter which is
197    derived form single-values level set function. The following script snippet shows the creation of the
198    \class{InversionCostFunction} for the case of a gravity inversion:
199    \begin{verbatim}
200    p=DensityMapping(...)
201    f=GravityModel(...)
202    J=InversionCostFunction(Regularization(...), \
203                            mappings=p, \
204                            forward_models=f)
205    \end{verbatim}
206    The argument \verb|...| refers to an appropriate argument list.
207    
208    If two forward models are coming into play using two different physical parameters
209    the \member{mappings} and \member{forward_models} are defined as lists in the following form:
210    \begin{verbatim}
211    p_rho=DensityMapping(...)
212    p_k=SusceptibilityMapping(...)
213    f_mag=MagneticModel(...)
214    f_grav=GravityModel(...)
215    
216    J=InversionCostFunction(Regularization(...), \
217                            mappings=[p_rho, p_k], \
218                            forward_models=[(f_mag, 1), (f_grav,0)])
219    \end{verbatim}
220    Here we define a joint inversion of gravity and magnetic data. \member{forward_models} is given as a list of
221    a tuple of a forward model and an index which referring to parameter in the \member{mappings} list to be used as an input.
222    The magnetic forward model \member{f_mag} is using the the second parameter (=\member{p_k}) in \member{mappings} list.
223    In this case the physical parameters are defined by a single-valued level set function. It is also possible
224    to link physical parameters to components of a level set function:
225    \begin{verbatim}
226    p_rho=DensityMapping(...)
227    p_k=SusceptibilityMapping(...)
228    f_mag=MagneticModel(...)
229    f_grav=GravityModel(...)
230    
231    J=InversionCostFunction(Regularization(numLevelSets=2,...), \
232                            mappings=[(p_rho,0), (p_k,1)], \
233                            forward_models=[[(f_mag, 1), (f_grav,0)])  
234    \end{verbatim}
235    The \member{mappings} argument is now a list of pairs where the first pair entry specifies the parameter mapping and
236    the second pair entry specifies the index of the component of the level set function to be used to evaluate the parameter.
237    In this case the level set function has two components, where the density mapping uses the first component of the level set function
238    while the susceptibility mapping uses the second component.
239    
240    The \class{InversionCostFunction} API is defined as follows:
241    
242    \begin{classdesc}{InversionCostFunction}{regularization, mappings, forward_models}
243    Constructor for the inversion cost function. \member{regularization} sets the regularization to be used, see Chapter~\ref{Chp:ref:regularization}.
244    \member{mappings} is a list of pairs where each pair comprises of a
245    physical parameter mapping (see Chapter~\ref{Chp:ref:mapping}) and an index which refers to the component of level set function
246    defined by the \member{regularization} to be used to calculate the corresponding physical parameter. If
247    the level set function has a single component the index can be omitted. If in addition there is a single physical parameter
248    the mapping can be given instead of a list. \member{forward_models} is a list of pairs where the
249    first pair component is a forward model ( see Chapter~\ref{Chp:ref:forward models}) and the second
250    pair component referring to the physical parameter in the \member{mappings} list providing the  physical parameter for the model.
251    If a single physical parameter is present the index can be omitted. If in addition a single  forward model is used this
252    forward model can be assigned to \member{forward_models} in replacement of a list.
253    \end{classdesc}
254    
255    
256  ===========  \begin{methoddesc}[InversionCostFunction]{getDomain}{}
257  \begin{classdesc}{SimpleCostFunction}{regularization, mapping, forwardmodel}          """
258      This is a simple cost function with a single continuous (mapped) variable.          returns the domain of the cost function
259      It is the sum of two weighted terms, a single forward model and a single          :rtype: 'Domain`
260      regularization term. This cost function is used by the provided gravity          """
261      and magnetic inversion implementations.          self.regularization.getDomain()
262  \end{classdesc}  \end{methoddesc}
263            
264    \begin{methoddesc}[InversionCostFunction]{getNumTradeOffFactors}{}
265            """
266            returns the number of trade-off factors being used including the
267            trade-off factors used in the regularization component.
268    
269            :rtype: ``int``
270            """
271            return self.__num_tradeoff_factors
272    \end{methoddesc}
273    
274     \begin{methoddesc}[InversionCostFunction]{getForwardModel}{idx=None}
275            """
276            returns the *idx*-th forward model.
277    
278            :param idx: model index. If cost function contains one model only `idx`
279                        can be omitted.
280            :type idx: ``int``
281            """
282            if idx==None: idx=0
283            f=self.forward_models[idx]
284            if isinstance(f, ForwardModel):
285                  F=f
286            else:
287                  F=f[0]
288            return F
289    \end{methoddesc}
290            
291    \begin{methoddesc}[InversionCostFunction]{getRegularization}{}
292            """
293            returns the regularization
294            """
295            return self.regularization
296    \end{methoddesc}
297    
298            
299    \begin{methoddesc}[InversionCostFunction]{setTradeOffFactorsModels}{mu=None}
300            """
301            sets the trade-off factors for the forward model components.
302            
303            :param mu: list of the trade-off factors. If not present ones are used.
304            :type mu: ``float`` in case of a single model or a ``list`` of ``float``
305                      with the length of the number of models.
306            """
307            if mu==None:
308                self.mu_model=np.ones((self.numModels, ))
309            else:
310                if self.numModels > 1:
311                   mu=np.asarray(mu)
312                   if min(mu) > 0:
313                      self.mu_model= mu
314                   else:
315                      raise ValueError("All value for trade-off factor mu must be positive.")
316                else:
317                  mu=float(mu)
318                  if mu > 0:
319                      self.mu_model= [mu, ]
320                  else:
321                      raise ValueError("Trade-off factor must be positive.")
322    \end{methoddesc}
323              
324    \begin{methoddesc}[InversionCostFunction]{setTradeOffFactorsRegularization}{mu=None, mu_c=None}
325            """
326            sets the trade of factors for the regularization component of the cost
327            function, see `Regularization` for details.
328            
329            :param mu: trade-off factors for the level-set variation part
330            :param mu_c: trade-off factors for the cross gradient variation part
331            """
332            self.regularization.setTradeOffFactorsForVariation(mu)
333            self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
334    \end{methoddesc}
335            
336    \begin{methoddesc}[InversionCostFunction]{setTradeOffFactors}{mu=None}
337            """
338            sets the trade-off factors for the forward model and regularization
339            terms.
340    
341            :param mu: list of trade-off factors.
342            :type mu: ``list`` of ``float``
343            """
344            if mu is None:
345                mu=np.ones((self.__num_tradeoff_factors,))
346            self.setTradeOffFactorsModels(mu[:self.numModels])
347            self.regularization.setTradeOffFactors(mu[self.numModels:])
348    \end{methoddesc}
349    
350    \begin{methoddesc}[InversionCostFunction]{createLevelSetFunction}{*props}
351            """
352            returns an instance of an object used to represent a level set function
353            initialized with zeros. Components can be overwritten by physical
354            properties 'props'. If present entries must correspond to the
355            `mappings` arguments in the constructor. Use `None` for properties for
356            which no value is given.
357            """
358            m=self.regularization.getPDE().createSolution()
359            if len(props) > 0:
360               for i in xrange(self.numMappings):
361                  if props[i]:
362                      mm=self.mappings[i]
363                      if isinstance(mm, Mapping):
364                          m=mm.getInverse(props[i])
365                      elif len(mm) == 1:
366                          m=mm[0].getInverse(props[i])
367                      else:
368                          m[mm[1]]=mm[0].getInverse(props[i])
369            return m
370    \end{methoddesc}
371        
372    \begin{methoddesc}[InversionCostFunction]{getProperties}{m, return_list=False}
373            """
374            returns a list of the physical properties from a given level set
375            function *m* using the mappings of the cost function.
376            
377            :param m: level set function
378            :type m: `Data`
379            :param return_list: if True a list is returned.
380                def _:type return_list: `bool`
381            :rtype: `list` of `Data`
382            """
383            props=[]
384            for i in xrange(self.numMappings):
385               mm=self.mappings[i]
386               if isinstance(mm, Mapping):
387                   p=mm.getValue(m)
388               elif len(mm) == 1:
389                   p=mm[0].getValue(m)
390               else:
391                   p=mm[0].getValue(m[mm[1]])
392               props.append(p)
393            if self.numMappings > 1 or return_list:
394               return props
395            else:
396               return props[0]
397    \end{methoddesc}
398              
399    \begin{methoddesc}[InversionCostFunction]{getDualProduct}{x, r}
400            """
401            Returns the dual product, see `Regularization.getDualProduct`
402    
403            :type x: `Data`
404            :type r: `ArithmeticTuple`            
405            :rtype: `float`
406            """
407            return self.regularization.getDualProduct(x, r)
408    \end{methoddesc}
409    
410    \begin{methoddesc}[InversionCostFunction]{getArguments}{m}
411            """
412            returns pre-computed values that are shared in the calculation of
413            *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
414            mapped value of ``m``, the arguments from the forward model and the
415            arguments from the regularization.
416            
417            :param m: current approximation of the level set function
418            :type m: `Data`
419            :return: tuple of of values of the parameters, pre-computed values for the forward model and
420                     pre-computed values for the regularization
421            :rtype: `tuple`
422            """
423            args_reg=self.regularization.getArguments(m)
424    
425            props=self.getProperties(m, return_list=True)
426            args_f=[]
427            for i in xrange(self.numModels):
428               f=self.forward_models[i]
429               if isinstance(f, ForwardModel):
430                  aa=f.getArguments(props[0])
431               elif len(f) == 1:
432                  aa=f[0].getArguments(props[0])
433               else:
434                  idx = f[1]
435                  f=f[0]
436                  if isinstance(idx, int):
437                     aa=f.getArguments(props[idx])
438                  else:
439                     pp=tuple( [ props[i] for i in idx] )
440                     aa=f.getArguments(*pp)
441               args_f.append(aa)
442              
443            return props, args_f, args_reg
444    \end{methoddesc}
445    
446    \begin{methoddesc}[InversionCostFunction]{getValue}{m, *args}
447            """
448            Returns the value *J(m)* of the cost function at *m*.
449            If the pre-computed values are not supplied `getArguments()` is called.
450    
451            :param m: current approximation of the level set function
452            :type m: `Data`
453            :param args: tuple of of values of the parameters, pre-computed values for the forward model and
454                     pre-computed values for the regularization
455            :rtype: `float`
456            """
457    
458            if len(args)==0:
459                args=self.getArguments(m)
460            
461            props=args[0]
462            args_f=args[1]
463            args_reg=args[2]
464            
465            J = self.regularization.getValue(m, *args_reg)
466            print "J_reg = %e"%J
467                    
468            for i in xrange(self.numModels):
469                    
470               f=self.forward_models[i]
471               if isinstance(f, ForwardModel):
472                  J_f = f.getValue(props[0],*args_f[i])
473               elif len(f) == 1:
474                  J_f=f[0].getValue(props[0],*args_f[i])
475               else:
476                  idx = f[1]
477                  f=f[0]
478                  if isinstance(idx, int):
479                     J_f = f.getValue(props[idx],*args_f[i])
480                  else:
481                     args=tuple( [ props[j] for j in idx] + args_f[i])
482                     J_f = f.getValue(*args)
483               print "J_f[%d] = %e"%(i, J_f)
484               print "mu_model[%d] = %e"%(i, self.mu_model[i])
485               J += self.mu_model[i] * J_f
486              
487            return   J
488    \end{methoddesc}
489    
490    \begin{methoddesc}[InversionCostFunction]{getGradient}{m, *args}
491            """
492            returns the gradient of the cost function  at *m*.
493            If the pre-computed values are not supplied `getArguments()` is called.
494    
495            :param m: current approximation of the level set function
496            :type m: `Data`
497            :param args: tuple of of values of the parameters, pre-computed values for the forward model and
498                     pre-computed values for the regularization
499                    
500            :rtype: `ArithmeticTuple`
501            """
502            if len(args)==0:
503                args = self.getArguments(m)
504            
505            props=args[0]
506            args_f=args[1]
507            args_reg=args[2]
508            
509            g_J = self.regularization.getGradient(m, *args_reg)
510            p_diffs=[]
511            for i in xrange(self.numMappings):
512               mm=self.mappings[i]
513               if isinstance(mm, Mapping):
514                   dpdm = mm.getDerivative(m)
515               elif len(mm) == 1:
516                   dpdm = mm[0].getDerivative(m)
517               else:
518                   dpdm = mm[0].getDerivative(m[mm[1]])
519               p_diffs.append(dpdm)
520              
521            Y=g_J[0]  
522            for i in xrange(self.numModels):
523               mu=self.mu_model[i]
524               f=self.forward_models[i]
525               if isinstance(f, ForwardModel):
526                  Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
527                  if self.numLevelSets == 1 :
528                     Y +=Ys
529                  else:
530                      Y[0] +=Ys
531               elif len(f) == 1:
532                  Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0]  * mu
533                  if self.numLevelSets == 1 :
534                     Y +=Ys
535                  else:
536                      Y[0] +=Ys
537               else:
538                  idx = f[1]
539                  f=f[0]
540                  if isinstance(idx, int):
541                     Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu
542                     if self.numLevelSets == 1 :
543                         if idx == 0:
544                             Y+=Ys
545                         else:
546                             raise IndexError("Illegal mapping index.")
547                     else:
548                         Y[idx] += Ys
549                  else:
550                     args=tuple( [ props[j] for j in idx] + args_f[i])
551                     Ys = f.getGradient(*args)
552                     for ii in xrange(len(idx)):
553                         Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]]  * mu
554    
555            return g_J
556    \end{methoddesc}
557    
558    
559    \begin{methoddesc}[InversionCostFunction]{getInverseHessianApproximation}{m, r, *args}
560            """
561            returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function
562            for a given gradient type *r* at a given location *m*: *H(m) p = r*
563    
564            :param m: level set approximation where to calculate Hessian inverse
565            :type m: `Data`
566            :param r: a given gradient
567            :type r: `ArithmeticTuple`
568            :param args: tuple of of values of the parameters, pre-computed values for the forward model and
569                     pre-computed values for the regularization
570            :rtype: `Data`
571            :note: in the current implementation only the regularization term is
572                   considered in the inverse Hessian approximation.
573    
574            """
575            m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
576            return m
577    
578    \end{methoddesc}
579            
580    \begin{methoddesc}[InversionCostFunction]{getNorm}{m}
581            """
582            returns the norm of ``m``
583    
584            :param m: level set function
585            :type m: `Data`
586            :rtype: ``float``
587            """
588    
589    \end{methoddesc}

Legend:
Removed from v.4132  
changed lines
  Added in v.4133

  ViewVC Help
Powered by ViewVC 1.1.26