/[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 4074 - (show annotations)
Thu Nov 15 03:30:59 2012 UTC (8 years ago) by gross
File MIME type: text/x-python
File size: 12998 byte(s)
this will break the downundermodule: major revision of the classes.
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 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-2012 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__ = ['InversionBase','SingleParameterInversionBase','GravityInversion','MagneticInversion']
26
27 import logging
28 from esys.escript import *
29 import esys.escript.unitsSI as U
30 from esys.weipa import createDataset
31
32 from .inversioncostfunctions import SimpleInversionCostFunction
33 from .forwardmodels import GravityModel, MagneticModel
34 from .mappings import *
35 from .minimizers import *
36 from .regularizations import Regularization
37 from .datasources import DataSource
38
39 class InversionBase(object):
40 """
41 Base class for inversions.
42 """
43 def __init__(self):
44 self.logger=logging.getLogger('inv.%s'%self.__class__.__name__)
45 self._solver_callback = None
46 self._solver_opts = {}
47 self._solver_tol = 1e-9
48 self._solver_maxiter = 200
49 # set default mapping and solver
50 self.setSolverClass()
51 self.__domain=None
52 self.__regularization=None
53 self.__forwardmodel=None
54 self.__mapping=None
55
56 def setDomain(self, domain):
57 """
58 sets the domain of the inversion
59
60 :param domain: domain of the inversion
61 :type domain: `Domain`
62 """
63 self.__domain=domain
64
65 def getDomain(self):
66 """
67 returns the domain of the inversion
68
69 :rtype: `Domain`
70 """
71 return self.__domain
72
73 def setMapping(self, mapping):
74 """
75 Sets the mapping object to map between model parameters and the data.
76 If no mapping is provided, an identity mapping is used
77 (`ScalingMapping` with constant 1).
78
79 :param mapping: parameter mapping
80 :type mapping: `Mapping`
81 """
82 self.__mapping=mapping
83
84 def getMapping(self):
85 """
86 return the mapping(s) used in the inversion
87
88 :rtype: `Mapping`
89 """
90 return self.__mapping
91
92
93 def setRegularization(self, regularization):
94 """
95 Sets the regularization for the inversion.
96
97 :param regularization: regularization
98 :type regularization: `Regularization`
99 """
100 self.__regularization=regularization
101
102 def getRegularization(self):
103 """
104 returns the regularization method(s)
105
106 :rtype: `Regularization`
107 """
108 return self.__regularization
109
110 def setForwardModel(self, forwardmodel):
111 """
112 Sets the forward model(s) for the inversion.
113
114 :param forwardmodel: forward model
115 :type forwardmodel: `ForwardModel`
116 """
117 self.__forwardmodel=forwardmodel
118
119 def getForwardModel(self):
120 """
121 returns the forward model
122
123 :rtype: `ForwardModel`
124 """
125 return self.__forwardmodel
126
127 def setSolverClass(self, solverclass=None):
128 """
129 The solver to be used in the inversion process. See the minimizers
130 module for available solvers.
131 By default, the L-BFGS minimizer is used.
132 """
133 if solverclass == None:
134 self.solverclass=MinimizerLBFGS
135 else:
136 self.solverclass=solverclass
137
138 def setSolverCallback(self, callback):
139 """
140 Sets the callback function which is called after every solver
141 iteration.
142 """
143 self._solver_callback=callback
144
145 def setSolverMaxIterations(self, maxiter):
146 """
147 Sets the maximum number of solver iterations to run.
148 """
149 if maxiter>0:
150 self._solver_maxiter=maxiter
151 else:
152 raise ValueError("maxiter must be positive.")
153
154 def setSolverOptions(self, **opts):
155 """
156 Sets additional solver options. The valid options depend on the
157 solver being used.
158 """
159 self._solver_opts.update(**opts)
160
161 def setSolverTolerance(self, tol):
162 """
163 Sets the error tolerance for the solver. An acceptable solution is
164 considered to be found once the tolerance is reached.
165 """
166 if tol>0:
167 self._solver_tol=tol
168 else:
169 raise ValueError("tolerance must be positive.")
170
171 def isSetup(self):
172 """
173 returns True if the inversion is set up and ready to run.
174 """
175 raise NotImplementedError
176
177 def run(self, *args):
178 """
179 This method starts the inversion process and must be overwritten.
180 Preferably, users should be able to set a starting point so
181 inversions can be continued after stopping.
182 """
183 raise NotImplementedError
184
185
186 class SingleParameterInversionBase(InversionBase):
187 """
188 Base class for inversions with a single parameter to be found.
189 """
190 def __init__(self):
191 super(SingleParameterInversionBase,self).__init__()
192 self.setWeights()
193
194 def setWeights(self, mu_reg=None, mu_model=1.):
195 """
196 Sets the weighting factors for the cost function.
197 """
198 self.logger.debug("Setting weighting factors...")
199 self.logger.debug("mu_reg = %s"%mu_reg)
200 self.logger.debug("mu_model = %s"%mu_model)
201 self._mu_reg=mu_reg
202 self._mu_model=mu_model
203
204 def siloWriterCallback(self, k, x, fx, gfx):
205 """
206 callback function that can be used to track the solution
207
208 :param k: iteration count
209 :param x: current m approximation
210 :param fx: value of cost function
211 :param gfx: gradient of f at x
212 """
213 fn='inv.%d'%k
214 ds=createDataset(rho=self.getMapping().getValue(x))
215 ds.setCycleAndTime(k,k)
216 ds.saveSilo(fn)
217 self.logger.debug("Jreg(m) = %e"%self.regularization.getValue(x))
218 self.logger.debug("f(m) = %e"%fx)
219
220
221 def isSetup(self):
222 if self.getRegularization() and self.getMapping() \
223 and self.getForwardModel() and self.getDomain():
224 return True
225 else:
226 return False
227
228 def run(self, initial_value=0.):
229 if not self.isSetup():
230 raise RuntimeError("Inversion is not setup properly.")
231
232 f=SimpleInversionCostFunction(self.getRegularization(), self.getMapping(), self.getForwardModel())
233 f.setWeights(mu_reg_1=self._mu_reg, mu_model=self._mu_model)
234
235 solver=self.solverclass(f)
236 solver.setCallback(self._solver_callback)
237 solver.setMaxIterations(self._solver_maxiter)
238 solver.setOptions(**self._solver_opts)
239 solver.setTolerance(self._solver_tol)
240 if not isinstance(initial_value, Data):
241 initial_value=Scalar(initial_value, ContinuousFunction(self.getDomain()))
242 m_init=self.getMapping().getInverse(initial_value)
243
244 self.logger.info("Starting solver...")
245 solver.run(m_init)
246 m_star=solver.getResult()
247 self.logger.info("m* = %s"%m_star)
248 value_star=self.getMapping().getValue(m_star)
249 self.logger.info("result* = %s"%value_star)
250 solver.logSummary()
251 return value_star
252
253
254 class GravityInversion(SingleParameterInversionBase):
255 """
256 Inversion of Gravity (Bouguer) anomaly data.
257 """
258 def setup(self, domainbuilder,
259 rho0=None, drho=None, z0=None, beta=None,
260 s0=None, s1=None):
261 """
262 Sets up the inversion parameters from a downunder.`DomainBuilder`.
263
264 :param domainbuilder: Domain builder object with gravity source(s)
265 :type domainbuilder: `DomainBuilder`
266
267
268 :param rho0: reference density. If not specified, zero is used.
269 :type rho0: ``float`` or `Scalar`
270 :param w0: weighting factor for L2 term in the regularization
271 :type w0: ``float`` or `Scalar`
272 :param w: weighting factor for H1 term in the regularization
273 :type w: ``list`` of float or `Vector`
274 """
275 self.logger.info('Retrieving domain...')
276 self.setDomain(domainbuilder.getDomain())
277 DIM=self.getDomain().getDim()
278 #========================
279 self.logger.info('Set mapping ...')
280 self.setMapping(DensityMapping(self.getDomain(), rho0=rho0, drho=drho, z0=z0, beta=beta))
281 #========================
282 self.logger.info("Setting up regularization...")
283 if s1 is None:
284 s1=[1.]*DIM
285 rho_mask = domainbuilder.getSetDensityMask()
286 self.setRegularization(Regularization(self.getDomain(), numLevelSets=1, \
287 s0=s0, s1=s1, location_of_set_m=rho_mask))
288 #====================================================================
289 self.logger.info("Retrieving gravity surveys...")
290 surveys=domainbuilder.getGravitySurveys()
291 g=[]
292 chi=[]
293 for g_i,sigma_i in surveys:
294 chi_i=safeDiv(1., sigma_i*sigma_i)
295 if g_i.getRank()==0:
296 g_i=g_i*kronecker(DIM)[DIM-1]
297 if chi_i.getRank()==0:
298 chi_i=chi_i*kronecker(DIM)[DIM-1]
299 g.append(g_i)
300 chi.append(chi_i)
301 self.logger.debug("Added gravity survey:")
302 self.logger.debug("g = %s"%g_i)
303 self.logger.debug("sigma = %s"%sigma_i)
304 self.logger.debug("chi = %s"%chi_i)
305 #====================================================================
306
307 self.logger.info("Setting up model...")
308 self.setForwardModel(GravityModel(self.getDomain(), chi, g))
309
310 # this is switched off for now:
311 if self._mu_reg is None and False:
312 print "SDSAD"
313 x=self.getDomain().getX()
314 l=0.
315 for i in range(DIM-1):
316 l=max(l, sup(x[i])-inf(x[i]))
317 G=U.Gravitational_Constant
318 mu_reg=0.5*(l*l*G)**2
319 self.setWeights(mu_reg=mu_reg)
320
321
322 class MagneticInversion(SingleParameterInversionBase):
323 """
324 Inversion of magnetic data.
325 """
326 def setup(self, domainbuilder,
327 k0=None, dk=None, z0=None, beta=None,
328 s0=None, s1=None):
329 """
330 Sets up the inversion parameters from a downunder.`DomainBuilder`.
331
332 :param domainbuilder: Domain builder object with magnetic source(s)
333 :type domainbuilder: `DomainBuilder`
334
335 :param k_ref: reference susceptibility. If not specified, zero is used.
336 :type k_ref: ``float`` or `Scalar`
337 :param w0: weighting factor for the L2 term in the regularization
338 :type w0: ``float`` or `Scalar`
339 :param w: weighting factor for the H1 term in the regularization
340 :type w: ``list`` of float or `Vector`
341 """
342 self.logger.info('Retrieving domain...')
343 self.setDomain(domainbuilder.getDomain())
344 DIM=self.getDomain().getDim()
345
346 #========================
347 self.logger.info('Set mapping ...')
348 self.setMapping(SusceptibilityMapping(self.getDomain(), k0=k0, dk=dk, z0=z0, beta=beta))
349
350 #========================
351 self.logger.info("Setting up regularization...")
352 if s1 is None:
353 s1=[1.]*DIM
354 k_mask = domainbuilder.getSetSusceptibilityMask()
355 self.setRegularization(Regularization(self.getDomain(), numLevelSets=1, \
356 s0=s0, s1=s1, location_of_set_m=k_mask))
357
358 #====================================================================
359 self.logger.info("Retrieving magnetic field surveys...")
360 surveys=domainbuilder.getMagneticSurveys()
361 B=[]
362 chi=[]
363 for B_i,sigma_i in surveys:
364 chi_i=safeDiv(1., sigma_i*sigma_i)
365 if B_i.getRank()==0:
366 B_i=B_i*kronecker(DIM)[DIM-1]
367 if chi_i.getRank()==0:
368 chi_i=chi_i*kronecker(DIM)[DIM-1]
369 B.append(B_i)
370 chi.append(chi_i)
371 self.logger.debug("Added magnetic survey:")
372 self.logger.debug("B = %s"%B_i)
373 self.logger.debug("sigma = %s"%sigma_i)
374 self.logger.debug("chi = %s"%chi_i)
375 #====================================================================
376 self.logger.info("Setting up model...")
377 self.setForwardModel(MagneticModel(self.getDomain(), chi, B, domainbuilder.getBackgroundMagneticField()))
378
379 # this is switched off for now:
380 if self._mu_reg is None and False:
381 x=self.getDomain().getX()
382 l=0.
383 for i in range(DIM-1):
384 l=max(l, sup(x[i])-inf(x[i]))
385 mu_reg=0.5*l**2
386 self.setWeights(mu_reg=mu_reg)
387

  ViewVC Help
Powered by ViewVC 1.1.26