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

Contents of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


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