/[escript]/trunk/downunder/test/python/run_inversioncostfunction.py
ViewVC logotype

Contents of /trunk/downunder/test/python/run_inversioncostfunction.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (19 months, 2 weeks ago) by jfenwick
File MIME type: text/x-python
File size: 10825 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 ##############################################################################
2 #
3 # Copyright (c) 2012-2018 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 from __future__ import print_function, division
17
18 __copyright__="""Copyright (c) 2012-2018 by The University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Apache License, version 2.0
22 http://www.apache.org/licenses/LICENSE-2.0"""
23 __url__="https://launchpad.net/escript-finley"
24
25 import esys.escriptcore.utestselect as unittest
26 from esys.escriptcore.testing import *
27 import sys
28 from esys.downunder import *
29 from esys.escript import *
30 import numpy as np
31
32 try:
33 import esys.ripley
34 HAVE_RIPLEY = True
35 except ImportError:
36 HAVE_RIPLEY = False
37
38 class LinearMappingX(Mapping):
39 def __init__(self, mat):
40 self.matrix=mat
41 self.det=mat[0,0]*mat[1,1]-mat[0,1]*mat[0,1]
42 self.inv=np.array([ [mat[1,1]/self.det, -mat[1,0]/self.det ], [-mat[0,1]/self.det, mat[0,0]/self.det] ])
43
44 def getDerivative(self, m):
45 return Data(self.matrix, m.getFunctionSpace())
46 def getValue(self, p):
47 return matrix_mult(self.matrix, p)
48 def getInverse(self, p):
49 return matrix_mult(self.inv, p)
50 def getTypicalDerivative(self):
51 return self.matrix[1,1]
52
53 class SimpleModel(ForwardModel):
54 def __init__(self, domain, coordinates, numComps=1):
55 self.domain = domain
56 self.trafo=makeTransformation(domain, coordinates)
57 self.numComps=numComps
58
59 def getCoordinateTransformation(self):
60 return self.trafo
61
62 def getArguments(self, x):
63 return tuple([ n+1 for n in range(self.numComps)] )
64
65 def getDefect(self, x, *args):
66 if self.numComps == 1:
67 out=(1.*x-2.)*(args[0]*x-2.)
68 else:
69 out=0
70 for n in range(self.numComps):
71 out=out+((n+1)*x[n]-(n+2))*(args[n]*x[n]-(n+2))
72 return integrate(out, Function(self.domain))/2
73
74 def getGradient(self, x, *args):
75 if self.numComps == 1:
76 Y=(1.*x-2.)*args[0]
77 else:
78 Y=Data(0.,(self.numComps,), Function(self.domain))
79 for n in range(self.numComps):
80 Y[n]=((n+1)*x[n]-(n+2))*args[n]
81 return Y
82
83 def getRealValue(self, s):
84 if self.numComps == 1:
85 return (s-2)**2 * 0.5
86 else:
87 return sum(((n+1)*s[n]-(n+2))**2 for n in range(self.numComps) )*.5
88
89 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
90 class TestInversionCostfunction(unittest.TestCase):
91 def setUp(self):
92 self.domain = esys.ripley.Rectangle(20*getMPISizeWorld()-1,20,d0=getMPISizeWorld()) # expected dimension 1mx1m
93 self.COORDINATES=CartesianReferenceSystem()
94
95 def tearDown(self):
96 del self.domain
97
98 # standart inversion case
99 def test_ConstantLevelSet1_Mapping1_Model1(self): # doesn't test the regularization
100 dom=self.domain
101 s0=10.
102
103 pm1=LinearMapping(s0, 0.) #p = a * m + p0
104 sm1=SimpleModel(dom, numComps=1, coordinates=self.COORDINATES)
105 reg=Regularization(dom, numLevelSets=1,
106 w1=[1,1.], # consider gradient terms
107 #wc=[[0,1],[0,0]], # and cross-gradient term
108 coordinates=self.COORDINATES)
109 #location_of_set_m=reg_mask)
110 cf=InversionCostFunction(reg, pm1, sm1)
111
112
113 m=Scalar(1., Solution(dom))
114 args=cf.getArguments(m)
115 df=cf.getValue(m, *args)
116 gf=cf.getGradient(m, *args)[0]
117
118 self.assertIsInstance(df, float)
119 self.assertAlmostEqual(df, sm1.getRealValue(s0))
120
121 self.assertIsInstance(gf, Data)
122 self.assertEqual(gf.getFunctionSpace(), Function(dom))
123 self.assertEqual(gf.getShape(), ())
124
125
126 STEP=0.001
127 m2=m+STEP
128 df2=cf.getValue(m2, *(cf.getArguments(m2)))
129 gf2=df2-df
130 self.assertTrue( Lsup(gf2-gf*STEP) < 1e-3 * Lsup(gf2) )
131
132 # strong magnetic-gravity coupling
133 def test_ConstantLevelSet1_Mapping2_Model1(self): # doesn't test the regularization
134 dom=self.domain
135 s=[10., 20]
136
137 pm0=LinearMapping(s[0], 0.) #p = a * m + p0
138 pm1=LinearMapping(s[1], 0.)
139
140 sm1=SimpleModel(dom, numComps=1, coordinates=self.COORDINATES)
141 sm2=SimpleModel(dom, numComps=1, coordinates=self.COORDINATES)
142
143 reg=Regularization(dom, numLevelSets=1,
144 w1=[1,1.], # consider gradient terms
145 #wc=[[0,1],[0,0]], # and cross-gradient term
146 coordinates=self.COORDINATES)
147 #location_of_set_m=reg_mask)
148 cf=InversionCostFunction(reg, [ pm0, pm1 ], [ (sm1,0), (sm2,1) ] )
149
150
151 m=Scalar(1., Solution(dom))
152 args=cf.getArguments(m)
153
154 df=cf.getValue(m, *args)
155 df_real=sm1.getRealValue(s[0])+sm2.getRealValue(s[1])
156 gf=cf.getGradient(m, *args)[0]
157
158 self.assertIsInstance(df, float)
159 self.assertAlmostEqual(df, df_real)
160
161 self.assertIsInstance(gf, Data)
162 self.assertEqual(gf.getFunctionSpace(), Function(dom))
163 self.assertEqual(gf.getShape(), ())
164
165 STEP=0.001
166 m2=m+STEP
167 df2=cf.getValue(m2, *(cf.getArguments(m2)))
168 gf2=df2-df
169 self.assertTrue( Lsup(gf2-gf*STEP) < 1e-3 * Lsup(gf2) )
170
171 # strong magnetic-gravity coupling
172 def test_ConstantLevelSet1_Mapping2_Model1_V2(self): # doesn't test the regularization
173 dom=self.domain
174 s=[10., 20]
175
176 pm0=LinearMapping(s[0], 0.) #p = a * m + p0
177 pm1=LinearMapping(s[1], 0.)
178
179 sm1=SimpleModel(dom, numComps=1, coordinates=self.COORDINATES)
180 sm2=SimpleModel(dom, numComps=1, coordinates=self.COORDINATES)
181
182 reg=Regularization(dom, numLevelSets=1,
183 w1=[1,1.], # consider gradient terms
184 #wc=[[0,1],[0,0]], # and cross-gradient term
185 coordinates=self.COORDINATES)
186 #location_of_set_m=reg_mask)
187 cf=InversionCostFunction(reg, [ (pm0,0), (pm1,0) ], [ (sm1,0), (sm2,1) ] )
188
189
190 m=Scalar(1., Solution(dom))
191 args=cf.getArguments(m)
192
193 df=cf.getValue(m, *args)
194 df_real=sm1.getRealValue(s[0])+sm2.getRealValue(s[1])
195 gf=cf.getGradient(m, *args)[0]
196
197
198 self.assertIsInstance(df, float)
199 self.assertAlmostEqual(df, df_real)
200
201 self.assertIsInstance(gf, Data)
202 self.assertEqual(gf.getFunctionSpace(), Function(dom))
203 self.assertEqual(gf.getShape(), ())
204
205
206 STEP=0.001
207 m2=m+STEP
208 df2=cf.getValue(m2, *(cf.getArguments(m2)))
209 gf2=df2-df
210 self.assertTrue( Lsup(gf2-gf*STEP) < 1e-3 * Lsup(gf2) )
211
212 # weak gravity- magnetic coupling case:
213 def test_ConstantLevelSet2_Mapping2_Model2(self): # doesn't test the regularization
214 dom=self.domain
215 s=[15., 7.]
216
217 pm0=LinearMapping(s[0], 0.) #p = a * m + p0
218 pm1=LinearMapping(s[1], 0.)
219
220 sm1=SimpleModel(dom, numComps=1, coordinates=self.COORDINATES)
221 sm2=SimpleModel(dom, numComps=1, coordinates=self.COORDINATES)
222
223 reg=Regularization(dom, numLevelSets=2,
224 w1=[[1,1.],[1.,1.]], # consider gradient terms
225 wc=[[0,0],[0,0]], # and cross-gradient term
226 coordinates=self.COORDINATES)
227 #location_of_set_m=reg_mask)
228 cf=InversionCostFunction(reg, [ (pm0,0), (pm1,1) ], [ (sm1,0), (sm2,1) ] )
229
230
231 m=Data([1.,2], Solution(dom))
232 args=cf.getArguments(m)
233
234 df=cf.getValue(m, *args)
235 df_real=sm1.getRealValue(s[0]*1.)+sm2.getRealValue(s[1]*2.)
236 gf=cf.getGradient(m, *args)[0]
237
238 self.assertIsInstance(df, float)
239 self.assertAlmostEqual(df, df_real)
240
241 self.assertIsInstance(gf, Data)
242 self.assertEqual(gf.getFunctionSpace(), Function(dom))
243 self.assertEqual(gf.getShape(), (2,))
244
245
246 STEP=0.001
247
248 m2=m+[STEP,0]
249 df2=cf.getValue(m2, *(cf.getArguments(m2)))
250 gf2=df2-df
251 self.assertTrue( Lsup(gf2-gf[0]*STEP) < 1e-3 * Lsup(gf2) )
252
253 m2=m+[0, STEP]
254 df2=cf.getValue(m2, *(cf.getArguments(m2)))
255 gf2=df2-df
256 self.assertTrue( Lsup(gf2-gf[1]*STEP) < 1e-3 * Lsup(gf2) )
257
258 # acoustic Wave form inversion case
259 def test_ConstantLevelSet2_Mapping1X2_Model1(self): # doesn't test the regularization
260 dom=self.domain
261 s=np.array([ [15., 7.] , [-4, 5.] ] )
262
263 pm0=LinearMappingX(s)
264
265 sm1=SimpleModel(dom, numComps=2, coordinates=self.COORDINATES)
266
267 reg=Regularization(dom, numLevelSets=2,
268 w1=[[1,1.],[1.,1.]], # consider gradient terms
269 wc=[[0,0],[0,0]], # and cross-gradient term
270 coordinates=self.COORDINATES)
271 #location_of_set_m=reg_mask)
272 cf=InversionCostFunction(reg, [ (pm0, (0,1))], [ (sm1,0) ] )
273
274
275 m=Data([1.,2], Solution(dom))
276 args=cf.getArguments(m)
277
278 df=cf.getValue(m, *args)
279 df_real=sm1.getRealValue(np.dot(pm0.matrix, np.array([1,2])))
280 gf=cf.getGradient(m, *args)[0]
281
282
283 self.assertIsInstance(df, float)
284 self.assertAlmostEqual(df, df_real)
285
286 self.assertIsInstance(gf, Data)
287 self.assertEqual(gf.getFunctionSpace(), Function(dom))
288 self.assertEqual(gf.getShape(), (2,))
289
290 STEP=0.001
291 m2=m+[STEP,0]
292 df2=cf.getValue(m2, *(cf.getArguments(m2)))
293 gf2=df2-df
294 self.assertTrue( Lsup(gf2-gf[0]*STEP) < 1e-3 * Lsup(gf2) )
295
296 m2=m+[0, STEP]
297 df2=cf.getValue(m2, *(cf.getArguments(m2)))
298 gf2=df2-df
299 self.assertTrue( Lsup(gf2-gf[1]*STEP) < 1e-3 * Lsup(gf2) )
300
301
302 if __name__ == '__main__':
303 run_tests(__name__, exit_on_failure=True)
304

  ViewVC Help
Powered by ViewVC 1.1.26