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

Annotation of /trunk/escript/py_src/generatediff

Parent Directory Parent Directory | Revision Log Revision Log


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