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

Contents of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2169 - (show annotations)
Wed Dec 17 03:08:58 2008 UTC (10 years, 10 months ago) by caltinay
File size: 165665 byte(s)
Assorted spelling, grammar, whitespace and copy/paste error fixes (Part 2).
All epydoc warnings for these files have been fixed.
This commit should be a no-op.

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 eigenvectors 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
1563 def minimumTEST(arg0,arg1):
1564 if isinstance(arg0,float):
1565 if isinstance(arg1,float):
1566 if arg0>arg1:
1567 return arg1
1568 else:
1569 return arg0
1570 else:
1571 arg0=numarray.ones(arg1.shape)*arg0
1572 else:
1573 if isinstance(arg1,float):
1574 arg1=numarray.ones(arg0.shape)*arg1
1575 out=numarray.zeros(arg0.shape,numarray.Float64)
1576 if arg0.rank==0:
1577 if arg0>arg1:
1578 out=arg1
1579 else:
1580 out=arg0
1581 elif arg0.rank==1:
1582 for i0 in range(arg0.shape[0]):
1583 if arg0[i0]>arg1[i0]:
1584 out[i0]=arg1[i0]
1585 else:
1586 out[i0]=arg0[i0]
1587 elif arg0.rank==2:
1588 for i0 in range(arg0.shape[0]):
1589 for i1 in range(arg0.shape[1]):
1590 if arg0[i0,i1]>arg1[i0,i1]:
1591 out[i0,i1]=arg1[i0,i1]
1592 else:
1593 out[i0,i1]=arg0[i0,i1]
1594 elif arg0.rank==3:
1595 for i0 in range(arg0.shape[0]):
1596 for i1 in range(arg0.shape[1]):
1597 for i2 in range(arg0.shape[2]):
1598 if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
1599 out[i0,i1,i2]=arg1[i0,i1,i2]
1600 else:
1601 out[i0,i1,i2]=arg0[i0,i1,i2]
1602 elif arg0.rank==4:
1603 for i0 in range(arg0.shape[0]):
1604 for i1 in range(arg0.shape[1]):
1605 for i2 in range(arg0.shape[2]):
1606 for i3 in range(arg0.shape[3]):
1607 if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
1608 out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
1609 else:
1610 out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
1611 return out
1612
1613 def unrollLoops(a,b,o,arg,tap="",x="x"):
1614 out=""
1615 if a.rank==1:
1616 z=""
1617 for i99 in range(a.shape[0]):
1618 if not z=="": z+="+"
1619 if o=="1":
1620 z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
1621 else:
1622 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
1623
1624 out+=tap+"%s=%s\n"%(arg,z)
1625
1626 elif a.rank==2:
1627 for i0 in range(a.shape[0]):
1628 z=""
1629 for i99 in range(a.shape[1]):
1630 if not z=="": z+="+"
1631 if o=="1":
1632 z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
1633 else:
1634 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
1635
1636 out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
1637 elif a.rank==3:
1638 for i0 in range(a.shape[0]):
1639 for i1 in range(a.shape[1]):
1640 z=""
1641 for i99 in range(a.shape[2]):
1642 if not z=="": z+="+"
1643 if o=="1":
1644 z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
1645 else:
1646 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
1647
1648 out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
1649 elif a.rank==4:
1650 for i0 in range(a.shape[0]):
1651 for i1 in range(a.shape[1]):
1652 for i2 in range(a.shape[2]):
1653 z=""
1654 for i99 in range(a.shape[3]):
1655 if not z=="": z+="+"
1656 if o=="1":
1657 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
1658 else:
1659 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
1660
1661 out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
1662 elif a.rank==5:
1663 for i0 in range(a.shape[0]):
1664 for i1 in range(a.shape[1]):
1665 for i2 in range(a.shape[2]):
1666 for i3 in range(a.shape[3]):
1667 z=""
1668 for i99 in range(a.shape[4]):
1669 if not z=="": z+="+"
1670 if o=="1":
1671 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
1672 else:
1673 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
1674
1675 out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
1676 return out
1677
1678 def unrollLoopsOfGrad(a,b,o,arg,tap=""):
1679 out=""
1680 if a.rank==1:
1681 for i99 in range(a.shape[0]):
1682 if o=="1":
1683 out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
1684 else:
1685 out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
1686
1687 elif a.rank==2:
1688 for i0 in range(a.shape[0]):
1689 for i99 in range(a.shape[1]):
1690 if o=="1":
1691 out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
1692 else:
1693 out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
1694
1695 elif a.rank==3:
1696 for i0 in range(a.shape[0]):
1697 for i1 in range(a.shape[1]):
1698 for i99 in range(a.shape[2]):
1699 if o=="1":
1700 out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
1701 else:
1702 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])
1703
1704 elif a.rank==4:
1705 for i0 in range(a.shape[0]):
1706 for i1 in range(a.shape[1]):
1707 for i2 in range(a.shape[2]):
1708 for i99 in range(a.shape[3]):
1709 if o=="1":
1710 out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
1711 else:
1712 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])
1713 return out
1714
1715 def unrollLoopsOfDiv(a,b,o,arg,tap=""):
1716 out=tap+arg+"="
1717 if o=="1":
1718 z=0.
1719 for i99 in range(a.shape[0]):
1720 z+=b[i99,i99]+a[i99,i99]
1721 out+="(%s)"%z
1722 else:
1723 z=0.
1724 for i99 in range(a.shape[0]):
1725 z+=b[i99,i99]
1726 if i99>0: out+="+"
1727 out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
1728 out+="+(%s)"%z
1729 return out
1730
1731 def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
1732 if where=="Function":
1733 xfac_o=1.
1734 xfac_op=0.
1735 z_fac=1./2.
1736 z_fac_s=""
1737 zo_fac_s="/(o+1.)"
1738 elif where=="FunctionOnBoundary":
1739 xfac_o=1.
1740 xfac_op=0.
1741 z_fac=1.
1742 z_fac_s="*dim"
1743 zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
1744 elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1745 xfac_o=0.
1746 xfac_op=1.
1747 z_fac=1./2.
1748 z_fac_s=""
1749 zo_fac_s="/(o+1.)"
1750 out=""
1751 if a.rank==1:
1752 zo=0.
1753 zop=0.
1754 z=0.
1755 for i99 in range(a.shape[0]):
1756 if i99==0:
1757 zo+= xfac_o*a[i99]
1758 zop+= xfac_op*a[i99]
1759 else:
1760 zo+=a[i99]
1761 z+=b[i99]
1762
1763 out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
1764 if zop==0.:
1765 out+="\n"
1766 else:
1767 out+="+(%s)*0.5**o\n"%zop
1768 elif a.rank==2:
1769 for i0 in range(a.shape[0]):
1770 zo=0.
1771 zop=0.
1772 z=0.
1773 for i99 in range(a.shape[1]):
1774 if i99==0:
1775 zo+= xfac_o*a[i0,i99]
1776 zop+= xfac_op*a[i0,i99]
1777 else:
1778 zo+=a[i0,i99]
1779 z+=b[i0,i99]
1780
1781 out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
1782 if zop==0.:
1783 out+="\n"
1784 else:
1785 out+="+(%s)*0.5**o\n"%zop
1786 elif a.rank==3:
1787 for i0 in range(a.shape[0]):
1788 for i1 in range(a.shape[1]):
1789 zo=0.
1790 zop=0.
1791 z=0.
1792 for i99 in range(a.shape[2]):
1793 if i99==0:
1794 zo+= xfac_o*a[i0,i1,i99]
1795 zop+= xfac_op*a[i0,i1,i99]
1796 else:
1797 zo+=a[i0,i1,i99]
1798 z+=b[i0,i1,i99]
1799
1800 out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
1801 if zop==0.:
1802 out+="\n"
1803 else:
1804 out+="+(%s)*0.5**o\n"%zop
1805 elif a.rank==4:
1806 for i0 in range(a.shape[0]):
1807 for i1 in range(a.shape[1]):
1808 for i2 in range(a.shape[2]):
1809 zo=0.
1810 zop=0.
1811 z=0.
1812 for i99 in range(a.shape[3]):
1813 if i99==0:
1814 zo+= xfac_o*a[i0,i1,i2,i99]
1815 zop+= xfac_op*a[i0,i1,i2,i99]
1816
1817 else:
1818 zo+=a[i0,i1,i2,i99]
1819 z+=b[i0,i1,i2,i99]
1820
1821 out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
1822 if zop==0.:
1823 out+="\n"
1824 else:
1825 out+="+(%s)*0.5**o\n"%zop
1826
1827 elif a.rank==5:
1828 for i0 in range(a.shape[0]):
1829 for i1 in range(a.shape[1]):
1830 for i2 in range(a.shape[2]):
1831 for i3 in range(a.shape[3]):
1832 zo=0.
1833 zop=0.
1834 z=0.
1835 for i99 in range(a.shape[4]):
1836 if i99==0:
1837 zo+= xfac_o*a[i0,i1,i2,i3,i99]
1838 zop+= xfac_op*a[i0,i1,i2,i3,i99]
1839
1840 else:
1841 zo+=a[i0,i1,i2,i3,i99]
1842 z+=b[i0,i1,i2,i3,i99]
1843 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)
1844 if zop==0.:
1845 out+="\n"
1846 else:
1847 out+="+(%s)*0.5**o\n"%zop
1848
1849 return out
1850
1851 def unrollLoopsSimplified(b,arg,tap=""):
1852 out=""
1853 if isinstance(b,float) or b.rank==0:
1854 out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
1855
1856 elif b.rank==1:
1857 for i0 in range(b.shape[0]):
1858 out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
1859 elif b.rank==2:
1860 for i0 in range(b.shape[0]):
1861 for i1 in range(b.shape[1]):
1862 out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
1863 elif b.rank==3:
1864 for i0 in range(b.shape[0]):
1865 for i1 in range(b.shape[1]):
1866 for i2 in range(b.shape[2]):
1867 out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
1868 elif b.rank==4:
1869 for i0 in range(b.shape[0]):
1870 for i1 in range(b.shape[1]):
1871 for i2 in range(b.shape[2]):
1872 for i3 in range(b.shape[3]):
1873 out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
1874 return out
1875
1876 def unrollLoopsOfL2(b,where,arg,tap=""):
1877 out=""
1878 z=[]
1879 if isinstance(b,float) or b.rank==0:
1880 z.append(b**2)
1881 elif b.rank==1:
1882 for i0 in range(b.shape[0]):
1883 z.append(b[i0]**2)
1884 elif b.rank==2:
1885 for i1 in range(b.shape[1]):
1886 s=0
1887 for i0 in range(b.shape[0]):
1888 s+=b[i0,i1]**2
1889 z.append(s)
1890 elif b.rank==3:
1891 for i2 in range(b.shape[2]):
1892 s=0
1893 for i0 in range(b.shape[0]):
1894 for i1 in range(b.shape[1]):
1895 s+=b[i0,i1,i2]**2
1896 z.append(s)
1897
1898 elif b.rank==4:
1899 for i3 in range(b.shape[3]):
1900 s=0
1901 for i0 in range(b.shape[0]):
1902 for i1 in range(b.shape[1]):
1903 for i2 in range(b.shape[2]):
1904 s+=b[i0,i1,i2,i3]**2
1905 z.append(s)
1906 if where=="Function":
1907 xfac_o=1.
1908 xfac_op=0.
1909 z_fac_s=""
1910 zo_fac_s=""
1911 zo_fac=1./3.
1912 elif where=="FunctionOnBoundary":
1913 xfac_o=1.
1914 xfac_op=0.
1915 z_fac_s="*dim"
1916 zo_fac_s="*(2.*dim+1.)/3."
1917 zo_fac=1.
1918 elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1919 xfac_o=0.
1920 xfac_op=1.
1921 z_fac_s=""
1922 zo_fac_s=""
1923 zo_fac=1./3.
1924 zo=0.
1925 zop=0.
1926 for i99 in range(len(z)):
1927 if i99==0:
1928 zo+=xfac_o*z[i99]
1929 zop+=xfac_op*z[i99]
1930 else:
1931 zo+=z[i99]
1932 out+=tap+"%s=sqrt((%s)%s"%(arg,zo*zo_fac,zo_fac_s)
1933 if zop==0.:
1934 out+=")\n"
1935 else:
1936 out+="+(%s))\n"%(zop*0.5**2)
1937 return out
1938 #==============================================================================
1939 # transpose
1940 #==============================================================================
1941 def transposeTest(r,offset):
1942 if isinstance(r,float): return r
1943 s=r.shape
1944 s1=1
1945 for i in s[:offset]: s1*=i
1946 s2=1
1947 for i in s[offset:]: s2*=i
1948 out=numarray.reshape(r,(s1,s2))
1949 out.transpose()
1950 return numarray.resize(out,s[offset:]+s[:offset])
1951
1952 name,tt="transpose",transposeTest
1953 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1954 for sh0 in [ (), (3,), (4,5), (6,2,2),(3,2,3,4)]:
1955 for offset in range(len(sh0)+1):
1956 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1957 tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1958 text+=" def %s(self):\n"%tname
1959 sh_t=sh0[offset:]+sh0[:offset]
1960
1961 # sh_t=list(sh0)
1962 # sh_t[offset+1]=sh_t[offset]
1963 # sh_t=tuple(sh_t)
1964 # sh_r=[]
1965 # for i in range(offset): sh_r.append(sh0[i])
1966 # for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1967 # sh_r=tuple(sh_r)
1968
1969 a_0=makeArray(sh0,[-1.,1])
1970 if case0 in ["taggedData", "expandedData"]:
1971 a1_0=makeArray(sh0,[-1.,1])
1972 else:
1973 a1_0=a_0
1974 r=tt(a_0,offset)
1975 r1=tt(a1_0,offset)
1976 text+=mkText(case0,"arg",a_0,a1_0)
1977 text+=" res=%s(arg,%s)\n"%(name,offset)
1978 if case0=="Symbol":
1979 text+=mkText("array","s",a_0,a1_0)
1980 text+=" sub=res.substitute({arg:s})\n"
1981 res="sub"
1982 text+=mkText("array","ref",r,r1)
1983 else:
1984 res="res"
1985 text+=mkText(case0,"ref",r,r1)
1986 text+=mkTypeAndShapeTest(case0,sh_t,"res")
1987 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1988
1989 if case0 == "taggedData" :
1990 t_prog_with_tags+=text
1991 else:
1992 t_prog+=text
1993
1994 print test_header
1995 # print t_prog
1996 print t_prog_with_tags
1997 print test_tail
1998 1/0
1999 #==============================================================================
2000 # L2
2001 #==============================================================================
2002 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2003 for data in ["Data","Symbol"]:
2004 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
2005 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2006 tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
2007 text+=" def %s(self):\n"%tname
2008 text+=" \"\"\"\n"
2009 text+=" tests L2-norm of %s on the %s\n\n"%(data,where)
2010 text+=" assumptions: self.domain supports integration on %s\n"%where
2011 text+=" \"\"\"\n"
2012 text+=" dim=self.domain.getDim()\n"
2013 text+=" w=%s(self.domain)\n"%where
2014 text+=" x=w.getX()\n"
2015 o="1"
2016 if len(sh)>0:
2017 sh_2=sh[:len(sh)-1]+(2,)
2018 sh_3=sh[:len(sh)-1]+(3,)
2019 b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
2020 b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
2021 else:
2022 sh_2=()
2023 sh_3=()
2024 b_2=makeArray(sh,[-1.,1])
2025 b_3=makeArray(sh,[-1.,1])
2026
2027 if data=="Symbol":
2028 val="s"
2029 res="sub"
2030 else:
2031 val="arg"
2032 res="res"
2033 text+=" if dim==2:\n"
2034 if data=="Symbol":
2035 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
2036
2037 text+=" %s=Data(0,%s,w)\n"%(val,sh_2)
2038 text+=unrollLoopsSimplified(b_2,val,tap=" ")
2039 text+=unrollLoopsOfL2(b_2,where,"ref",tap=" ")
2040 text+="\n else:\n"
2041 if data=="Symbol":
2042 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
2043 text+=" %s=Data(0,%s,w)\n"%(val,sh_3)
2044 text+=unrollLoopsSimplified(b_3,val,tap=" ")
2045 text+=unrollLoopsOfL2(b_3,where,"ref",tap=" ")
2046 text+="\n res=L2(arg)\n"
2047 if data=="Symbol":
2048 text+=" sub=res.substitute({arg:s})\n"
2049 text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
2050 text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
2051 else:
2052 text+=" self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
2053 text+=" self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
2054 t_prog+=text
2055 print t_prog
2056 1/0
2057
2058 #==============================================================================
2059 # div
2060 #==============================================================================
2061 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2062 for data in ["Data","Symbol"]:
2063 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
2064 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2065 tname="test_div_on%s_from%s_%s"%(where,data,case)
2066 text+=" def %s(self):\n"%tname
2067 text+=" \"\"\"\n"
2068 text+=" tests divergence of %s on the %s\n\n"%(data,where)
2069 text+=" assumptions: %s(self.domain) exists\n"%case
2070 text+=" self.domain supports div on %s\n"%where
2071 text+=" \"\"\"\n"
2072 if case=="ReducedSolution":
2073 text+=" o=1\n"
2074 o="1"
2075 else:
2076 text+=" o=self.order\n"
2077 o="o"
2078 text+=" dim=self.domain.getDim()\n"
2079 text+=" w_ref=%s(self.domain)\n"%where
2080 text+=" x_ref=w_ref.getX()\n"
2081 text+=" w=%s(self.domain)\n"%case
2082 text+=" x=w.getX()\n"
2083 a_2=makeArray((2,2),[-1.,1])
2084 b_2=makeArray((2,2),[-1.,1])
2085 a_3=makeArray((3,3),[-1.,1])
2086 b_3=makeArray((3,3),[-1.,1])
2087 if data=="Symbol":
2088 text+=" arg=Symbol(shape=(dim,),dim=dim)\n"
2089 val="s"
2090 res="sub"
2091 else:
2092 val="arg"
2093 res="res"
2094 text+=" %s=Vector(0,w)\n"%val
2095 text+=" if dim==2:\n"
2096 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
2097 text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap=" ")
2098 text+="\n else:\n"
2099
2100 text+=unrollLoops(a_3,b_3,o,val,tap=" ")
2101 text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap=" ")
2102 text+="\n res=div(arg,where=w_ref)\n"
2103 if data=="Symbol":
2104 text+=" sub=res.substitute({arg:s})\n"
2105 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
2106 text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
2107 text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
2108 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
2109
2110
2111 t_prog+=text
2112 print t_prog
2113 1/0
2114
2115 #==============================================================================
2116 # interpolation
2117 #==============================================================================
2118 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
2119 for data in ["Data","Symbol"]:
2120 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
2121 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
2122 if where==case or \
2123 ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
2124 ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
2125 (case=="ContinuousFunction" and where in ["Solution","ReducedSolution"]) or \
2126 (case=="Solution" and where=="ReducedSolution") :
2127
2128
2129 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2130 tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
2131 text+=" def %s(self):\n"%tname
2132 text+=" \"\"\"\n"
2133 text+=" tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
2134 text+=" assumptions: self.domain supports interpolation from %s onto %s\n"%(case,where)
2135 text+=" \"\"\"\n"
2136 if case=="ReducedSolution" or where=="ReducedSolution":
2137 text+=" o=1\n"
2138 o="1"
2139 else:
2140 text+=" o=self.order\n"
2141 o="o"
2142 text+=" dim=self.domain.getDim()\n"
2143 text+=" w_ref=%s(self.domain)\n"%where
2144 text+=" x_ref=w_ref.getX()\n"
2145 text+=" w=%s(self.domain)\n"%case
2146 text+=" x=w.getX()\n"
2147 a_2=makeArray(sh+(2,),[-1.,1])
2148 b_2=makeArray(sh+(2,),[-1.,1])
2149 a_3=makeArray(sh+(3,),[-1.,1])
2150 b_3=makeArray(sh+(3,),[-1.,1])
2151 if data=="Symbol":
2152 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
2153 val="s"
2154 res="sub"
2155 else:
2156 val="arg"
2157 res="res"
2158 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
2159 text+=" ref=Data(0,%s,w_ref)\n"%str(sh)
2160 text+=" if dim==2:\n"
2161 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
2162 text+=unrollLoops(a_2,b_2,o,"ref",tap=" ",x="x_ref")
2163 text+=" else:\n"
2164
2165 text+=unrollLoops(a_3,b_3,o,val