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

Annotation of /trunk/escript/py_src/generatediff

Parent Directory Parent Directory | Revision Log Revision Log


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