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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4777 - (hide annotations)
Wed Mar 19 23:52:36 2014 UTC (5 years, 6 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 gross 4737
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 sshaw 4777 return tuple([ n+1 for n in range(self.numComps)] )
43 gross 4737
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 sshaw 4777 for n in range(self.numComps):
50 gross 4737 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 sshaw 4777 for n in range(self.numComps):
59 gross 4737 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 sshaw 4777 return sum(((n+1)*s[n]-(n+2))**2 for n in range(self.numComps) )*.5
67 gross 4737
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 sshaw 4748 self.domain = esys.ripley.Rectangle(20*getMPISizeWorld()-1,20,d0=getMPISizeWorld()) # expected dimension 1mx1m
87 gross 4737 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