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

Contents of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


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