/[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 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (19 months, 1 week 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 gross 4737 ##############################################################################
2     #
3 jfenwick 6651 # Copyright (c) 2012-2018 by The University of Queensland
4 gross 4737 # http://www.uq.edu.au
5     #
6     # Primary Business: Queensland, Australia
7 jfenwick 6112 # Licensed under the Apache License, version 2.0
8     # http://www.apache.org/licenses/LICENSE-2.0
9 gross 4737 #
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 sshaw 5706 from __future__ import print_function, division
17    
18 jfenwick 6651 __copyright__="""Copyright (c) 2012-2018 by The University of Queensland
19 gross 4737 http://www.uq.edu.au
20     Primary Business: Queensland, Australia"""
21 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
22     http://www.apache.org/licenses/LICENSE-2.0"""
23 gross 4737 __url__="https://launchpad.net/escript-finley"
24    
25 jfenwick 4938 import esys.escriptcore.utestselect as unittest
26 sshaw 4984 from esys.escriptcore.testing import *
27 gross 4737 import sys
28     from esys.downunder import *
29     from esys.escript import *
30     import numpy as np
31    
32 sshaw 5288 try:
33     import esys.ripley
34     HAVE_RIPLEY = True
35     except ImportError:
36     HAVE_RIPLEY = False
37 gross 4737
38 sshaw 5031 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 gross 4737 class SimpleModel(ForwardModel):
54     def __init__(self, domain, coordinates, numComps=1):
55     self.domain = domain
56 sshaw 5780 self.trafo=makeTransformation(domain, coordinates)
57 gross 4737 self.numComps=numComps
58    
59     def getCoordinateTransformation(self):
60     return self.trafo
61    
62     def getArguments(self, x):
63 sshaw 4777 return tuple([ n+1 for n in range(self.numComps)] )
64 gross 4737
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 sshaw 4777 for n in range(self.numComps):
71 gross 4737 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 sshaw 4777 for n in range(self.numComps):
80 gross 4737 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 sshaw 4777 return sum(((n+1)*s[n]-(n+2))**2 for n in range(self.numComps) )*.5
88 gross 5004
89 sshaw 5288 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
90 gross 4737 class TestInversionCostfunction(unittest.TestCase):
91     def setUp(self):
92 sshaw 4748 self.domain = esys.ripley.Rectangle(20*getMPISizeWorld()-1,20,d0=getMPISizeWorld()) # expected dimension 1mx1m
93 gross 4737 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 sshaw 4984 if __name__ == '__main__':
303     run_tests(__name__, exit_on_failure=True)
304 gross 4737

  ViewVC Help
Powered by ViewVC 1.1.26