/[escript]/trunk/escript/py_src/generateutil
ViewVC logotype

Contents of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


Revision 785 - (show annotations)
Tue Jul 25 03:48:10 2006 UTC (13 years, 2 months ago) by gross
File size: 164408 byte(s)
some new functions around tensor product added.
1 #!/usr/bin/python
2 # $Id$
3
4 """
5 program generates parts of the util.py and the test_util.py script
6 """
7 test_header=""
8 test_header+="import unittest\n"
9 test_header+="import numarray\n"
10 test_header+="from esys.escript import *\n"
11 test_header+="from esys.finley import Rectangle\n"
12 test_header+="class Test_util2(unittest.TestCase):\n"
13 test_header+=" RES_TOL=1.e-7\n"
14 test_header+=" def setUp(self):\n"
15 test_header+=" self.__dom =Rectangle(11,11,2)\n"
16 test_header+=" self.functionspace = FunctionOnBoundary(self.__dom)\n"
17 test_tail=""
18 test_tail+="suite = unittest.TestSuite()\n"
19 test_tail+="suite.addTest(unittest.makeSuite(Test_util2))\n"
20 test_tail+="unittest.TextTestRunner(verbosity=2).run(suite)\n"
21
22 case_set=["float","array","constData","taggedData","expandedData","Symbol"]
23 shape_set=[ (),(2,), (4,5), (6,2,2),(3,2,3,4)]
24
25 t_prog=""
26 t_prog_with_tags=""
27 t_prog_failing=""
28 u_prog=""
29
30 def wherepos(arg):
31 if arg>0.:
32 return 1.
33 else:
34 return 0.
35
36
37 class OPERATOR:
38 def __init__(self,nickname,rng=[-1000.,1000],test_expr="",math_expr=None,
39 numarray_expr="",symbol_expr=None,diff=None,name=""):
40 self.nickname=nickname
41 self.rng=rng
42 self.test_expr=test_expr
43 if math_expr==None:
44 self.math_expr=test_expr.replace("%a1%","arg")
45 else:
46 self.math_expr=math_expr
47 self.numarray_expr=numarray_expr
48 self.symbol_expr=symbol_expr
49 self.diff=diff
50 self.name=name
51
52 import random
53 import numarray
54 import math
55 finc=1.e-6
56
57 def getResultCaseForBin(case0,case1):
58 if case0=="float":
59 if case1=="float":
60 case="float"
61 elif case1=="array":
62 case="array"
63 elif case1=="constData":
64 case="constData"
65 elif case1=="taggedData":
66 case="taggedData"
67 elif case1=="expandedData":
68 case="expandedData"
69 elif case1=="Symbol":
70 case="Symbol"
71 else:
72 raise ValueError,"unknown case1=%s"%case1
73 elif case0=="array":
74 if case1=="float":
75 case="array"
76 elif case1=="array":
77 case="array"
78 elif case1=="constData":
79 case="constData"
80 elif case1=="taggedData":
81 case="taggedData"
82 elif case1=="expandedData":
83 case="expandedData"
84 elif case1=="Symbol":
85 case="Symbol"
86 else:
87 raise ValueError,"unknown case1=%s"%case1
88 elif case0=="constData":
89 if case1=="float":
90 case="constData"
91 elif case1=="array":
92 case="constData"
93 elif case1=="constData":
94 case="constData"
95 elif case1=="taggedData":
96 case="taggedData"
97 elif case1=="expandedData":
98 case="expandedData"
99 elif case1=="Symbol":
100 case="Symbol"
101 else:
102 raise ValueError,"unknown case1=%s"%case1
103 elif case0=="taggedData":
104 if case1=="float":
105 case="taggedData"
106 elif case1=="array":
107 case="taggedData"
108 elif case1=="constData":
109 case="taggedData"
110 elif case1=="taggedData":
111 case="taggedData"
112 elif case1=="expandedData":
113 case="expandedData"
114 elif case1=="Symbol":
115 case="Symbol"
116 else:
117 raise ValueError,"unknown case1=%s"%case1
118 elif case0=="expandedData":
119 if case1=="float":
120 case="expandedData"
121 elif case1=="array":
122 case="expandedData"
123 elif case1=="constData":
124 case="expandedData"
125 elif case1=="taggedData":
126 case="expandedData"
127 elif case1=="expandedData":
128 case="expandedData"
129 elif case1=="Symbol":
130 case="Symbol"
131 else:
132 raise ValueError,"unknown case1=%s"%case1
133 elif case0=="Symbol":
134 if case1=="float":
135 case="Symbol"
136 elif case1=="array":
137 case="Symbol"
138 elif case1=="constData":
139 case="Symbol"
140 elif case1=="taggedData":
141 case="Symbol"
142 elif case1=="expandedData":
143 case="Symbol"
144 elif case1=="Symbol":
145 case="Symbol"
146 else:
147 raise ValueError,"unknown case1=%s"%case1
148 else:
149 raise ValueError,"unknown case0=%s"%case0
150 return case
151
152
153 def makeArray(shape,rng):
154 l=rng[1]-rng[0]
155 out=numarray.zeros(shape,numarray.Float64)
156 if isinstance(l,int):
157 if len(shape)==0:
158 out=int(l*random.random()+rng[0])*1.
159 elif len(shape)==1:
160 for i0 in range(shape[0]):
161 out[i0]=int(l*random.random()+rng[0])
162 elif len(shape)==2:
163 for i0 in range(shape[0]):
164 for i1 in range(shape[1]):
165 out[i0,i1]=int(l*random.random()+rng[0])
166 elif len(shape)==3:
167 for i0 in range(shape[0]):
168 for i1 in range(shape[1]):
169 for i2 in range(shape[2]):
170 out[i0,i1,i2]=int(l*random.random()+rng[0])
171 elif len(shape)==4:
172 for i0 in range(shape[0]):
173 for i1 in range(shape[1]):
174 for i2 in range(shape[2]):
175 for i3 in range(shape[3]):
176 out[i0,i1,i2,i3]=int(l*random.random()+rng[0])
177 elif len(shape)==5:
178 for i0 in range(shape[0]):
179 for i1 in range(shape[1]):
180 for i2 in range(shape[2]):
181 for i3 in range(shape[3]):
182 for i4 in range(shape[4]):
183 out[i0,i1,i2,i3,i4]=int(l*ranm.random()+rng[0])
184 else:
185 raise SystemError,"rank is restricted to 5"
186 else:
187 if len(shape)==0:
188 out=l*random.random()+rng[0]
189 elif len(shape)==1:
190 for i0 in range(shape[0]):
191 out[i0]=l*random.random()+rng[0]
192 elif len(shape)==2:
193 for i0 in range(shape[0]):
194 for i1 in range(shape[1]):
195 out[i0,i1]=l*random.random()+rng[0]
196 elif len(shape)==3:
197 for i0 in range(shape[0]):
198 for i1 in range(shape[1]):
199 for i2 in range(shape[2]):
200 out[i0,i1,i2]=l*random.random()+rng[0]
201 elif len(shape)==4:
202 for i0 in range(shape[0]):
203 for i1 in range(shape[1]):
204 for i2 in range(shape[2]):
205 for i3 in range(shape[3]):
206 out[i0,i1,i2,i3]=l*random.random()+rng[0]
207 elif len(shape)==5:
208 for i0 in range(shape[0]):
209 for i1 in range(shape[1]):
210 for i2 in range(shape[2]):
211 for i3 in range(shape[3]):
212 for i4 in range(shape[4]):
213 out[i0,i1,i2,i3,i4]=l*ranm.random()+rng[0]
214 else:
215 raise SystemError,"rank is restricted to 5"
216 return out
217
218 def makeNumberedArray(shape,s=1.):
219 out=numarray.zeros(shape,numarray.Float64)
220 if len(shape)==0:
221 out=s*1.
222 elif len(shape)==1:
223 for i0 in range(shape[0]):
224 out[i0]=s*int(8*random.random()+1)
225 elif len(shape)==2:
226 for i0 in range(shape[0]):
227 for i1 in range(shape[1]):
228 out[i0,i1]=s*int(8*random.random()+1)
229 elif len(shape)==3:
230 for i0 in range(shape[0]):
231 for i1 in range(shape[1]):
232 for i2 in range(shape[2]):
233 out[i0,i1,i2]=s*int(8*random.random()+1)
234 elif len(shape)==4:
235 for i0 in range(shape[0]):
236 for i1 in range(shape[1]):
237 for i2 in range(shape[2]):
238 for i3 in range(shape[3]):
239 out[i0,i1,i2,i3]=s*int(8*random.random()+1)
240 else:
241 raise SystemError,"rank is restricted to 4"
242 return out
243
244 def makeResult(val,test_expr):
245 if isinstance(val,float):
246 out=eval(test_expr.replace("%a1%","val"))
247 elif len(val.shape)==0:
248 out=eval(test_expr.replace("%a1%","val"))
249 elif len(val.shape)==1:
250 out=numarray.zeros(val.shape,numarray.Float64)
251 for i0 in range(val.shape[0]):
252 out[i0]=eval(test_expr.replace("%a1%","val[i0]"))
253 elif len(val.shape)==2:
254 out=numarray.zeros(val.shape,numarray.Float64)
255 for i0 in range(val.shape[0]):
256 for i1 in range(val.shape[1]):
257 out[i0,i1]=eval(test_expr.replace("%a1%","val[i0,i1]"))
258 elif len(val.shape)==3:
259 out=numarray.zeros(val.shape,numarray.Float64)
260 for i0 in range(val.shape[0]):
261 for i1 in range(val.shape[1]):
262 for i2 in range(val.shape[2]):
263 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val[i0,i1,i2]"))
264 elif len(val.shape)==4:
265 out=numarray.zeros(val.shape,numarray.Float64)
266 for i0 in range(val.shape[0]):
267 for i1 in range(val.shape[1]):
268 for i2 in range(val.shape[2]):
269 for i3 in range(val.shape[3]):
270 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val[i0,i1,i2,i3]"))
271 else:
272 raise SystemError,"rank is restricted to 4"
273 return out
274
275 def makeResult2(val0,val1,test_expr):
276 if isinstance(val0,float):
277 if isinstance(val1,float):
278 out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
279 elif len(val1.shape)==0:
280 out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
281 elif len(val1.shape)==1:
282 out=numarray.zeros(val1.shape,numarray.Float64)
283 for i0 in range(val1.shape[0]):
284 out[i0]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0]"))
285 elif len(val1.shape)==2:
286 out=numarray.zeros(val1.shape,numarray.Float64)
287 for i0 in range(val1.shape[0]):
288 for i1 in range(val1.shape[1]):
289 out[i0,i1]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1]"))
290 elif len(val1.shape)==3:
291 out=numarray.zeros(val1.shape,numarray.Float64)
292 for i0 in range(val1.shape[0]):
293 for i1 in range(val1.shape[1]):
294 for i2 in range(val1.shape[2]):
295 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2]"))
296 elif len(val1.shape)==4:
297 out=numarray.zeros(val1.shape,numarray.Float64)
298 for i0 in range(val1.shape[0]):
299 for i1 in range(val1.shape[1]):
300 for i2 in range(val1.shape[2]):
301 for i3 in range(val1.shape[3]):
302 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2,i3]"))
303 else:
304 raise SystemError,"rank of val1 is restricted to 4"
305 elif len(val0.shape)==0:
306 if isinstance(val1,float):
307 out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
308 elif len(val1.shape)==0:
309 out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
310 elif len(val1.shape)==1:
311 out=numarray.zeros(val1.shape,numarray.Float64)
312 for i0 in range(val1.shape[0]):
313 out[i0]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0]"))
314 elif len(val1.shape)==2:
315 out=numarray.zeros(val1.shape,numarray.Float64)
316 for i0 in range(val1.shape[0]):
317 for i1 in range(val1.shape[1]):
318 out[i0,i1]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1]"))
319 elif len(val1.shape)==3:
320 out=numarray.zeros(val1.shape,numarray.Float64)
321 for i0 in range(val1.shape[0]):
322 for i1 in range(val1.shape[1]):
323 for i2 in range(val1.shape[2]):
324 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2]"))
325 elif len(val1.shape)==4:
326 out=numarray.zeros(val1.shape,numarray.Float64)
327 for i0 in range(val1.shape[0]):
328 for i1 in range(val1.shape[1]):
329 for i2 in range(val1.shape[2]):
330 for i3 in range(val1.shape[3]):
331 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2,i3]"))
332 else:
333 raise SystemError,"rank of val1 is restricted to 4"
334 elif len(val0.shape)==1:
335 if isinstance(val1,float):
336 out=numarray.zeros(val0.shape,numarray.Float64)
337 for i0 in range(val0.shape[0]):
338 out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1"))
339 elif len(val1.shape)==0:
340 out=numarray.zeros(val0.shape,numarray.Float64)
341 for i0 in range(val0.shape[0]):
342 out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1"))
343 elif len(val1.shape)==1:
344 out=numarray.zeros(val0.shape,numarray.Float64)
345 for i0 in range(val0.shape[0]):
346 out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0]"))
347 elif len(val1.shape)==2:
348 out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
349 for i0 in range(val0.shape[0]):
350 for j1 in range(val1.shape[1]):
351 out[i0,j1]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1]"))
352 elif len(val1.shape)==3:
353 out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
354 for i0 in range(val0.shape[0]):
355 for j1 in range(val1.shape[1]):
356 for j2 in range(val1.shape[2]):
357 out[i0,j1,j2]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1,j2]"))
358 elif len(val1.shape)==4:
359 out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
360 for i0 in range(val0.shape[0]):
361 for j1 in range(val1.shape[1]):
362 for j2 in range(val1.shape[2]):
363 for j3 in range(val1.shape[3]):
364 out[i0,j1,j2,j3]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1,j2,j3]"))
365 else:
366 raise SystemError,"rank of val1 is restricted to 4"
367 elif len(val0.shape)==2:
368 if isinstance(val1,float):
369 out=numarray.zeros(val0.shape,numarray.Float64)
370 for i0 in range(val0.shape[0]):
371 for i1 in range(val0.shape[1]):
372 out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1"))
373 elif len(val1.shape)==0:
374 out=numarray.zeros(val0.shape,numarray.Float64)
375 for i0 in range(val0.shape[0]):
376 for i1 in range(val0.shape[1]):
377 out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1"))
378 elif len(val1.shape)==1:
379 out=numarray.zeros(val0.shape,numarray.Float64)
380 for i0 in range(val0.shape[0]):
381 for i1 in range(val0.shape[1]):
382 out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0]"))
383 elif len(val1.shape)==2:
384 out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
385 for i0 in range(val0.shape[0]):
386 for i1 in range(val0.shape[1]):
387 out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1]"))
388 elif len(val1.shape)==3:
389 out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
390 for i0 in range(val0.shape[0]):
391 for i1 in range(val0.shape[1]):
392 for j2 in range(val1.shape[2]):
393 out[i0,i1,j2]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1,j2]"))
394 elif len(val1.shape)==4:
395 out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
396 for i0 in range(val0.shape[0]):
397 for i1 in range(val0.shape[1]):
398 for j2 in range(val1.shape[2]):
399 for j3 in range(val1.shape[3]):
400 out[i0,i1,j2,j3]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1,j2,j3]"))
401 else:
402 raise SystemError,"rank of val1 is restricted to 4"
403 elif len(val0.shape)==3:
404 if isinstance(val1,float):
405 out=numarray.zeros(val0.shape,numarray.Float64)
406 for i0 in range(val0.shape[0]):
407 for i1 in range(val0.shape[1]):
408 for i2 in range(val0.shape[2]):
409 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1"))
410 elif len(val1.shape)==0:
411 out=numarray.zeros(val0.shape,numarray.Float64)
412 for i0 in range(val0.shape[0]):
413 for i1 in range(val0.shape[1]):
414 for i2 in range(val0.shape[2]):
415 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1"))
416 elif len(val1.shape)==1:
417 out=numarray.zeros(val0.shape,numarray.Float64)
418 for i0 in range(val0.shape[0]):
419 for i1 in range(val0.shape[1]):
420 for i2 in range(val0.shape[2]):
421 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0]"))
422 elif len(val1.shape)==2:
423 out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
424 for i0 in range(val0.shape[0]):
425 for i1 in range(val0.shape[1]):
426 for i2 in range(val0.shape[2]):
427 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1]"))
428 elif len(val1.shape)==3:
429 out=numarray.zeros(val0.shape,numarray.Float64)
430 for i0 in range(val0.shape[0]):
431 for i1 in range(val0.shape[1]):
432 for i2 in range(val0.shape[2]):
433 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1,i2]"))
434 elif len(val1.shape)==4:
435 out=numarray.zeros(val0.shape+val1.shape[3:],numarray.Float64)
436 for i0 in range(val0.shape[0]):
437 for i1 in range(val0.shape[1]):
438 for i2 in range(val0.shape[2]):
439 for j3 in range(val1.shape[3]):
440 out[i0,i1,i2,j3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1,i2,j3]"))
441 else:
442 raise SystemError,"rank of val1 is rargs[1]estricted to 4"
443 elif len(val0.shape)==4:
444 if isinstance(val1,float):
445 out=numarray.zeros(val0.shape,numarray.Float64)
446 for i0 in range(val0.shape[0]):
447 for i1 in range(val0.shape[1]):
448 for i2 in range(val0.shape[2]):
449 for i3 in range(val0.shape[3]):
450 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1"))
451 elif len(val1.shape)==0:
452 out=numarray.zeros(val0.shape,numarray.Float64)
453 for i0 in range(val0.shape[0]):
454 for i1 in range(val0.shape[1]):
455 for i2 in range(val0.shape[2]):
456 for i3 in range(val0.shape[3]):
457 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1"))
458 elif len(val1.shape)==1:
459 out=numarray.zeros(val0.shape,numarray.Float64)
460 for i0 in range(val0.shape[0]):
461 for i1 in range(val0.shape[1]):
462 for i2 in range(val0.shape[2]):
463 for i3 in range(val0.shape[3]):
464 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0]"))
465 elif len(val1.shape)==2:
466 out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
467 for i0 in range(val0.shape[0]):
468 for i1 in range(val0.shape[1]):
469 for i2 in range(val0.shape[2]):
470 for i3 in range(val0.shape[3]):
471 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0,i1]"))
472 elif len(val1.shape)==3:
473 out=numarray.zeros(val0.shape,numarray.Float64)
474 for i0 in range(val0.shape[0]):
475 for i1 in range(val0.shape[1]):
476 for i2 in range(val0.shape[2]):
477 for i3 in range(val0.shape[3]):
478 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0,i1,i2]"))
479 elif len(val1.shape)==4:
480 out=numarray.zeros(val0.shape,numarray.Float64)
481 for i0 in range(val0.shape[0]):
482 for i1 in range(val0.shape[1]):
483 for i2 in range(val0.shape[2]):
484 for i3 in range(val0.shape[3]):
485 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0,i1,i2,i3]"))
486 else:
487 raise SystemError,"rank of val1 is restricted to 4"
488 else:
489 raise SystemError,"rank is restricted to 4"
490 return out
491
492
493 def mkText(case,name,a,a1=None,use_tagging_for_expanded_data=False):
494 t_out=""
495 if case=="float":
496 if isinstance(a,float):
497 t_out+=" %s=%s\n"%(name,a)
498 elif a.rank==0:
499 t_out+=" %s=%s\n"%(name,a)
500 else:
501 t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
502 elif case=="array":
503 if isinstance(a,float):
504 t_out+=" %s=numarray.array(%s)\n"%(name,a)
505 elif a.rank==0:
506 t_out+=" %s=numarray.array(%s)\n"%(name,a)
507 else:
508 t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
509 elif case=="constData":
510 if isinstance(a,float):
511 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
512 elif a.rank==0:
513 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
514 else:
515 t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
516 elif case=="taggedData":
517 if isinstance(a,float):
518 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
519 t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
520 elif a.rank==0:
521 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
522 t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
523 else:
524 t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
525 t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
526 elif case=="expandedData":
527 if use_tagging_for_expanded_data:
528 if isinstance(a,float):
529 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
530 t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
531 elif a.rank==0:
532 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
533 t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
534 else:
535 t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
536 t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
537 t_out+=" %s.expand()\n"%name
538 else:
539 t_out+=" msk_%s=1-whereZero(self.functionspace.getX()[0],1.e-8)\n"%name
540 if isinstance(a,float):
541 t_out+=" %s=msk_%s*(%s)+(1-msk_%s)*(%s)\n"%(name,name,a,name,a1)
542 elif a.rank==0:
543 t_out+=" %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a,name,a1)
544 else:
545 t_out+=" %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a.tolist(),name,a1.tolist())
546 elif case=="Symbol":
547 if isinstance(a,float):
548 t_out+=" %s=Symbol(shape=())\n"%(name)
549 elif a.rank==0:
550 t_out+=" %s=Symbol(shape=())\n"%(name)
551 else:
552 t_out+=" %s=Symbol(shape=%s)\n"%(name,str(a.shape))
553
554 return t_out
555
556 def mkTypeAndShapeTest(case,sh,argstr,name=""):
557 text=""
558 if case=="float":
559 text+=" self.failUnless(isinstance(%s,float),\"wrong type of result%s.\")\n"%(argstr,name)
560 elif case=="array":
561 text+=" self.failUnless(isinstance(%s,numarray.NumArray),\"wrong type of result%s.\")\n"%(argstr,name)
562 text+=" self.failUnlessEqual(%s.shape,%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
563 elif case in ["constData","taggedData","expandedData"]:
564 text+=" self.failUnless(isinstance(%s,Data),\"wrong type of result%s.\")\n"%(argstr,name)
565 text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
566 elif case=="Symbol":
567 text+=" self.failUnless(isinstance(%s,Symbol),\"wrong type of result%s.\")\n"%(argstr,name)
568 text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
569 return text
570
571 def mkCode(txt,args=[],intend=""):
572 s=txt.split("\n")
573 if len(s)>1:
574 out=""
575 for l in s:
576 out+=intend+l+"\n"
577 else:
578 out="%sreturn %s\n"%(intend,txt)
579 c=1
580 for r in args:
581 out=out.replace("%%a%s%%"%c,r)
582 return out
583 #=================================================
584 def testMatrixMult(arg0,arg1,sh_s):
585 return numarray.matrixmultiply(arg0,arg1)
586
587 def shapeTestMatrixMult(sh0,sh1,sh_s):
588 if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2):
589 return sh0+sh_s,sh_s+sh1
590
591 def testTensorMult(arg0,arg1,sh_s):
592 if arg0.rank==2:
593 return numarray.matrixmultiply(arg0,arg1)
594 else:
595 if arg1.rank==4:
596 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
597 for i0 in range(arg0.shape[0]):
598 for i1 in range(arg0.shape[1]):
599 for i2 in range(arg0.shape[2]):
600 for i3 in range(arg0.shape[3]):
601 for j2 in range(arg1.shape[2]):
602 for j3 in range(arg1.shape[3]):
603 out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
604 elif arg1.rank==3:
605 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
606 for i0 in range(arg0.shape[0]):
607 for i1 in range(arg0.shape[1]):
608 for i2 in range(arg0.shape[2]):
609 for i3 in range(arg0.shape[3]):
610 for j2 in range(arg1.shape[2]):
611 out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
612 elif arg1.rank==2:
613 out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
614 for i0 in range(arg0.shape[0]):
615 for i1 in range(arg0.shape[1]):
616 for i2 in range(arg0.shape[2]):
617 for i3 in range(arg0.shape[3]):
618 out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
619 return out
620
621 def shapeTestTensorMult(sh0,sh1,sh_s):
622 if ( len(sh_s)==1 and len(sh0+sh_s)==2 and ( len(sh1+sh_s)==2 or len(sh1+sh_s)==1 )) or (len(sh_s)==2 and len(sh0+sh_s)==4 and (len(sh1+sh_s)==2 or len(sh1+sh_s)==3 or len(sh1+sh_s)==4)):
623 return sh0+sh_s,sh_s+sh1
624
625 def tensorProductTest(arg0,arg1,sh_s):
626 if isinstance(arg0,float):
627 out=numarray.array(arg0*arg1)
628 elif isinstance(arg1,float):
629 out=numarray.array(arg0*arg1)
630 elif len(sh_s)==0:
631 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
632 else:
633 l=len(sh_s)
634 sh0=arg0.shape[:arg0.rank-l]
635 sh1=arg1.shape[l:]
636 ls,l0,l1=1,1,1
637 for i in sh_s: ls*=i
638 for i in sh0: l0*=i
639 for i in sh1: l1*=i
640 out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
641 out2=numarray.zeros((l0,l1),numarray.Float)
642 for i0 in range(l0):
643 for i1 in range(l1):
644 for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
645 out=out2.resize(sh0+sh1)
646 return out
647 def tensorProductShape(sh0,sh1,sh_s):
648 return sh0+sh_s,sh_s+sh1
649
650 def transpose(r,offset):
651 if isinstance(r,float): return r
652 s=r.shape
653 s1=1
654 for i in s[:offset]: s1*=i
655 s2=1
656 for i in s[offset:]: s2*=i
657 out=numarray.reshape(r,(s1,s2))
658 out.transpose()
659 return numarray.resize(out,s[offset:]+s[:offset])
660
661 def shapeTestTMatrixMult(sh0,sh1,sh_s):
662 if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2):
663 return sh_s+sh0,sh_s+sh1
664
665 def testTMatrixMult(arg0,arg1,sh_s):
666 return numarray.matrixmultiply(transpose(arg0,1),arg1)
667
668 def shapeTestTTensorMult(sh0,sh1,sh_s):
669 if ( len(sh_s)==1 and len(sh0+sh_s)==2 and ( len(sh1+sh_s)==2 or len(sh1+sh_s)==1 )) or (len(sh_s)==2 and len(sh0+sh_s)==4 and (len(sh1+sh_s)==2 or len(sh1+sh_s)==3 or len(sh1+sh_s)==4)):
670 return sh_s+sh0,sh_s+sh1
671
672 def testTTensorMult(arg0,arg1,sh_s):
673 if arg0.rank==2:
674 return numarray.matrixmultiply(transpose(arg0,1),arg1)
675 else:
676 if arg1.rank==4:
677 out=numarray.zeros((arg0.shape[2],arg0.shape[3],arg1.shape[2],arg1.shape[3]),numarray.Float)
678 for i0 in range(arg0.shape[2]):
679 for i1 in range(arg0.shape[3]):
680 for i2 in range(arg0.shape[0]):
681 for i3 in range(arg0.shape[1]):
682 for j2 in range(arg1.shape[2]):
683 for j3 in range(arg1.shape[3]):
684 out[i0,i1,j2,j3]+=arg0[i2,i3,i0,i1]*arg1[i2,i3,j2,j3]
685 elif arg1.rank==3:
686 out=numarray.zeros((arg0.shape[2],arg0.shape[3],arg1.shape[2]),numarray.Float)
687 for i0 in range(arg0.shape[2]):
688 for i1 in range(arg0.shape[3]):
689 for i2 in range(arg0.shape[0]):
690 for i3 in range(arg0.shape[1]):
691 for j2 in range(arg1.shape[2]):
692 out[i0,i1,j2]+=arg0[i2,i3,i0,i1]*arg1[i2,i3,j2]
693 elif arg1.rank==2:
694 out=numarray.zeros((arg0.shape[2],arg0.shape[3]),numarray.Float)
695 for i0 in range(arg0.shape[2]):
696 for i1 in range(arg0.shape[3]):
697 for i2 in range(arg0.shape[0]):
698 for i3 in range(arg0.shape[1]):
699 out[i0,i1]+=arg0[i2,i3,i0,i1]*arg1[i2,i3]
700 return out
701
702 def tensorTProductShape(sh0,sh1,sh_s):
703 return sh_s+sh0,sh_s+sh1
704
705 def tensorTProductTest(arg0,arg1,sh_s):
706 if isinstance(arg0,float):
707 out=numarray.array(arg0*arg1)
708 elif isinstance(arg1,float):
709 out=numarray.array(arg0*arg1)
710 elif len(sh_s)==0:
711 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
712 else:
713 l=len(sh_s)
714 sh0=arg0.shape[l:]
715 sh1=arg1.shape[l:]
716 ls,l0,l1=1,1,1
717 for i in sh_s: ls*=i
718 for i in sh0: l0*=i
719 for i in sh1: l1*=i
720 out1=numarray.outerproduct(arg0,arg1).resize((ls,l0,ls,l1))
721 out2=numarray.zeros((l0,l1),numarray.Float)
722 for i0 in range(l0):
723 for i1 in range(l1):
724 for i in range(ls): out2[i0,i1]+=out1[i,i0,i,i1]
725 out=out2.resize(sh0+sh1)
726 return out
727
728 def shapeTestMatrixTMult(sh0,sh1,sh_s):
729 if len(sh_s)==1 and len(sh0+sh_s)==2 and len(sh_s+sh1)==2:
730 return sh0+sh_s,sh1+sh_s
731
732 def testMatrixTMult(arg0,arg1,sh_s):
733 return numarray.matrixmultiply(arg0,transpose(arg1,1))
734
735 def shapeTestTensorTMult(sh0,sh1,sh_s):
736 if ( len(sh_s)==1 and len(sh0+sh_s)==2 and ( len(sh1+sh_s)==2 or len(sh1+sh_s)==1 )) or (len(sh_s)==2 and len(sh0+sh_s)==4 and (len(sh1+sh_s)==2 or len(sh1+sh_s)==3 or len(sh1+sh_s)==4)):
737 return sh0+sh_s,sh1+sh_s
738
739 def testTensorTMult(arg0,arg1,sh_s):
740 if arg0.rank==2:
741 return numarray.matrixmultiply(arg0,transpose(arg1,1))
742 else:
743 if arg1.rank==4:
744 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[0],arg1.shape[1]),numarray.Float)
745 for i0 in range(arg0.shape[0]):
746 for i1 in range(arg0.shape[1]):
747 for i2 in range(arg0.shape[2]):
748 for i3 in range(arg0.shape[3]):
749 for j0 in range(arg1.shape[0]):
750 for j1 in range(arg1.shape[1]):
751 out[i0,i1,j0,j1]+=arg0[i0,i1,i2,i3]*arg1[j0,j1,i2,i3]
752 elif arg1.rank==3:
753 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[0]),numarray.Float)
754 for i0 in range(arg0.shape[0]):
755 for i1 in range(arg0.shape[1]):
756 for i2 in range(arg0.shape[2]):
757 for i3 in range(arg0.shape[3]):
758 for j0 in range(arg1.shape[0]):
759 out[i0,i1,j0]+=arg0[i0,i1,i2,i3]*arg1[j0,i2,i3]
760 elif arg1.rank==2:
761 out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
762 for i0 in range(arg0.shape[0]):
763 for i1 in range(arg0.shape[1]):
764 for i2 in range(arg0.shape[2]):
765 for i3 in range(arg0.shape[3]):
766 out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
767 return out
768
769 def tensorProductTTest(arg0,arg1,sh_s):
770 if isinstance(arg0,float):
771 out=numarray.array(arg0*arg1)
772 elif isinstance(arg1,float):
773 out=numarray.array(arg0*arg1)
774 elif len(sh_s)==0:
775 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
776 else:
777 l=len(sh_s)
778 sh0=arg0.shape[:arg0.rank-l]
779 sh1=arg1.shape[:arg1.rank-l]
780 ls,l0,l1=1,1,1
781 for i in sh_s: ls*=i
782 for i in sh0: l0*=i
783 for i in sh1: l1*=i
784 out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,l1,ls))
785 out2=numarray.zeros((l0,l1),numarray.Float)
786 for i0 in range(l0):
787 for i1 in range(l1):
788 for i in range(ls): out2[i0,i1]+=out1[i0,i,i1,i]
789 out=out2.resize(sh0+sh1)
790 return out
791
792 def tensorProductTShape(sh0,sh1,sh_s):
793 return sh0+sh_s,sh1+sh_s
794
795 #=======================================================================================================
796 # tensor multiply
797 #=======================================================================================================
798 # oper=
799 # oper=
800 for oper in [ ["tensor_mult",testTensorMult,shapeTestTensorMult], \
801 ["generalTensorProduct",tensorProductTest,tensorProductShape], \
802 ["matrix_mult",testMatrixMult,shapeTestMatrixMult], \
803 ["transposed_matrix_mult",testTMatrixMult,shapeTestTMatrixMult], \
804 ["transposed_tensor_mult",testTTensorMult,shapeTestTTensorMult], \
805 ["generalTransposedTensorProduct",tensorTProductTest,tensorTProductShape], \
806 ["matrix_transposed_mult",testMatrixTMult,shapeTestMatrixTMult], \
807 ["tensor_transposed_mult",testTensorTMult,shapeTestTensorTMult], \
808 ["generalTensorTransposedProduct",tensorProductTTest,tensorProductTShape] ]:
809
810 for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
811 for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
812 for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
813 for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
814 for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
815 if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
816 and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
817 arg_shape=oper[2](sh0,sh1,sh_s)
818 if not arg_shape==None:
819 case=getResultCaseForBin(case0,case1)
820 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
821 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
822 if oper[0] in [ "generalTensorProduct", "generalTransposedTensorProduct", "generalTensorTransposedProduct"]:
823 tname="test_%s_%s_rank%s_%s_rank%s_offset%s"%(oper[0],case0,len(sh0+sh_s),case1,len(sh_s+sh1),len(sh_s))
824
825 res_text=" res=%s(arg0,arg1,axis_offset=%s)\n"%(oper[0],len(sh_s))
826 else:
827 tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0+sh_s),case1,len(sh_s+sh1))
828
829 res_text=" res=%s(arg0,arg1)\n"%oper[0]
830
831 text+=" def %s(self):\n"%tname
832 a_0=makeArray(arg_shape[0],[-8,8])
833 if case0 in ["taggedData", "expandedData"]:
834 a1_0=makeArray(arg_shape[0],[-8,8])
835 else:
836 a1_0=a_0
837
838 a_1=makeArray(arg_shape[1],[-8,8])
839 if case1 in ["taggedData", "expandedData"]:
840 a1_1=makeArray(arg_shape[1],[-8,8])
841 else:
842 a1_1=a_1
843 r=oper[1](a_0,a_1,sh_s)
844 r1=oper[1](a1_0,a1_1,sh_s)
845 text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
846 text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
847 text+=res_text
848 if case=="Symbol":
849 c0_res,c1_res=case0,case1
850 subs="{"
851 if case0=="Symbol":
852 text+=mkText("array","s0",a_0,a1_0)
853 subs+="arg0:s0"
854 c0_res="array"
855 if case1=="Symbol":
856 text+=mkText("array","s1",a_1,a1_1)
857 if not subs.endswith("{"): subs+=","
858 subs+="arg1:s1"
859 c1_res="array"
860 subs+="}"
861 text+=" sub=res.substitute(%s)\n"%subs
862 res="sub"
863 text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
864 else:
865 res="res"
866 text+=mkText(case,"ref",r,r1)
867 text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
868 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
869 if case0 == "taggedData" or case1 == "taggedData":
870 t_prog_with_tags+=text
871 else:
872 t_prog+=text
873 print test_header
874 print t_prog
875 # print t_prog_with_tags
876 print test_tail
877 1/0
878
879
880
881 #=======================================================================================================
882 # eigenvalues and eigen vectors 2D:
883 #=======================================================================================================
884 alltests= \
885 [ ("case0",[[0.0, 0.0], [0.0, 0.0]],[0.0, 0.0]) \
886 , ("case3",[[-1.0, 0.0], [0.0, -1.0]],[-1.0, -1.0]) \
887 , ("case5",[[-0.99999999999999967, -6.4606252205695602e-16], [-6.4606252205695602e-16, -0.99999999999999967]],[-1.0, -1.0]) \
888 , ("case6",[[0.0, 0.0], [0.0, 0.0001]],[0.0, 0.0001]) \
889 , ("case7",[[0.0001, 0.0], [0.0, 0.0]],[0.0, 0.0001]) \
890 , ("case8",[[6.0598371831785722e-06, 2.3859213977648625e-05], [2.3859213977648629e-05, 9.3940162816821425e-05]],[0.0, 0.0001]) \
891 , ("case9",[[1.0, 0.0], [0.0, 2.0]],[1.0, 2.0]) \
892 , ("case10",[[2.0, 0.0], [0.0, 1.0]],[1.0, 2.0]) \
893 , ("case11",[[1.0605983718317855, 0.23859213977648688], [0.23859213977648688, 1.9394016281682138]],[1.0, 2.0]) \
894 , ("case12",[[1.0, 0.0], [0.0, 1000000.0]],[1.0, 1000000.0]) \
895 , ("case13",[[1000000.0, 0.0], [0.0, 1.0]],[1.0, 1000000.0]) \
896 , ("case14",[[60599.311233413886, 238591.90118434647], [238591.90118434647, 939401.68876658613]],[1.0, 1000000.0]) \
897 ]
898 dim=2
899
900
901 alltests= \
902 [ ("case0",[[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],[0.0, 0.0, 0.0]) \
903 , ("case5",[[10.0, 0.0, 0.0], [0.0, 10.0, 0.0], [0.0, 0.0, 10.0]],[10.0, 10.0, 10.0]) \
904 , ("case10",[[0.9, 0.0, 0.0], [0.0, 0.9, 0.0], [0.0, 0.0, 1.0]],[0.9, 0.9, 1.0]) \
905 , ("case11",[[0.9, 0.0, 0.0], [0.0, 0.97060899725040983, -0.045555123008643325], [0.0, -0.045555123008643339, 0.92939100274959041]],[0.9, 0.9, 1.0]) \
906 , ("case12",[[0.92694799760252555, 0.0, 0.044368966468320177], [0.0, 0.9, 0.0], [0.044368966468320184, 0.0, 0.97305200239747425]],[0.9, 0.9, 1.0]) \
907 , ("case13",[[1.0, 0.0, 0.0], [0.0, 0.9, 0.], [0.0, 0., 0.9]],[0.9, 0.9, 1.0]) \
908 , ("case14",[[0.92379770619813639, 0.041031106298491521, -0.011396846732439278], [0.041031106298491535, 0.97074428392640366, -0.019650012730342326], [-0.011396846732439236, -0.019650012730342337, 0.90545800987545966]],[0.9, 0.9, 1.0]) \
909 , ("case15",[[1.0, 0.0, 0.0], [0.0, 1.1, 0.0], [0.0, 0.0, 1.1]],[1.0, 1.1, 1.1]) \
910 , ("case17",[[1.0269479976025255, 0.0, 0.044368966468320309], [0.0, 1.1, 0.0], [0.044368966468320295, 0.0, 1.0730520023974743]],[1.0, 1.1, 1.1]) \
911 , ("case18",[[1.1, 0.0, 0.0], [0.0, 1.0153410887977139, -0.036038311201720394], [0.0, -0.036038311201720373, 1.084658911202286]],[1.0, 1.1, 1.1]) \
912 , ("case19",[[1.035487967756175, 0.026317079185831614, -0.039960133424212368], [0.026317079185831618, 1.0892641940924184, 0.016301362071911414], [-0.039960133424212355, 0.016301362071911431, 1.0752478381514063]],[1.0, 1.1, 1.1]) \
913 , ("case20",[[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 3.0]],[1.0, 2.0, 3.0]) \
914 , ("case21",[[1.0, 0.0, 0.0], [0.0, 2.7060899725040968, -0.45555123008643206], [0.0, -0.45555123008643228, 2.2939100274959037]],[1.0, 2.0, 3.0]) \
915 , ("case22",[[1.5389599520505153, 0.0, 0.88737932936638753], [0.0, 2.0, 0.0], [0.88737932936638753, 0.0, 2.4610400479494858]],[1.0, 2.0, 3.0]) \
916 , ("case23",[[3.0, 0.0, 0.0], [0.0, 1.153410887977139, -0.36038311201720391], [0.0, -0.36038311201720391, 1.8465891120228608]],[1.0, 2.0, 3.0]) \
917 , ("case24",[[1.5928567395431172, 0.67348185484323142, -0.51356980156651744], [0.67348185484323153, 2.6000847801882254, -0.033486506584313548], [-0.51356980156651744, -0.033486506584313541, 1.8070584802686565]],[1.0, 2.0, 3.0]) \
918 , ("case25",[[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 30000.0]],[1.0, 2.0, 30000.0]) \
919 , ("case26",[[1.0, 0.0, 0.0], [0.0, 21183.286995177881, -13665.625800132779], [0.0, -13665.625800132779, 8818.7130048221279]],[1.0, 2.0, 30000.0]) \
920 , ("case27",[[8085.1298007817086, 0.0, 13310.246250831115], [0.0, 2.0, 0.0], [13310.246250831115, 0.0, 21915.870199218316]],[1.0, 2.0, 30000.0]) \
921 , ("case28",[[30000.0, 0.0, 0.0], [0.0, 1.153410887977139, -0.36038311201720391], [0.0, -0.36038311201720391, 1.8465891120228608]],[1.0, 2.0, 30000.0]) \
922 , ("case29",[[7140.1907849945546, 12308.774438213351, -3419.2256841313947], [12308.774438213351, 21223.762934183575, -5894.4478052274408], [-3419.2256841313947, -5894.4478052274408, 1639.0462808218595]],[1.0, 2.0, 30000.0]) \
923 ]
924
925 dim=3
926
927 alltests= \
928 [ ("case0",[[0.0]],[0.0]) \
929 , ("case1",[[1.0]],[1.0]) \
930 ]
931 dim=1
932
933 for case in ["constData","taggedData","expandedData"]:
934 n=0
935 while n<len(alltests):
936 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
937 tname="test_eigenvalues_and_eigenvectors_%s_dim%s_%s"%(case,dim,alltests[n][0])
938 text+=" def %s(self):\n"%tname
939 a_0=numarray.array(alltests[n][1],numarray.Float64)
940 ev_0=numarray.array(alltests[n][2],numarray.Float64)
941 if case in ["taggedData", "expandedData"]:
942 if n+1<len(alltests):
943 a1_0=numarray.array(alltests[n+1][1],numarray.Float64)
944 ev1_0=numarray.array(alltests[n+1][2],numarray.Float64)
945 else:
946 a1_0=numarray.array(alltests[0][1],numarray.Float64)
947 ev1_0=numarray.array(alltests[0][2],numarray.Float64)
948 n+=2
949 else:
950 a1_0=a_0
951 ev1_0=ev_0
952 n+=1
953 text+=mkText(case,"arg",a_0,a1_0)
954 text+=" res=eigenvalues_and_eigenvectors(arg)\n"
955 text+=mkText(case,"ref_ev",ev_0,ev1_0)
956 text+=mkTypeAndShapeTest(case,(dim,),"res[0]"," for eigenvalues")
957 text+=mkTypeAndShapeTest(case,(dim,dim),"res[1]"," for eigenvectors")
958 text+=" self.failUnless(Lsup(res[0]-ref_ev)<=self.RES_TOL*Lsup(ref_ev),\"wrong eigenvalues\")\n"
959 for i in range(dim):
960 text+=" self.failUnless(Lsup(matrixmult(arg,res[1][:,%s])-res[0][%s]*res[1][:,%s])<=self.RES_TOL*Lsup(res[0]),\"wrong eigenvector %s\")\n"%(i,i,i,i)
961 text+=" self.failUnless(Lsup(matrixmult(transpose(res[1]),res[1])-kronecker(%s))<=self.RES_TOL,\"eigenvectors are not orthonormal\")\n"%dim
962 if case == "taggedData" :
963 t_prog_with_tags+=text
964 else:
965 t_prog+=text
966 print test_header
967 print t_prog
968 print t_prog_with_tags
969 print test_tail
970 1/0
971 #=======================================================================================================
972 # get slices
973 #=======================================================================================================
974 from esys.escript import *
975 for case0 in ["constData","taggedData","expandedData"]:
976 for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,3,3)]:
977 # get perm:
978 if len(sh0)==2:
979 check=[[1,0]]
980 elif len(sh0)==3:
981 check=[[1,0,2],
982 [1,2,0],
983 [2,1,0],
984 [2,0,2],
985 [0,2,1]]
986 elif len(sh0)==4:
987 check=[[0,1,3,2],
988 [0,2,1,3],
989 [0,2,3,1],
990 [0,3,2,1],
991 [0,3,1,2] ,
992 [1,0,2,3],
993 [1,0,3,2],
994 [1,2,0,3],
995 [1,2,3,0],
996 [1,3,2,0],
997 [1,3,0,2],
998 [2,0,1,3],
999 [2,0,3,1],
1000 [2,1,0,3],
1001 [2,1,3,0],
1002 [2,3,1,0],
1003 [2,3,0,1],
1004 [3,0,1,2],
1005 [3,0,2,1],
1006 [3,1,0,2],
1007 [3,1,2,0],
1008 [3,2,1,0],
1009 [3,2,0,1]]
1010 else:
1011 check=[]
1012
1013 # create the test cases:
1014 processed=[]
1015 l=["R","U","L","P","C","N"]
1016 c=[""]
1017 for i in range(len(sh0)):
1018 tmp=[]
1019 for ci in c:
1020 tmp+=[ci+li for li in l]
1021 c=tmp
1022 # SHUFFLE
1023 c2=[]
1024 while len(c)>0:
1025 i=int(random.random()*len(c))
1026 c2.append(c[i])
1027 del c[i]
1028 c=c2
1029 for ci in c:
1030 t=""
1031 sh=()
1032 sl=()
1033 for i in range(len(ci)):
1034 if ci[i]=="R":
1035 s="%s:%s"%(1,sh0[i]-1)
1036 sl=sl+(slice(1,sh0[i]-1),)
1037 sh=sh+(sh0[i]-2,)
1038 if ci[i]=="U":
1039 s=":%s"%(sh0[i]-1)
1040 sh=sh+(sh0[i]-1,)
1041 sl=sl+(slice(0,sh0[i]-1),)
1042 if ci[i]=="L":
1043 s="2:"
1044 sh=sh+(sh0[i]-2,)
1045 sl=sl+(slice(2,sh0[i]),)
1046 if ci[i]=="P":
1047 s="%s"%(int(sh0[i]/2))
1048 sl=sl+(int(sh0[i]/2),)
1049 if ci[i]=="C":
1050 s=":"
1051 sh=sh+(sh0[i],)
1052 sl=sl+(slice(0,sh0[i]),)
1053 if ci[i]=="N":
1054 s=""
1055 sh=sh+(sh0[i],)
1056 if len(s)>0:
1057 if not t=="": t+=","
1058 t+=s
1059 if len(sl)==1: sl=sl[0]
1060 N_found=False
1061 noN_found=False
1062 process=len(t)>0
1063 for i in ci:
1064 if i=="N":
1065 if not noN_found and N_found: process=False
1066 N_found=True
1067 else:
1068 if N_found: process=False
1069 noNfound=True
1070 # is there a similar one processed allready
1071 if process and ci.find("N")==-1:
1072 for ci2 in processed:
1073 for chi in check:
1074 is_perm=True
1075 for i in range(len(chi)):
1076 if not ci[i]==ci2[chi[i]]: is_perm=False
1077 if is_perm: process=False
1078 # if not process: print ci," rejected"
1079 if process:
1080 processed.append(ci)
1081 for case1 in ["array","constData","taggedData","expandedData"]:
1082 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1083 tname="test_setslice_%s_rank%s_%s_%s"%(case0,len(sh0),case1,ci)
1084 text+=" def %s(self):\n"%tname
1085 a_0=makeNumberedArray(sh0)
1086 if case0 in ["taggedData", "expandedData"]:
1087 a1_0=makeNumberedArray(sh0)
1088 else:
1089 a1_0=a_0*1.
1090
1091 a_1=makeNumberedArray(sh)
1092 if case1 in ["taggedData", "expandedData"]:
1093 a1_1=makeNumberedArray(sh)
1094 else:
1095 a1_1=a_1*1.
1096
1097 text+=mkText(case0,"arg",a_0,a1_0)
1098 text+=mkText(case1,"val",a_1,a1_1)
1099 text+=" arg[%s]=val\n"%t
1100 a_0.__setitem__(sl,a_1)
1101 a1_0.__setitem__(sl,a1_1)
1102 if Lsup(a_0-a1_0)==0.:
1103 text+=mkText("constData","ref",a_0,a1_0)
1104 else:
1105 text+=mkText("expandedData","ref",a_0,a1_0)
1106 text+=" self.failUnless(Lsup(arg-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"
1107
1108 if case0 == "taggedData" or case1 == "taggedData":
1109 t_prog_with_tags+=text
1110 else:
1111 t_prog+=text
1112
1113 print test_header
1114 # print t_prog
1115 print t_prog_with_tags
1116 print test_tail
1117 1/0
1118
1119 #=======================================================================================================
1120 # (non)symmetric part
1121 #=======================================================================================================
1122 from esys.escript import *
1123 for name in ["symmetric", "nonsymmetric"]:
1124 f=1.
1125 if name=="nonsymmetric": f=-1
1126 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1127 for sh0 in [ (3,3), (2,3,2,3)]:
1128 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1129 tname="test_%s_%s_rank%s"%(name,case0,len(sh0))
1130 text+=" def %s(self):\n"%tname
1131 a_0=makeNumberedArray(sh0,s=1.)
1132 r_0=(a_0+f*transpose(a_0))/2.
1133 if case0 in ["taggedData", "expandedData"]:
1134 a1_0=makeNumberedArray(sh0,s=-1.)
1135 r1_0=(a1_0+f*transpose(a1_0))/2.
1136 else:
1137 a1_0=a_0
1138 r1_0=r_0
1139 text+=mkText(case0,"arg",a_0,a1_0)
1140 text+=" res=%s(arg)\n"%name
1141 if case0=="Symbol":
1142 text+=mkText("array","s",a_0,a1_0)
1143 text+=" sub=res.substitute({arg:s})\n"
1144 res="sub"
1145 text+=mkText("array","ref",r_0,r1_0)
1146 else:
1147 res="res"
1148 text+=mkText(case0,"ref",r_0,r1_0)
1149 text+=mkTypeAndShapeTest(case0,sh0,"res")
1150 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1151
1152 if case0 == "taggedData" :
1153 t_prog_with_tags+=text
1154 else:
1155 t_prog+=text
1156 print test_header
1157 print t_prog
1158 # print t_prog_with_tags
1159 print test_tail
1160 1/0
1161
1162 #=======================================================================================================
1163 # eigenvalues
1164 #=======================================================================================================
1165 import numarray.linear_algebra
1166 name="eigenvalues"
1167 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1168 for sh0 in [ (1,1), (2,2), (3,3)]:
1169 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1170 tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1171 text+=" def %s(self):\n"%tname
1172 a_0=makeArray(sh0,[-1.,1])
1173 a_0=(a_0+numarray.transpose(a_0))/2.
1174 ev=numarray.linear_algebra.eigenvalues(a_0)
1175 ev.sort()
1176 if case0 in ["taggedData", "expandedData"]:
1177 a1_0=makeArray(sh0,[-1.,1])
1178 a1_0=(a1_0+numarray.transpose(a1_0))/2.
1179 ev1=numarray.linear_algebra.eigenvalues(a1_0)
1180 ev1.sort()
1181 else:
1182 a1_0=a_0
1183 ev1=ev
1184 text+=mkText(case0,"arg",a_0,a1_0)
1185 text+=" res=%s(arg)\n"%name
1186 if case0=="Symbol":
1187 text+=mkText("array","s",a_0,a1_0)
1188 text+=" sub=res.substitute({arg:s})\n"
1189 res="sub"
1190 text+=mkText("array","ref",ev,ev1)
1191 else:
1192 res="res"
1193 text+=mkText(case0,"ref",ev,ev1)
1194 text+=mkTypeAndShapeTest(case0,(sh0[0],),"res")
1195 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1196
1197 if case0 == "taggedData" :
1198 t_prog_with_tags+=text
1199 else:
1200 t_prog+=text
1201 print test_header
1202 # print t_prog
1203 print t_prog_with_tags
1204 print test_tail
1205 1/0
1206
1207 #=======================================================================================================
1208 # get slices
1209 #=======================================================================================================
1210 for case0 in ["constData","taggedData","expandedData","Symbol"]:
1211 for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,5,3)]:
1212 # get perm:
1213 if len(sh0)==2:
1214 check=[[1,0]]
1215 elif len(sh0)==3:
1216 check=[[1,0,2],
1217 [1,2,0],
1218 [2,1,0],
1219 [2,0,2],
1220 [0,2,1]]
1221 elif len(sh0)==4:
1222 check=[[0,1,3,2],
1223 [0,2,1,3],
1224 [0,2,3,1],
1225 [0,3,2,1],
1226 [0,3,1,2] ,
1227 [1,0,2,3],
1228 [1,0,3,2],
1229 [1,2,0,3],
1230 [1,2,3,0],
1231 [1,3,2,0],
1232 [1,3,0,2],
1233 [2,0,1,3],
1234 [2,0,3,1],
1235 [2,1,0,3],
1236 [2,1,3,0],
1237 [2,3,1,0],
1238 [2,3,0,1],
1239 [3,0,1,2],
1240 [3,0,2,1],
1241 [3,1,0,2],
1242 [3,1,2,0],
1243 [3,2,1,0],
1244 [3,2,0,1]]
1245 else:
1246 check=[]
1247
1248 # create the test cases:
1249 processed=[]
1250 l=["R","U","L","P","C","N"]
1251 c=[""]
1252 for i in range(len(sh0)):
1253 tmp=[]
1254 for ci in c:
1255 tmp+=[ci+li for li in l]
1256 c=tmp
1257 # SHUFFLE
1258 c2=[]
1259 while len(c)>0:
1260 i=int(random.random()*len(c))
1261 c2.append(c[i])
1262 del c[i]
1263 c=c2
1264 for ci in c:
1265 t=""
1266 sh=()
1267 for i in range(len(ci)):
1268 if ci[i]=="R":
1269 s="%s:%s"%(1,sh0[i]-1)
1270 sh=sh+(sh0[i]-2,)
1271 if ci[i]=="U":
1272 s=":%s"%(sh0[i]-1)
1273 sh=sh+(sh0[i]-1,)
1274 if ci[i]=="L":
1275 s="2:"
1276 sh=sh+(sh0[i]-2,)
1277 if ci[i]=="P":
1278 s="%s"%(int(sh0[i]/2))
1279 if ci[i]=="C":
1280 s=":"
1281 sh=sh+(sh0[i],)
1282 if ci[i]=="N":
1283 s=""
1284 sh=sh+(sh0[i],)
1285 if len(s)>0:
1286 if not t=="": t+=","
1287 t+=s
1288 N_found=False
1289 noN_found=False
1290 process=len(t)>0
1291 for i in ci:
1292 if i=="N":
1293 if not noN_found and N_found: process=False
1294 N_found=True
1295 else:
1296 if N_found: process=False
1297 noNfound=True
1298 # is there a similar one processed allready
1299 if process and ci.find("N")==-1:
1300 for ci2 in processed:
1301 for chi in check:
1302 is_perm=True
1303 for i in range(len(chi)):
1304 if not ci[i]==ci2[chi[i]]: is_perm=False
1305 if is_perm: process=False
1306 # if not process: print ci," rejected"
1307 if process:
1308 processed.append(ci)
1309 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1310 tname="test_getslice_%s_rank%s_%s"%(case0,len(sh0),ci)
1311 text+=" def %s(self):\n"%tname
1312 a_0=makeNumberedArray(sh0,s=1)
1313 if case0 in ["taggedData", "expandedData"]:
1314 a1_0=makeNumberedArray(sh0,s=-1.)
1315 else:
1316 a1_0=a_0
1317 r=eval("a_0[%s]"%t)
1318 r1=eval("a1_0[%s]"%t)
1319 text+=mkText(case0,"arg",a_0,a1_0)
1320 text+=" res=arg[%s]\n"%t
1321 if case0=="Symbol":
1322 text+=mkText("array","s",a_0,a1_0)
1323 text+=" sub=res.substitute({arg:s})\n"
1324 res="sub"
1325 text+=mkText("array","ref",r,r1)
1326 else:
1327 res="res"
1328 text+=mkText(case0,"ref",r,r1)
1329 text+=mkTypeAndShapeTest(case0,sh,"res")
1330 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1331
1332 if case0 == "taggedData" :
1333 t_prog_with_tags+=text
1334 else:
1335 t_prog+=text
1336
1337 print test_header
1338 # print t_prog
1339 print t_prog_with_tags
1340 print test_tail
1341 1/0
1342 #============================================================================================
1343 def innerTEST(arg0,arg1):
1344 if isinstance(arg0,float):
1345 out=numarray.array(arg0*arg1)
1346 else:
1347 out=(arg0*arg1).sum()
1348 return out
1349
1350 def outerTEST(arg0,arg1):
1351 if isinstance(arg0,float):
1352 out=numarray.array(arg0*arg1)
1353 elif isinstance(arg1,float):
1354 out=numarray.array(arg0*arg1)
1355 else:
1356 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1357 return out
1358
1359 def tensorProductTest(arg0,arg1,sh_s):
1360 if isinstance(arg0,float):
1361 out=numarray.array(arg0*arg1)
1362 elif isinstance(arg1,float):
1363 out=numarray.array(arg0*arg1)
1364 elif len(sh_s)==0:
1365 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1366 else:
1367 l=len(sh_s)
1368 sh0=arg0.shape[:arg0.rank-l]
1369 sh1=arg1.shape[l:]
1370 ls,l0,l1=1,1,1
1371 for i in sh_s: ls*=i
1372 for i in sh0: l0*=i
1373 for i in sh1: l1*=i
1374 out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
1375 out2=numarray.zeros((l0,l1),numarray.Float)
1376 for i0 in range(l0):
1377 for i1 in range(l1):
1378 for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
1379 out=out2.resize(sh0+sh1)
1380 return out
1381
1382 def testMatrixMult(arg0,arg1,sh_s):
1383 return numarray.matrixmultiply(arg0,arg1)
1384
1385
1386 def testTensorMult(arg0,arg1,sh_s):
1387 if len(arg0)==2:
1388 return numarray.matrixmultiply(arg0,arg1)
1389 else:
1390 if arg1.rank==4:
1391 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
1392 for i0 in range(arg0.shape[0]):
1393 for i1 in range(arg0.shape[1]):
1394 for i2 in range(arg0.shape[2]):
1395 for i3 in range(arg0.shape[3]):
1396 for j2 in range(arg1.shape[2]):
1397 for j3 in range(arg1.shape[3]):
1398 out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
1399 elif arg1.rank==3:
1400 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
1401 for i0 in range(arg0.shape[0]):
1402 for i1 in range(arg0.shape[1]):
1403 for i2 in range(arg0.shape[2]):
1404 for i3 in range(arg0.shape[3]):
1405 for j2 in range(arg1.shape[2]):
1406 out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
1407 elif arg1.rank==2:
1408 out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
1409 for i0 in range(arg0.shape[0]):
1410 for i1 in range(arg0.shape[1]):
1411 for i2 in range(arg0.shape[2]):
1412 for i3 in range(arg0.shape[3]):
1413 out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
1414 return out
1415
1416 def testReduce(arg0,init_val,test_expr,post_expr):
1417 out=init_val
1418 if isinstance(arg0,float):
1419 out=eval(test_expr.replace("%a1%","arg0"))
1420 elif arg0.rank==0:
1421 out=eval(test_expr.replace("%a1%","arg0"))
1422 elif arg0.rank==1:
1423 for i0 in range(arg0.shape[0]):
1424 out=eval(test_expr.replace("%a1%","arg0[i0]"))
1425 elif arg0.rank==2:
1426 for i0 in range(arg0.shape[0]):
1427 for i1 in range(arg0.shape[1]):
1428 out=eval(test_expr.replace("%a1%","arg0[i0,i1]"))
1429 elif arg0.rank==3:
1430 for i0 in range(arg0.shape[0]):
1431 for i1 in range(arg0.shape[1]):
1432 for i2 in range(arg0.shape[2]):
1433 out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2]"))
1434 elif arg0.rank==4:
1435 for i0 in range(arg0.shape[0]):
1436 for i1 in range(arg0.shape[1]):
1437 for i2 in range(arg0.shape[2]):
1438 for i3 in range(arg0.shape[3]):
1439 out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))
1440 return eval(post_expr)
1441
1442 def clipTEST(arg0,mn,mx):
1443 if isinstance(arg0,float):
1444 return max(min(arg0,mx),mn)
1445 out=numarray.zeros(arg0.shape,numarray.Float64)
1446 if arg0.rank==1:
1447 for i0 in range(arg0.shape[0]):
1448 out[i0]=max(min(arg0[i0],mx),mn)
1449 elif arg0.rank==2:
1450 for i0 in range(arg0.shape[0]):
1451 for i1 in range(arg0.shape[1]):
1452 out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
1453 elif arg0.rank==3:
1454 for i0 in range(arg0.shape[0]):
1455 for i1 in range(arg0.shape[1]):
1456 for i2 in range(arg0.shape[2]):
1457 out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
1458 elif arg0.rank==4:
1459 for i0 in range(arg0.shape[0]):
1460 for i1 in range(arg0.shape[1]):
1461 for i2 in range(arg0.shape[2]):
1462 for i3 in range(arg0.shape[3]):
1463 out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
1464 return out
1465 def minimumTEST(arg0,arg1):
1466 if isinstance(arg0,float):
1467 if isinstance(arg1,float):
1468 if arg0>arg1:
1469 return arg1
1470 else:
1471 return arg0
1472 else:
1473 arg0=numarray.ones(arg1.shape)*arg0
1474 else:
1475 if isinstance(arg1,float):
1476 arg1=numarray.ones(arg0.shape)*arg1
1477 out=numarray.zeros(arg0.shape,numarray.Float64)
1478 if arg0.rank==0:
1479 if arg0>arg1:
1480 out=arg1
1481 else:
1482 out=arg0
1483 elif arg0.rank==1:
1484 for i0 in range(arg0.shape[0]):
1485 if arg0[i0]>arg1[i0]:
1486 out[i0]=arg1[i0]
1487 else:
1488 out[i0]=arg0[i0]
1489 elif arg0.rank==2:
1490 for i0 in range(arg0.shape[0]):
1491 for i1 in range(arg0.shape[1]):
1492 if arg0[i0,i1]>arg1[i0,i1]:
1493 out[i0,i1]=arg1[i0,i1]
1494 else:
1495 out[i0,i1]=arg0[i0,i1]
1496 elif arg0.rank==3:
1497 for i0 in range(arg0.shape[0]):
1498 for i1 in range(arg0.shape[1]):
1499 for i2 in range(arg0.shape[2]):
1500 if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
1501 out[i0,i1,i2]=arg1[i0,i1,i2]
1502 else:
1503 out[i0,i1,i2]=arg0[i0,i1,i2]
1504 elif arg0.rank==4:
1505 for i0 in range(arg0.shape[0]):
1506 for i1 in range(arg0.shape[1]):
1507 for i2 in range(arg0.shape[2]):
1508 for i3 in range(arg0.shape[3]):
1509 if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
1510 out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
1511 else:
1512 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
1513 return out
1514
1515 def unrollLoops(a,b,o,arg,tap="",x="x"):
1516 out=""
1517 if a.rank==1:
1518 z=""
1519 for i99 in range(a.shape[0]):
1520 if not z=="": z+="+"
1521 if o=="1":
1522 z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
1523 else:
1524 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
1525
1526 out+=tap+"%s=%s\n"%(arg,z)
1527
1528 elif a.rank==2:
1529 for i0 in range(a.shape[0]):
1530 z=""
1531 for i99 in range(a.shape[1]):
1532 if not z=="": z+="+"
1533 if o=="1":
1534 z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
1535 else:
1536 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
1537
1538 out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
1539 elif a.rank==3:
1540 for i0 in range(a.shape[0]):
1541 for i1 in range(a.shape[1]):
1542 z=""
1543 for i99 in range(a.shape[2]):
1544 if not z=="": z+="+"
1545 if o=="1":
1546 z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
1547 else:
1548 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
1549
1550 out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
1551 elif a.rank==4:
1552 for i0 in range(a.shape[0]):
1553 for i1 in range(a.shape[1]):
1554 for i2 in range(a.shape[2]):
1555 z=""
1556 for i99 in range(a.shape[3]):
1557 if not z=="": z+="+"
1558 if o=="1":
1559 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
1560 else:
1561 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
1562
1563 out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
1564 elif a.rank==5:
1565 for i0 in range(a.shape[0]):
1566 for i1 in range(a.shape[1]):
1567 for i2 in range(a.shape[2]):
1568 for i3 in range(a.shape[3]):
1569 z=""
1570 for i99 in range(a.shape[4]):
1571 if not z=="": z+="+"
1572 if o=="1":
1573 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
1574 else:
1575 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
1576
1577 out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
1578 return out
1579
1580 def unrollLoopsOfGrad(a,b,o,arg,tap=""):
1581 out=""
1582 if a.rank==1:
1583 for i99 in range(a.shape[0]):
1584 if o=="1":
1585 out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
1586 else:
1587 out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
1588
1589 elif a.rank==2:
1590 for i0 in range(a.shape[0]):
1591 for i99 in range(a.shape[1]):
1592 if o=="1":
1593 out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
1594 else:
1595 out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
1596
1597 elif a.rank==3:
1598 for i0 in range(a.shape[0]):
1599 for i1 in range(a.shape[1]):
1600 for i99 in range(a.shape[2]):
1601 if o=="1":
1602 out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
1603 else:
1604 out+=tap+"%s[%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99],i99,b[i0,i1,i99])
1605
1606 elif a.rank==4:
1607 for i0 in range(a.shape[0]):
1608 for i1 in range(a.shape[1]):
1609 for i2 in range(a.shape[2]):
1610 for i99 in range(a.shape[3]):
1611 if o=="1":
1612 out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
1613 else:
1614 out+=tap+"%s[%s,%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99],i99,b[i0,i1,i2,i99])
1615 return out
1616 def unrollLoopsOfDiv(a,b,o,arg,tap=""):
1617 out=tap+arg+"="
1618 if o=="1":
1619 z=0.
1620 for i99 in range(a.shape[0]):
1621 z+=b[i99,i99]+a[i99,i99]
1622 out+="(%s)"%z
1623 else:
1624 z=0.
1625 for i99 in range(a.shape[0]):
1626 z+=b[i99,i99]
1627 if i99>0: out+="+"
1628 out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
1629 out+="+(%s)"%z
1630 return out
1631
1632 def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
1633 if where=="Function":
1634 xfac_o=1.
1635 xfac_op=0.
1636 z_fac=1./2.
1637 z_fac_s=""
1638 zo_fac_s="/(o+1.)"
1639 elif where=="FunctionOnBoundary":
1640 xfac_o=1.
1641 xfac_op=0.
1642 z_fac=1.
1643 z_fac_s="*dim"
1644 zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
1645 elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1646 xfac_o=0.
1647 xfac_op=1.
1648 z_fac=1./2.
1649 z_fac_s=""
1650 zo_fac_s="/(o+1.)"
1651 out=""
1652 if a.rank==1:
1653 zo=0.
1654 zop=0.
1655 z=0.
1656 for i99 in range(a.shape[0]):
1657 if i99==0:
1658 zo+= xfac_o*a[i99]
1659 zop+= xfac_op*a[i99]
1660 else:
1661 zo+=a[i99]
1662 z+=b[i99]
1663
1664 out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
1665 if zop==0.:
1666 out+="\n"
1667 else:
1668 out+="+(%s)*0.5**o\n"%zop
1669 elif a.rank==2:
1670 for i0 in range(a.shape[0]):
1671 zo=0.
1672 zop=0.
1673 z=0.
1674 for i99 in range(a.shape[1]):
1675 if i99==0:
1676 zo+= xfac_o*a[i0,i99]
1677 zop+= xfac_op*a[i0,i99]
1678 else:
1679 zo+=a[i0,i99]
1680 z+=b[i0,i99]
1681
1682 out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
1683 if zop==0.:
1684 out+="\n"
1685 else:
1686 out+="+(%s)*0.5**o\n"%zop
1687 elif a.rank==3:
1688 for i0 in range(a.shape[0]):
1689 for i1 in range(a.shape[1]):
1690 zo=0.
1691 zop=0.
1692 z=0.
1693 for i99 in range(a.shape[2]):
1694 if i99==0:
1695 zo+= xfac_o*a[i0,i1,i99]
1696 zop+= xfac_op*a[i0,i1,i99]
1697 else:
1698 zo+=a[i0,i1,i99]
1699 z+=b[i0,i1,i99]
1700
1701 out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
1702 if zop==0.:
1703 out+="\n"
1704 else:
1705 out+="+(%s)*0.5**o\n"%zop
1706 elif a.rank==4:
1707 for i0 in range(a.shape[0]):
1708 for i1 in range(a.shape[1]):
1709 for i2 in range(a.shape[2]):
1710 zo=0.
1711 zop=0.
1712 z=0.
1713 for i99 in range(a.shape[3]):
1714 if i99==0:
1715 zo+= xfac_o*a[i0,i1,i2,i99]
1716 zop+= xfac_op*a[i0,i1,i2,i99]
1717
1718 else:
1719 zo+=a[i0,i1,i2,i99]
1720 z+=b[i0,i1,i2,i99]
1721
1722 out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
1723 if zop==0.:
1724 out+="\n"
1725 else:
1726 out+="+(%s)*0.5**o\n"%zop
1727
1728 elif a.rank==5:
1729 for i0 in range(a.shape[0]):
1730 for i1 in range(a.shape[1]):
1731 for i2 in range(a.shape[2]):
1732 for i3 in range(a.shape[3]):
1733 zo=0.
1734 zop=0.
1735 z=0.
1736 for i99 in range(a.shape[4]):
1737 if i99==0:
1738 zo+= xfac_o*a[i0,i1,i2,i3,i99]
1739 zop+= xfac_op*a[i0,i1,i2,i3,i99]
1740
1741 else:
1742 zo+=a[i0,i1,i2,i3,i99]
1743 z+=b[i0,i1,i2,i3,i99]
1744 out+=tap+"%s[%s,%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,i3,zo,zo_fac_s,z*z_fac,z_fac_s)
1745 if zop==0.:
1746 out+="\n"
1747 else:
1748 out+="+(%s)*0.5**o\n"%zop
1749
1750 return out
1751 def unrollLoopsSimplified(b,arg,tap=""):
1752 out=""
1753 if isinstance(b,float) or b.rank==0:
1754 out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
1755
1756 elif b.rank==1:
1757 for i0 in range(b.shape[0]):
1758 out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
1759 elif b.rank==2:
1760 for i0 in range(b.shape[0]):
1761 for i1 in range(b.shape[1]):
1762 out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
1763 elif b.rank==3:
1764 for i0 in range(b.shape[0]):
1765 for i1 in range(b.shape[1]):
1766 for i2 in range(b.shape[2]):
1767 out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
1768 elif b.rank==4:
1769 for i0 in range(b.shape[0]):
1770 for i1 in range(b.shape[1]):
1771 for i2 in range(b.shape[2]):
1772 for i3 in range(b.shape[3]):
1773 out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
1774 return out
1775
1776 def unrollLoopsOfL2(b,where,arg,tap=""):
1777 out=""
1778 z=[]
1779 if isinstance(b,float) or b.rank==0:
1780 z.append(b**2)
1781 elif b.rank==1:
1782 for i0 in range(b.shape[0]):
1783 z.append(b[i0]**2)
1784 elif b.rank==2:
1785 for i1 in range(b.shape[1]):
1786 s=0
1787 for i0 in range(b.shape[0]):
1788 s+=b[i0,i1]**2
1789 z.append(s)
1790 elif b.rank==3:
1791 for i2 in range(b.shape[2]):
1792 s=0
1793 for i0 in range(b.shape[0]):
1794 for i1 in range(b.shape[1]):
1795 s+=b[i0,i1,i2]**2
1796 z.append(s)
1797
1798 elif b.rank==4:
1799 for i3 in range(b.shape[3]):
1800 s=0
1801 for i0 in range(b.shape[0]):
1802 for i1 in range(b.shape[1]):
1803 for i2 in range(b.shape[2]):
1804 s+=b[i0,i1,i2,i3]**2
1805 z.append(s)
1806 if where=="Function":
1807 xfac_o=1.
1808 xfac_op=0.
1809 z_fac_s=""
1810 zo_fac_s=""
1811 zo_fac=1./3.
1812 elif where=="FunctionOnBoundary":
1813 xfac_o=1.
1814 xfac_op=0.
1815 z_fac_s="*dim"
1816 zo_fac_s="*(2.*dim+1.)/3."
1817 zo_fac=1.
1818 elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1819 xfac_o=0.
1820 xfac_op=1.
1821 z_fac_s=""
1822 zo_fac_s=""
1823 zo_fac=1./3.
1824 zo=0.
1825 zop=0.
1826 for i99 in range(len(z)):
1827 if i99==0:
1828 zo+=xfac_o*z[i99]
1829 zop+=xfac_op*z[i99]
1830 else:
1831 zo+=z[i99]
1832 out+=tap+"%s=sqrt((%s)%s"%(arg,zo*zo_fac,zo_fac_s)
1833 if zop==0.:
1834 out+=")\n"
1835 else:
1836 out+="+(%s))\n"%(zop*0.5**2)
1837 return out
1838 #=======================================================================================================
1839 # transpose
1840 #=======================================================================================================
1841 def transposeTest(r,offset):
1842 if isinstance(r,float): return r
1843 s=r.shape
1844 s1=1
1845 for i in s[:offset]: s1*=i
1846 s2=1
1847 for i in s[offset:]: s2*=i
1848 out=numarray.reshape(r,(s1,s2))
1849 out.transpose()
1850 return numarray.resize(out,s[offset:]+s[:offset])
1851
1852 name,tt="transpose",transposeTest
1853 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1854 for sh0 in [ (), (3,), (4,5), (6,2,2),(3,2,3,4)]:
1855 for offset in range(len(sh0)+1):
1856 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1857 tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1858 text+=" def %s(self):\n"%tname
1859 sh_t=sh0[offset:]+sh0[:offset]
1860
1861 # sh_t=list(sh0)
1862 # sh_t[offset+1]=sh_t[offset]
1863 # sh_t=tuple(sh_t)
1864 # sh_r=[]
1865 # for i in range(offset): sh_r.append(sh0[i])
1866 # for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1867 # sh_r=tuple(sh_r)
1868
1869 a_0=makeArray(sh0,[-1.,1])
1870 if case0 in ["taggedData", "expandedData"]:
1871 a1_0=makeArray(sh0,[-1.,1])
1872 else:
1873 a1_0=a_0
1874 r=tt(a_0,offset)
1875 r1=tt(a1_0,offset)
1876 text+=mkText(case0,"arg",a_0,a1_0)
1877 text+=" res=%s(arg,%s)\n"%(name,offset)
1878 if case0=="Symbol":
1879 text+=mkText("array","s",a_0,a1_0)
1880 text+=" sub=res.substitute({arg:s})\n"
1881 res="sub"
1882 text+=mkText("array","ref",r,r1)
1883 else:
1884 res="res"
1885 text+=mkText(case0,"ref",r,r1)
1886 text+=mkTypeAndShapeTest(case0,sh_t,"res")
1887 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1888
1889 if case0 == "taggedData" :
1890 t_prog_with_tags+=text
1891 else:
1892 t_prog+=text
1893
1894 print test_header
1895 # print t_prog
1896 print t_prog_with_tags
1897 print test_tail
1898 1/0
1899 #=======================================================================================================
1900 # L2
1901 #=======================================================================================================
1902 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1903 for data in ["Data","Symbol"]:
1904 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1905 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1906 tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
1907 text+=" def %s(self):\n"%tname
1908 text+=" \"\"\"\n"
1909 text+=" tests L2-norm of %s on the %s\n\n"%(data,where)
1910 text+=" assumptions: self.domain supports integration on %s\n"%where
1911 text+=" \"\"\"\n"
1912 text+=" dim=self.domain.getDim()\n"
1913 text+=" w=%s(self.domain)\n"%where
1914 text+=" x=w.getX()\n"
1915 o="1"
1916 if len(sh)>0:
1917 sh_2=sh[:len(sh)-1]+(2,)
1918 sh_3=sh[:len(sh)-1]+(3,)
1919 b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
1920 b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
1921 else:
1922 sh_2=()
1923 sh_3=()
1924 b_2=makeArray(sh,[-1.,1])
1925 b_3=makeArray(sh,[-1.,1])
1926
1927 if data=="Symbol":
1928 val="s"
1929 res="sub"
1930 else:
1931 val="arg"
1932 res="res"
1933 text+=" if dim==2:\n"
1934 if data=="Symbol":
1935 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
1936
1937 text+=" %s=Data(0,%s,w)\n"%(val,sh_2)
1938 text+=unrollLoopsSimplified(b_2,val,tap=" ")
1939 text+=unrollLoopsOfL2(b_2,where,"ref",tap=" ")
1940 text+="\n else:\n"
1941 if data=="Symbol":
1942 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
1943 text+=" %s=Data(0,%s,w)\n"%(val,sh_3)
1944 text+=unrollLoopsSimplified(b_3,val,tap=" ")
1945 text+=unrollLoopsOfL2(b_3,where,"ref",tap=" ")
1946 text+="\n res=L2(arg)\n"
1947 if data=="Symbol":
1948 text+=" sub=res.substitute({arg:s})\n"
1949 text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1950 text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1951 else:
1952 text+=" self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1953 text+=" self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
1954 t_prog+=text
1955 print t_prog
1956 1/0
1957
1958 #=======================================================================================================
1959 # div
1960 #=======================================================================================================
1961 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1962 for data in ["Data","Symbol"]:
1963 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1964 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1965 tname="test_div_on%s_from%s_%s"%(where,data,case)
1966 text+=" def %s(self):\n"%tname
1967 text+=" \"\"\"\n"
1968 text+=" tests divergence of %s on the %s\n\n"%(data,where)
1969 text+=" assumptions: %s(self.domain) exists\n"%case
1970 text+=" self.domain supports div on %s\n"%where
1971 text+=" \"\"\"\n"
1972 if case=="ReducedSolution":
1973 text+=" o=1\n"
1974 o="1"
1975 else:
1976 text+=" o=self.order\n"
1977 o="o"
1978 text+=" dim=self.domain.getDim()\n"
1979 text+=" w_ref=%s(self.domain)\n"%where
1980 text+=" x_ref=w_ref.getX()\n"
1981 text+=" w=%s(self.domain)\n"%case
1982 text+=" x=w.getX()\n"
1983 a_2=makeArray((2,2),[-1.,1])
1984 b_2=makeArray((2,2),[-1.,1])
1985 a_3=makeArray((3,3),[-1.,1])
1986 b_3=makeArray((3,3),[-1.,1])
1987 if data=="Symbol":
1988 text+=" arg=Symbol(shape=(dim,),dim=dim)\n"
1989 val="s"
1990 res="sub"
1991 else:
1992 val="arg"
1993 res="res"
1994 text+=" %s=Vector(0,w)\n"%val
1995 text+=" if dim==2:\n"
1996 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1997 text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap=" ")
1998 text+="\n else:\n"
1999
2000 text+=unrollLoops(a_3,b_3,o,val,tap=" ")
2001 text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap=" ")
2002 text+="\n res=div(arg,where=w_ref)\n"
2003 if data=="Symbol":
2004 text+=" sub=res.substitute({arg:s})\n"
2005 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
2006 text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
2007 text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
2008 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
2009
2010
2011 t_prog+=text
2012 print t_prog
2013 1/0
2014
2015 #=======================================================================================================
2016 # interpolation
2017 #=======================================================================================================
2018 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
2019 for data in ["Data","Symbol"]:
2020 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2021 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
2022 if where==case or \
2023 ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
2024 ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
2025 (case=="ContinuousFunction" and where in ["Solution","ReducedSolution"]) or \
2026 (case=="Solution" and where=="ReducedSolution") :
2027
2028
2029 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2030 tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
2031 text+=" def %s(self):\n"%tname
2032 text+=" \"\"\"\n"
2033 text+=" tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
2034 text+=" assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
2035 text+=" \"\"\"\n"
2036 if case=="ReducedSolution" or where=="ReducedSolution":
2037 text+=" o=1\n"
2038 o="1"
2039 else:
2040 text+=" o=self.order\n"
2041 o="o"
2042 text+=" dim=self.domain.getDim()\n"
2043 text+=" w_ref=%s(self.domain)\n"%where
2044 text+=" x_ref=w_ref.getX()\n"
2045 text+=" w=%s(self.domain)\n"%case
2046 text+=" x=w.getX()\n"
2047 a_2=makeArray(sh+(2,),[-1.,1])
2048 b_2=makeArray(sh+(2,),[-1.,1])
2049 a_3=makeArray(sh+(3,),[-1.,1])
2050 b_3=makeArray(sh+(3,),[-1.,1])
2051 if data=="Symbol":
2052 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
2053 val="s"
2054 res="sub"
2055 else:
2056 val="arg"
2057 res="res"
2058 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
2059 text+=" ref=Data(0,%s,w_ref)\n"%str(sh)
2060 text+=" if dim==2:\n"
2061 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
2062 text+=unrollLoops(a_2,b_2,o,"ref",tap=" ",x="x_ref")
2063 text+=" else:\n"
2064
2065 text+=unrollLoops(a_3,b_3,o,val,tap=" ")
2066 text+=unrollLoops(a_3,b_3,o,"ref",tap=" ",x="x_ref")
2067 text+=" res=interpolate(arg,where=w_ref)\n"
2068 if data=="Symbol":
2069 text+=" sub=res.substitute({arg:s})\n"
2070 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
2071 text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
2072 text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
2073 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
2074 t_prog+=text
2075 print test_header
2076 print t_prog
2077 print test_tail
2078 1/0
2079
2080 #=======================================================================================================
2081 # grad
2082 #=======================================================================================================
2083 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2084 for data in ["Data","Symbol"]:
2085 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
2086 for sh in [ (),(2,), (4,5), (6,2,2)]:
2087 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2088 tname="test_grad_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
2089 text+=" def %s(self):\n"%tname
2090 text+=" \"\"\"\n"
2091 text+=" tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
2092 text+=" assumptions: %s(self.domain) exists\n"%case
2093 text+=" self.domain supports gardient on %s\n"%where
2094 text+=" \"\"\"\n"
2095 if case=="ReducedSolution":
2096 text+=" o=1\n"
2097 o="1"
2098 else:
2099 text+=" o=self.order\n"
2100 o="o"
2101 text+=" dim=self.domain.getDim()\n"
2102 text+=" w_ref=%s(self.domain)\n"%where
2103 text+=" x_ref=w_ref.getX()\n"
2104 text+=" w=%s(self.domain)\n"%case
2105 text+=" x=w.getX()\n"
2106 a_2=makeArray(sh+(2,),[-1.,1])
2107 b_2=makeArray(sh+(2,),[-1.,1])
2108 a_3=makeArray(sh+(3,),[-1.,1])
2109 b_3=makeArray(sh+(3,),[-1.,1])
2110 if data=="Symbol":
2111 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
2112 val="s"
2113 res="sub"
2114 else:
2115 val="arg"
2116 res="res"
2117 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
2118 text+=" ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
2119 text+=" if dim==2:\n"
2120 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
2121 text+=unrollLoopsOfGrad(a_2,b_2,o,"ref",tap=" ")
2122 text+=" else:\n"
2123
2124 text+=unrollLoops(a_3,b_3,o,val,tap=" ")
2125 text+=unrollLoopsOfGrad(a_3,b_3,o,"ref",tap=" ")
2126 text+=" res=grad(arg,where=w_ref)\n"
2127 if data=="Symbol":
2128 text+=" sub=res.substitute({arg:s})\n"
2129 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
2130 text+=" self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
2131 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
2132
2133
2134 t_prog+=text
2135 print test_header
2136 print t_prog
2137 print test_tail
2138 1/0
2139
2140
2141 #=======================================================================================================
2142 # integrate
2143 #=======================================================================================================
2144 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2145 for data in ["Data","Symbol"]:
2146 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2147 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
2148 if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:
2149 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2150 tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
2151 text+=" def %s(self):\n"%tname
2152 text+=" \"\"\"\n"
2153 text+=" tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
2154 text+=" assumptions: %s(self.domain) exists\n"%case
2155 text+=" self.domain supports integral on %s\n"%where
2156
2157 text+=" \"\"\"\n"
2158 if case=="ReducedSolution":
2159 text+=" o=1\n"
2160 o="1"
2161 else:
2162 text+=" o=self.order\n"
2163 o="o"
2164 text+=" dim