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

Annotation of /trunk/escript/py_src/generatediff

Parent Directory Parent Directory | Revision Log Revision Log


Revision 517 - (hide annotations)
Tue Feb 14 02:25:02 2006 UTC (13 years, 8 months ago) by gross
Original Path: trunk/escript/py_src/generateutil
File size: 135913 byte(s)
tests for __getitem__ added
1 jgs 154 #!/usr/bin/python
2 gross 283 # $Id$
3 jgs 154
4     """
5     program generates parts of the util.py and the test_util.py script
6     """
7 gross 157 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 jgs 154 t_prog=""
26 gross 157 t_prog_with_tags=""
27     t_prog_failing=""
28     u_prog=""
29 jgs 154
30 gross 157 def wherepos(arg):
31     if arg>0.:
32     return 1.
33     else:
34     return 0.
35    
36 gross 313
37 gross 157 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 jgs 154 import random
53     import numarray
54     import math
55     finc=1.e-6
56    
57 gross 157 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 jgs 154
152 gross 157
153 jgs 154 def makeArray(shape,rng):
154     l=rng[1]-rng[0]
155     out=numarray.zeros(shape,numarray.Float64)
156     if len(shape)==0:
157 gross 157 out=l*random.random()+rng[0]
158 jgs 154 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 gross 437 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*random.random()+rng[0]
183 jgs 154 else:
184 gross 437 raise SystemError,"rank is restricted to 5"
185 jgs 154 return out
186    
187 gross 517 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*i0
194     elif len(shape)==2:
195     for i0 in range(shape[0]):
196     for i1 in range(shape[1]):
197     out[i0,i1]=s*(i1+shape[1]*i0)
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*(i2+shape[2]*i1+shape[2]*shape[1]*i0)
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*(i3+shape[3]*i2+shape[3]*shape[2]*i1+shape[3]*shape[2]*shape[1]*i0)
209     else:
210     raise SystemError,"rank is restricted to 4"
211     return out
212 jgs 154
213 gross 157 def makeResult(val,test_expr):
214 jgs 154 if isinstance(val,float):
215 gross 157 out=eval(test_expr.replace("%a1%","val"))
216 jgs 154 elif len(val.shape)==0:
217 gross 157 out=eval(test_expr.replace("%a1%","val"))
218 jgs 154 elif len(val.shape)==1:
219     out=numarray.zeros(val.shape,numarray.Float64)
220     for i0 in range(val.shape[0]):
221 gross 157 out[i0]=eval(test_expr.replace("%a1%","val[i0]"))
222 jgs 154 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 gross 157 out[i0,i1]=eval(test_expr.replace("%a1%","val[i0,i1]"))
227 jgs 154 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 gross 157 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val[i0,i1,i2]"))
233 jgs 154 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 gross 157 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val[i0,i1,i2,i3]"))
240 jgs 154 else:
241     raise SystemError,"rank is restricted to 4"
242     return out
243    
244 gross 157 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 jgs 154
461 gross 157
462     def mkText(case,name,a,a1=None,use_tagging_for_expanded_data=False):
463 jgs 154 t_out=""
464 gross 157 if case=="float":
465 jgs 154 if isinstance(a,float):
466     t_out+=" %s=%s\n"%(name,a)
467 gross 291 elif a.rank==0:
468 jgs 154 t_out+=" %s=%s\n"%(name,a)
469     else:
470     t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
471 gross 157 elif case=="array":
472     if isinstance(a,float):
473     t_out+=" %s=numarray.array(%s)\n"%(name,a)
474 gross 291 elif a.rank==0:
475 gross 157 t_out+=" %s=numarray.array(%s)\n"%(name,a)
476     else:
477     t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
478 jgs 154 elif case=="constData":
479     if isinstance(a,float):
480     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
481 gross 291 elif a.rank==0:
482 jgs 154 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 gross 291 elif a.rank==0:
490 jgs 154 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 gross 157 t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
495 jgs 154 elif case=="expandedData":
496 gross 157 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 gross 291 elif a.rank==0:
501 gross 157 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=whereNegative(self.functionspace.getX()[0]-0.5)\n"%name
509     if isinstance(a,float):
510     t_out+=" %s=msk_%s*(%s)+(1.-msk_%s)*(%s)\n"%(name,name,a,name,a1)
511 gross 291 elif a.rank==0:
512 gross 157 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 jgs 154 elif case=="Symbol":
516     if isinstance(a,float):
517     t_out+=" %s=Symbol(shape=())\n"%(name)
518 gross 291 elif a.rank==0:
519 jgs 154 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 gross 157 def mkTypeAndShapeTest(case,sh,argstr):
526     text=""
527     if case=="float":
528     text+=" self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%argstr
529     elif case=="array":
530     text+=" self.failUnless(isinstance(%s,numarray.NumArray),\"wrong type of result.\")\n"%argstr
531     text+=" self.failUnlessEqual(%s.shape,%s,\"wrong shape of result.\")\n"%(argstr,str(sh))
532     elif case in ["constData","taggedData","expandedData"]:
533     text+=" self.failUnless(isinstance(%s,Data),\"wrong type of result.\")\n"%argstr
534     text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result.\")\n"%(argstr,str(sh))
535     elif case=="Symbol":
536     text+=" self.failUnless(isinstance(%s,Symbol),\"wrong type of result.\")\n"%argstr
537     text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result.\")\n"%(argstr,str(sh))
538     return text
539 jgs 154
540 gross 157 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 jgs 154
553 gross 517 #=======================================================================================================
554     # slicing
555     #=======================================================================================================
556     for case0 in ["constData","taggedData","expandedData","Symbol"]:
557     for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,5,3)]:
558     # get perm:
559     if len(sh0)==2:
560     check=[[1,0]]
561     elif len(sh0)==3:
562     check=[[1,0,2],
563     [1,2,0],
564     [2,1,0],
565     [2,0,2],
566     [0,2,1]]
567     elif len(sh0)==4:
568     check=[[0,1,3,2],
569     [0,2,1,3],
570     [0,2,3,1],
571     [0,3,2,1],
572     [0,3,1,2] ,
573     [1,0,2,3],
574     [1,0,3,2],
575     [1,2,0,3],
576     [1,2,3,0],
577     [1,3,2,0],
578     [1,3,0,2],
579     [2,0,1,3],
580     [2,0,3,1],
581     [2,1,0,3],
582     [2,1,3,0],
583     [2,3,1,0],
584     [2,3,0,1],
585     [3,0,1,2],
586     [3,0,2,1],
587     [3,1,0,2],
588     [3,1,2,0],
589     [3,2,1,0],
590     [3,2,0,1]]
591     else:
592     check=[]
593    
594     # create the test cases:
595     processed=[]
596     l=["R","U","L","P","C","N"]
597     c=[""]
598     for i in range(len(sh0)):
599     tmp=[]
600     for ci in c:
601     tmp+=[ci+li for li in l]
602     c=tmp
603     # SHUFFLE
604     c2=[]
605     while len(c)>0:
606     i=int(random.random()*len(c))
607     c2.append(c[i])
608     del c[i]
609     c=c2
610     for ci in c:
611     t=""
612     sh=()
613     for i in range(len(ci)):
614     if ci[i]=="R":
615     s="%s:%s"%(1,sh0[i]-1)
616     sh=sh+(sh0[i]-2,)
617     if ci[i]=="U":
618     s=":%s"%(sh0[i]-1)
619     sh=sh+(sh0[i]-1,)
620     if ci[i]=="L":
621     s="2:"
622     sh=sh+(sh0[i]-2,)
623     if ci[i]=="P":
624     s="%s"%(int(sh0[i]/2))
625     if ci[i]=="C":
626     s=":"
627     sh=sh+(sh0[i],)
628     if ci[i]=="N":
629     s=""
630     sh=sh+(sh0[i],)
631     if len(s)>0:
632     if not t=="": t+=","
633     t+=s
634     N_found=False
635     noN_found=False
636     process=len(t)>0
637     for i in ci:
638     if i=="N":
639     if not noN_found and N_found: process=False
640     N_found=True
641     else:
642     if N_found: process=False
643     noNfound=True
644     # is there a similar one processed allready
645     if process and ci.find("N")==-1:
646     for ci2 in processed:
647     for chi in check:
648     is_perm=True
649     for i in range(len(chi)):
650     if not ci[i]==ci2[chi[i]]: is_perm=False
651     if is_perm: process=False
652     # if not process: print ci," rejected"
653     if process:
654     processed.append(ci)
655     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
656     tname="test_getslice_%s_rank%s_%s"%(case0,len(sh0),ci)
657     text+=" def %s(self):\n"%tname
658     a_0=makeNumberedArray(sh0,s=1)
659     if case0 in ["taggedData", "expandedData"]:
660     a1_0=makeNumberedArray(sh0,s=-1.)
661     else:
662     a1_0=a_0
663     r=eval("a_0[%s]"%t)
664     r1=eval("a1_0[%s]"%t)
665     text+=mkText(case0,"arg",a_0,a1_0)
666     text+=" res=arg[%s]\n"%t
667     if case0=="Symbol":
668     text+=mkText("array","s",a_0,a1_0)
669     text+=" sub=res.substitute({arg:s})\n"
670     res="sub"
671     text+=mkText("array","ref",r,r1)
672     else:
673     res="res"
674     text+=mkText(case0,"ref",r,r1)
675     text+=mkTypeAndShapeTest(case0,sh,"res")
676     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
677    
678     if case0 == "taggedData" :
679     t_prog_with_tags+=text
680     else:
681     t_prog+=text
682    
683     print test_header
684     # print t_prog
685     print t_prog_with_tags
686     print test_tail
687     1/0
688     #============================================================================================
689 gross 291 def innerTEST(arg0,arg1):
690     if isinstance(arg0,float):
691     out=numarray.array(arg0*arg1)
692     else:
693     out=(arg0*arg1).sum()
694     return out
695    
696     def outerTEST(arg0,arg1):
697     if isinstance(arg0,float):
698     out=numarray.array(arg0*arg1)
699     elif isinstance(arg1,float):
700     out=numarray.array(arg0*arg1)
701     else:
702     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
703     return out
704    
705     def tensorProductTest(arg0,arg1,sh_s):
706     if isinstance(arg0,float):
707     out=numarray.array(arg0*arg1)
708     elif isinstance(arg1,float):
709     out=numarray.array(arg0*arg1)
710     elif len(sh_s)==0:
711     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
712     else:
713     l=len(sh_s)
714     sh0=arg0.shape[:arg0.rank-l]
715     sh1=arg1.shape[l:]
716     ls,l0,l1=1,1,1
717     for i in sh_s: ls*=i
718     for i in sh0: l0*=i
719     for i in sh1: l1*=i
720     out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
721     out2=numarray.zeros((l0,l1),numarray.Float)
722     for i0 in range(l0):
723     for i1 in range(l1):
724     for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
725     out=out2.resize(sh0+sh1)
726     return out
727    
728     def testMatrixMult(arg0,arg1,sh_s):
729     return numarray.matrixmultiply(arg0,arg1)
730    
731    
732     def testTensorMult(arg0,arg1,sh_s):
733     if len(arg0)==2:
734     return numarray.matrixmultiply(arg0,arg1)
735     else:
736     if arg1.rank==4:
737     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
738     for i0 in range(arg0.shape[0]):
739     for i1 in range(arg0.shape[1]):
740     for i2 in range(arg0.shape[2]):
741     for i3 in range(arg0.shape[3]):
742     for j2 in range(arg1.shape[2]):
743     for j3 in range(arg1.shape[3]):
744     out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
745     elif arg1.rank==3:
746     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
747     for i0 in range(arg0.shape[0]):
748     for i1 in range(arg0.shape[1]):
749     for i2 in range(arg0.shape[2]):
750     for i3 in range(arg0.shape[3]):
751     for j2 in range(arg1.shape[2]):
752     out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
753     elif arg1.rank==2:
754     out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
755     for i0 in range(arg0.shape[0]):
756     for i1 in range(arg0.shape[1]):
757     for i2 in range(arg0.shape[2]):
758     for i3 in range(arg0.shape[3]):
759     out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
760     return out
761 gross 313
762     def testReduce(arg0,init_val,test_expr,post_expr):
763     out=init_val
764     if isinstance(arg0,float):
765     out=eval(test_expr.replace("%a1%","arg0"))
766     elif arg0.rank==0:
767     out=eval(test_expr.replace("%a1%","arg0"))
768     elif arg0.rank==1:
769     for i0 in range(arg0.shape[0]):
770     out=eval(test_expr.replace("%a1%","arg0[i0]"))
771     elif arg0.rank==2:
772     for i0 in range(arg0.shape[0]):
773     for i1 in range(arg0.shape[1]):
774     out=eval(test_expr.replace("%a1%","arg0[i0,i1]"))
775     elif arg0.rank==3:
776     for i0 in range(arg0.shape[0]):
777     for i1 in range(arg0.shape[1]):
778     for i2 in range(arg0.shape[2]):
779     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2]"))
780     elif arg0.rank==4:
781     for i0 in range(arg0.shape[0]):
782     for i1 in range(arg0.shape[1]):
783     for i2 in range(arg0.shape[2]):
784     for i3 in range(arg0.shape[3]):
785     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))
786     return eval(post_expr)
787 gross 396
788     def clipTEST(arg0,mn,mx):
789     if isinstance(arg0,float):
790     return max(min(arg0,mx),mn)
791     out=numarray.zeros(arg0.shape,numarray.Float64)
792     if arg0.rank==1:
793     for i0 in range(arg0.shape[0]):
794     out[i0]=max(min(arg0[i0],mx),mn)
795     elif arg0.rank==2:
796     for i0 in range(arg0.shape[0]):
797     for i1 in range(arg0.shape[1]):
798     out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
799     elif arg0.rank==3:
800     for i0 in range(arg0.shape[0]):
801     for i1 in range(arg0.shape[1]):
802     for i2 in range(arg0.shape[2]):
803     out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
804     elif arg0.rank==4:
805     for i0 in range(arg0.shape[0]):
806     for i1 in range(arg0.shape[1]):
807     for i2 in range(arg0.shape[2]):
808     for i3 in range(arg0.shape[3]):
809     out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
810     return out
811     def minimumTEST(arg0,arg1):
812     if isinstance(arg0,float):
813     if isinstance(arg1,float):
814     if arg0>arg1:
815     return arg1
816     else:
817     return arg0
818     else:
819     arg0=numarray.ones(arg1.shape)*arg0
820     else:
821     if isinstance(arg1,float):
822     arg1=numarray.ones(arg0.shape)*arg1
823     out=numarray.zeros(arg0.shape,numarray.Float64)
824     if arg0.rank==0:
825     if arg0>arg1:
826     out=arg1
827     else:
828     out=arg0
829     elif arg0.rank==1:
830     for i0 in range(arg0.shape[0]):
831     if arg0[i0]>arg1[i0]:
832     out[i0]=arg1[i0]
833     else:
834     out[i0]=arg0[i0]
835     elif arg0.rank==2:
836     for i0 in range(arg0.shape[0]):
837     for i1 in range(arg0.shape[1]):
838     if arg0[i0,i1]>arg1[i0,i1]:
839     out[i0,i1]=arg1[i0,i1]
840     else:
841     out[i0,i1]=arg0[i0,i1]
842     elif arg0.rank==3:
843     for i0 in range(arg0.shape[0]):
844     for i1 in range(arg0.shape[1]):
845     for i2 in range(arg0.shape[2]):
846     if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
847     out[i0,i1,i2]=arg1[i0,i1,i2]
848     else:
849     out[i0,i1,i2]=arg0[i0,i1,i2]
850     elif arg0.rank==4:
851     for i0 in range(arg0.shape[0]):
852     for i1 in range(arg0.shape[1]):
853     for i2 in range(arg0.shape[2]):
854     for i3 in range(arg0.shape[3]):
855     if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
856     out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
857     else:
858     out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
859     return out
860 gross 443
861 gross 439 def unrollLoops(a,b,o,arg,tap="",x="x"):
862 gross 437 out=""
863     if a.rank==1:
864     z=""
865     for i99 in range(a.shape[0]):
866     if not z=="": z+="+"
867     if o=="1":
868 gross 439 z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
869 gross 437 else:
870 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
871 gross 437
872     out+=tap+"%s=%s\n"%(arg,z)
873    
874     elif a.rank==2:
875     for i0 in range(a.shape[0]):
876     z=""
877     for i99 in range(a.shape[1]):
878     if not z=="": z+="+"
879     if o=="1":
880     z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
881     else:
882 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
883 gross 437
884     out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
885     elif a.rank==3:
886     for i0 in range(a.shape[0]):
887     for i1 in range(a.shape[1]):
888     z=""
889     for i99 in range(a.shape[2]):
890     if not z=="": z+="+"
891     if o=="1":
892 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
893 gross 437 else:
894 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
895 gross 437
896     out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
897     elif a.rank==4:
898     for i0 in range(a.shape[0]):
899     for i1 in range(a.shape[1]):
900     for i2 in range(a.shape[2]):
901     z=""
902     for i99 in range(a.shape[3]):
903     if not z=="": z+="+"
904     if o=="1":
905 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
906 gross 437 else:
907 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
908 gross 437
909     out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
910     elif a.rank==5:
911     for i0 in range(a.shape[0]):
912     for i1 in range(a.shape[1]):
913     for i2 in range(a.shape[2]):
914     for i3 in range(a.shape[3]):
915     z=""
916     for i99 in range(a.shape[4]):
917     if not z=="": z+="+"
918     if o=="1":
919 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
920 gross 437 else:
921 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
922 gross 437
923     out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
924     return out
925    
926     def unrollLoopsOfGrad(a,b,o,arg,tap=""):
927     out=""
928     if a.rank==1:
929     for i99 in range(a.shape[0]):
930     if o=="1":
931     out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
932     else:
933     out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
934    
935     elif a.rank==2:
936     for i0 in range(a.shape[0]):
937     for i99 in range(a.shape[1]):
938     if o=="1":
939     out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
940     else:
941     out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
942    
943     elif a.rank==3:
944     for i0 in range(a.shape[0]):
945     for i1 in range(a.shape[1]):
946     for i99 in range(a.shape[2]):
947     if o=="1":
948     out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
949     else:
950     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])
951    
952     elif a.rank==4:
953     for i0 in range(a.shape[0]):
954     for i1 in range(a.shape[1]):
955     for i2 in range(a.shape[2]):
956     for i99 in range(a.shape[3]):
957     if o=="1":
958     out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
959     else:
960     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])
961 gross 438 return out
962 gross 441 def unrollLoopsOfDiv(a,b,o,arg,tap=""):
963     out=tap+arg+"="
964     if o=="1":
965     z=0.
966     for i99 in range(a.shape[0]):
967     z+=b[i99,i99]+a[i99,i99]
968     out+="(%s)"%z
969     else:
970     z=0.
971     for i99 in range(a.shape[0]):
972     z+=b[i99,i99]
973     if i99>0: out+="+"
974     out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
975     out+="+(%s)"%z
976     return out
977    
978 gross 437 def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
979     if where=="Function":
980     xfac_o=1.
981     xfac_op=0.
982     z_fac=1./2.
983     z_fac_s=""
984     zo_fac_s="/(o+1.)"
985     elif where=="FunctionOnBoundary":
986     xfac_o=1.
987     xfac_op=0.
988     z_fac=1.
989     z_fac_s="*dim"
990     zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
991     elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
992     xfac_o=0.
993     xfac_op=1.
994     z_fac=1./2.
995     z_fac_s=""
996     zo_fac_s="/(o+1.)"
997     out=""
998     if a.rank==1:
999     zo=0.
1000     zop=0.
1001     z=0.
1002     for i99 in range(a.shape[0]):
1003     if i99==0:
1004     zo+= xfac_o*a[i99]
1005     zop+= xfac_op*a[i99]
1006     else:
1007     zo+=a[i99]
1008     z+=b[i99]
1009    
1010     out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
1011     if zop==0.:
1012     out+="\n"
1013     else:
1014     out+="+(%s)*0.5**o\n"%zop
1015     elif a.rank==2:
1016     for i0 in range(a.shape[0]):
1017     zo=0.
1018     zop=0.
1019     z=0.
1020     for i99 in range(a.shape[1]):
1021     if i99==0:
1022     zo+= xfac_o*a[i0,i99]
1023     zop+= xfac_op*a[i0,i99]
1024     else:
1025     zo+=a[i0,i99]
1026     z+=b[i0,i99]
1027    
1028     out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
1029     if zop==0.:
1030     out+="\n"
1031     else:
1032     out+="+(%s)*0.5**o\n"%zop
1033     elif a.rank==3:
1034     for i0 in range(a.shape[0]):
1035     for i1 in range(a.shape[1]):
1036     zo=0.
1037     zop=0.
1038     z=0.
1039     for i99 in range(a.shape[2]):
1040     if i99==0:
1041     zo+= xfac_o*a[i0,i1,i99]
1042     zop+= xfac_op*a[i0,i1,i99]
1043     else:
1044     zo+=a[i0,i1,i99]
1045     z+=b[i0,i1,i99]
1046    
1047     out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
1048     if zop==0.:
1049     out+="\n"
1050     else:
1051     out+="+(%s)*0.5**o\n"%zop
1052     elif a.rank==4:
1053     for i0 in range(a.shape[0]):
1054     for i1 in range(a.shape[1]):
1055     for i2 in range(a.shape[2]):
1056     zo=0.
1057     zop=0.
1058     z=0.
1059     for i99 in range(a.shape[3]):
1060     if i99==0:
1061     zo+= xfac_o*a[i0,i1,i2,i99]
1062     zop+= xfac_op*a[i0,i1,i2,i99]
1063    
1064     else:
1065     zo+=a[i0,i1,i2,i99]
1066     z+=b[i0,i1,i2,i99]
1067    
1068     out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
1069     if zop==0.:
1070     out+="\n"
1071     else:
1072     out+="+(%s)*0.5**o\n"%zop
1073    
1074     elif a.rank==5:
1075     for i0 in range(a.shape[0]):
1076     for i1 in range(a.shape[1]):
1077     for i2 in range(a.shape[2]):
1078     for i3 in range(a.shape[3]):
1079     zo=0.
1080     zop=0.
1081     z=0.
1082     for i99 in range(a.shape[4]):
1083     if i99==0:
1084     zo+= xfac_o*a[i0,i1,i2,i3,i99]
1085     zop+= xfac_op*a[i0,i1,i2,i3,i99]
1086    
1087     else:
1088     zo+=a[i0,i1,i2,i3,i99]
1089     z+=b[i0,i1,i2,i3,i99]
1090     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)
1091     if zop==0.:
1092     out+="\n"
1093     else:
1094     out+="+(%s)*0.5**o\n"%zop
1095    
1096     return out
1097 gross 443 def unrollLoopsSimplified(b,arg,tap=""):
1098     out=""
1099     if isinstance(b,float) or b.rank==0:
1100     out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
1101    
1102     elif b.rank==1:
1103     for i0 in range(b.shape[0]):
1104     out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
1105     elif b.rank==2:
1106     for i0 in range(b.shape[0]):
1107     for i1 in range(b.shape[1]):
1108     out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
1109     elif b.rank==3:
1110     for i0 in range(b.shape[0]):
1111     for i1 in range(b.shape[1]):
1112     for i2 in range(b.shape[2]):
1113     out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
1114     elif b.rank==4:
1115     for i0 in range(b.shape[0]):
1116     for i1 in range(b.shape[1]):
1117     for i2 in range(b.shape[2]):
1118     for i3 in range(b.shape[3]):
1119     out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
1120     return out
1121    
1122     def unrollLoopsOfL2(b,where,arg,tap=""):
1123     out=""
1124     z=[]
1125     if isinstance(b,float) or b.rank==0:
1126     z.append(b**2)
1127     elif b.rank==1:
1128     for i0 in range(b.shape[0]):
1129     z.append(b[i0]**2)
1130     elif b.rank==2:
1131     for i1 in range(b.shape[1]):
1132     s=0
1133     for i0 in range(b.shape[0]):
1134     s+=b[i0,i1]**2
1135     z.append(s)
1136     elif b.rank==3:
1137     for i2 in range(b.shape[2]):
1138     s=0
1139     for i0 in range(b.shape[0]):
1140     for i1 in range(b.shape[1]):
1141     s+=b[i0,i1,i2]**2
1142     z.append(s)
1143    
1144     elif b.rank==4:
1145     for i3 in range(b.shape[3]):
1146     s=0
1147     for i0 in range(b.shape[0]):
1148     for i1 in range(b.shape[1]):
1149     for i2 in range(b.shape[2]):
1150     s+=b[i0,i1,i2,i3]**2
1151     z.append(s)
1152     if where=="Function":
1153     xfac_o=1.
1154     xfac_op=0.
1155     z_fac_s=""
1156     zo_fac_s=""
1157     zo_fac=1./3.
1158     elif where=="FunctionOnBoundary":
1159     xfac_o=1.
1160     xfac_op=0.
1161     z_fac_s="*dim"
1162     zo_fac_s="*(2.*dim+1.)/3."
1163     zo_fac=1.
1164     elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1165     xfac_o=0.
1166     xfac_op=1.
1167     z_fac_s=""
1168     zo_fac_s=""
1169     zo_fac=1./3.
1170     zo=0.
1171     zop=0.
1172     for i99 in range(len(z)):
1173     if i99==0:
1174     zo+=xfac_o*z[i99]
1175     zop+=xfac_op*z[i99]
1176     else:
1177     zo+=z[i99]
1178     out+=tap+"%s=sqrt((%s)%s"%(arg,zo*zo_fac,zo_fac_s)
1179     if zop==0.:
1180     out+=")\n"
1181     else:
1182     out+="+(%s))\n"%(zop*0.5**2)
1183     return out
1184 gross 493 #=======================================================================================================
1185     # transpose
1186     #=======================================================================================================
1187     def transposeTest(r,offset):
1188     if isinstance(r,float): return r
1189     s=r.shape
1190     s1=1
1191     for i in s[:offset]: s1*=i
1192     s2=1
1193     for i in s[offset:]: s2*=i
1194     out=numarray.reshape(r,(s1,s2))
1195     out.transpose()
1196     return numarray.resize(out,s[offset:]+s[:offset])
1197 gross 443
1198 gross 493 name,tt="transpose",transposeTest
1199     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1200     for sh0 in [ (), (3,), (4,5), (6,2,2),(3,2,3,4)]:
1201     for offset in range(len(sh0)+1):
1202     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1203     tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1204     text+=" def %s(self):\n"%tname
1205     sh_t=sh0[offset:]+sh0[:offset]
1206    
1207     # sh_t=list(sh0)
1208     # sh_t[offset+1]=sh_t[offset]
1209     # sh_t=tuple(sh_t)
1210     # sh_r=[]
1211     # for i in range(offset): sh_r.append(sh0[i])
1212     # for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1213     # sh_r=tuple(sh_r)
1214    
1215     a_0=makeArray(sh0,[-1.,1])
1216     if case0 in ["taggedData", "expandedData"]:
1217     a1_0=makeArray(sh0,[-1.,1])
1218     else:
1219     a1_0=a_0
1220     r=tt(a_0,offset)
1221     r1=tt(a1_0,offset)
1222     text+=mkText(case0,"arg",a_0,a1_0)
1223     text+=" res=%s(arg,%s)\n"%(name,offset)
1224     if case0=="Symbol":
1225     text+=mkText("array","s",a_0,a1_0)
1226     text+=" sub=res.substitute({arg:s})\n"
1227     res="sub"
1228     text+=mkText("array","ref",r,r1)
1229     else:
1230     res="res"
1231     text+=mkText(case0,"ref",r,r1)
1232     text+=mkTypeAndShapeTest(case0,sh_t,"res")
1233     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1234    
1235     if case0 == "taggedData" :
1236     t_prog_with_tags+=text
1237     else:
1238     t_prog+=text
1239    
1240     print test_header
1241     # print t_prog
1242     print t_prog_with_tags
1243     print test_tail
1244     1/0
1245 gross 441 #=======================================================================================================
1246 gross 443 # L2
1247     #=======================================================================================================
1248     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1249     for data in ["Data","Symbol"]:
1250     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1251     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1252     tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
1253     text+=" def %s(self):\n"%tname
1254     text+=" \"\"\"\n"
1255     text+=" tests L2-norm of %s on the %s\n\n"%(data,where)
1256     text+=" assumptions: self.domain supports integration on %s\n"%where
1257     text+=" \"\"\"\n"
1258     text+=" dim=self.domain.getDim()\n"
1259     text+=" w=%s(self.domain)\n"%where
1260     text+=" x=w.getX()\n"
1261     o="1"
1262     if len(sh)>0:
1263     sh_2=sh[:len(sh)-1]+(2,)
1264     sh_3=sh[:len(sh)-1]+(3,)
1265     b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
1266     b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
1267     else:
1268     sh_2=()
1269     sh_3=()
1270     b_2=makeArray(sh,[-1.,1])
1271     b_3=makeArray(sh,[-1.,1])
1272    
1273     if data=="Symbol":
1274     val="s"
1275     res="sub"
1276     else:
1277     val="arg"
1278     res="res"
1279     text+=" if dim==2:\n"
1280     if data=="Symbol":
1281     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
1282    
1283     text+=" %s=Data(0,%s,w)\n"%(val,sh_2)
1284     text+=unrollLoopsSimplified(b_2,val,tap=" ")
1285     text+=unrollLoopsOfL2(b_2,where,"ref",tap=" ")
1286     text+="\n else:\n"
1287     if data=="Symbol":
1288     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
1289     text+=" %s=Data(0,%s,w)\n"%(val,sh_3)
1290     text+=unrollLoopsSimplified(b_3,val,tap=" ")
1291     text+=unrollLoopsOfL2(b_3,where,"ref",tap=" ")
1292     text+="\n res=L2(arg)\n"
1293     if data=="Symbol":
1294     text+=" sub=res.substitute({arg:s})\n"
1295     text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1296     text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1297     else:
1298     text+=" self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1299     text+=" self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
1300     t_prog+=text
1301     print t_prog
1302     1/0
1303    
1304     #=======================================================================================================
1305 gross 441 # div
1306     #=======================================================================================================
1307     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1308     for data in ["Data","Symbol"]:
1309     for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1310     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1311     tname="test_div_on%s_from%s_%s"%(where,data,case)
1312     text+=" def %s(self):\n"%tname
1313     text+=" \"\"\"\n"
1314     text+=" tests divergence of %s on the %s\n\n"%(data,where)
1315     text+=" assumptions: %s(self.domain) exists\n"%case
1316     text+=" self.domain supports div on %s\n"%where
1317     text+=" \"\"\"\n"
1318     if case=="ReducedSolution":
1319     text+=" o=1\n"
1320     o="1"
1321     else:
1322     text+=" o=self.order\n"
1323     o="o"
1324     text+=" dim=self.domain.getDim()\n"
1325     text+=" w_ref=%s(self.domain)\n"%where
1326     text+=" x_ref=w_ref.getX()\n"
1327     text+=" w=%s(self.domain)\n"%case
1328     text+=" x=w.getX()\n"
1329     a_2=makeArray((2,2),[-1.,1])
1330     b_2=makeArray((2,2),[-1.,1])
1331     a_3=makeArray((3,3),[-1.,1])
1332     b_3=makeArray((3,3),[-1.,1])
1333     if data=="Symbol":
1334     text+=" arg=Symbol(shape=(dim,),dim=dim)\n"
1335     val="s"
1336     res="sub"
1337     else:
1338     val="arg"
1339     res="res"
1340     text+=" %s=Vector(0,w)\n"%val
1341     text+=" if dim==2:\n"
1342     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1343     text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap=" ")
1344     text+="\n else:\n"
1345    
1346     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1347     text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap=" ")
1348     text+="\n res=div(arg,where=w_ref)\n"
1349     if data=="Symbol":
1350     text+=" sub=res.substitute({arg:s})\n"
1351     text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1352     text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1353     text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
1354     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1355 gross 437
1356 gross 441
1357     t_prog+=text
1358     print t_prog
1359     1/0
1360    
1361 gross 439 #=======================================================================================================
1362     # interpolation
1363     #=======================================================================================================
1364     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
1365     for data in ["Data","Symbol"]:
1366     for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1367     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1368     if where==case or \
1369     ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
1370     ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
1371     (case=="ContinuousFunction" and where in ["Solution","ReducedSolution"]) or \
1372     (case=="Solution" and where=="ReducedSolution") :
1373    
1374    
1375     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1376     tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1377     text+=" def %s(self):\n"%tname
1378     text+=" \"\"\"\n"
1379     text+=" tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
1380     text+=" assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
1381     text+=" \"\"\"\n"
1382     if case=="ReducedSolution" or where=="ReducedSolution":
1383     text+=" o=1\n"
1384     o="1"
1385     else:
1386     text+=" o=self.order\n"
1387     o="o"
1388     text+=" dim=self.domain.getDim()\n"
1389     text+=" w_ref=%s(self.domain)\n"%where
1390     text+=" x_ref=w_ref.getX()\n"
1391     text+=" w=%s(self.domain)\n"%case
1392     text+=" x=w.getX()\n"
1393     a_2=makeArray(sh+(2,),[-1.,1])
1394     b_2=makeArray(sh+(2,),[-1.,1])
1395     a_3=makeArray(sh+(3,),[-1.,1])
1396     b_3=makeArray(sh+(3,),[-1.,1])
1397     if data=="Symbol":
1398     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1399     val="s"
1400     res="sub"
1401     else:
1402     val="arg"
1403     res="res"
1404     text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1405     text+=" ref=Data(0,%s,w_ref)\n"%str(sh)
1406     text+=" if dim==2:\n"
1407     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1408     text+=unrollLoops(a_2,b_2,o,"ref",tap=" ",x="x_ref")
1409     text+=" else:\n"
1410    
1411     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1412     text+=unrollLoops(a_3,b_3,o,"ref",tap=" ",x="x_ref")
1413     text+=" res=interpolate(arg,where=w_ref)\n"
1414     if data=="Symbol":
1415     text+=" sub=res.substitute({arg:s})\n"
1416     text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1417     text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
1418     text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1419     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1420     t_prog+=text
1421     print test_header
1422     print t_prog
1423     print test_tail
1424     1/0
1425 gross 437
1426 gross 429 #=======================================================================================================
1427 gross 438 # grad
1428 gross 437 #=======================================================================================================
1429     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1430     for data in ["Data","Symbol"]:
1431 gross 438 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1432     for sh in [ (),(2,), (4,5), (6,2,2)]:
1433 gross 437 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1434 gross 438 tname="test_grad_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1435 gross 437 text+=" def %s(self):\n"%tname
1436     text+=" \"\"\"\n"
1437 gross 438 text+=" tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
1438 gross 437 text+=" assumptions: %s(self.domain) exists\n"%case
1439 gross 438 text+=" self.domain supports gardient on %s\n"%where
1440 gross 437 text+=" \"\"\"\n"
1441     if case=="ReducedSolution":
1442     text+=" o=1\n"
1443     o="1"
1444     else:
1445     text+=" o=self.order\n"
1446     o="o"
1447     text+=" dim=self.domain.getDim()\n"
1448     text+=" w_ref=%s(self.domain)\n"%where
1449 gross 438 text+=" x_ref=w_ref.getX()\n"
1450 gross 437 text+=" w=%s(self.domain)\n"%case
1451     text+=" x=w.getX()\n"
1452     a_2=makeArray(sh+(2,),[-1.,1])
1453     b_2=makeArray(sh+(2,),[-1.,1])
1454     a_3=makeArray(sh+(3,),[-1.,1])
1455     b_3=makeArray(sh+(3,),[-1.,1])
1456     if data=="Symbol":
1457 gross 438 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1458 gross 437 val="s"
1459     res="sub"
1460     else:
1461     val="arg"
1462     res="res"
1463     text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1464 gross 438 text+=" ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1465 gross 437 text+=" if dim==2:\n"
1466     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1467 gross 438 text+=unrollLoopsOfGrad(a_2,b_2,o,"ref",tap=" ")
1468 gross 437 text+=" else:\n"
1469    
1470     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1471 gross 438 text+=unrollLoopsOfGrad(a_3,b_3,o,"ref",tap=" ")
1472     text+=" res=grad(arg,where=w_ref)\n"
1473 gross 437 if data=="Symbol":
1474     text+=" sub=res.substitute({arg:s})\n"
1475 gross 438 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1476     text+=" self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1477 gross 437 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1478    
1479    
1480     t_prog+=text
1481     print test_header
1482     print t_prog
1483     print test_tail
1484     1/0
1485 gross 438
1486    
1487 gross 437 #=======================================================================================================
1488 gross 438 # integrate
1489 gross 437 #=======================================================================================================
1490 gross 438 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1491 gross 437 for data in ["Data","Symbol"]:
1492 gross 438 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1493     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1494     if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:
1495 gross 437 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1496 gross 438 tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1497 gross 437 text+=" def %s(self):\n"%tname
1498     text+=" \"\"\"\n"
1499 gross 438 text+=" tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
1500     text+=" assumptions: %s(self.domain) exists\n"%case
1501     text+=" self.domain supports integral on %s\n"%where
1502    
1503 gross 437 text+=" \"\"\"\n"
1504     if case=="ReducedSolution":
1505     text+=" o=1\n"
1506     o="1"
1507     else:
1508     text+=" o=self.order\n"
1509     o="o"
1510     text+=" dim=self.domain.getDim()\n"
1511     text+=" w_ref=%s(self.domain)\n"%where
1512     text+=" w=%s(self.domain)\n"%case
1513     text+=" x=w.getX()\n"
1514     a_2=makeArray(sh+(2,),[-1.,1])
1515     b_2=makeArray(sh+(2,),[-1.,1])
1516     a_3=makeArray(sh+(3,),[-1.,1])
1517     b_3=makeArray(sh+(3,),[-1.,1])
1518     if data=="Symbol":
1519 gross 438 text+=" arg=Symbol(shape=%s)\n"%str(sh)
1520 gross 437 val="s"
1521     res="sub"
1522     else:
1523     val="arg"
1524     res="res"
1525 gross 438
1526 gross 437 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1527 gross 438 if not len(sh)==0:
1528     text+=" ref=numarray.zeros(%s,numarray.Float)\n"%str(sh)
1529 gross 437 text+=" if dim==2:\n"
1530     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1531 gross 438 text+=unrollLoopsOfInteriorIntegral(a_2,b_2,where,"ref",tap=" ")
1532 gross 437 text+=" else:\n"
1533    
1534     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1535 gross 438 text+=unrollLoopsOfInteriorIntegral(a_3,b_3,where,"ref",tap=" ")
1536     if case in ["ContinuousFunction","Solution","ReducedSolution"]:
1537     text+=" res=integrate(arg,where=w_ref)\n"
1538     else:
1539     text+=" res=integrate(arg)\n"
1540    
1541 gross 437 if data=="Symbol":
1542     text+=" sub=res.substitute({arg:s})\n"
1543 gross 438 if len(sh)==0 and data=="Data":
1544     text+=" self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%res
1545     else:
1546     if data=="Symbol":
1547     text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1548     text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1549     else:
1550     text+=" self.failUnless(isinstance(res,numarray.NumArray),\"wrong type of result.\")\n"
1551     text+=" self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
1552 gross 437 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1553    
1554    
1555     t_prog+=text
1556     print test_header
1557     print t_prog
1558     print test_tail
1559     1/0
1560     #=======================================================================================================
1561 gross 433 # inverse
1562     #=======================================================================================================
1563     name="inverse"
1564     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1565     for sh0 in [ (1,1), (2,2), (3,3)]:
1566     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1567     tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1568     text+=" def %s(self):\n"%tname
1569     a_0=makeArray(sh0,[-1.,1])
1570     for i in range(sh0[0]): a_0[i,i]+=2
1571     if case0 in ["taggedData", "expandedData"]:
1572     a1_0=makeArray(sh0,[-1.,1])
1573     for i in range(sh0[0]): a1_0[i,i]+=3
1574     else:
1575     a1_0=a_0
1576    
1577     text+=mkText(case0,"arg",a_0,a1_0)
1578     text+=" res=%s(arg)\n"%name
1579     if case0=="Symbol":
1580     text+=mkText("array","s",a_0,a1_0)
1581     text+=" sub=res.substitute({arg:s})\n"
1582     res="sub"
1583     ref="s"
1584     else:
1585     ref="arg"
1586     res="res"
1587     text+=mkTypeAndShapeTest(case0,sh0,"res")
1588     text+=" self.failUnless(Lsup(matrixmult(%s,%s)-kronecker(%s))<=self.RES_TOL,\"wrong result\")\n"%(res,ref,sh0[0])
1589    
1590     if case0 == "taggedData" :
1591     t_prog_with_tags+=text
1592     else:
1593     t_prog+=text
1594    
1595     print test_header
1596     # print t_prog
1597     print t_prog_with_tags
1598     print test_tail
1599     1/0
1600    
1601     #=======================================================================================================
1602 gross 429 # trace
1603     #=======================================================================================================
1604     def traceTest(r,offset):
1605     sh=r.shape
1606     r1=1
1607     for i in range(offset): r1*=sh[i]
1608     r2=1
1609     for i in range(offset+2,len(sh)): r2*=sh[i]
1610     r_s=numarray.reshape(r,(r1,sh[offset],sh[offset],r2))
1611     s=numarray.zeros([r1,r2],numarray.Float)
1612     for i1 in range(r1):
1613     for i2 in range(r2):
1614     for j in range(sh[offset]): s[i1,i2]+=r_s[i1,j,j,i2]
1615     return s.resize(sh[:offset]+sh[offset+2:])
1616     name,tt="trace",traceTest
1617     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1618     for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
1619     for offset in range(len(sh0)-1):
1620     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1621     tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1622     text+=" def %s(self):\n"%tname
1623     sh_t=list(sh0)
1624     sh_t[offset+1]=sh_t[offset]
1625     sh_t=tuple(sh_t)
1626     sh_r=[]
1627     for i in range(offset): sh_r.append(sh0[i])
1628     for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1629     sh_r=tuple(sh_r)
1630     a_0=makeArray(sh_t,[-1.,1])
1631     if case0 in ["taggedData", "expandedData"]:
1632     a1_0=makeArray(sh_t,[-1.,1])
1633     else:
1634     a1_0=a_0
1635     r=tt(a_0,offset)
1636     r1=tt(a1_0,offset)
1637     text+=mkText(case0,"arg",a_0,a1_0)
1638     text+=" res=%s(arg,%s)\n"%(name,offset)
1639     if case0=="Symbol":
1640     text+=mkText("array","s",a_0,a1_0)
1641     text+=" sub=res.substitute({arg:s})\n"
1642     res="sub"
1643     text+=mkText("array","ref",r,r1)
1644     else:
1645     res="res"
1646     text+=mkText(case0,"ref",r,r1)
1647     text+=mkTypeAndShapeTest(case0,sh_r,"res")
1648     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1649    
1650     if case0 == "taggedData" :
1651     t_prog_with_tags+=text
1652     else:
1653     t_prog+=text
1654 gross 396
1655 gross 429 print test_header
1656     # print t_prog
1657     print t_prog_with_tags
1658     print test_tail
1659     1/0
1660 gross 396
1661 gross 157 #=======================================================================================================
1662 gross 396 # clip
1663     #=======================================================================================================
1664     oper_L=[["clip",clipTEST]]
1665     for oper in oper_L:
1666     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1667     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1668     if len(sh0)==0 or not case0=="float":
1669     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1670     tname="test_%s_%s_rank%s"%(oper[0],case0,len(sh0))
1671     text+=" def %s(self):\n"%tname
1672     a_0=makeArray(sh0,[-1.,1])
1673     if case0 in ["taggedData", "expandedData"]:
1674     a1_0=makeArray(sh0,[-1.,1])
1675     else:
1676     a1_0=a_0
1677    
1678     r=oper[1](a_0,-0.3,0.5)
1679     r1=oper[1](a1_0,-0.3,0.5)
1680     text+=mkText(case0,"arg",a_0,a1_0)
1681     text+=" res=%s(arg,-0.3,0.5)\n"%oper[0]
1682     if case0=="Symbol":
1683     text+=mkText("array","s",a_0,a1_0)
1684     text+=" sub=res.substitute({arg:s})\n"
1685     res="sub"
1686     text+=mkText("array","ref",r,r1)
1687     else:
1688     res="res"
1689     text+=mkText(case0,"ref",r,r1)
1690     text+=mkTypeAndShapeTest(case0,sh0,"res")
1691     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1692    
1693     if case0 == "taggedData" :
1694     t_prog_with_tags+=text
1695     else:
1696     t_prog+=text
1697    
1698     print test_header
1699     # print t_prog
1700     print t_prog_with_tags
1701     print test_tail
1702     1/0
1703 gross 429
1704 gross 396 #=======================================================================================================
1705     # maximum, minimum, clipping
1706     #=======================================================================================================
1707     oper_L=[ ["maximum",maximumTEST],
1708     ["minimum",minimumTEST]]
1709     for oper in oper_L:
1710     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1711     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1712     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1713     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1714     if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1715     and (sh0==sh1 or len(sh0)==0 or len(sh1)==0) :
1716     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1717    
1718     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1719     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1720     text+=" def %s(self):\n"%tname
1721     a_0=makeArray(sh0,[-1.,1])
1722     if case0 in ["taggedData", "expandedData"]:
1723     a1_0=makeArray(sh0,[-1.,1])
1724     else:
1725     a1_0=a_0
1726    
1727     a_1=makeArray(sh1,[-1.,1])
1728     if case1 in ["taggedData", "expandedData"]:
1729     a1_1=makeArray(sh1,[-1.,1])
1730     else:
1731     a1_1=a_1
1732     r=oper[1](a_0,a_1)
1733     r1=oper[1](a1_0,a1_1)
1734     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1735     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1736     text+=" res=%s(arg0,arg1)\n"%oper[0]
1737     case=getResultCaseForBin(case0,case1)
1738     if case=="Symbol":
1739     c0_res,c1_res=case0,case1
1740     subs="{"
1741     if case0=="Symbol":
1742     text+=mkText("array","s0",a_0,a1_0)
1743     subs+="arg0:s0"
1744     c0_res="array"
1745     if case1=="Symbol":
1746     text+=mkText("array","s1",a_1,a1_1)
1747     if not subs.endswith("{"): subs+=","
1748     subs+="arg1:s1"
1749     c1_res="array"
1750     subs+="}"
1751     text+=" sub=res.substitute(%s)\n"%subs
1752     res="sub"
1753     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1754     else:
1755     res="res"
1756     text+=mkText(case,"ref",r,r1)
1757     if len(sh0)>len(sh1):
1758     text+=mkTypeAndShapeTest(case,sh0,"res")
1759     else:
1760     text+=mkTypeAndShapeTest(case,sh1,"res")
1761     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1762    
1763     if case0 == "taggedData" or case1 == "taggedData":
1764     t_prog_with_tags+=text
1765     else:
1766     t_prog+=text
1767    
1768     print test_header
1769     # print t_prog
1770     print t_prog_with_tags
1771     print test_tail
1772     1/0
1773    
1774    
1775     #=======================================================================================================
1776     # outer inner
1777     #=======================================================================================================
1778     oper=["outer",outerTEST]
1779     # oper=["inner",innerTEST]
1780     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1781     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1782     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1783     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1784     if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1785     and len(sh0+sh1)<5:
1786     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1787    
1788     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1789     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1790     text+=" def %s(self):\n"%tname
1791     a_0=makeArray(sh0,[-1.,1])
1792     if case0 in ["taggedData", "expandedData"]:
1793     a1_0=makeArray(sh0,[-1.,1])
1794     else:
1795     a1_0=a_0
1796    
1797     a_1=makeArray(sh1,[-1.,1])
1798     if case1 in ["taggedData", "expandedData"]:
1799     a1_1=makeArray(sh1,[-1.,1])
1800     else:
1801     a1_1=a_1
1802     r=oper[1](a_0,a_1)
1803     r1=oper[1](a1_0,a1_1)
1804     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1805     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1806     text+=" res=%s(arg0,arg1)\n"%oper[0]
1807     case=getResultCaseForBin(case0,case1)
1808     if case=="Symbol":
1809     c0_res,c1_res=case0,case1
1810     subs="{"
1811     if case0=="Symbol":
1812     text+=mkText("array","s0",a_0,a1_0)
1813     subs+="arg0:s0"
1814     c0_res="array"
1815     if case1=="Symbol":
1816     text+=mkText("array","s1",a_1,a1_1)
1817     if not subs.endswith("{"): subs+=","
1818     subs+="arg1:s1"
1819     c1_res="array"
1820     subs+="}"
1821     text+=" sub=res.substitute(%s)\n"%subs
1822     res="sub"
1823     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1824     else:
1825     res="res"
1826     text+=mkText(case,"ref",r,r1)
1827     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1828     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1829    
1830     if case0 == "taggedData" or case1 == "taggedData":
1831     t_prog_with_tags+=text
1832     else:
1833     t_prog+=text
1834    
1835     print test_header
1836     # print t_prog
1837     print t_prog_with_tags
1838     print test_tail
1839     1/0
1840    
1841     #=======================================================================================================
1842 gross 313 # local reduction
1843     #=======================================================================================================
1844     for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],
1845     ["maxval",-1.e99,"max(out,%a1%)","out"],
1846     ["minval",1.e99,"min(out,%a1%)","out"] ]:
1847     for case in case_set:
1848     for sh in shape_set:
1849     if not case=="float" or len(sh)==0:
1850     text=""
1851     text+=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1852     tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))
1853     text+=" %s(self):\n"%tname
1854     a=makeArray(sh,[-1.,1.])
1855     a1=makeArray(sh,[-1.,1.])
1856     r1=testReduce(a1,oper[1],oper[2],oper[3])
1857     r=testReduce(a,oper[1],oper[2],oper[3])
1858    
1859     text+=mkText(case,"arg",a,a1)
1860     text+=" res=%s(arg)\n"%oper[0]
1861     if case=="Symbol":
1862     text+=mkText("array","s",a,a1)
1863