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

Annotation of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


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