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

Annotation of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


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