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

  ViewVC Help
Powered by ViewVC 1.1.26