/[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 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 7263 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
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-2015 by The University of Queensland
21 http://www.uq.edu.au
22 Primary Business: Queensland, Australia"""
23 __license__="""Licensed under the Open Software License version 3.0
24 http://www.opensource.org/licenses/osl-3.0.php"""
25 __url__="https://launchpad.net/escript-finley"
26
27 __all__ = ['ForwardModel','ForwardModelWithPotential']
28
29 from esys.downunder.coordinates import makeTranformation
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 = makeTranformation(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 A += integrate(abs(inner(self.__weight[s], self.__data[s]) * inner(self.__weight[s], 1/self.edge_lengths) * fetch_factor))
143 if A > 0:
144 A=sqrt(scale/A)/self.diameter
145 if not self.__trafo.isCartesian():
146 A*=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
147 for s in range(len(self.__weight)):
148 self.__weight[s]*=A
149 else:
150 raise ValueError("Rescaling of weights failed.")
151
152 def getDomain(self):
153 """
154 Returns the domain of the forward model.
155
156 :rtype: `Domain`
157 """
158 return self.__domain
159
160 def getCoordinateTransformation(self):
161 """
162 returns the coordinate transformation being used
163
164 :rtype: ``CoordinateTransformation``
165 """
166 return self.__trafo
167
168 def getPDE(self):
169 """
170 Return the underlying PDE.
171
172 :rtype: `LinearPDE`
173 """
174 return self.__pde
175
176 def _getDefect(self, result):
177 """
178 Returns the defect value.
179
180 :param result: a result vector
181 :type result: `Vector`
182 :rtype: ``float``
183 """
184 A=0.
185 for s in range(len(self.__weight)):
186 A += integrate( inner(self.__weight[s], self.__data[s]-result)**2 )
187 return A/2
188
189 def getDefectGradient(self, result):
190 Y=0.
191 for s in range(len(self.__weight)):
192 Y = inner(self.__weight[s], self.__data[s]-result) * self.__weight[s] + Y
193 return Y
194
195 def getSurvey(self, index=None):
196 """
197 Returns the pair (data_index, weight_index), where data_i is the data
198 of survey i, weight_i is the weighting factor for survey i.
199 If index is None, all surveys will be returned in a pair of lists.
200 """
201 if index is None:
202 return self.__data, self.__weight
203 if index>=len(self.__data):
204 raise IndexError("Forward model only has %d surveys"%len(self.__data))
205 return self.__data[index], self.__weight[index]
206

  ViewVC Help
Powered by ViewVC 1.1.26