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

  ViewVC Help
Powered by ViewVC 1.1.26