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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5288 - (hide annotations)
Tue Dec 2 23:18:40 2014 UTC (4 years, 10 months ago) by sshaw
File MIME type: text/x-python
File size: 10862 byte(s)
fixing tests for cases where required domains not built
1 sshaw 5288 from __future__ import print_function
2 gross 5001 ##############################################################################
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 sshaw 5288 try:
32     import esys.ripley
33     HAVE_RIPLEY = True
34     except ImportError:
35     HAVE_RIPLEY = False
36    
37     @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
38 sshaw 5002 class Test_Regularizaton2D(unittest.TestCase):
39 gross 5001 def setUp(self):
40 sshaw 5017 self.domain = esys.ripley.Rectangle(20*getMPISizeWorld()-1,20*getMPISizeWorld()-1, d0=getMPISizeWorld()) # expected dimension 1mx1m
41 gross 5001 self.COORDINATES=CartesianReferenceSystem()
42    
43     def tearDown(self):
44     del self.domain
45 gross 5005
46     def test_ConstantLevelSet1(self):
47 gross 5004 a=[1,2]
48    
49     reg=Regularization(self.domain, numLevelSets=1,
50 gross 5001 w1=[1,1.], # consider gradient terms
51     #wc=[[0,1],[0,0]], # and cross-gradient term
52     coordinates=self.COORDINATES)
53     #location_of_set_m=reg_mask)
54    
55    
56 gross 5004 m=Scalar(1., Solution(self.domain))
57     args=reg.getArguments(m)
58     df=reg.getValue(m, *args)
59 gross 5001 self.assertIsInstance(df, float)
60 gross 5004 self.assertAlmostEqual(df,0.)
61 gross 5001
62 gross 5004 gf=reg.getGradient(m, *args)
63     self.assertIsInstance(gf[0], Data)
64     self.assertIsInstance(gf[1], Data)
65     self.assertEqual(gf[0].getFunctionSpace(), Function(self.domain))
66     self.assertEqual(gf[1].getFunctionSpace(), Function(self.domain))
67     self.assertEqual(gf[0].getShape(), ())
68     self.assertEqual(gf[1].getShape(), (2,))
69     self.assertAlmostEqual(Lsup(gf[0]),0.)
70     self.assertAlmostEqual(Lsup(gf[1]),0.)
71 gross 5001
72 gross 5004 # lets try another m
73     x=Solution(self.domain).getX()
74     m=a[0]*x[0]+a[1]*x[1]
75    
76     args=reg.getArguments(m)
77     df=reg.getValue(m, *args)
78     self.assertAlmostEqual(df,(a[0]**2+a[1]**2)/2./2.)
79     gf=reg.getGradient(m, *args)
80     # and now the derivatives:
81 gross 5001 STEP=0.001
82 gross 5004 dm=STEP*x[0]**2
83     m2=m+dm
84     df2=reg.getValue(m2, *(reg.getArguments(m2)))
85 gross 5001 gf2=df2-df
86 gross 5004 self.assertTrue( abs(gf2-reg.getDualProduct(dm, gf)) < 1e-3 * abs(gf2) )
87 gross 5005
88     def test_ConstantLevelSet1Hessian(self):
89     x=Solution(self.domain).getX()
90    
91    
92    
93     reg=Regularization(self.domain, numLevelSets=1,
94     w1=[1,1.], # consider gradient terms
95     #wc=[[0,1],[0,0]], # and cross-gradient term
96     coordinates=self.COORDINATES,
97     location_of_set_m=whereZero(x[0]-inf(x[0])) )
98    
99     x=Solution(self.domain).getX()
100     m=x[0]
101    
102     gf=reg.getGradient(m, *(reg.getArguments(m)))
103    
104     STEP=0.001
105     dm=STEP*x[0]*(x[0]-1)
106     m2=m+dm
107     gf2=reg.getGradient(m2, *(reg.getArguments(m2)))
108    
109     p=reg.getInverseHessianApproximation(m, gf2-gf, *(reg.getArguments(m)))
110     self.assertAlmostEqual(Lsup(p-dm)/Lsup(dm),0.)
111    
112 gross 5001
113 gross 5005 def test_ConstantLevelSet2_noCrossGradient(self):
114 gross 5004 a=[1,2]
115    
116     reg=Regularization(self.domain, numLevelSets=2,
117 gross 5001 w1=[[1,1.],[1.,1.]], # consider gradient terms
118     wc=[[0,0],[0,0]], # and cross-gradient term
119     coordinates=self.COORDINATES)
120     #location_of_set_m=reg_mask)
121    
122    
123 gross 5004 m=Data([1.,2], Solution(self.domain))
124     args=reg.getArguments(m)
125 gross 5001
126 gross 5004 df=reg.getValue(m, *args)
127    
128 gross 5001 self.assertIsInstance(df, float)
129 gross 5004 self.assertAlmostEqual(df,0.)
130 gross 5001
131 gross 5004 gf=reg.getGradient(m, *args)
132     self.assertIsInstance(gf[0], Data)
133     self.assertIsInstance(gf[1], Data)
134     self.assertEqual(gf[0].getFunctionSpace(), Function(self.domain))
135     self.assertEqual(gf[1].getFunctionSpace(), Function(self.domain))
136     self.assertEqual(gf[0].getShape(), (2,))
137     self.assertEqual(gf[1].getShape(), (2,2))
138     self.assertAlmostEqual(Lsup(gf[0]),0.)
139     self.assertAlmostEqual(Lsup(gf[1]),0.)
140 gross 5001
141 gross 5004 # lets try another m
142     x=Solution(self.domain).getX()
143     m=x*a
144 gross 5001
145 gross 5004 args=reg.getArguments(m)
146     df=reg.getValue(m, *args)
147     self.assertAlmostEqual(df,(a[0]**2+a[1]**2)/2./2.)
148     gf=reg.getGradient(m, *args)
149    
150     # and now the derivatives:
151     STEP=0.0002
152     dm=STEP*x[0]*x[1]
153    
154     m2=m+dm*[1,0]
155     df2=reg.getValue(m2, *(reg.getArguments(m2)))
156     gf2=df2-df
157     self.assertTrue( abs(gf2-reg.getDualProduct(dm*[1,0], gf)) < 1e-3 * abs(gf2) )
158 gross 5001
159 gross 5004 m2=m+dm*[0,1]
160     df2=reg.getValue(m2, *(reg.getArguments(m2)))
161 gross 5001 gf2=df2-df
162 gross 5004 self.assertTrue( abs(gf2-reg.getDualProduct(dm*[0,1], gf)) < 1e-3 * abs(gf2) )
163 gross 5001
164 gross 5005 def test_ConstantLevelSet2_noCrossGradientHessian(self):
165     x=Solution(self.domain).getX()
166    
167     reg=Regularization(self.domain, numLevelSets=2,
168     w1=[[1,1.],[1.,1.]], # consider gradient terms
169     wc=[[0,0],[0,0]], # and cross-gradient term
170     coordinates=self.COORDINATES,
171     location_of_set_m=whereZero(x[0]-inf(x[0]))*[1,1] )
172    
173     m=x[0]*[1,0]+x[1]*x[1]*[0,1]
174    
175     gf=reg.getGradient(m, *(reg.getArguments(m)))
176     STEP=0.001
177     dm=STEP*x[0]*(x[0]-1)
178 gross 5006
179 gross 5005 m2=m+dm*[1,0]
180     gf2=reg.getGradient(m2, *(reg.getArguments(m2)))
181     p=reg.getInverseHessianApproximation(m, gf2-gf, *(reg.getArguments(m)))
182     self.assertAlmostEqual(Lsup(p[0]-dm)/Lsup(dm),0.)
183     self.assertAlmostEqual(Lsup(p[1])/Lsup(dm),0.)
184    
185 gross 5006 m2=m+dm*[0,1]
186     gf2=reg.getGradient(m2, *(reg.getArguments(m2)))
187     p=reg.getInverseHessianApproximation(m, gf2-gf, *(reg.getArguments(m)))
188     self.assertAlmostEqual(Lsup(p[1]-dm)/Lsup(dm),0.)
189     self.assertAlmostEqual(Lsup(p[0])/Lsup(dm),0.)
190    
191 gross 5005
192 gross 5004 def test_ConstantLevelSet2_WithCrossGradient(self): # doesn't test the regularization
193     reg=Regularization(self.domain, numLevelSets=2,
194     w1=[[1,1.],[1.,1.]], # consider gradient terms
195     wc=[[0,1],[0,0]], # and cross-gradient term
196     coordinates=self.COORDINATES)
197     #location_of_set_m=reg_mask)
198 gross 5006 reg.setTradeOffFactorsForVariation([1.e-8,1.e-8])
199 gross 5004
200     m=Data([1.,2], Solution(self.domain))
201     args=reg.getArguments(m)
202    
203     df=reg.getValue(m, *args)
204    
205     self.assertIsInstance(df, float)
206     self.assertAlmostEqual(df,0.)
207    
208     gf=reg.getGradient(m, *args)
209     self.assertIsInstance(gf[0], Data)
210     self.assertIsInstance(gf[1], Data)
211     self.assertEqual(gf[0].getFunctionSpace(), Function(self.domain))
212     self.assertEqual(gf[1].getFunctionSpace(), Function(self.domain))
213     self.assertEqual(gf[0].getShape(), (2,))
214     self.assertEqual(gf[1].getShape(), (2,2))
215     self.assertAlmostEqual(Lsup(gf[0]),0.)
216     self.assertAlmostEqual(Lsup(gf[1]),0.)
217    
218     # lets try another m:
219     x=Solution(self.domain).getX()
220     a1=1.
221     a2=2.
222     # for this one cross gradient should not impact on cost function
223     f=0.1
224     m=(a1*x[0]+a2*x[1]) * [1,0] + (f*a1*x[0]+f*a2*x[1]) * [0,1] # cross gradient term is zero!
225     args=reg.getArguments(m)
226     df=reg.getValue(m, *args)
227 gross 5006 #self.assertAlmostEqual(df,(a1**2+a2**2)*(1+f**2)/2./2.)
228     self.assertAlmostEqual(df,0.)
229 gross 5004 gf=reg.getGradient(m, *args)
230 gross 5006 self.assertAlmostEqual(Lsup(gf[0]),0.)
231     self.assertAlmostEqual(Lsup(gf[1]),0.)
232 gross 5004
233     # for this gives maximum impact on on cost function
234     m=(a1*x[0]+a2*x[1]) * [1,0] + (a2*x[0]-a1*x[1]) * [0,1] # cross gradient term is zero!
235     args=reg.getArguments(m)
236     df=reg.getValue(m, *args)
237 gross 5006 self.assertAlmostEqual(df,((a1**2+a2**2)**2/4.)/2.)
238 gross 5004 gf=reg.getGradient(m, *args)
239     # and now the derivatives:
240     STEP=0.002
241     dm=STEP*x[0]*x[1]
242    
243     m2=m+dm*[1,0]
244     df2=reg.getValue(m2, *(reg.getArguments(m2)))
245     gf2=df2-df
246     self.assertTrue( abs(gf2-reg.getDualProduct(dm*[1,0], gf)) < 1e-3 * abs(gf2) )
247 gross 5001
248 gross 5004 m2=m+dm*[0,1]
249     df2=reg.getValue(m2, *(reg.getArguments(m2)))
250     gf2=df2-df
251     self.assertTrue( abs(gf2-reg.getDualProduct(dm*[0,1], gf)) < 1e-3 * abs(gf2) )
252 sshaw 5017
253 gross 5006 def test_ConstantLevelSet2_WithCrossGradientHessian(self):
254     x=Solution(self.domain).getX()
255 gross 5001
256 gross 5006 reg=Regularization(self.domain, numLevelSets=2,
257     w1=[[1,1.],[1.,1.]], # consider gradient terms
258     wc=[[0,1],[0,0]], # and cross-gradient term
259     coordinates=self.COORDINATES,
260     location_of_set_m=whereZero(x[0]-inf(x[0]))*[1,1] )
261 sshaw 5017 reg.setTradeOffFactorsForVariation([1e-10,1e-10])
262 gross 5006 m=x[0]*[1,0]+x[1]*x[1]*[0,1]
263 gross 5004
264 gross 5006 gf=reg.getGradient(m, *(reg.getArguments(m)))
265     STEP=0.001
266     dm=STEP*x[0]*x[1]
267    
268 sshaw 5017 for vector in [[1,0], [0,1]]:
269     m2=m+dm*vector
270     gf2=reg.getGradient(m2, *(reg.getArguments(m2)))
271 gross 5006
272 sshaw 5017 p=reg.getInverseHessianApproximation(m, gf2-gf, *(reg.getArguments(m)),
273     solve=False)
274     X = (gf2-gf)[1]
275     A = p.getCoefficient("A")
276     res = generalTensorProduct(A, grad(dm*vector), axis_offset=2)
277     self.assertLessEqual(Lsup(res-X), 5e-3 * Lsup(X))
278 gross 5006
279 gross 5001 if __name__ == '__main__':
280     run_tests(__name__, exit_on_failure=True)
281    

  ViewVC Help
Powered by ViewVC 1.1.26