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

Contents of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


Revision 804 - (show annotations)
Thu Aug 10 01:12:16 2006 UTC (13 years, 3 months ago) by gross
File size: 168742 byte(s)
the new function swap_axes + tests added. (It replaces swap).


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
585 #=======================================================================================================
586 # swap axes
587 #=======================================================================================================
588 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
589 for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
590 for axis0 in range(len(sh0)):
591 for axis1 in range(axis0+1,len(sh0)):
592 tname="test_%s_%s_rank%s_axes_%s%s"%("swapaxes",case0,len(sh0),axis0,axis1)
593 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
594 text+=" def %s(self):\n"%tname
595 a_0=makeArray(sh0,[-8,8])
596 if case0 in ["taggedData", "expandedData"]:
597 a1_0=makeArray(sh0,[-8,8])
598 else:
599 a1_0=a_0
600 r=numarray.swapaxes(a_0,axis0,axis1)
601 r1=numarray.swapaxes(a1_0,axis0,axis1)
602 text+=mkText(case0,"arg",a_0,a1_0)
603 text+=" res=swap_axes(arg,axis0=%s,axis1=%s)\n"%(axis0,axis1)
604 if case0=="Symbol":
605 text+=mkText("array","s",a_0,a1_0)
606 text+=" sub=res.substitute({arg:s})\n"
607 res="sub"
608 text+=mkText("array","ref",r,r1)
609 else:
610 res="res"
611 text+=mkText(case0,"ref",r,r1)
612 text+=mkTypeAndShapeTest(case0,r.shape,"res")
613 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
614
615 if case0 == "taggedData":
616 t_prog_with_tags+=text
617 else:
618 t_prog+=text
619 print test_header
620 # print t_prog
621 print t_prog_with_tags
622 print test_tail
623 1/0
624 while True:
625 for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
626 if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
627 and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
628 arg_shape=oper[2](sh0,sh1,sh_s)
629 if not arg_shape==None:
630 case=getResultCaseForBin(case0,case1)
631 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
632 if oper[0] in [ "generalTensorProduct", "generalTransposedTensorProduct", "generalTensorTransposedProduct"]:
633 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))
634
635 res_text=" res=%s(arg0,arg1,axis_offset=%s)\n"%(oper[0],len(sh_s))
636 else:
637 tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0+sh_s),case1,len(sh_s+sh1))
638
639 res_text=" res=%s(arg0,arg1)\n"%oper[0]
640
641 a_0=makeArray(arg_shape[0],[-8,8])
642 if case0 in ["taggedData", "expandedData"]:
643 a1_0=makeArray(arg_shape[0],[-8,8])
644 else:
645 a1_0=a_0
646
647 a_1=makeArray(arg_shape[1],[-8,8])
648 if case1 in ["taggedData", "expandedData"]:
649 a1_1=makeArray(arg_shape[1],[-8,8])
650 else:
651 a1_1=a_1
652 r=oper[1](a_0,a_1,sh_s)
653 r1=oper[1](a1_0,a1_1,sh_s)
654 text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
655 text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
656 text+=res_text
657 if case=="Symbol":
658 c0_res,c1_res=case0,case1
659 subs="{"
660 if case0=="Symbol":
661 text+=mkText("array","s0",a_0,a1_0)
662 subs+="arg0:s0"
663 c0_res="array"
664 if case1=="Symbol":
665 text+=mkText("array","s1",a_1,a1_1)
666 if not subs.endswith("{"): subs+=","
667 subs+="arg1:s1"
668 c1_res="array"
669 subs+="}"
670 text+=" sub=res.substitute(%s)\n"%subs
671 res="sub"
672 text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
673 else:
674 res="res"
675 text+=mkText(case,"ref",r,r1)
676 text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
677 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
678 1/0
679
680 #=================================================
681 def testMatrixMult(arg0,arg1,sh_s):
682 return numarray.matrixmultiply(arg0,arg1)
683
684 def shapeTestMatrixMult(sh0,sh1,sh_s):
685 if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2):
686 return sh0+sh_s,sh_s+sh1
687
688 def testTensorMult(arg0,arg1,sh_s):
689 if arg0.rank==2:
690 return numarray.matrixmultiply(arg0,arg1)
691 else:
692 if arg1.rank==4:
693 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
694 for i0 in range(arg0.shape[0]):
695 for i1 in range(arg0.shape[1]):
696 for i2 in range(arg0.shape[2]):
697 for i3 in range(arg0.shape[3]):
698 for j2 in range(arg1.shape[2]):
699 for j3 in range(arg1.shape[3]):
700 out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
701 elif arg1.rank==3:
702 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
703 for i0 in range(arg0.shape[0]):
704 for i1 in range(arg0.shape[1]):
705 for i2 in range(arg0.shape[2]):
706 for i3 in range(arg0.shape[3]):
707 for j2 in range(arg1.shape[2]):
708 out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
709 elif arg1.rank==2:
710 out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
711 for i0 in range(arg0.shape[0]):
712 for i1 in range(arg0.shape[1]):
713 for i2 in range(arg0.shape[2]):
714 for i3 in range(arg0.shape[3]):
715 out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
716 return out
717
718 def shapeTestTensorMult(sh0,sh1,sh_s):
719 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)):
720 return sh0+sh_s,sh_s+sh1
721
722 def tensorProductTest(arg0,arg1,sh_s):
723 if isinstance(arg0,float):
724 out=numarray.array(arg0*arg1)
725 elif isinstance(arg1,float):
726 out=numarray.array(arg0*arg1)
727 elif len(sh_s)==0:
728 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
729 else:
730 l=len(sh_s)
731 sh0=arg0.shape[:arg0.rank-l]
732 sh1=arg1.shape[l:]
733 ls,l0,l1=1,1,1
734 for i in sh_s: ls*=i
735 for i in sh0: l0*=i
736 for i in sh1: l1*=i
737 out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
738 out2=numarray.zeros((l0,l1),numarray.Float)
739 for i0 in range(l0):
740 for i1 in range(l1):
741 for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
742 out=out2.resize(sh0+sh1)
743 return out
744 def tensorProductShape(sh0,sh1,sh_s):
745 return sh0+sh_s,sh_s+sh1
746
747 def transpose(r,offset):
748 if isinstance(r,float): return r
749 s=r.shape
750 s1=1
751 for i in s[:offset]: s1*=i
752 s2=1
753 for i in s[offset:]: s2*=i
754 out=numarray.reshape(r,(s1,s2))
755 out.transpose()
756 return numarray.resize(out,s[offset:]+s[:offset])
757
758 def shapeTestTMatrixMult(sh0,sh1,sh_s):
759 if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2):
760 return sh_s+sh0,sh_s+sh1
761
762 def testTMatrixMult(arg0,arg1,sh_s):
763 return numarray.matrixmultiply(transpose(arg0,1),arg1)
764
765 def shapeTestTTensorMult(sh0,sh1,sh_s):
766 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)):
767 return sh_s+sh0,sh_s+sh1
768
769 def testTTensorMult(arg0,arg1,sh_s):
770 if arg0.rank==2:
771 return numarray.matrixmultiply(transpose(arg0,1),arg1)
772 else:
773 if arg1.rank==4:
774 out=numarray.zeros((arg0.shape[2],arg0.shape[3],arg1.shape[2],arg1.shape[3]),numarray.Float)
775 for i0 in range(arg0.shape[2]):
776 for i1 in range(arg0.shape[3]):
777 for i2 in range(arg0.shape[0]):
778 for i3 in range(arg0.shape[1]):
779 for j2 in range(arg1.shape[2]):
780 for j3 in range(arg1.shape[3]):
781 out[i0,i1,j2,j3]+=arg0[i2,i3,i0,i1]*arg1[i2,i3,j2,j3]
782 elif arg1.rank==3:
783 out=numarray.zeros((arg0.shape[2],arg0.shape[3],arg1.shape[2]),numarray.Float)
784 for i0 in range(arg0.shape[2]):
785 for i1 in range(arg0.shape[3]):
786 for i2 in range(arg0.shape[0]):
787 for i3 in range(arg0.shape[1]):
788 for j2 in range(arg1.shape[2]):
789 out[i0,i1,j2]+=arg0[i2,i3,i0,i1]*arg1[i2,i3,j2]
790 elif arg1.rank==2:
791 out=numarray.zeros((arg0.shape[2],arg0.shape[3]),numarray.Float)
792 for i0 in range(arg0.shape[2]):
793 for i1 in range(arg0.shape[3]):
794 for i2 in range(arg0.shape[0]):
795 for i3 in range(arg0.shape[1]):
796 out[i0,i1]+=arg0[i2,i3,i0,i1]*arg1[i2,i3]
797 return out
798
799 def tensorTProductShape(sh0,sh1,sh_s):
800 return sh_s+sh0,sh_s+sh1
801
802 def tensorTProductTest(arg0,arg1,sh_s):
803 if isinstance(arg0,float):
804 out=numarray.array(arg0*arg1)
805 elif isinstance(arg1,float):
806 out=numarray.array(arg0*arg1)
807 elif len(sh_s)==0:
808 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
809 else:
810 l=len(sh_s)
811 sh0=arg0.shape[l:]
812 sh1=arg1.shape[l:]
813 ls,l0,l1=1,1,1
814 for i in sh_s: ls*=i
815 for i in sh0: l0*=i
816 for i in sh1: l1*=i
817 out1=numarray.outerproduct(arg0,arg1).resize((ls,l0,ls,l1))
818 out2=numarray.zeros((l0,l1),numarray.Float)
819 for i0 in range(l0):
820 for i1 in range(l1):
821 for i in range(ls): out2[i0,i1]+=out1[i,i0,i,i1]
822 out=out2.resize(sh0+sh1)
823 return out
824
825 def shapeTestMatrixTMult(sh0,sh1,sh_s):
826 if len(sh_s)==1 and len(sh0+sh_s)==2 and len(sh_s+sh1)==2:
827 return sh0+sh_s,sh1+sh_s
828
829 def testMatrixTMult(arg0,arg1,sh_s):
830 return numarray.matrixmultiply(arg0,transpose(arg1,1))
831
832 def shapeTestTensorTMult(sh0,sh1,sh_s):
833 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)):
834 return sh0+sh_s,sh1+sh_s
835
836 def testTensorTMult(arg0,arg1,sh_s):
837 if arg0.rank==2:
838 return numarray.matrixmultiply(arg0,transpose(arg1,1))
839 else:
840 if arg1.rank==4:
841 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[0],arg1.shape[1]),numarray.Float)
842 for i0 in range(arg0.shape[0]):
843 for i1 in range(arg0.shape[1]):
844 for i2 in range(arg0.shape[2]):
845 for i3 in range(arg0.shape[3]):
846 for j0 in range(arg1.shape[0]):
847 for j1 in range(arg1.shape[1]):
848 out[i0,i1,j0,j1]+=arg0[i0,i1,i2,i3]*arg1[j0,j1,i2,i3]
849 elif arg1.rank==3:
850 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[0]),numarray.Float)
851 for i0 in range(arg0.shape[0]):
852 for i1 in range(arg0.shape[1]):
853 for i2 in range(arg0.shape[2]):
854 for i3 in range(arg0.shape[3]):
855 for j0 in range(arg1.shape[0]):
856 out[i0,i1,j0]+=arg0[i0,i1,i2,i3]*arg1[j0,i2,i3]
857 elif arg1.rank==2:
858 out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
859 for i0 in range(arg0.shape[0]):
860 for i1 in range(arg0.shape[1]):
861 for i2 in range(arg0.shape[2]):
862 for i3 in range(arg0.shape[3]):
863 out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
864 return out
865
866 def tensorProductTTest(arg0,arg1,sh_s):
867 if isinstance(arg0,float):
868 out=numarray.array(arg0*arg1)
869 elif isinstance(arg1,float):
870 out=numarray.array(arg0*arg1)
871 elif len(sh_s)==0:
872 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
873 else:
874 l=len(sh_s)
875 sh0=arg0.shape[:arg0.rank-l]
876 sh1=arg1.shape[:arg1.rank-l]
877 ls,l0,l1=1,1,1
878 for i in sh_s: ls*=i
879 for i in sh0: l0*=i
880 for i in sh1: l1*=i
881 out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,l1,ls))
882 out2=numarray.zeros((l0,l1),numarray.Float)
883 for i0 in range(l0):
884 for i1 in range(l1):
885 for i in range(ls): out2[i0,i1]+=out1[i0,i,i1,i]
886 out=out2.resize(sh0+sh1)
887 return out
888
889 def tensorProductTShape(sh0,sh1,sh_s):
890 return sh0+sh_s,sh1+sh_s
891
892 #=======================================================================================================
893 # tensor multiply
894 #=======================================================================================================
895 # oper=
896 # oper=
897 for oper in [ ["tensor_mult",testTensorMult,shapeTestTensorMult], \
898 ["generalTensorProduct",tensorProductTest,tensorProductShape], \
899 ["matrix_mult",testMatrixMult,shapeTestMatrixMult], \
900 ["transposed_matrix_mult",testTMatrixMult,shapeTestTMatrixMult], \
901 ["transposed_tensor_mult",testTTensorMult,shapeTestTTensorMult], \
902 ["generalTransposedTensorProduct",tensorTProductTest,tensorTProductShape], \
903 ["matrix_transposed_mult",testMatrixTMult,shapeTestMatrixTMult], \
904 ["tensor_transposed_mult",testTensorTMult,shapeTestTensorTMult], \
905 ["generalTensorTransposedProduct",tensorProductTTest,tensorProductTShape] ]:
906
907 for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
908 for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
909 for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
910 for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
911 for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
912 if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
913 and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
914 arg_shape=oper[2](sh0,sh1,sh_s)
915 if not arg_shape==None:
916 case=getResultCaseForBin(case0,case1)
917 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
918 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
919 if oper[0] in [ "generalTensorProduct", "generalTransposedTensorProduct", "generalTensorTransposedProduct"]:
920 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))
921
922 res_text=" res=%s(arg0,arg1,axis_offset=%s)\n"%(oper[0],len(sh_s))
923 else:
924 tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0+sh_s),case1,len(sh_s+sh1))
925
926 res_text=" res=%s(arg0,arg1)\n"%oper[0]
927
928 text+=" def %s(self):\n"%tname
929 a_0=makeArray(arg_shape[0],[-8,8])
930 if case0 in ["taggedData", "expandedData"]:
931 a1_0=makeArray(arg_shape[0],[-8,8])
932 else:
933 a1_0=a_0
934
935 a_1=makeArray(arg_shape[1],[-8,8])
936 if case1 in ["taggedData", "expandedData"]:
937 a1_1=makeArray(arg_shape[1],[-8,8])
938 else:
939 a1_1=a_1
940 r=oper[1](a_0,a_1,sh_s)
941 r1=oper[1](a1_0,a1_1,sh_s)
942 text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
943 text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
944 text+=res_text
945 if case=="Symbol":
946 c0_res,c1_res=case0,case1
947 subs="{"
948 if case0=="Symbol":
949 text+=mkText("array","s0",a_0,a1_0)
950 subs+="arg0:s0"
951 c0_res="array"
952 if case1=="Symbol":
953 text+=mkText("array","s1",a_1,a1_1)
954 if not subs.endswith("{"): subs+=","
955 subs+="arg1:s1"
956 c1_res="array"
957 subs+="}"
958 text+=" sub=res.substitute(%s)\n"%subs
959 res="sub"
960 text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
961 else:
962 res="res"
963 text+=mkText(case,"ref",r,r1)
964 text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
965 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
966 if case0 == "taggedData" or case1 == "taggedData":
967 t_prog_with_tags+=text
968 else:
969 t_prog+=text
970 print test_header
971 print t_prog
972 # print t_prog_with_tags
973 print test_tail
974 1/0
975
976
977
978 #=======================================================================================================
979 # eigenvalues and eigen vectors 2D:
980 #=======================================================================================================
981 alltests= \
982 [ ("case0",[[0.0, 0.0], [0.0, 0.0]],[0.0, 0.0]) \
983 , ("case3",[[-1.0, 0.0], [0.0, -1.0]],[-1.0, -1.0]) \
984 , ("case5",[[-0.99999999999999967, -6.4606252205695602e-16], [-6.4606252205695602e-16, -0.99999999999999967]],[-1.0, -1.0]) \
985 , ("case6",[[0.0, 0.0], [0.0, 0.0001]],[0.0, 0.0001]) \
986 , ("case7",[[0.0001, 0.0], [0.0, 0.0]],[0.0, 0.0001]) \
987 , ("case8",[[6.0598371831785722e-06, 2.3859213977648625e-05], [2.3859213977648629e-05, 9.3940162816821425e-05]],[0.0, 0.0001]) \
988 , ("case9",[[1.0, 0.0], [0.0, 2.0]],[1.0, 2.0]) \
989 , ("case10",[[2.0, 0.0], [0.0, 1.0]],[1.0, 2.0]) \
990 , ("case11",[[1.0605983718317855, 0.23859213977648688], [0.23859213977648688, 1.9394016281682138]],[1.0, 2.0]) \
991 , ("case12",[[1.0, 0.0], [0.0, 1000000.0]],[1.0, 1000000.0]) \
992 , ("case13",[[1000000.0, 0.0], [0.0, 1.0]],[1.0, 1000000.0]) \
993 , ("case14",[[60599.311233413886, 238591.90118434647], [238591.90118434647, 939401.68876658613]],[1.0, 1000000.0]) \
994 ]
995 dim=2
996
997
998 alltests= \
999 [ ("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]) \
1000 , ("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]) \
1001 , ("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]) \
1002 , ("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]) \
1003 , ("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]) \
1004 , ("case13",[[1.0, 0.0, 0.0], [0.0, 0.9, 0.], [0.0, 0., 0.9]],[0.9, 0.9, 1.0]) \
1005 , ("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]) \
1006 , ("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]) \
1007 , ("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]) \
1008 , ("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]) \
1009 , ("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]) \
1010 , ("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]) \
1011 , ("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]) \
1012 , ("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]) \
1013 , ("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]) \
1014 , ("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]) \
1015 , ("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]) \
1016 , ("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]) \
1017 , ("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]) \
1018 , ("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]) \
1019 , ("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]) \
1020 ]
1021
1022 dim=3
1023
1024 alltests= \
1025 [ ("case0",[[0.0]],[0.0]) \
1026 , ("case1",[[1.0]],[1.0]) \
1027 ]
1028 dim=1
1029
1030 for case in ["constData","taggedData","expandedData"]:
1031 n=0
1032 while n<len(alltests):
1033 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1034 tname="test_eigenvalues_and_eigenvectors_%s_dim%s_%s"%(case,dim,alltests[n][0])
1035 text+=" def %s(self):\n"%tname
1036 a_0=numarray.array(alltests[n][1],numarray.Float64)
1037 ev_0=numarray.array(alltests[n][2],numarray.Float64)
1038 if case in ["taggedData", "expandedData"]:
1039 if n+1<len(alltests):
1040 a1_0=numarray.array(alltests[n+1][1],numarray.Float64)
1041 ev1_0=numarray.array(alltests[n+1][2],numarray.Float64)
1042 else:
1043 a1_0=numarray.array(alltests[0][1],numarray.Float64)
1044 ev1_0=numarray.array(alltests[0][2],numarray.Float64)
1045 n+=2
1046 else:
1047 a1_0=a_0
1048 ev1_0=ev_0
1049 n+=1
1050 text+=mkText(case,"arg",a_0,a1_0)
1051 text+=" res=eigenvalues_and_eigenvectors(arg)\n"
1052 text+=mkText(case,"ref_ev",ev_0,ev1_0)
1053 text+=mkTypeAndShapeTest(case,(dim,),"res[0]"," for eigenvalues")
1054 text+=mkTypeAndShapeTest(case,(dim,dim),"res[1]"," for eigenvectors")
1055 text+=" self.failUnless(Lsup(res[0]-ref_ev)<=self.RES_TOL*Lsup(ref_ev),\"wrong eigenvalues\")\n"
1056 for i in range(dim):
1057 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)
1058 text+=" self.failUnless(Lsup(matrixmult(transpose(res[1]),res[1])-kronecker(%s))<=self.RES_TOL,\"eigenvectors are not orthonormal\")\n"%dim
1059 if case == "taggedData" :
1060 t_prog_with_tags+=text
1061 else:
1062 t_prog+=text
1063 print test_header
1064 print t_prog
1065 print t_prog_with_tags
1066 print test_tail
1067 1/0
1068 #=======================================================================================================
1069 # get slices
1070 #=======================================================================================================
1071 from esys.escript import *
1072 for case0 in ["constData","taggedData","expandedData"]:
1073 for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,3,3)]:
1074 # get perm:
1075 if len(sh0)==2:
1076 check=[[1,0]]
1077 elif len(sh0)==3:
1078 check=[[1,0,2],
1079 [1,2,0],
1080 [2,1,0],
1081 [2,0,2],
1082 [0,2,1]]
1083 elif len(sh0)==4:
1084 check=[[0,1,3,2],
1085 [0,2,1,3],
1086 [0,2,3,1],
1087 [0,3,2,1],
1088 [0,3,1,2] ,
1089 [1,0,2,3],
1090 [1,0,3,2],
1091 [1,2,0,3],
1092 [1,2,3,0],
1093 [1,3,2,0],
1094 [1,3,0,2],
1095 [2,0,1,3],
1096 [2,0,3,1],
1097 [2,1,0,3],
1098 [2,1,3,0],
1099 [2,3,1,0],
1100 [2,3,0,1],
1101 [3,0,1,2],
1102 [3,0,2,1],
1103 [3,1,0,2],
1104 [3,1,2,0],
1105 [3,2,1,0],
1106 [3,2,0,1]]
1107 else:
1108 check=[]
1109
1110 # create the test cases:
1111 processed=[]
1112 l=["R","U","L","P","C","N"]
1113 c=[""]
1114 for i in range(len(sh0)):
1115 tmp=[]
1116 for ci in c:
1117 tmp+=[ci+li for li in l]
1118 c=tmp
1119 # SHUFFLE
1120 c2=[]
1121 while len(c)>0:
1122 i=int(random.random()*len(c))
1123 c2.append(c[i])
1124 del c[i]
1125 c=c2
1126 for ci in c:
1127 t=""
1128 sh=()
1129 sl=()
1130 for i in range(len(ci)):
1131 if ci[i]=="R":
1132 s="%s:%s"%(1,sh0[i]-1)
1133 sl=sl+(slice(1,sh0[i]-1),)
1134 sh=sh+(sh0[i]-2,)
1135 if ci[i]=="U":
1136 s=":%s"%(sh0[i]-1)
1137 sh=sh+(sh0[i]-1,)
1138 sl=sl+(slice(0,sh0[i]-1),)
1139 if ci[i]=="L":
1140 s="2:"
1141 sh=sh+(sh0[i]-2,)
1142 sl=sl+(slice(2,sh0[i]),)
1143 if ci[i]=="P":
1144 s="%s"%(int(sh0[i]/2))
1145 sl=sl+(int(sh0[i]/2),)
1146 if ci[i]=="C":
1147 s=":"
1148 sh=sh+(sh0[i],)
1149 sl=sl+(slice(0,sh0[i]),)
1150 if ci[i]=="N":
1151 s=""
1152 sh=sh+(sh0[i],)
1153 if len(s)>0:
1154 if not t=="": t+=","
1155 t+=s
1156 if len(sl)==1: sl=sl[0]
1157 N_found=False
1158 noN_found=False
1159 process=len(t)>0
1160 for i in ci:
1161 if i=="N":
1162 if not noN_found and N_found: process=False
1163 N_found=True
1164 else:
1165 if N_found: process=False
1166 noNfound=True
1167 # is there a similar one processed allready
1168 if process and ci.find("N")==-1:
1169 for ci2 in processed:
1170 for chi in check:
1171 is_perm=True
1172 for i in range(len(chi)):
1173 if not ci[i]==ci2[chi[i]]: is_perm=False
1174 if is_perm: process=False
1175 # if not process: print ci," rejected"
1176 if process:
1177 processed.append(ci)
1178 for case1 in ["array","constData","taggedData","expandedData"]:
1179 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1180 tname="test_setslice_%s_rank%s_%s_%s"%(case0,len(sh0),case1,ci)
1181 text+=" def %s(self):\n"%tname
1182 a_0=makeNumberedArray(sh0)
1183 if case0 in ["taggedData", "expandedData"]:
1184 a1_0=makeNumberedArray(sh0)
1185 else:
1186 a1_0=a_0*1.
1187
1188 a_1=makeNumberedArray(sh)
1189 if case1 in ["taggedData", "expandedData"]:
1190 a1_1=makeNumberedArray(sh)
1191 else:
1192 a1_1=a_1*1.
1193
1194 text+=mkText(case0,"arg",a_0,a1_0)
1195 text+=mkText(case1,"val",a_1,a1_1)
1196 text+=" arg[%s]=val\n"%t
1197 a_0.__setitem__(sl,a_1)
1198 a1_0.__setitem__(sl,a1_1)
1199 if Lsup(a_0-a1_0)==0.:
1200 text+=mkText("constData","ref",a_0,a1_0)
1201 else:
1202 text+=mkText("expandedData","ref",a_0,a1_0)
1203 text+=" self.failUnless(Lsup(arg-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"
1204
1205 if case0 == "taggedData" or case1 == "taggedData":
1206 t_prog_with_tags+=text
1207 else:
1208 t_prog+=text
1209
1210 print test_header
1211 # print t_prog
1212 print t_prog_with_tags
1213 print test_tail
1214 1/0
1215
1216 #=======================================================================================================
1217 # (non)symmetric part
1218 #=======================================================================================================
1219 from esys.escript import *
1220 for name in ["symmetric", "nonsymmetric"]:
1221 f=1.
1222 if name=="nonsymmetric": f=-1
1223 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1224 for sh0 in [ (3,3), (2,3,2,3)]:
1225 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1226 tname="test_%s_%s_rank%s"%(name,case0,len(sh0))
1227 text+=" def %s(self):\n"%tname
1228 a_0=makeNumberedArray(sh0,s=1.)
1229 r_0=(a_0+f*transpose(a_0))/2.
1230 if case0 in ["taggedData", "expandedData"]:
1231 a1_0=makeNumberedArray(sh0,s=-1.)
1232 r1_0=(a1_0+f*transpose(a1_0))/2.
1233 else:
1234 a1_0=a_0
1235 r1_0=r_0
1236 text+=mkText(case0,"arg",a_0,a1_0)
1237 text+=" res=%s(arg)\n"%name
1238 if case0=="Symbol":
1239 text+=mkText("array","s",a_0,a1_0)
1240 text+=" sub=res.substitute({arg:s})\n"
1241 res="sub"
1242 text+=mkText("array","ref",r_0,r1_0)
1243 else:
1244 res="res"
1245 text+=mkText(case0,"ref",r_0,r1_0)
1246 text+=mkTypeAndShapeTest(case0,sh0,"res")
1247 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1248
1249 if case0 == "taggedData" :
1250 t_prog_with_tags+=text
1251 else:
1252 t_prog+=text
1253 print test_header
1254 print t_prog
1255 # print t_prog_with_tags
1256 print test_tail
1257 1/0
1258
1259 #=======================================================================================================
1260 # eigenvalues
1261 #=======================================================================================================
1262 import numarray.linear_algebra
1263 name="eigenvalues"
1264 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1265 for sh0 in [ (1,1), (2,2), (3,3)]:
1266 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1267 tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1268 text+=" def %s(self):\n"%tname
1269 a_0=makeArray(sh0,[-1.,1])
1270 a_0=(a_0+numarray.transpose(a_0))/2.
1271 ev=numarray.linear_algebra.eigenvalues(a_0)
1272 ev.sort()
1273 if case0 in ["taggedData", "expandedData"]:
1274 a1_0=makeArray(sh0,[-1.,1])
1275 a1_0=(a1_0+numarray.transpose(a1_0))/2.
1276 ev1=numarray.linear_algebra.eigenvalues(a1_0)
1277 ev1.sort()
1278 else:
1279 a1_0=a_0
1280 ev1=ev
1281 text+=mkText(case0,"arg",a_0,a1_0)
1282 text+=" res=%s(arg)\n"%name
1283 if case0=="Symbol":
1284 text+=mkText("array","s",a_0,a1_0)
1285 text+=" sub=res.substitute({arg:s})\n"
1286 res="sub"
1287 text+=mkText("array","ref",ev,ev1)
1288 else:
1289 res="res"
1290 text+=mkText(case0,"ref",ev,ev1)
1291 text+=mkTypeAndShapeTest(case0,(sh0[0],),"res")
1292 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1293
1294 if case0 == "taggedData" :
1295 t_prog_with_tags+=text
1296 else:
1297 t_prog+=text
1298 print test_header
1299 # print t_prog
1300 print t_prog_with_tags
1301 print test_tail
1302 1/0
1303
1304 #=======================================================================================================
1305 # get slices
1306 #=======================================================================================================
1307 for case0 in ["constData","taggedData","expandedData","Symbol"]:
1308 for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,5,3)]:
1309 # get perm:
1310 if len(sh0)==2:
1311 check=[[1,0]]
1312 elif len(sh0)==3:
1313 check=[[1,0,2],
1314 [1,2,0],
1315 [2,1,0],
1316 [2,0,2],
1317 [0,2,1]]
1318 elif len(sh0)==4:
1319 check=[[0,1,3,2],
1320 [0,2,1,3],
1321 [0,2,3,1],
1322 [0,3,2,1],
1323 [0,3,1,2] ,
1324 [1,0,2,3],
1325 [1,0,3,2],
1326 [1,2,0,3],
1327 [1,2,3,0],
1328 [1,3,2,0],
1329 [1,3,0,2],
1330 [2,0,1,3],
1331 [2,0,3,1],
1332 [2,1,0,3],
1333 [2,1,3,0],
1334 [2,3,1,0],
1335 [2,3,0,1],
1336 [3,0,1,2],
1337 [3,0,2,1],
1338 [3,1,0,2],
1339 [3,1,2,0],
1340 [3,2,1,0],
1341 [3,2,0,1]]
1342 else:
1343 check=[]
1344
1345 # create the test cases:
1346 processed=[]
1347 l=["R","U","L","P","C","N"]
1348 c=[""]
1349 for i in range(len(sh0)):
1350 tmp=[]
1351 for ci in c:
1352 tmp+=[ci+li for li in l]
1353 c=tmp
1354 # SHUFFLE
1355 c2=[]
1356 while len(c)>0:
1357 i=int(random.random()*len(c))
1358 c2.append(c[i])
1359 del c[i]
1360 c=c2
1361 for ci in c:
1362 t=""
1363 sh=()
1364 for i in range(len(ci)):
1365 if ci[i]=="R":
1366 s="%s:%s"%(1,sh0[i]-1)
1367 sh=sh+(sh0[i]-2,)
1368 if ci[i]=="U":
1369 s=":%s"%(sh0[i]-1)
1370 sh=sh+(sh0[i]-1,)
1371 if ci[i]=="L":
1372 s="2:"
1373 sh=sh+(sh0[i]-2,)
1374 if ci[i]=="P":
1375 s="%s"%(int(sh0[i]/2))
1376 if ci[i]=="C":
1377 s=":"
1378 sh=sh+(sh0[i],)
1379 if ci[i]=="N":
1380 s=""
1381 sh=sh+(sh0[i],)
1382 if len(s)>0:
1383 if not t=="": t+=","
1384 t+=s
1385 N_found=False
1386 noN_found=False
1387 process=len(t)>0
1388 for i in ci:
1389 if i=="N":
1390 if not noN_found and N_found: process=False
1391 N_found=True
1392 else:
1393 if N_found: process=False
1394 noNfound=True
1395 # is there a similar one processed allready
1396 if process and ci.find("N")==-1:
1397 for ci2 in processed:
1398 for chi in check:
1399 is_perm=True
1400 for i in range(len(chi)):
1401 if not ci[i]==ci2[chi[i]]: is_perm=False
1402 if is_perm: process=False
1403 # if not process: print ci," rejected"
1404 if process:
1405 processed.append(ci)
1406 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1407 tname="test_getslice_%s_rank%s_%s"%(case0,len(sh0),ci)
1408 text+=" def %s(self):\n"%tname
1409 a_0=makeNumberedArray(sh0,s=1)
1410 if case0 in ["taggedData", "expandedData"]:
1411 a1_0=makeNumberedArray(sh0,s=-1.)
1412 else:
1413 a1_0=a_0
1414 r=eval("a_0[%s]"%t)
1415 r1=eval("a1_0[%s]"%t)
1416 text+=mkText(case0,"arg",a_0,a1_0)
1417 text+=" res=arg[%s]\n"%t
1418 if case0=="Symbol":
1419 text+=mkText("array","s",a_0,a1_0)
1420 text+=" sub=res.substitute({arg:s})\n"
1421 res="sub"
1422 text+=mkText("array","ref",r,r1)
1423 else:
1424 res="res"
1425 text+=mkText(case0,"ref",r,r1)
1426 text+=mkTypeAndShapeTest(case0,sh,"res")
1427 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1428
1429 if case0 == "taggedData" :
1430 t_prog_with_tags+=text
1431 else:
1432 t_prog+=text
1433
1434 print test_header
1435 # print t_prog
1436 print t_prog_with_tags
1437 print test_tail
1438 1/0
1439 #============================================================================================
1440 def innerTEST(arg0,arg1):
1441 if isinstance(arg0,float):
1442 out=numarray.array(arg0*arg1)
1443 else:
1444 out=(arg0*arg1).sum()
1445 return out
1446
1447 def outerTEST(arg0,arg1):
1448 if isinstance(arg0,float):
1449 out=numarray.array(arg0*arg1)
1450 elif isinstance(arg1,float):
1451 out=numarray.array(arg0*arg1)
1452 else:
1453 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1454 return out
1455
1456 def tensorProductTest(arg0,arg1,sh_s):
1457 if isinstance(arg0,float):
1458 out=numarray.array(arg0*arg1)
1459 elif isinstance(arg1,float):
1460 out=numarray.array(arg0*arg1)
1461 elif len(sh_s)==0:
1462 out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1463 else:
1464 l=len(sh_s)
1465 sh0=arg0.shape[:arg0.rank-l]
1466 sh1=arg1.shape[l:]
1467 ls,l0,l1=1,1,1
1468 for i in sh_s: ls*=i
1469 for i in sh0: l0*=i
1470 for i in sh1: l1*=i
1471 out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
1472 out2=numarray.zeros((l0,l1),numarray.Float)
1473 for i0 in range(l0):
1474 for i1 in range(l1):
1475 for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
1476 out=out2.resize(sh0+sh1)
1477 return out
1478
1479 def testMatrixMult(arg0,arg1,sh_s):
1480 return numarray.matrixmultiply(arg0,arg1)
1481
1482
1483 def testTensorMult(arg0,arg1,sh_s):
1484 if len(arg0)==2:
1485 return numarray.matrixmultiply(arg0,arg1)
1486 else:
1487 if arg1.rank==4:
1488 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
1489 for i0 in range(arg0.shape[0]):
1490 for i1 in range(arg0.shape[1]):
1491 for i2 in range(arg0.shape[2]):
1492 for i3 in range(arg0.shape[3]):
1493 for j2 in range(arg1.shape[2]):
1494 for j3 in range(arg1.shape[3]):
1495 out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
1496 elif arg1.rank==3:
1497 out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
1498 for i0 in range(arg0.shape[0]):
1499 for i1 in range(arg0.shape[1]):
1500 for i2 in range(arg0.shape[2]):
1501 for i3 in range(arg0.shape[3]):
1502 for j2 in range(arg1.shape[2]):
1503 out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
1504 elif arg1.rank==2:
1505 out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
1506 for i0 in range(arg0.shape[0]):
1507 for i1 in range(arg0.shape[1]):
1508 for i2 in range(arg0.shape[2]):
1509 for i3 in range(arg0.shape[3]):
1510 out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
1511 return out
1512
1513 def testReduce(arg0,init_val,test_expr,post_expr):
1514 out=init_val
1515 if isinstance(arg0,float):
1516 out=eval(test_expr.replace("%a1%","arg0"))
1517 elif arg0.rank==0:
1518 out=eval(test_expr.replace("%a1%","arg0"))
1519 elif arg0.rank==1:
1520 for i0 in range(arg0.shape[0]):
1521 out=eval(test_expr.replace("%a1%","arg0[i0]"))
1522 elif arg0.rank==2:
1523 for i0 in range(arg0.shape[0]):
1524 for i1 in range(arg0.shape[1]):
1525 out=eval(test_expr.replace("%a1%","arg0[i0,i1]"))
1526 elif arg0.rank==3:
1527 for i0 in range(arg0.shape[0]):
1528 for i1 in range(arg0.shape[1]):
1529 for i2 in range(arg0.shape[2]):
1530 out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2]"))
1531 elif arg0.rank==4:
1532 for i0 in range(arg0.shape[0]):
1533 for i1 in range(arg0.shape[1]):
1534 for i2 in range(arg0.shape[2]):
1535 for i3 in range(arg0.shape[3]):
1536 out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))
1537 return eval(post_expr)
1538
1539 def clipTEST(arg0,mn,mx):
1540 if isinstance(arg0,float):
1541 return max(min(arg0,mx),mn)
1542 out=numarray.zeros(arg0.shape,numarray.Float64)
1543 if arg0.rank==1:
1544 for i0 in range(arg0.shape[0]):
1545 out[i0]=max(min(arg0[i0],mx),mn)
1546 elif arg0.rank==2:
1547 for i0 in range(arg0.shape[0]):
1548 for i1 in range(arg0.shape[1]):
1549 out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
1550 elif arg0.rank==3:
1551 for i0 in range(arg0.shape[0]):
1552 for i1 in range(arg0.shape[1]):
1553 for i2 in range(arg0.shape[2]):
1554 out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
1555 elif arg0.rank==4:
1556 for i0 in range(arg0.shape[0]):
1557 for i1 in range(arg0.shape[1]):
1558 for i2 in range(arg0.shape[2]):
1559 for i3 in range(arg0.shape[3]):
1560 out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
1561 return out
1562 def minimumTEST(arg0,arg1):
1563 if isinstance(arg0,float):
1564 if isinstance(arg1,float):
1565 if arg0>arg1:
1566 return arg1
1567 else:
1568 return arg0
1569 else:
1570 arg0=numarray.ones(arg1.shape)*arg0
1571 else:
1572 if isinstance(arg1,float):
1573 arg1=numarray.ones(arg0.shape)*arg1
1574 out=numarray.zeros(arg0.shape,numarray.Float64)
1575 if arg0.rank==0:
1576 if arg0>arg1:
1577 out=arg1
1578 else:
1579 out=arg0
1580 elif arg0.rank==1:
1581 for i0 in range(arg0.shape[0]):
1582 if arg0[i0]>arg1[i0]:
1583 out[i0]=arg1[i0]
1584 else:
1585 out[i0]=arg0[i0]
1586 elif arg0.rank==2:
1587 for i0 in range(arg0.shape[0]):
1588 for i1 in range(arg0.shape[1]):
1589 if arg0[i0,i1]>arg1[i0,i1]:
1590 out[i0,i1]=arg1[i0,i1]
1591 else:
1592 out[i0,i1]=arg0[i0,i1]
1593 elif arg0.rank==3:
1594 for i0 in range(arg0.shape[0]):
1595 for i1 in range(arg0.shape[1]):
1596 for i2 in range(arg0.shape[2]):
1597 if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
1598 out[i0,i1,i2]=arg1[i0,i1,i2]
1599 else:
1600 out[i0,i1,i2]=arg0[i0,i1,i2]
1601 elif arg0.rank==4:
1602 for i0 in range(arg0.shape[0]):
1603 for i1 in range(arg0.shape[1]):
1604 for i2 in range(arg0.shape[2]):
1605 for i3 in range(arg0.shape[3]):
1606 if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
1607 out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
1608 else:
1609 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
1610 return out
1611
1612 def unrollLoops(a,b,o,arg,tap="",x="x"):
1613 out=""
1614 if a.rank==1:
1615 z=""
1616 for i99 in range(a.shape[0]):
1617 if not z=="": z+="+"
1618 if o=="1":
1619 z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
1620 else:
1621 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
1622
1623 out+=tap+"%s=%s\n"%(arg,z)
1624
1625 elif a.rank==2:
1626 for i0 in range(a.shape[0]):
1627 z=""
1628 for i99 in range(a.shape[1]):
1629 if not z=="": z+="+"
1630 if o=="1":
1631 z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
1632 else:
1633 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
1634
1635 out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
1636 elif a.rank==3:
1637 for i0 in range(a.shape[0]):
1638 for i1 in range(a.shape[1]):
1639 z=""
1640 for i99 in range(a.shape[2]):
1641 if not z=="": z+="+"
1642 if o=="1":
1643 z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
1644 else:
1645 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
1646
1647 out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
1648 elif a.rank==4:
1649 for i0 in range(a.shape[0]):
1650 for i1 in range(a.shape[1]):
1651 for i2 in range(a.shape[2]):
1652 z=""
1653 for i99 in range(a.shape[3]):
1654 if not z=="": z+="+"
1655 if o=="1":
1656 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
1657 else:
1658 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
1659
1660 out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
1661 elif a.rank==5:
1662 for i0 in range(a.shape[0]):
1663 for i1 in range(a.shape[1]):
1664 for i2 in range(a.shape[2]):
1665 for i3 in range(a.shape[3]):
1666 z=""
1667 for i99 in range(a.shape[4]):
1668 if not z=="": z+="+"
1669 if o=="1":
1670 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
1671 else:
1672 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
1673
1674 out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
1675 return out
1676
1677 def unrollLoopsOfGrad(a,b,o,arg,tap=""):
1678 out=""
1679 if a.rank==1:
1680 for i99 in range(a.shape[0]):
1681 if o=="1":
1682 out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
1683 else:
1684 out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
1685
1686 elif a.rank==2:
1687 for i0 in range(a.shape[0]):
1688 for i99 in range(a.shape[1]):
1689 if o=="1":
1690 out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
1691 else:
1692 out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
1693
1694 elif a.rank==3:
1695 for i0 in range(a.shape[0]):
1696 for i1 in range(a.shape[1]):
1697 for i99 in range(a.shape[2]):
1698 if o=="1":
1699 out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
1700 else:
1701 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])
1702
1703 elif a.rank==4:
1704 for i0 in range(a.shape[0]):
1705 for i1 in range(a.shape[1]):
1706 for i2 in range(a.shape[2]):
1707 for i99 in range(a.shape[3]):
1708 if o=="1":
1709 out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
1710 else:
1711 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])
1712 return out
1713 def unrollLoopsOfDiv(a,b,o,arg,tap=""):
1714 out=tap+arg+"="
1715 if o=="1":
1716 z=0.
1717 for i99 in range(a.shape[0]):
1718 z+=b[i99,i99]+a[i99,i99]
1719 out+="(%s)"%z
1720 else:
1721 z=0.
1722 for i99 in range(a.shape[0]):
1723 z+=b[i99,i99]
1724 if i99>0: out+="+"
1725 out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
1726 out+="+(%s)"%z
1727 return out
1728
1729 def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
1730 if where=="Function":
1731 xfac_o=1.
1732 xfac_op=0.
1733 z_fac=1./2.
1734 z_fac_s=""
1735 zo_fac_s="/(o+1.)"
1736 elif where=="FunctionOnBoundary":
1737 xfac_o=1.
1738 xfac_op=0.
1739 z_fac=1.
1740 z_fac_s="*dim"
1741 zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
1742 elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1743 xfac_o=0.
1744 xfac_op=1.
1745 z_fac=1./2.
1746 z_fac_s=""
1747 zo_fac_s="/(o+1.)"
1748 out=""
1749 if a.rank==1:
1750 zo=0.
1751 zop=0.
1752 z=0.
1753 for i99 in range(a.shape[0]):
1754 if i99==0:
1755 zo+= xfac_o*a[i99]
1756 zop+= xfac_op*a[i99]
1757 else:
1758 zo+=a[i99]
1759 z+=b[i99]
1760
1761 out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
1762 if zop==0.:
1763 out+="\n"
1764 else:
1765 out+="+(%s)*0.5**o\n"%zop
1766 elif a.rank==2:
1767 for i0 in range(a.shape[0]):
1768 zo=0.
1769 zop=0.
1770 z=0.
1771 for i99 in range(a.shape[1]):
1772 if i99==0:
1773 zo+= xfac_o*a[i0,i99]
1774 zop+= xfac_op*a[i0,i99]
1775 else:
1776 zo+=a[i0,i99]
1777 z+=b[i0,i99]
1778
1779 out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
1780 if zop==0.:
1781 out+="\n"
1782 else:
1783 out+="+(%s)*0.5**o\n"%zop
1784 elif a.rank==3:
1785 for i0 in range(a.shape[0]):
1786 for i1 in range(a.shape[1]):
1787 zo=0.
1788 zop=0.
1789 z=0.
1790 for i99 in range(a.shape[2]):
1791 if i99==0:
1792 zo+= xfac_o*a[i0,i1,i99]
1793 zop+= xfac_op*a[i0,i1,i99]
1794 else:
1795 zo+=a[i0,i1,i99]
1796 z+=b[i0,i1,i99]
1797
1798 out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
1799 if zop==0.:
1800 out+="\n"
1801 else:
1802 out+="+(%s)*0.5**o\n"%zop
1803 elif a.rank==4:
1804 for i0 in range(a.shape[0]):
1805 for i1 in range(a.shape[1]):
1806 for i2 in range(a.shape[2]):
1807 zo=0.
1808 zop=0.
1809 z=0.
1810 for i99 in range(a.shape[3]):
1811 if i99==0:
1812 zo+= xfac_o*a[i0,i1,i2,i99]
1813 zop+= xfac_op*a[i0,i1,i2,i99]
1814
1815 else:
1816 zo+=a[i0,i1,i2,i99]
1817 z+=b[i0,i1,i2,i99]
1818
1819 out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
1820 if zop==0.:
1821 out+="\n"
1822 else:
1823 out+="+(%s)*0.5**o\n"%zop
1824
1825 elif a.rank==5:
1826 for i0 in range(a.shape[0]):
1827 for i1 in range(a.shape[1]):
1828 for i2 in range(a.shape[2]):
1829 for i3 in range(a.shape[3]):
1830 zo=0.
1831 zop=0.
1832 z=0.
1833 for i99 in range(a.shape[4]):
1834 if i99==0:
1835 zo+= xfac_o*a[i0,i1,i2,i3,i99]
1836 zop+= xfac_op*a[i0,i1,i2,i3,i99]
1837
1838 else:
1839 zo+=a[i0,i1,i2,i3,i99]
1840 z+=b[i0,i1,i2,i3,i99]
1841 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)
1842 if zop==0.:
1843 out+="\n"
1844 else:
1845 out+="+(%s)*0.5**o\n"%zop
1846
1847 return out
1848 def unrollLoopsSimplified(b,arg,tap=""):
1849 out=""
1850 if isinstance(b,float) or b.rank==0:
1851 out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
1852
1853 elif b.rank==1:
1854 for i0 in range(b.shape[0]):
1855 out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
1856 elif b.rank==2:
1857 for i0 in range(b.shape[0]):
1858 for i1 in range(b.shape[1]):
1859 out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
1860 elif b.rank==3:
1861 for i0 in range(b.shape[0]):
1862 for i1 in range(b.shape[1]):
1863 for i2 in range(b.shape[2]):
1864 out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
1865 elif b.rank==4:
1866 for i0 in range(b.shape[0]):
1867 for i1 in range(b.shape[1]):
1868 for i2 in range(b.shape[2]):
1869 for i3 in range(b.shape[3]):
1870 out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
1871 return out
1872
1873 def unrollLoopsOfL2(b,where,arg,tap=""):
1874 out=""
1875 z=[]
1876 if isinstance(b,float) or b.rank==0:
1877 z.append(b**2)
1878 elif b.rank==1:
1879 for i0 in range(b.shape[0]):
1880 z.append(b[i0]**2)
1881 elif b.rank==2:
1882 for i1 in range(b.shape[1]):
1883 s=0
1884 for i0 in range(b.shape[0]):
1885 s+=b[i0,i1]**2
1886 z.append(s)
1887 elif b.rank==3:
1888 for i2 in range(b.shape[2]):
1889 s=0
1890 for i0 in range(b.shape[0]):
1891 for i1 in range(b.shape[1]):
1892 s+=b[i0,i1,i2]**2
1893 z.append(s)
1894
1895 elif b.rank==4:
1896 for i3 in range(b.shape[3]):
1897 s=0
1898 for i0 in range(b.shape[0]):
1899 for i1 in range(b.shape[1]):
1900 for i2 in range(b.shape[2]):
1901 s+=b[i0,i1,i2,i3]**2
1902 z.append(s)
1903 if where=="Function":
1904 xfac_o=1.
1905 xfac_op=0.
1906 z_fac_s=""
1907 zo_fac_s=""
1908 zo_fac=1./3.
1909 elif where=="FunctionOnBoundary":
1910 xfac_o=1.
1911 xfac_op=0.
1912 z_fac_s="*dim"
1913 zo_fac_s="*(2.*dim+1.)/3."
1914 zo_fac=1.
1915 elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1916 xfac_o=0.
1917 xfac_op=1.
1918 z_fac_s=""
1919 zo_fac_s=""
1920 zo_fac=1./3.
1921 zo=0.
1922 zop=0.
1923 for i99 in range(len(z)):
1924 if i99==0:
1925 zo+=xfac_o*z[i99]
1926 zop+=xfac_op*z[i99]
1927 else:
1928 zo+=z[i99]
1929 out+=tap+"%s=sqrt((%s)%s"%(arg,zo*zo_fac,zo_fac_s)
1930 if zop==0.:
1931 out+=")\n"
1932 else:
1933 out+="+(%s))\n"%(zop*0.5**2)
1934 return out
1935 #=======================================================================================================
1936 # transpose
1937 #=======================================================================================================
1938 def transposeTest(r,offset):
1939 if isinstance(r,float): return r
1940 s=r.shape
1941 s1=1
1942 for i in s[:offset]: s1*=i
1943 s2=1
1944 for i in s[offset:]: s2*=i
1945 out=numarray.reshape(r,(s1,s2))
1946 out.transpose()
1947 return numarray.resize(out,s[offset:]+s[:offset])
1948
1949 name,tt="transpose",transposeTest
1950 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1951 for sh0 in [ (), (3,), (4,5), (6,2,2),(3,2,3,4)]:
1952 for offset in range(len(sh0)+1):
1953 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1954 tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1955 text+=" def %s(self):\n"%tname
1956 sh_t=sh0[offset:]+sh0[:offset]
1957
1958 # sh_t=list(sh0)
1959 # sh_t[offset+1]=sh_t[offset]
1960 # sh_t=tuple(sh_t)
1961 # sh_r=[]
1962 # for i in range(offset): sh_r.append(sh0[i])
1963 # for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1964 # sh_r=tuple(sh_r)
1965
1966 a_0=makeArray(sh0,[-1.,1])
1967 if case0 in ["taggedData", "expandedData"]:
1968 a1_0=makeArray(sh0,[-1.,1])
1969 else:
1970 a1_0=a_0
1971 r=tt(a_0,offset)
1972 r1=tt(a1_0,offset)
1973 text+=mkText(case0,"arg",a_0,a1_0)
1974 text+=" res=%s(arg,%s)\n"%(name,offset)
1975 if case0=="Symbol":
1976 text+=mkText("array","s",a_0,a1_0)
1977 text+=" sub=res.substitute({arg:s})\n"
1978 res="sub"
1979 text+=mkText("array","ref",r,r1)
1980 else:
1981 res="res"
1982 text+=mkText(case0,"ref",r,r1)
1983 text+=mkTypeAndShapeTest(case0,sh_t,"res")
1984 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1985
1986 if case0 == "taggedData" :
1987 t_prog_with_tags+=text
1988 else:
1989 t_prog+=text
1990
1991 print test_header
1992 # print t_prog
1993 print t_prog_with_tags
1994 print test_tail
1995 1/0
1996 #=======================================================================================================
1997 # L2
1998 #=======================================================================================================
1999 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2000 for data in ["Data","Symbol"]:
2001 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
2002 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2003 tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
2004 text+=" def %s(self):\n"%tname
2005 text+=" \"\"\"\n"
2006 text+=" tests L2-norm of %s on the %s\n\n"%(data,where)
2007 text+=" assumptions: self.domain supports integration on %s\n"%where
2008 text+=" \"\"\"\n"
2009 text+=" dim=self.domain.getDim()\n"
2010 text+=" w=%s(self.domain)\n"%where
2011 text+=" x=w.getX()\n"
2012 o="1"
2013 if len(sh)>0:
2014 sh_2=sh[:len(sh)-1]+(2,)
2015 sh_3=sh[:len(sh)-1]+(3,)
2016 b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
2017 b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
2018 else:
2019 sh_2=()
2020 sh_3=()
2021 b_2=makeArray(sh,[-1.,1])
2022 b_3=makeArray(sh,[-1.,1])
2023
2024 if data=="Symbol":
2025 val="s"
2026 res="sub"
2027 else:
2028 val="arg"
2029 res="res"
2030 text+=" if dim==2:\n"
2031 if data=="Symbol":
2032 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
2033
2034 text+=" %s=Data(0,%s,w)\n"%(val,sh_2)
2035 text+=unrollLoopsSimplified(b_2,val,tap=" ")
2036 text+=unrollLoopsOfL2(b_2,where,"ref",tap=" ")
2037 text+="\n else:\n"
2038 if data=="Symbol":
2039 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
2040 text+=" %s=Data(0,%s,w)\n"%(val,sh_3)
2041 text+=unrollLoopsSimplified(b_3,val,tap=" ")
2042 text+=unrollLoopsOfL2(b_3,where,"ref",tap=" ")
2043 text+="\n res=L2(arg)\n"
2044 if data=="Symbol":
2045 text+=" sub=res.substitute({arg:s})\n"
2046 text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
2047 text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
2048 else:
2049 text+=" self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
2050 text+=" self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
2051 t_prog+=text
2052 print t_prog
2053 1/0
2054
2055 #=======================================================================================================
2056 # div
2057 #=======================================================================================================
2058 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2059 for data in ["Data","Symbol"]:
2060 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
2061 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2062 tname="test_div_on%s_from%s_%s"%(where,data,case)
2063 text+=" def %s(self):\n"%tname
2064 text+=" \"\"\"\n"
2065 text+=" tests divergence of %s on the %s\n\n"%(data,where)
2066 text+=" assumptions: %s(self.domain) exists\n"%case
2067 text+=" self.domain supports div on %s\n"%where
2068 text+=" \"\"\"\n"
2069 if case=="ReducedSolution":
2070 text+=" o=1\n"
2071 o="1"
2072 else:
2073 text+=" o=self.order\n"
2074 o="o"
2075 text+=" dim=self.domain.getDim()\n"
2076 text+=" w_ref=%s(self.domain)\n"%where
2077 text+=" x_ref=w_ref.getX()\n"
2078 text+=" w=%s(self.domain)\n"%case
2079 text+=" x=w.getX()\n"
2080 a_2=makeArray((2,2),[-1.,1])
2081 b_2=makeArray((2,2),[-1.,1])
2082 a_3=makeArray((3,3),[-1.,1])
2083 b_3=makeArray((3,3),[-1.,1])
2084 if data=="Symbol":
2085 text+=" arg=Symbol(shape=(dim,),dim=dim)\n"
2086 val="s"
2087 res="sub"
2088 else:
2089 val="arg"
2090 res="res"
2091 text+=" %s=Vector(0,w)\n"%val
2092 text+=" if dim==2:\n"
2093 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
2094 text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap=" ")
2095 text+="\n else:\n"
2096
2097 text+=unrollLoops(a_3,b_3,o,val,tap=" ")
2098 text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap=" ")
2099 text+="\n res=div(arg,where=w_ref)\n"
2100 if data=="Symbol":
2101 text+=" sub=res.substitute({arg:s})\n"
2102 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
2103 text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
2104 text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
2105 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
2106
2107
2108 t_prog+=text
2109 print t_prog
2110 1/0
2111
2112 #=======================================================================================================
2113 # interpolation
2114 #=======================================================================================================
2115 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
2116 for data in ["Data","Symbol"]:
2117 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2118 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
2119 if where==case or \
2120 ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
2121 ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
2122 (case=="ContinuousFunction" and where in ["Solution","ReducedSolution"]) or \
2123 (case=="Solution" and where=="ReducedSolution") :
2124
2125
2126 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2127 tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
2128 text+=" def %s(self):\n"%tname
2129 text+=" \"\"\"\n"
2130 text+=" tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
2131 text+=" assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
2132 text+=" \"\"\"\n"
2133 if case=="ReducedSolution" or where=="ReducedSolution":
2134 text+=" o=1\n"
2135 o="1"
2136 else:
2137 text+=" o=self.order\n"
2138 o="o"
2139 text+=" dim=self.domain.getDim()\n"
2140 text+=" w_ref=%s(self.domain)\n"%where
2141 text+=" x_ref=w_ref.getX()\n"
2142 text+=" w=%s(self.domain)\n"%case
2143 text+=" x=w.getX()\n"
2144 a_2=makeArray(sh+(2,),[-1.,1])
2145 b_2=makeArray(sh+(2,),[-1.,1])
2146 a_3=makeArray(sh+(3,),[-1.,1])
2147 b_3=makeArray(sh+(3,),[-1.,1])
2148 if data=="Symbol":
2149 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
2150 val="s"
2151 res="sub"
2152 else:
2153 val="arg"
2154 res="res"
2155 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
2156 text+=" ref=Data(0,%s,w_ref)\n"%str(sh)
2157 text+=" if dim==2:\n"
2158 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
2159 text+=unrollLoops(a_2,b_2,o,"ref",tap=" ",x="x_ref"