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

Revision 4138 - (show 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 \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{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 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}