/[escript]/trunk/downunder/py_src/forwardmodels/base.py
ViewVC logotype

Contents of /trunk/downunder/py_src/forwardmodels/base.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6496 - (show annotations)
Tue Feb 14 00:38:08 2017 UTC (2 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 7939 byte(s)
bug to obtain function space fixed.
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2016 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Apache License, version 2.0
8 # http://www.apache.org/licenses/LICENSE-2.0
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15
16 """Base classes for forward models"""
17
18 from __future__ import division, print_function
19
20 __copyright__="""Copyright (c) 2003-2016 by The University of Queensland
21 http://www.uq.edu.au
22 Primary Business: Queensland, Australia"""
23 __license__="""Licensed under the Apache License, version 2.0
24 http://www.apache.org/licenses/LICENSE-2.0"""
25 __url__="https://launchpad.net/escript-finley"
26
27 __all__ = ['ForwardModel','ForwardModelWithPotential']
28
29 from esys.downunder.coordinates import makeTransformation
30 from esys.escript.linearPDEs import LinearSinglePDE
31 from esys.escript.util import *
32 import numpy as np
33
34 class ForwardModel(object):
35 """
36 An abstract forward model that can be plugged into a cost function.
37 Subclasses need to implement `getDefect()`, `getGradient()`, and possibly
38 `getArguments()` and 'getCoordinateTransformation'.
39 """
40 def __init__(self):
41 pass
42
43 def getArguments(self, x):
44 return ()
45
46 def getCoordinateTransformation(self):
47 return None
48
49 def getDefect(self, x, *args):
50 raise NotImplementedError
51
52 def getGradient(self, x, *args):
53 raise NotImplementedError
54
55
56 class ForwardModelWithPotential(ForwardModel):
57 """
58 Base class for a forward model using a potential such as magnetic or
59 gravity. It defines a cost function:
60
61 defect = 1/2 sum_s integrate( ( weight_i[s] * ( r_i - data_i[s] ) )**2 )
62
63 where s runs over the survey, weight_i are weighting factors, data_i are
64 the data, and r_i are the results produced by the forward model.
65 It is assumed that the forward model is produced through postprocessing
66 of the solution of a potential PDE.
67 """
68 def __init__(self, domain, w, data, coordinates=None,
69 fixPotentialAtBottom=False,
70 tol=1e-8):
71 """
72 initializes a new forward model with potential.
73
74 :param domain: domain of the model
75 :type domain: `Domain`
76 :param w: data weighting factors
77 :type w: ``Vector`` or list of ``Vector``
78 :param data: data
79 :type data: ``Vector`` or list of ``Vector``
80 :param coordinates: defines coordinate system to be used
81 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
82 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
83 in addition to the top.
84 :type fixPotentialAtBottom: ``bool``
85 :param tol: tolerance of underlying PDE
86 :type tol: positive ``float``
87 """
88 super(ForwardModelWithPotential, self).__init__()
89 self.__domain = domain
90 self.__trafo = makeTransformation(domain, coordinates)
91
92 try:
93 n=len(w)
94 m=len(data)
95 if not m == n:
96 raise ValueError("Length of weight and data must be the same.")
97 self.__weight = w
98 self.__data = data
99 except TypeError:
100 self.__weight = [w]
101 self.__data = [data]
102
103 BX = boundingBox(domain)
104 DIM = domain.getDim()
105 x = domain.getX()
106 self.__pde=LinearSinglePDE(domain)
107 self.__pde.getSolverOptions().setTolerance(tol)
108 self.__pde.setSymmetryOn()
109 z=x[DIM-1]
110 q0=whereZero(z-BX[DIM-1][1])
111 if fixPotentialAtBottom: q0+=whereZero(z-BX[DIM-1][0])
112 self.__pde.setValue(q=q0)
113
114 self.edge_lengths=np.asarray(boundingBoxEdgeLengths(domain))
115 self.diameter=1./sqrt(sum(1./self.edge_lengths**2))
116
117 self.__origweight=[]
118 for s in range(len(self.__weight)):
119 # save a copy of the original weights in case of rescaling
120 self.__origweight.append(1.*self.__weight[s])
121
122 if not self.__trafo.isCartesian():
123 fd=1./self.__trafo.getScalingFactors()
124 fw=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
125 for s in range(len(self.__weight)):
126 self.__weight[s] = fw * self.__weight[s]
127 self.__data[s] = fd * self.__data[s]
128
129 def _rescaleWeights(self, scale=1., fetch_factor=1.):
130 """
131 rescales the weights such that
132
133 *sum_s integrate( ( weight_i[s] *data_i[s]) (weight_j[s]*1/L_j) * L**2 * fetch_factor )=scale*
134 """
135 if not scale > 0:
136 raise ValueError("Value for scale must be positive.")
137 A=0
138 # copy back original weights before rescaling
139 self.__weight=[1.*ow for ow in self.__origweight]
140
141 for s in range(len(self.__weight)):
142 if self.__data[s].getShape() == ():
143 ff=self.__weight[s]**2*self.__data[s]/length(self.edge_lengths)
144 else:
145 ff=inner(self.__weight[s], self.__data[s]) * inner(self.__weight[s], 1/self.edge_lengths)
146 A += integrate(abs(ff * fetch_factor))
147 if A > 0:
148 A=sqrt(scale/A)/self.diameter
149 if not self.__trafo.isCartesian():
150 A*=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
151 for s in range(len(self.__weight)):
152 self.__weight[s]*=A
153 else:
154 raise ValueError("Rescaling of weights failed.")
155
156 def getDomain(self):
157 """
158 Returns the domain of the forward model.
159
160 :rtype: `Domain`
161 """
162 return self.__domain
163
164 def getMisfitWeights(self):
165 """
166 Returns the weights of the misfit function
167
168 :rtype: ``list`` of ``Data``
169 """
170 return self.__weight
171
172 def getData(self):
173 """
174 Returns the data
175
176 :rtype: ``list`` of ``Data``
177 """
178 return self.__data
179 def getDataFunctionSpace(self):
180 """
181 Returns the ``FunctionSpace`` of the data
182
183 :rtype: ``FunctionSpace``
184 """
185 return self.getData()[0].getFunctionSpace()
186
187 def getCoordinateTransformation(self):
188 """
189 returns the coordinate transformation being used
190
191 :rtype: ``CoordinateTransformation``
192 """
193 return self.__trafo
194
195 def getPDE(self):
196 """
197 Return the underlying PDE.
198
199 :rtype: `LinearPDE`
200 """
201 return self.__pde
202
203 def _getDefect(self, result):
204 """
205 Returns the defect value.
206
207 :param result: a result vector
208 :type result: `Vector`
209 :rtype: ``float``
210 """
211 A=0.
212 for s in range(len(self.__weight)):
213 A += integrate( inner(self.__weight[s], self.__data[s]-result)**2 )
214 return A/2
215
216 def getDefectGradient(self, result):
217 Y=0.
218 for s in range(len(self.__weight)):
219 Y = inner(self.__weight[s], self.__data[s]-result) * self.__weight[s] + Y
220 return Y
221
222 def getSurvey(self, index=None):
223 """
224 Returns the pair (data_index, weight_index), where data_i is the data
225 of survey i, weight_i is the weighting factor for survey i.
226 If index is None, all surveys will be returned in a pair of lists.
227 """
228 if index is None:
229 return self.__data, self.__weight
230 if index>=len(self.__data):
231 raise IndexError("Forward model only has %d surveys"%len(self.__data))
232 return self.__data[index], self.__weight[index]
233

  ViewVC Help
Powered by ViewVC 1.1.26