/[escript]/trunk/downunder/py_src/inversions.py
ViewVC logotype

Contents of /trunk/downunder/py_src/inversions.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4397 - (show annotations)
Thu May 9 06:17:15 2013 UTC (6 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 25636 byte(s)
Added some doco to make it a little easier to work out what is going on.
Cleaned up some reST ugliness

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 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 """Higher-level classes that allow running inversions with minimal set up"""
17
18 __copyright__="""Copyright (c) 2003-2013 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__ = ['InversionDriver', 'GravityInversion','MagneticInversion', 'JointGravityMagneticInversion']
26
27 import logging
28 import numpy as np
29
30 from esys.escript import *
31 from esys.escript import unitsSI as U
32 from esys.weipa import createDataset
33
34 from .inversioncostfunctions import InversionCostFunction
35 from .forwardmodels import GravityModel, MagneticModel
36 from .mappings import *
37 from .minimizers import *
38 from .regularizations import Regularization
39 from .datasources import DataSource
40 from .coordinates import makeTranformation
41
42 class InversionDriver(object):
43 """
44 Base class for running an inversion
45 """
46 def __init__(self, solverclass=None, debug=False):
47 """
48 creates an instance of an inversion driver.
49
50 :param solverclass: class of the solver to be used.
51 :type solverclass: 'AbstractMinimizer'.
52 """
53 # This configures basic logging to get some meaningful output
54 if debug:
55 level=logging.DEBUG
56 else:
57 level=logging.INFO
58 logging.getLogger('inv').setLevel(level)
59 self.logger=logging.getLogger('inv.%s'%self.__class__.__name__)
60 if solverclass is None: solverclass=MinimizerLBFGS
61 self.__solver=solverclass()
62 self.initial_value = None
63 self.m=None
64 self.fixGravityPotentialAtBottom()
65 self.fixMagneticPotentialAtBottom()
66
67 def fixGravityPotentialAtBottom(self, status=False):
68 """
69 indicates to fix the gravity potential at the bottom to zero (in addition to the top)
70
71 :param status: if True gravity potential at the bottom is set to zero
72 :type status: ``bool``
73 """
74 self._fixGravityPotentialAtBottom=status
75
76 def fixMagneticPotentialAtBottom(self, status=True):
77 """
78 indicates to fix the magnetic potential at the bottom to zero (in addition to the top)
79
80 :param status: if True magnetic potential at the bottom is set to zero
81 :type status: ``bool``
82 """
83 self._fixMagneticPotentialAtBottom=status
84
85 def setCostFunction(self, costfunction):
86 """
87 sets the cost function of the inversion. This function needs to be called
88 before the inversion iteration can be started.
89
90 :param costfunction: domain of the inversion
91 :type costfunction: 'InversionCostFunction'
92 """
93 self.__solver.setCostFunction(costfunction)
94
95 def getCostFunction(self):
96 """
97 returns the domain of the inversion
98
99 :rtype: 'InversionCostFunction'
100 """
101 if self.isSetUp():
102 return self.getSolver().getCostFunction()
103 else:
104 raise RuntimeError("Inversion is not set up.")
105
106 def getSolver(self):
107 """
108 The solver to be used in the inversion process. See the minimizers
109 module for available solvers. By default, the L-BFGS minimizer is used.
110
111 :rtype: 'AbstractMinimizer'.
112 """
113 return self.__solver
114
115 def isSetUp(self):
116 """
117 returns True if the inversion is set up and is ready to run.
118
119 :rtype: `bool`
120 """
121 if self.getSolver().getCostFunction():
122 return True
123 else:
124 return False
125
126 def getDomain(self):
127 """
128 returns the domain of the inversion
129
130 :rtype: `Domain`
131 """
132 self.getCostFunction().getDomain()
133
134 def setSolverCallback(self, callback=None):
135 """
136 Sets the callback function which is called after every solver
137 iteration.
138 """
139 self.getSolver().setCallback(callback)
140
141 def setSolverMaxIterations(self, maxiter=None):
142 """
143 Sets the maximum number of solver iterations to run.
144 If `maxiter` is reached the iteration is terminated and
145 `MinimizerMaxIterReached` is thrown.
146
147 :param maxiter: maximum number of iteration steps.
148 :type maxiter: positive ``int``
149 """
150 if maxiter == None: maxiter = 200
151 if maxiter>0:
152 self.getSolver().setMaxIterations(maxiter)
153 else:
154 raise ValueError("maxiter must be positive.")
155
156
157 def setSolverTolerance(self, m_tol=None, J_tol=None):
158 """
159 Sets the error tolerance for the solver. An acceptable solution is
160 considered to be found once the tolerance is reached.
161
162 :param m_tol: tolerance for changes to level set function. If `None` changes to the
163 level set function are not checked for convergence during iteration.
164 :param J_tol: tolerance for changes to cost function. If `None` changes to the
165 cost function are not checked for convergence during iteration.
166 :type m_tol: `float` or `None`
167 :type J_tol: `float` or `None`
168 :note: if both arguments are equal to `None` the default setting m_tol=1e-4, J_tol=None is used.
169
170 """
171 if m_tol == None and J_tol==None:
172 m_tol=1e-4
173 self.getSolver().setTolerance(m_tol=m_tol, J_tol=J_tol)
174
175 def setInitialGuess(self, *props):
176 """
177 Sets the initial guess for the inversion iteration.
178 By default zero is used.
179 """
180 self.initial_value=self.getCostFunction().createLevelSetFunction(*props)
181
182 def run(self):
183 """
184 This function runs the inversion.
185
186 :return: physical parameters as result of the inversion
187 :rtype: ``list`` of physical parameters or a physical parameter
188 """
189 if not self.isSetUp():
190 raise RuntimeError("Inversion is not setup.")
191
192 if self.initial_value == None: self.setInitialGuess()
193 self.logger.info("Starting solver, please stand by...")
194 try:
195 self.getSolver().run(self.initial_value)
196 self.m=self.getSolver().getResult()
197 self.p=self.getCostFunction().getProperties(self.m)
198 except MinimizerException as e:
199 self.m=self.getSolver().getResult()
200 self.p=self.getCostFunction().getProperties(self.m)
201 self.logger.info("Solver iteration failed!")
202 raise e
203 self.logger.info("Solver completed successfully!")
204 self.getSolver().logSummary()
205 return self.p
206
207 def getLevelSetFunction(self):
208 """
209 returns the level set function as solution of the optimization problem
210
211 :rtype: `Data`
212 """
213 return self.m
214
215 def setup(self, *args, **k_args):
216 """
217 returns True if the inversion is set up and ready to run.
218 """
219 pass
220
221 class GravityInversion(InversionDriver):
222 """
223 Driver class to perform an inversion of Gravity (Bouguer) anomaly data.
224 The class uses the standard `Regularization` class for a single level set
225 function, `DensityMapping` mapping, and the gravity forward model
226 `GravityModel`.
227 """
228 def setup(self, domainbuilder,
229 rho0=None, drho=None, z0=None, beta=None,
230 w0=None, w1=None, rho_at_depth=None):
231 """
232 Sets up the inversion parameters from a `DomainBuilder`.
233
234 :param domainbuilder: Domain builder object with gravity source(s)
235 :type domainbuilder: `DomainBuilder`
236 :param rho0: reference density, see `DensityMapping`. If not specified, zero is used.
237 :type rho0: ``float`` or ``Scalar``
238 :param drho: density scale, see `DensityMapping`. If not specified, 2750kg/m^3 is used.
239 :type drho: ``float`` or ``Scalar``
240 :param z0: reference depth for depth weighting, see `DensityMapping`. If not specified, zero is used.
241 :type z0: ``float`` or ``Scalar``
242 :param beta: exponent for depth weighting, see `DensityMapping`. If not specified, zero is used.
243 :type beta: ``float`` or ``Scalar``
244 :param w0: weighting factor for level set term regularization. If not set zero is assumed.
245 :type w0: ``Scalar`` or ``float``
246 :param w1: weighting factor for the gradient term in the regularization. If not set zero is assumed
247 :type w1: ``Vector`` or list of ``float``
248 :param rho_at_depth: value for density at depth, see `DomainBuilder`.
249 :type rho_at_depth: ``float`` or ``None``
250 """
251 self.logger.info('Retrieving domain...')
252 dom=domainbuilder.getDomain()
253 trafo=makeTranformation(dom, domainbuilder.getReferenceSystem())
254 DIM=dom.getDim()
255 rho_mask = domainbuilder.getSetDensityMask()
256 #========================
257 self.logger.info('Creating mapping...')
258 if rho_at_depth:
259 rho2 = rho_mask * rho_at_depth + (1-rho_mask) * rho0
260 elif rho0:
261 rho2 = (1-rho_mask) * rho0
262 else:
263 rho2 = 0.
264
265 rho_mapping=DensityMapping(dom, rho0=rho2, drho=drho, z0=z0, beta=beta)
266 scale_mapping=rho_mapping.getTypicalDerivative()
267 #========================
268 self.logger.info("Setting up regularization...")
269 if w1 is None:
270 w1=[1.]*DIM
271
272 regularization=Regularization(dom, numLevelSets=1,\
273 w0=w0, w1=w1, location_of_set_m=rho_mask,
274 coordinates=trafo)
275 #====================================================================
276 self.logger.info("Retrieving gravity surveys...")
277 surveys=domainbuilder.getGravitySurveys()
278 g=[]
279 w=[]
280 for g_i,sigma_i in surveys:
281 w_i=safeDiv(1., sigma_i)
282 if g_i.getRank()==0:
283 g_i=g_i*kronecker(DIM)[DIM-1]
284 if w_i.getRank()==0:
285 w_i=w_i*kronecker(DIM)[DIM-1]
286 g.append(g_i)
287 w.append(w_i)
288 self.logger.debug("Added gravity survey:")
289 self.logger.debug("g = %s"%g_i)
290 self.logger.debug("sigma = %s"%sigma_i)
291 self.logger.debug("w = %s"%w_i)
292 #====================================================================
293
294 self.logger.info("Setting up model...")
295 forward_model=GravityModel(dom, w, g, fixPotentialAtBottom=self._fixGravityPotentialAtBottom, coordinates=trafo)
296 forward_model.rescaleWeights(rho_scale=scale_mapping)
297
298 #====================================================================
299 self.logger.info("Setting cost function...")
300 self.setCostFunction(InversionCostFunction(regularization, rho_mapping, forward_model))
301
302 def setInitialGuess(self, rho=None):
303 """
304 set the initial guess *rho* for density the inversion iteration. If no *rho* present
305 then an appropriate initial guess is chosen.
306
307 :param rho: initial value for the density anomaly.
308 :type rho: `Scalar`
309 """
310 if rho:
311 super(GravityInversion,self).setInitialGuess(rho)
312 else:
313 super(GravityInversion,self).setInitialGuess()
314
315 def siloWriterCallback(self, k, m, Jm, g_Jm):
316 """
317 callback function that can be used to track the solution
318
319 :param k: iteration count
320 :param m: current m approximation
321 :param Jm: value of cost function
322 :param g_Jm: gradient of f at x
323 """
324 fn='inv.%d'%k
325 ds=createDataset(density=self.getCostFunction().getProperties(m))
326 ds.setCycleAndTime(k,k)
327 ds.saveSilo(fn)
328 self.logger.debug("J(m) = %e"%Jm)
329
330 class MagneticInversion(InversionDriver):
331 """
332 Driver class to perform an inversion of magnetic anomaly data. The class
333 uses the standard `Regularization` class for a single level set function,
334 `SusceptibilityMapping` mapping and the linear magnetic forward model
335 `MagneticModel`.
336 """
337 def setup(self, domainbuilder,
338 k0=None, dk=None, z0=None, beta=None,
339 w0=None, w1=None,k_at_depth=None ):
340 """
341 Sets up the inversion from a `DomainBuilder`.
342 If magnetic data are given as scalar it is assumed that values are
343 collected in direction of the background magnetic field.
344
345 :param domainbuilder: Domain builder object with gravity source(s)
346 :type domainbuilder: `DomainBuilder`
347 :param k0: reference susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
348 :type k0: ``float`` or ``Scalar``
349 :param dk: susceptibility scale, see `SusceptibilityMapping`. If not specified, 2750kg/m^3 is used.
350 :type dk: ``float`` or ``Scalar``
351 :param z0: reference depth for depth weighting, see `SusceptibilityMapping`. If not specified, zero is used.
352 :type z0: ``float`` or ``Scalar``
353 :param beta: exponent for depth weighting, see `SusceptibilityMapping`. If not specified, zero is used.
354 :type beta: ``float`` or ``Scalar``
355 :param w0: weighting factor for level set term regularization. If not set zero is assumed.
356 :type w0: ``Scalar`` or ``float``
357 :param w1: weighting factor for the gradient term in the regularization. If not set zero is assumed
358 :type w1: ``Vector`` or list of ``float``
359 :param k_at_depth: value for susceptibility at depth, see `DomainBuilder`.
360 :type k_at_depth: ``float`` or ``None``
361 """
362 self.logger.info('Retrieving domain...')
363 dom=domainbuilder.getDomain()
364 DIM=dom.getDim()
365 trafo=makeTranformation(dom, domainbuilder.getReferenceSystem())
366 #========================
367 self.logger.info('Creating mapping ...')
368 k_mask = domainbuilder.getSetSusceptibilityMask()
369 if k_at_depth:
370 k2= k_mask * k_at_depth + (1-k_mask) * k0
371 elif k0:
372 k2= (1-k_mask) * k0
373 else:
374 k2=0
375
376 k_mapping=SusceptibilityMapping(dom, k0=k2, dk=dk, z0=z0, beta=beta)
377 scale_mapping=k_mapping.getTypicalDerivative()
378 #========================
379 self.logger.info("Setting up regularization...")
380 if w1 is None:
381 w1=[1.]*DIM
382
383 regularization=Regularization(dom, numLevelSets=1,w0=w0, w1=w1, location_of_set_m=k_mask, coordinates=trafo)
384
385 #====================================================================
386 self.logger.info("Retrieving magnetic field surveys...")
387 d_b=normalize(domainbuilder.getBackgroundMagneticFluxDensity())
388 surveys=domainbuilder.getMagneticSurveys()
389 B=[]
390 w=[]
391 for B_i,sigma_i in surveys:
392 w_i=safeDiv(1., sigma_i)
393 if B_i.getRank()==0:
394 B_i=B_i*d_b
395 if w_i.getRank()==0:
396 w_i=w_i*d_b
397 B.append(B_i)
398 w.append(w_i)
399 self.logger.debug("Added magnetic survey:")
400 self.logger.debug("B = %s"%B_i)
401 self.logger.debug("sigma = %s"%sigma_i)
402 self.logger.debug("w = %s"%w_i)
403 #====================================================================
404 self.logger.info("Setting up model...")
405 forward_model=MagneticModel(dom, w, B, domainbuilder.getBackgroundMagneticFluxDensity(), fixPotentialAtBottom=self._fixMagneticPotentialAtBottom, coordinates=trafo)
406 forward_model.rescaleWeights(k_scale=scale_mapping)
407
408 #====================================================================
409 self.logger.info("Setting cost function...")
410 self.setCostFunction(InversionCostFunction(regularization, k_mapping, forward_model))
411
412 def setInitialGuess(self, k=None):
413 """
414 set the initial guess *k* for susceptibility for the inversion iteration. If no *k* present
415 then an appropriate initial guess is chosen.
416
417 :param k: initial value for the susceptibility anomaly.
418 :type k: `Scalar`
419 """
420 if k:
421 super(MagneticInversion,self).setInitialGuess(k)
422 else:
423 super(MagneticInversion,self).setInitialGuess()
424
425 def siloWriterCallback(self, k, m, Jm, g_Jm):
426 """
427 callback function that can be used to track the solution
428
429 :param k: iteration count
430 :param m: current m approximation
431 :param Jm: value of cost function
432 :param g_Jm: gradient of f at x
433 """
434 fn='inv.%d'%k
435 ds=createDataset(susceptibility=self.getCostFunction().getProperties(m))
436 ds.setCycleAndTime(k,k)
437 ds.saveSilo(fn)
438 self.logger.debug("J(m) = %e"%Jm)
439
440 class JointGravityMagneticInversion(InversionDriver):
441 """
442 Driver class to perform a joint inversion of Gravity (Bouguer) and
443 magnetic anomaly data. The class uses the standard `Regularization` class
444 for two level set functions with cross-gradient correlation,
445 `DensityMapping` and `SusceptibilityMapping` mappings, the gravity forward
446 model `GravityModel` and the linear magnetic forward model `MagneticModel`.
447 """
448 DENSITY=0
449 SUSCEPTIBILITY=1
450
451 def setup(self, domainbuilder,
452 rho0=None, drho=None, rho_z0=None, rho_beta=None,
453 k0=None, dk=None, k_z0=None, k_beta=None, w0=None, w1=None,
454 w_gc=None,rho_at_depth=None, k_at_depth=None):
455 """
456 Sets up the inversion from an instance ``domainbuilder`` of a
457 `DomainBuilder`. Gravity and magnetic data attached to the
458 ``domainbuilder`` are considered in the inversion.
459 If magnetic data are given as scalar it is assumed that values are
460 collected in direction of the background magnetic field.
461
462 :param domainbuilder: Domain builder object with gravity source(s)
463 :type domainbuilder: `DomainBuilder`
464 :param rho0: reference density, see `DensityMapping`. If not specified, zero is used.
465 :type rho0: ``float`` or `Scalar`
466 :param drho: density scale, see `DensityMapping`. If not specified, 2750kg/m^3 is used.
467 :type drho: ``float`` or `Scalar`
468 :param rho_z0: reference depth for depth weighting for density, see `DensityMapping`. If not specified, zero is used.
469 :type rho_z0: ``float`` or `Scalar`
470 :param rho_beta: exponent for depth weighting for density, see `DensityMapping`. If not specified, zero is used.
471 :type rho_beta: ``float`` or `Scalar`
472 :param k0: reference susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
473 :type k0: ``float`` or `Scalar`
474 :param dk: susceptibility scale, see `SusceptibilityMapping`. If not specified, 2750kg/m^3 is used.
475 :type dk: ``float`` or `Scalar`
476 :param k_z0: reference depth for depth weighting for susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
477 :type k_z0: ``float`` or `Scalar`
478 :param k_beta: exponent for depth weighting for susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
479 :type k_beta: ``float`` or `Scalar`
480 :param w0: weighting factors for level set term regularization, see `Regularization`. If not set zero is assumed.
481 :type w0: `Data` or ``ndarray`` of shape (2,)
482 :param w1: weighting factor for the gradient term in the regularization see `Regularization`. If not set zero is assumed
483 :type w1: `Data` or ``ndarray`` of shape (2,DIM)
484 :param w_gc: weighting factor for the cross gradient term in the regularization, see `Regularization`. If not set one is assumed
485 :type w_gc: `Scalar` or `float`
486 :param k_at_depth: value for susceptibility at depth, see `DomainBuilder`.
487 :type k_at_depth: ``float`` or ``None``
488 :param rho_at_depth: value for density at depth, see `DomainBuilder`.
489 :type rho_at_depth: ``float`` or ``None``
490 """
491 self.logger.info('Retrieving domain...')
492 dom=domainbuilder.getDomain()
493 DIM=dom.getDim()
494 trafo=makeTranformation(dom, domainbuilder.getReferenceSystem())
495 #========================
496 self.logger.info('Creating mappings ...')
497 rho_mask=domainbuilder.getSetDensityMask()
498 if rho_at_depth:
499 rho2= rho_mask * rho_at_depth + (1-rho_mask) * rho0
500 elif rho0:
501 rho2= (1-rho_mask) * rho0
502 else:
503 rho2=0
504
505 k_mask = domainbuilder.getSetSusceptibilityMask()
506 if k_at_depth:
507 k2= k_mask * k_at_depth + (1-k_mask) * k0
508 elif k0:
509 k2= (1-k_mask) * k0
510 else:
511 k2=0
512
513 rho_mapping=DensityMapping(dom, rho0=rho2, drho=drho, z0=rho_z0, beta=rho_beta)
514 rho_scale_mapping=rho_mapping.getTypicalDerivative()
515 self.logger.debug("rho_scale_mapping = %s"%rho_scale_mapping)
516 k_mapping=SusceptibilityMapping(dom, k0=k2, dk=dk, z0=k_z0, beta=k_beta)
517 k_scale_mapping=k_mapping.getTypicalDerivative()
518 self.logger.debug("k_scale_mapping = %s"%k_scale_mapping)
519 #========================
520 self.logger.info("Setting up regularization...")
521
522 if w1 is None:
523 w1=np.ones((2,DIM))
524
525 wc=Data(0.,(2,2), Function(dom))
526 if w_gc is None:
527 wc[0,1]=1
528 else:
529 wc[0,1]=w_gc
530
531 reg_mask=Data(0.,(2,), Solution(dom))
532 reg_mask[self.DENSITY] = rho_mask
533 reg_mask[self.SUSCEPTIBILITY] = k_mask
534 regularization=Regularization(dom, numLevelSets=2,\
535 w0=w0, w1=w1,wc=wc, location_of_set_m=reg_mask, coordinates=trafo)
536 #====================================================================
537 self.logger.info("Retrieving gravity surveys...")
538 surveys=domainbuilder.getGravitySurveys()
539 g=[]
540 w=[]
541 for g_i,sigma_i in surveys:
542 w_i=safeDiv(1., sigma_i)
543 if g_i.getRank()==0:
544 g_i=g_i*kronecker(DIM)[DIM-1]
545 if w_i.getRank()==0:
546 w_i=w_i*kronecker(DIM)[DIM-1]
547 g.append(g_i)
548 w.append(w_i)
549 self.logger.debug("Added gravity survey:")
550 self.logger.debug("g = %s"%g_i)
551 self.logger.debug("sigma = %s"%sigma_i)
552 self.logger.debug("w = %s"%w_i)
553
554 self.logger.info("Setting up gravity model...")
555 gravity_model=GravityModel(dom, w, g, fixPotentialAtBottom=self._fixGravityPotentialAtBottom, coordinates=trafo)
556 gravity_model.rescaleWeights(rho_scale=rho_scale_mapping)
557 #====================================================================
558 self.logger.info("Retrieving magnetic field surveys...")
559 d_b=normalize(domainbuilder.getBackgroundMagneticFluxDensity())
560 surveys=domainbuilder.getMagneticSurveys()
561 B=[]
562 w=[]
563 for B_i,sigma_i in surveys:
564 w_i=safeDiv(1., sigma_i)
565 if B_i.getRank()==0:
566 B_i=B_i*d_b
567 if w_i.getRank()==0:
568 w_i=w_i*d_b
569 B.append(B_i)
570 w.append(w_i)
571 self.logger.debug("Added magnetic survey:")
572 self.logger.debug("B = %s"%B_i)
573 self.logger.debug("sigma = %s"%sigma_i)
574 self.logger.debug("w = %s"%w_i)
575
576 self.logger.info("Setting up magnetic model...")
577 magnetic_model=MagneticModel(dom, w, B, domainbuilder.getBackgroundMagneticFluxDensity(), fixPotentialAtBottom=self._fixMagneticPotentialAtBottom, coordinates=trafo)
578 magnetic_model.rescaleWeights(k_scale=k_scale_mapping)
579 #====================================================================
580 self.logger.info("Setting cost function...")
581 self.setCostFunction(InversionCostFunction(regularization,
582 ((rho_mapping,self.DENSITY), (k_mapping, self.SUSCEPTIBILITY)),
583 ((gravity_model,0), (magnetic_model,1)) ))
584
585 def setInitialGuess(self, rho=None, k=None):
586 """
587 set the initial guess *rho* for density and *k* for susceptibility for the inversion iteration.
588
589 :param rho: initial value for the density anomaly.
590 :type rho: `Scalar`
591 :param k: initial value for the susceptibility anomaly.
592 :type k: `Scalar`
593 """
594 super(JointGravityMagneticInversion,self).setInitialGuess(rho, k)
595
596 def siloWriterCallback(self, k, m, Jm, g_Jm):
597 """
598 callback function that can be used to track the solution
599
600 :param k: iteration count
601 :param m: current m approximation
602 :param Jm: value of cost function
603 :param g_Jm: gradient of f at x
604 """
605 fn='inv.%d'%k
606 p=self.getCostFunction().getProperties(m)
607 ds=createDataset(density=p[self.DENSITY], susceptibility=p[self.SUSCEPTIBILITY])
608 ds.setCycleAndTime(k,k)
609 ds.saveSilo(fn)
610 self.logger.debug("J(m) = %e"%Jm)
611

  ViewVC Help
Powered by ViewVC 1.1.26