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

Contents of /trunk/downunder/test/python/run_regularization.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: 10862 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 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
38 class Test_Regularizaton2D(unittest.TestCase):
39 def setUp(self):
40 self.domain = esys.ripley.Rectangle(20*getMPISizeWorld()-1,20*getMPISizeWorld()-1, d0=getMPISizeWorld()) # expected dimension 1mx1m
41 self.COORDINATES=CartesianReferenceSystem()
42
43 def tearDown(self):
44 del self.domain
45
46 def test_ConstantLevelSet1(self):
47 a=[1,2]
48
49 reg=Regularization(self.domain, numLevelSets=1,
50 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 m=Scalar(1., Solution(self.domain))
57 args=reg.getArguments(m)
58 df=reg.getValue(m, *args)
59 self.assertIsInstance(df, float)
60 self.assertAlmostEqual(df,0.)
61
62 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
72 # 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 STEP=0.001
82 dm=STEP*x[0]**2
83 m2=m+dm
84 df2=reg.getValue(m2, *(reg.getArguments(m2)))
85 gf2=df2-df
86 self.assertTrue( abs(gf2-reg.getDualProduct(dm, gf)) < 1e-3 * abs(gf2) )
87
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
113 def test_ConstantLevelSet2_noCrossGradient(self):
114 a=[1,2]
115
116 reg=Regularization(self.domain, numLevelSets=2,
117 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 m=Data([1.,2], Solution(self.domain))
124 args=reg.getArguments(m)
125
126 df=reg.getValue(m, *args)
127
128 self.assertIsInstance(df, float)
129 self.assertAlmostEqual(df,0.)
130
131 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
141 # lets try another m
142 x=Solution(self.domain).getX()
143 m=x*a
144
145 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
159 m2=m+dm*[0,1]
160 df2=reg.getValue(m2, *(reg.getArguments(m2)))
161 gf2=df2-df
162 self.assertTrue( abs(gf2-reg.getDualProduct(dm*[0,1], gf)) < 1e-3 * abs(gf2) )
163
164 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
179 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 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
192 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 reg.setTradeOffFactorsForVariation([1.e-8,1.e-8])
199
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 #self.assertAlmostEqual(df,(a1**2+a2**2)*(1+f**2)/2./2.)
228 self.assertAlmostEqual(df,0.)
229 gf=reg.getGradient(m, *args)
230 self.assertAlmostEqual(Lsup(gf[0]),0.)
231 self.assertAlmostEqual(Lsup(gf[1]),0.)
232
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 self.assertAlmostEqual(df,((a1**2+a2**2)**2/4.)/2.)
238 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
248 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
253 def test_ConstantLevelSet2_WithCrossGradientHessian(self):
254 x=Solution(self.domain).getX()
255
256 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 reg.setTradeOffFactorsForVariation([1e-10,1e-10])
262 m=x[0]*[1,0]+x[1]*x[1]*[0,1]
263
264 gf=reg.getGradient(m, *(reg.getArguments(m)))
265 STEP=0.001
266 dm=STEP*x[0]*x[1]
267
268 for vector in [[1,0], [0,1]]:
269 m2=m+dm*vector
270 gf2=reg.getGradient(m2, *(reg.getArguments(m2)))
271
272 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
279 if __name__ == '__main__':
280 run_tests(__name__, exit_on_failure=True)
281

  ViewVC Help
Powered by ViewVC 1.1.26