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

  ViewVC Help
Powered by ViewVC 1.1.26