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

  ViewVC Help
Powered by ViewVC 1.1.26