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

Annotation of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1387 - (hide annotations)
Fri Jan 11 07:45:26 2008 UTC (11 years, 9 months ago) by trankine
Original Path: temp/escript/py_src/generateutil
File size: 168742 byte(s)
Restore the trunk that existed before the windows changes were committed to the (now moved to branches) old trunk.
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 gross 785 if isinstance(l,int):
157     if len(shape)==0:
158     out=int(l*random.random()+rng[0])*1.
159     elif len(shape)==1:
160     for i0 in range(shape[0]):
161     out[i0]=int(l*random.random()+rng[0])
162     elif len(shape)==2:
163     for i0 in range(shape[0]):
164     for i1 in range(shape[1]):
165     out[i0,i1]=int(l*random.random()+rng[0])
166     elif len(shape)==3:
167     for i0 in range(shape[0]):
168     for i1 in range(shape[1]):
169     for i2 in range(shape[2]):
170     out[i0,i1,i2]=int(l*random.random()+rng[0])
171     elif len(shape)==4:
172     for i0 in range(shape[0]):
173     for i1 in range(shape[1]):
174     for i2 in range(shape[2]):
175     for i3 in range(shape[3]):
176     out[i0,i1,i2,i3]=int(l*random.random()+rng[0])
177     elif len(shape)==5:
178     for i0 in range(shape[0]):
179     for i1 in range(shape[1]):
180     for i2 in range(shape[2]):
181     for i3 in range(shape[3]):
182     for i4 in range(shape[4]):
183     out[i0,i1,i2,i3,i4]=int(l*ranm.random()+rng[0])
184     else:
185     raise SystemError,"rank is restricted to 5"
186 jgs 154 else:
187 gross 785 if len(shape)==0:
188     out=l*random.random()+rng[0]
189     elif len(shape)==1:
190     for i0 in range(shape[0]):
191     out[i0]=l*random.random()+rng[0]
192     elif len(shape)==2:
193     for i0 in range(shape[0]):
194     for i1 in range(shape[1]):
195     out[i0,i1]=l*random.random()+rng[0]
196     elif len(shape)==3:
197     for i0 in range(shape[0]):
198     for i1 in range(shape[1]):
199     for i2 in range(shape[2]):
200     out[i0,i1,i2]=l*random.random()+rng[0]
201     elif len(shape)==4:
202     for i0 in range(shape[0]):
203     for i1 in range(shape[1]):
204     for i2 in range(shape[2]):
205     for i3 in range(shape[3]):
206     out[i0,i1,i2,i3]=l*random.random()+rng[0]
207     elif len(shape)==5:
208     for i0 in range(shape[0]):
209     for i1 in range(shape[1]):
210     for i2 in range(shape[2]):
211     for i3 in range(shape[3]):
212     for i4 in range(shape[4]):
213     out[i0,i1,i2,i3,i4]=l*ranm.random()+rng[0]
214     else:
215     raise SystemError,"rank is restricted to 5"
216 jgs 154 return out
217    
218 gross 517 def makeNumberedArray(shape,s=1.):
219     out=numarray.zeros(shape,numarray.Float64)
220     if len(shape)==0:
221     out=s*1.
222     elif len(shape)==1:
223     for i0 in range(shape[0]):
224 gross 550 out[i0]=s*int(8*random.random()+1)
225 gross 517 elif len(shape)==2:
226     for i0 in range(shape[0]):
227     for i1 in range(shape[1]):
228 gross 550 out[i0,i1]=s*int(8*random.random()+1)
229 gross 517 elif len(shape)==3:
230     for i0 in range(shape[0]):
231     for i1 in range(shape[1]):
232     for i2 in range(shape[2]):
233 gross 550 out[i0,i1,i2]=s*int(8*random.random()+1)
234 gross 517 elif len(shape)==4:
235     for i0 in range(shape[0]):
236     for i1 in range(shape[1]):
237     for i2 in range(shape[2]):
238     for i3 in range(shape[3]):
239 gross 550 out[i0,i1,i2,i3]=s*int(8*random.random()+1)
240 gross 517 else:
241     raise SystemError,"rank is restricted to 4"
242     return out
243 jgs 154
244 gross 157 def makeResult(val,test_expr):
245 jgs 154 if isinstance(val,float):
246 gross 157 out=eval(test_expr.replace("%a1%","val"))
247 jgs 154 elif len(val.shape)==0:
248 gross 157 out=eval(test_expr.replace("%a1%","val"))
249 jgs 154 elif len(val.shape)==1:
250     out=numarray.zeros(val.shape,numarray.Float64)
251     for i0 in range(val.shape[0]):
252 gross 157 out[i0]=eval(test_expr.replace("%a1%","val[i0]"))
253 jgs 154 elif len(val.shape)==2:
254     out=numarray.zeros(val.shape,numarray.Float64)
255     for i0 in range(val.shape[0]):
256     for i1 in range(val.shape[1]):
257 gross 157 out[i0,i1]=eval(test_expr.replace("%a1%","val[i0,i1]"))
258 jgs 154 elif len(val.shape)==3:
259     out=numarray.zeros(val.shape,numarray.Float64)
260     for i0 in range(val.shape[0]):
261     for i1 in range(val.shape[1]):
262     for i2 in range(val.shape[2]):
263 gross 157 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val[i0,i1,i2]"))
264 jgs 154 elif len(val.shape)==4:
265     out=numarray.zeros(val.shape,numarray.Float64)
266     for i0 in range(val.shape[0]):
267     for i1 in range(val.shape[1]):
268     for i2 in range(val.shape[2]):
269     for i3 in range(val.shape[3]):
270 gross 157 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val[i0,i1,i2,i3]"))
271 jgs 154 else:
272     raise SystemError,"rank is restricted to 4"
273     return out
274    
275 gross 157 def makeResult2(val0,val1,test_expr):
276     if isinstance(val0,float):
277     if isinstance(val1,float):
278     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
279     elif len(val1.shape)==0:
280     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
281     elif len(val1.shape)==1:
282     out=numarray.zeros(val1.shape,numarray.Float64)
283     for i0 in range(val1.shape[0]):
284     out[i0]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0]"))
285     elif len(val1.shape)==2:
286     out=numarray.zeros(val1.shape,numarray.Float64)
287     for i0 in range(val1.shape[0]):
288     for i1 in range(val1.shape[1]):
289     out[i0,i1]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1]"))
290     elif len(val1.shape)==3:
291     out=numarray.zeros(val1.shape,numarray.Float64)
292     for i0 in range(val1.shape[0]):
293     for i1 in range(val1.shape[1]):
294     for i2 in range(val1.shape[2]):
295     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2]"))
296     elif len(val1.shape)==4:
297     out=numarray.zeros(val1.shape,numarray.Float64)
298     for i0 in range(val1.shape[0]):
299     for i1 in range(val1.shape[1]):
300     for i2 in range(val1.shape[2]):
301     for i3 in range(val1.shape[3]):
302     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2,i3]"))
303     else:
304     raise SystemError,"rank of val1 is restricted to 4"
305     elif len(val0.shape)==0:
306     if isinstance(val1,float):
307     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
308     elif len(val1.shape)==0:
309     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
310     elif len(val1.shape)==1:
311     out=numarray.zeros(val1.shape,numarray.Float64)
312     for i0 in range(val1.shape[0]):
313     out[i0]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0]"))
314     elif len(val1.shape)==2:
315     out=numarray.zeros(val1.shape,numarray.Float64)
316     for i0 in range(val1.shape[0]):
317     for i1 in range(val1.shape[1]):
318     out[i0,i1]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1]"))
319     elif len(val1.shape)==3:
320     out=numarray.zeros(val1.shape,numarray.Float64)
321     for i0 in range(val1.shape[0]):
322     for i1 in range(val1.shape[1]):
323     for i2 in range(val1.shape[2]):
324     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2]"))
325     elif len(val1.shape)==4:
326     out=numarray.zeros(val1.shape,numarray.Float64)
327     for i0 in range(val1.shape[0]):
328     for i1 in range(val1.shape[1]):
329     for i2 in range(val1.shape[2]):
330     for i3 in range(val1.shape[3]):
331     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2,i3]"))
332     else:
333     raise SystemError,"rank of val1 is restricted to 4"
334     elif len(val0.shape)==1:
335     if isinstance(val1,float):
336     out=numarray.zeros(val0.shape,numarray.Float64)
337     for i0 in range(val0.shape[0]):
338     out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1"))
339     elif len(val1.shape)==0:
340     out=numarray.zeros(val0.shape,numarray.Float64)
341     for i0 in range(val0.shape[0]):
342     out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1"))
343     elif len(val1.shape)==1:
344     out=numarray.zeros(val0.shape,numarray.Float64)
345     for i0 in range(val0.shape[0]):
346     out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0]"))
347     elif len(val1.shape)==2:
348     out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
349     for i0 in range(val0.shape[0]):
350     for j1 in range(val1.shape[1]):
351     out[i0,j1]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1]"))
352     elif len(val1.shape)==3:
353     out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
354     for i0 in range(val0.shape[0]):
355     for j1 in range(val1.shape[1]):
356     for j2 in range(val1.shape[2]):
357     out[i0,j1,j2]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1,j2]"))
358     elif len(val1.shape)==4:
359     out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
360     for i0 in range(val0.shape[0]):
361     for j1 in range(val1.shape[1]):
362     for j2 in range(val1.shape[2]):
363     for j3 in range(val1.shape[3]):
364     out[i0,j1,j2,j3]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1,j2,j3]"))
365     else:
366     raise SystemError,"rank of val1 is restricted to 4"
367     elif len(val0.shape)==2:
368     if isinstance(val1,float):
369     out=numarray.zeros(val0.shape,numarray.Float64)
370     for i0 in range(val0.shape[0]):
371     for i1 in range(val0.shape[1]):
372     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1"))
373     elif len(val1.shape)==0:
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     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1"))
378     elif len(val1.shape)==1:
379     out=numarray.zeros(val0.shape,numarray.Float64)
380     for i0 in range(val0.shape[0]):
381     for i1 in range(val0.shape[1]):
382     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0]"))
383     elif len(val1.shape)==2:
384     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
385     for i0 in range(val0.shape[0]):
386     for i1 in range(val0.shape[1]):
387     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1]"))
388     elif len(val1.shape)==3:
389     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
390     for i0 in range(val0.shape[0]):
391     for i1 in range(val0.shape[1]):
392     for j2 in range(val1.shape[2]):
393     out[i0,i1,j2]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1,j2]"))
394     elif len(val1.shape)==4:
395     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
396     for i0 in range(val0.shape[0]):
397     for i1 in range(val0.shape[1]):
398     for j2 in range(val1.shape[2]):
399     for j3 in range(val1.shape[3]):
400     out[i0,i1,j2,j3]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1,j2,j3]"))
401     else:
402     raise SystemError,"rank of val1 is restricted to 4"
403     elif len(val0.shape)==3:
404     if isinstance(val1,float):
405     out=numarray.zeros(val0.shape,numarray.Float64)
406     for i0 in range(val0.shape[0]):
407     for i1 in range(val0.shape[1]):
408     for i2 in range(val0.shape[2]):
409     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1"))
410     elif len(val1.shape)==0:
411     out=numarray.zeros(val0.shape,numarray.Float64)
412     for i0 in range(val0.shape[0]):
413     for i1 in range(val0.shape[1]):
414     for i2 in range(val0.shape[2]):
415     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1"))
416     elif len(val1.shape)==1:
417     out=numarray.zeros(val0.shape,numarray.Float64)
418     for i0 in range(val0.shape[0]):
419     for i1 in range(val0.shape[1]):
420     for i2 in range(val0.shape[2]):
421     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0]"))
422     elif len(val1.shape)==2:
423     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
424     for i0 in range(val0.shape[0]):
425     for i1 in range(val0.shape[1]):
426     for i2 in range(val0.shape[2]):
427     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1]"))
428     elif len(val1.shape)==3:
429     out=numarray.zeros(val0.shape,numarray.Float64)
430     for i0 in range(val0.shape[0]):
431     for i1 in range(val0.shape[1]):
432     for i2 in range(val0.shape[2]):
433     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1,i2]"))
434     elif len(val1.shape)==4:
435     out=numarray.zeros(val0.shape+val1.shape[3:],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 j3 in range(val1.shape[3]):
440     out[i0,i1,i2,j3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1,i2,j3]"))
441     else:
442     raise SystemError,"rank of val1 is rargs[1]estricted to 4"
443     elif len(val0.shape)==4:
444     if isinstance(val1,float):
445     out=numarray.zeros(val0.shape,numarray.Float64)
446     for i0 in range(val0.shape[0]):
447     for i1 in range(val0.shape[1]):
448     for i2 in range(val0.shape[2]):
449     for i3 in range(val0.shape[3]):
450     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1"))
451     elif len(val1.shape)==0:
452     out=numarray.zeros(val0.shape,numarray.Float64)
453     for i0 in range(val0.shape[0]):
454     for i1 in range(val0.shape[1]):
455     for i2 in range(val0.shape[2]):
456     for i3 in range(val0.shape[3]):
457     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1"))
458     elif len(val1.shape)==1:
459     out=numarray.zeros(val0.shape,numarray.Float64)
460     for i0 in range(val0.shape[0]):
461     for i1 in range(val0.shape[1]):
462     for i2 in range(val0.shape[2]):
463     for i3 in range(val0.shape[3]):
464     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0]"))
465     elif len(val1.shape)==2:
466     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
467     for i0 in range(val0.shape[0]):
468     for i1 in range(val0.shape[1]):
469     for i2 in range(val0.shape[2]):
470     for i3 in range(val0.shape[3]):
471     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0,i1]"))
472     elif len(val1.shape)==3:
473     out=numarray.zeros(val0.shape,numarray.Float64)
474     for i0 in range(val0.shape[0]):
475     for i1 in range(val0.shape[1]):
476     for i2 in range(val0.shape[2]):
477     for i3 in range(val0.shape[3]):
478     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0,i1,i2]"))
479     elif len(val1.shape)==4:
480     out=numarray.zeros(val0.shape,numarray.Float64)
481     for i0 in range(val0.shape[0]):
482     for i1 in range(val0.shape[1]):
483     for i2 in range(val0.shape[2]):
484     for i3 in range(val0.shape[3]):
485     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0,i1,i2,i3]"))
486     else:
487     raise SystemError,"rank of val1 is restricted to 4"
488     else:
489     raise SystemError,"rank is restricted to 4"
490     return out
491 jgs 154
492 gross 157
493     def mkText(case,name,a,a1=None,use_tagging_for_expanded_data=False):
494 jgs 154 t_out=""
495 gross 157 if case=="float":
496 jgs 154 if isinstance(a,float):
497     t_out+=" %s=%s\n"%(name,a)
498 gross 291 elif a.rank==0:
499 jgs 154 t_out+=" %s=%s\n"%(name,a)
500     else:
501     t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
502 gross 157 elif case=="array":
503     if isinstance(a,float):
504     t_out+=" %s=numarray.array(%s)\n"%(name,a)
505 gross 291 elif a.rank==0:
506 gross 157 t_out+=" %s=numarray.array(%s)\n"%(name,a)
507     else:
508     t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
509 jgs 154 elif case=="constData":
510     if isinstance(a,float):
511     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
512 gross 291 elif a.rank==0:
513 jgs 154 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
514     else:
515     t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
516     elif case=="taggedData":
517     if isinstance(a,float):
518     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
519     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
520 gross 291 elif a.rank==0:
521 jgs 154 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
522     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
523     else:
524     t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
525 gross 157 t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
526 jgs 154 elif case=="expandedData":
527 gross 157 if use_tagging_for_expanded_data:
528     if isinstance(a,float):
529     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
530     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
531 gross 291 elif a.rank==0:
532 gross 157 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
533     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
534     else:
535     t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
536     t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
537     t_out+=" %s.expand()\n"%name
538     else:
539 gross 785 t_out+=" msk_%s=1-whereZero(self.functionspace.getX()[0],1.e-8)\n"%name
540 gross 157 if isinstance(a,float):
541 gross 785 t_out+=" %s=msk_%s*(%s)+(1-msk_%s)*(%s)\n"%(name,name,a,name,a1)
542 gross 291 elif a.rank==0:
543 gross 157 t_out+=" %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a,name,a1)
544     else:
545     t_out+=" %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a.tolist(),name,a1.tolist())
546 jgs 154 elif case=="Symbol":
547     if isinstance(a,float):
548     t_out+=" %s=Symbol(shape=())\n"%(name)
549 gross 291 elif a.rank==0:
550 jgs 154 t_out+=" %s=Symbol(shape=())\n"%(name)
551     else:
552     t_out+=" %s=Symbol(shape=%s)\n"%(name,str(a.shape))
553    
554     return t_out
555    
556 gross 587 def mkTypeAndShapeTest(case,sh,argstr,name=""):
557 gross 157 text=""
558     if case=="float":
559 gross 587 text+=" self.failUnless(isinstance(%s,float),\"wrong type of result%s.\")\n"%(argstr,name)
560 gross 157 elif case=="array":
561 gross 587 text+=" self.failUnless(isinstance(%s,numarray.NumArray),\"wrong type of result%s.\")\n"%(argstr,name)
562     text+=" self.failUnlessEqual(%s.shape,%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
563 gross 157 elif case in ["constData","taggedData","expandedData"]:
564 gross 587 text+=" self.failUnless(isinstance(%s,Data),\"wrong type of result%s.\")\n"%(argstr,name)
565     text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
566 gross 157 elif case=="Symbol":
567 gross 587 text+=" self.failUnless(isinstance(%s,Symbol),\"wrong type of result%s.\")\n"%(argstr,name)
568     text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
569 gross 157 return text
570 jgs 154
571 gross 157 def mkCode(txt,args=[],intend=""):
572     s=txt.split("\n")
573     if len(s)>1:
574     out=""
575     for l in s:
576     out+=intend+l+"\n"
577     else:
578     out="%sreturn %s\n"%(intend,txt)
579     c=1
580     for r in args:
581     out=out.replace("%%a%s%%"%c,r)
582     return out
583 gross 804
584    
585     #=======================================================================================================
586     # swap axes
587     #=======================================================================================================
588     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
589     for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
590     for axis0 in range(len(sh0)):
591     for axis1 in range(axis0+1,len(sh0)):
592     tname="test_%s_%s_rank%s_axes_%s%s"%("swapaxes",case0,len(sh0),axis0,axis1)
593     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
594     text+=" def %s(self):\n"%tname
595     a_0=makeArray(sh0,[-8,8])
596     if case0 in ["taggedData", "expandedData"]:
597     a1_0=makeArray(sh0,[-8,8])
598     else:
599     a1_0=a_0
600     r=numarray.swapaxes(a_0,axis0,axis1)
601     r1=numarray.swapaxes(a1_0,axis0,axis1)
602     text+=mkText(case0,"arg",a_0,a1_0)
603     text+=" res=swap_axes(arg,axis0=%s,axis1=%s)\n"%(axis0,axis1)
604     if case0=="Symbol":
605     text+=mkText("array","s",a_0,a1_0)
606     text+=" sub=res.substitute({arg:s})\n"
607     res="sub"
608     text+=mkText("array","ref",r,r1)
609     else:
610     res="res"
611     text+=mkText(case0,"ref",r,r1)
612     text+=mkTypeAndShapeTest(case0,r.shape,"res")
613     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
614    
615     if case0 == "taggedData":
616     t_prog_with_tags+=text
617     else:
618     t_prog+=text
619     print test_header
620     # print t_prog
621     print t_prog_with_tags
622     print test_tail
623     1/0
624     while True:
625     for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
626     if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
627     and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
628     arg_shape=oper[2](sh0,sh1,sh_s)
629     if not arg_shape==None:
630     case=getResultCaseForBin(case0,case1)
631     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
632     if oper[0] in [ "generalTensorProduct", "generalTransposedTensorProduct", "generalTensorTransposedProduct"]:
633     tname="test_%s_%s_rank%s_%s_rank%s_offset%s"%(oper[0],case0,len(sh0+sh_s),case1,len(sh_s+sh1),len(sh_s))
634    
635     res_text=" res=%s(arg0,arg1,axis_offset=%s)\n"%(oper[0],len(sh_s))
636     else:
637     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0+sh_s),case1,len(sh_s+sh1))
638    
639     res_text=" res=%s(arg0,arg1)\n"%oper[0]
640    
641     a_0=makeArray(arg_shape[0],[-8,8])
642     if case0 in ["taggedData", "expandedData"]:
643     a1_0=makeArray(arg_shape[0],[-8,8])
644     else:
645     a1_0=a_0
646    
647     a_1=makeArray(arg_shape[1],[-8,8])
648     if case1 in ["taggedData", "expandedData"]:
649     a1_1=makeArray(arg_shape[1],[-8,8])
650     else:
651     a1_1=a_1
652     r=oper[1](a_0,a_1,sh_s)
653     r1=oper[1](a1_0,a1_1,sh_s)
654     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
655     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
656     text+=res_text
657     if case=="Symbol":
658     c0_res,c1_res=case0,case1
659     subs="{"
660     if case0=="Symbol":
661     text+=mkText("array","s0",a_0,a1_0)
662     subs+="arg0:s0"
663     c0_res="array"
664     if case1=="Symbol":
665     text+=mkText("array","s1",a_1,a1_1)
666     if not subs.endswith("{"): subs+=","
667     subs+="arg1:s1"
668     c1_res="array"
669     subs+="}"
670     text+=" sub=res.substitute(%s)\n"%subs
671     res="sub"
672     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
673     else:
674     res="res"
675     text+=mkText(case,"ref",r,r1)
676     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
677     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
678     1/0
679    
680 gross 785 #=================================================
681     def testMatrixMult(arg0,arg1,sh_s):
682     return numarray.matrixmultiply(arg0,arg1)
683 gross 587
684 gross 785 def shapeTestMatrixMult(sh0,sh1,sh_s):
685     if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2):
686     return sh0+sh_s,sh_s+sh1
687 gross 587
688 gross 785 def testTensorMult(arg0,arg1,sh_s):
689     if arg0.rank==2:
690     return numarray.matrixmultiply(arg0,arg1)
691     else:
692     if arg1.rank==4:
693     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
694     for i0 in range(arg0.shape[0]):
695     for i1 in range(arg0.shape[1]):
696     for i2 in range(arg0.shape[2]):
697     for i3 in range(arg0.shape[3]):
698     for j2 in range(arg1.shape[2]):
699     for j3 in range(arg1.shape[3]):
700     out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
701     elif arg1.rank==3:
702     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
703     for i0 in range(arg0.shape[0]):
704     for i1 in range(arg0.shape[1]):
705     for i2 in range(arg0.shape[2]):
706     for i3 in range(arg0.shape[3]):
707     for j2 in range(arg1.shape[2]):
708     out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
709     elif arg1.rank==2:
710     out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
711     for i0 in range(arg0.shape[0]):
712     for i1 in range(arg0.shape[1]):
713     for i2 in range(arg0.shape[2]):
714     for i3 in range(arg0.shape[3]):
715     out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
716     return out
717    
718     def shapeTestTensorMult(sh0,sh1,sh_s):
719     if ( len(sh_s)==1 and len(sh0+sh_s)==2 and ( len(sh1+sh_s)==2 or len(sh1+sh_s)==1 )) or (len(sh_s)==2 and len(sh0+sh_s)==4 and (len(sh1+sh_s)==2 or len(sh1+sh_s)==3 or len(sh1+sh_s)==4)):
720     return sh0+sh_s,sh_s+sh1
721    
722     def tensorProductTest(arg0,arg1,sh_s):
723     if isinstance(arg0,float):
724     out=numarray.array(arg0*arg1)
725     elif isinstance(arg1,float):
726     out=numarray.array(arg0*arg1)
727     elif len(sh_s)==0:
728     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
729     else:
730     l=len(sh_s)
731     sh0=arg0.shape[:arg0.rank-l]
732     sh1=arg1.shape[l:]
733     ls,l0,l1=1,1,1
734     for i in sh_s: ls*=i
735     for i in sh0: l0*=i
736     for i in sh1: l1*=i
737     out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
738     out2=numarray.zeros((l0,l1),numarray.Float)
739     for i0 in range(l0):
740     for i1 in range(l1):
741     for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
742     out=out2.resize(sh0+sh1)
743     return out
744     def tensorProductShape(sh0,sh1,sh_s):
745     return sh0+sh_s,sh_s+sh1
746    
747     def transpose(r,offset):
748     if isinstance(r,float): return r
749     s=r.shape
750     s1=1
751     for i in s[:offset]: s1*=i
752     s2=1
753     for i in s[offset:]: s2*=i
754     out=numarray.reshape(r,(s1,s2))
755     out.transpose()
756     return numarray.resize(out,s[offset:]+s[:offset])
757    
758     def shapeTestTMatrixMult(sh0,sh1,sh_s):
759     if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2):
760     return sh_s+sh0,sh_s+sh1
761    
762     def testTMatrixMult(arg0,arg1,sh_s):
763     return numarray.matrixmultiply(transpose(arg0,1),arg1)
764    
765     def shapeTestTTensorMult(sh0,sh1,sh_s):
766     if ( len(sh_s)==1 and len(sh0+sh_s)==2 and ( len(sh1+sh_s)==2 or len(sh1+sh_s)==1 )) or (len(sh_s)==2 and len(sh0+sh_s)==4 and (len(sh1+sh_s)==2 or len(sh1+sh_s)==3 or len(sh1+sh_s)==4)):
767     return sh_s+sh0,sh_s+sh1
768    
769     def testTTensorMult(arg0,arg1,sh_s):
770     if arg0.rank==2:
771     return numarray.matrixmultiply(transpose(arg0,1),arg1)
772     else:
773     if arg1.rank==4:
774     out=numarray.zeros((arg0.shape[2],arg0.shape[3],arg1.shape[2],arg1.shape[3]),numarray.Float)
775     for i0 in range(arg0.shape[2]):
776     for i1 in range(arg0.shape[3]):
777     for i2 in range(arg0.shape[0]):
778     for i3 in range(arg0.shape[1]):
779     for j2 in range(arg1.shape[2]):
780     for j3 in range(arg1.shape[3]):
781     out[i0,i1,j2,j3]+=arg0[i2,i3,i0,i1]*arg1[i2,i3,j2,j3]
782     elif arg1.rank==3:
783     out=numarray.zeros((arg0.shape[2],arg0.shape[3],arg1.shape[2]),numarray.Float)
784     for i0 in range(arg0.shape[2]):
785     for i1 in range(arg0.shape[3]):
786     for i2 in range(arg0.shape[0]):
787     for i3 in range(arg0.shape[1]):
788     for j2 in range(arg1.shape[2]):
789     out[i0,i1,j2]+=arg0[i2,i3,i0,i1]*arg1[i2,i3,j2]
790     elif arg1.rank==2:
791     out=numarray.zeros((arg0.shape[2],arg0.shape[3]),numarray.Float)
792     for i0 in range(arg0.shape[2]):
793     for i1 in range(arg0.shape[3]):
794     for i2 in range(arg0.shape[0]):
795     for i3 in range(arg0.shape[1]):
796     out[i0,i1]+=arg0[i2,i3,i0,i1]*arg1[i2,i3]
797     return out
798    
799     def tensorTProductShape(sh0,sh1,sh_s):
800     return sh_s+sh0,sh_s+sh1
801    
802     def tensorTProductTest(arg0,arg1,sh_s):
803     if isinstance(arg0,float):
804     out=numarray.array(arg0*arg1)
805     elif isinstance(arg1,float):
806     out=numarray.array(arg0*arg1)
807     elif len(sh_s)==0:
808     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
809     else:
810     l=len(sh_s)
811     sh0=arg0.shape[l:]
812     sh1=arg1.shape[l:]
813     ls,l0,l1=1,1,1
814     for i in sh_s: ls*=i
815     for i in sh0: l0*=i
816     for i in sh1: l1*=i
817     out1=numarray.outerproduct(arg0,arg1).resize((ls,l0,ls,l1))
818     out2=numarray.zeros((l0,l1),numarray.Float)
819     for i0 in range(l0):
820     for i1 in range(l1):
821     for i in range(ls): out2[i0,i1]+=out1[i,i0,i,i1]
822     out=out2.resize(sh0+sh1)
823     return out
824    
825     def shapeTestMatrixTMult(sh0,sh1,sh_s):
826     if len(sh_s)==1 and len(sh0+sh_s)==2 and len(sh_s+sh1)==2:
827     return sh0+sh_s,sh1+sh_s
828    
829     def testMatrixTMult(arg0,arg1,sh_s):
830     return numarray.matrixmultiply(arg0,transpose(arg1,1))
831    
832     def shapeTestTensorTMult(sh0,sh1,sh_s):
833     if ( len(sh_s)==1 and len(sh0+sh_s)==2 and ( len(sh1+sh_s)==2 or len(sh1+sh_s)==1 )) or (len(sh_s)==2 and len(sh0+sh_s)==4 and (len(sh1+sh_s)==2 or len(sh1+sh_s)==3 or len(sh1+sh_s)==4)):
834     return sh0+sh_s,sh1+sh_s
835    
836     def testTensorTMult(arg0,arg1,sh_s):
837     if arg0.rank==2:
838     return numarray.matrixmultiply(arg0,transpose(arg1,1))
839     else:
840     if arg1.rank==4:
841     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[0],arg1.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     for j0 in range(arg1.shape[0]):
847     for j1 in range(arg1.shape[1]):
848     out[i0,i1,j0,j1]+=arg0[i0,i1,i2,i3]*arg1[j0,j1,i2,i3]
849     elif arg1.rank==3:
850     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[0]),numarray.Float)
851     for i0 in range(arg0.shape[0]):
852     for i1 in range(arg0.shape[1]):
853     for i2 in range(arg0.shape[2]):
854     for i3 in range(arg0.shape[3]):
855     for j0 in range(arg1.shape[0]):
856     out[i0,i1,j0]+=arg0[i0,i1,i2,i3]*arg1[j0,i2,i3]
857     elif arg1.rank==2:
858     out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
859     for i0 in range(arg0.shape[0]):
860     for i1 in range(arg0.shape[1]):
861     for i2 in range(arg0.shape[2]):
862     for i3 in range(arg0.shape[3]):
863     out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
864     return out
865    
866     def tensorProductTTest(arg0,arg1,sh_s):
867     if isinstance(arg0,float):
868     out=numarray.array(arg0*arg1)
869     elif isinstance(arg1,float):
870     out=numarray.array(arg0*arg1)
871     elif len(sh_s)==0:
872     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
873     else:
874     l=len(sh_s)
875     sh0=arg0.shape[:arg0.rank-l]
876     sh1=arg1.shape[:arg1.rank-l]
877     ls,l0,l1=1,1,1
878     for i in sh_s: ls*=i
879     for i in sh0: l0*=i
880     for i in sh1: l1*=i
881     out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,l1,ls))
882     out2=numarray.zeros((l0,l1),numarray.Float)
883     for i0 in range(l0):
884     for i1 in range(l1):
885     for i in range(ls): out2[i0,i1]+=out1[i0,i,i1,i]
886     out=out2.resize(sh0+sh1)
887     return out
888    
889     def tensorProductTShape(sh0,sh1,sh_s):
890     return sh0+sh_s,sh1+sh_s
891    
892 gross 530 #=======================================================================================================
893 gross 785 # tensor multiply
894     #=======================================================================================================
895     # oper=
896     # oper=
897     for oper in [ ["tensor_mult",testTensorMult,shapeTestTensorMult], \
898     ["generalTensorProduct",tensorProductTest,tensorProductShape], \
899     ["matrix_mult",testMatrixMult,shapeTestMatrixMult], \
900     ["transposed_matrix_mult",testTMatrixMult,shapeTestTMatrixMult], \
901     ["transposed_tensor_mult",testTTensorMult,shapeTestTTensorMult], \
902     ["generalTransposedTensorProduct",tensorTProductTest,tensorTProductShape], \
903     ["matrix_transposed_mult",testMatrixTMult,shapeTestMatrixTMult], \
904     ["tensor_transposed_mult",testTensorTMult,shapeTestTensorTMult], \
905     ["generalTensorTransposedProduct",tensorProductTTest,tensorProductTShape] ]:
906    
907     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
908     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
909     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
910     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
911     for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
912     if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
913     and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
914     arg_shape=oper[2](sh0,sh1,sh_s)
915     if not arg_shape==None:
916     case=getResultCaseForBin(case0,case1)
917     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
918     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
919     if oper[0] in [ "generalTensorProduct", "generalTransposedTensorProduct", "generalTensorTransposedProduct"]:
920     tname="test_%s_%s_rank%s_%s_rank%s_offset%s"%(oper[0],case0,len(sh0+sh_s),case1,len(sh_s+sh1),len(sh_s))
921    
922     res_text=" res=%s(arg0,arg1,axis_offset=%s)\n"%(oper[0],len(sh_s))
923     else:
924     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0+sh_s),case1,len(sh_s+sh1))
925    
926     res_text=" res=%s(arg0,arg1)\n"%oper[0]
927    
928     text+=" def %s(self):\n"%tname
929     a_0=makeArray(arg_shape[0],[-8,8])
930     if case0 in ["taggedData", "expandedData"]:
931     a1_0=makeArray(arg_shape[0],[-8,8])
932     else:
933     a1_0=a_0
934    
935     a_1=makeArray(arg_shape[1],[-8,8])
936     if case1 in ["taggedData", "expandedData"]:
937     a1_1=makeArray(arg_shape[1],[-8,8])
938     else:
939     a1_1=a_1
940     r=oper[1](a_0,a_1,sh_s)
941     r1=oper[1](a1_0,a1_1,sh_s)
942     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
943     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
944     text+=res_text
945     if case=="Symbol":
946     c0_res,c1_res=case0,case1
947     subs="{"
948     if case0=="Symbol":
949     text+=mkText("array","s0",a_0,a1_0)
950     subs+="arg0:s0"
951     c0_res="array"
952     if case1=="Symbol":
953     text+=mkText("array","s1",a_1,a1_1)
954     if not subs.endswith("{"): subs+=","
955     subs+="arg1:s1"
956     c1_res="array"
957     subs+="}"
958     text+=" sub=res.substitute(%s)\n"%subs
959     res="sub"
960     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
961     else:
962     res="res"
963     text+=mkText(case,"ref",r,r1)
964     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
965     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
966     if case0 == "taggedData" or case1 == "taggedData":
967     t_prog_with_tags+=text
968     else:
969     t_prog+=text
970     print test_header
971     print t_prog
972     # print t_prog_with_tags
973     print test_tail
974     1/0
975    
976    
977    
978     #=======================================================================================================
979 gross 587 # eigenvalues and eigen vectors 2D:
980     #=======================================================================================================
981     alltests= \
982     [ ("case0",[[0.0, 0.0], [0.0, 0.0]],[0.0, 0.0]) \
983     , ("case3",[[-1.0, 0.0], [0.0, -1.0]],[-1.0, -1.0]) \
984     , ("case5",[[-0.99999999999999967, -6.4606252205695602e-16], [-6.4606252205695602e-16, -0.99999999999999967]],[-1.0, -1.0]) \
985     , ("case6",[[0.0, 0.0], [0.0, 0.0001]],[0.0, 0.0001]) \
986     , ("case7",[[0.0001, 0.0], [0.0, 0.0]],[0.0, 0.0001]) \
987     , ("case8",[[6.0598371831785722e-06, 2.3859213977648625e-05], [2.3859213977648629e-05, 9.3940162816821425e-05]],[0.0, 0.0001]) \
988     , ("case9",[[1.0, 0.0], [0.0, 2.0]],[1.0, 2.0]) \
989     , ("case10",[[2.0, 0.0], [0.0, 1.0]],[1.0, 2.0]) \
990     , ("case11",[[1.0605983718317855, 0.23859213977648688], [0.23859213977648688, 1.9394016281682138]],[1.0, 2.0]) \
991     , ("case12",[[1.0, 0.0], [0.0, 1000000.0]],[1.0, 1000000.0]) \
992     , ("case13",[[1000000.0, 0.0], [0.0, 1.0]],[1.0, 1000000.0]) \
993     , ("case14",[[60599.311233413886, 238591.90118434647], [238591.90118434647, 939401.68876658613]],[1.0, 1000000.0]) \
994     ]
995     dim=2
996 gross 588
997    
998     alltests= \
999     [ ("case0",[[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],[0.0, 0.0, 0.0]) \
1000     , ("case5",[[10.0, 0.0, 0.0], [0.0, 10.0, 0.0], [0.0, 0.0, 10.0]],[10.0, 10.0, 10.0]) \
1001     , ("case10",[[0.9, 0.0, 0.0], [0.0, 0.9, 0.0], [0.0, 0.0, 1.0]],[0.9, 0.9, 1.0]) \
1002     , ("case11",[[0.9, 0.0, 0.0], [0.0, 0.97060899725040983, -0.045555123008643325], [0.0, -0.045555123008643339, 0.92939100274959041]],[0.9, 0.9, 1.0]) \
1003     , ("case12",[[0.92694799760252555, 0.0, 0.044368966468320177], [0.0, 0.9, 0.0], [0.044368966468320184, 0.0, 0.97305200239747425]],[0.9, 0.9, 1.0]) \
1004     , ("case13",[[1.0, 0.0, 0.0], [0.0, 0.9, 0.], [0.0, 0., 0.9]],[0.9, 0.9, 1.0]) \
1005     , ("case14",[[0.92379770619813639, 0.041031106298491521, -0.011396846732439278], [0.041031106298491535, 0.97074428392640366, -0.019650012730342326], [-0.011396846732439236, -0.019650012730342337, 0.90545800987545966]],[0.9, 0.9, 1.0]) \
1006     , ("case15",[[1.0, 0.0, 0.0], [0.0, 1.1, 0.0], [0.0, 0.0, 1.1]],[1.0, 1.1, 1.1]) \
1007     , ("case17",[[1.0269479976025255, 0.0, 0.044368966468320309], [0.0, 1.1, 0.0], [0.044368966468320295, 0.0, 1.0730520023974743]],[1.0, 1.1, 1.1]) \
1008     , ("case18",[[1.1, 0.0, 0.0], [0.0, 1.0153410887977139, -0.036038311201720394], [0.0, -0.036038311201720373, 1.084658911202286]],[1.0, 1.1, 1.1]) \
1009     , ("case19",[[1.035487967756175, 0.026317079185831614, -0.039960133424212368], [0.026317079185831618, 1.0892641940924184, 0.016301362071911414], [-0.039960133424212355, 0.016301362071911431, 1.0752478381514063]],[1.0, 1.1, 1.1]) \
1010     , ("case20",[[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 3.0]],[1.0, 2.0, 3.0]) \
1011     , ("case21",[[1.0, 0.0, 0.0], [0.0, 2.7060899725040968, -0.45555123008643206], [0.0, -0.45555123008643228, 2.2939100274959037]],[1.0, 2.0, 3.0]) \
1012     , ("case22",[[1.5389599520505153, 0.0, 0.88737932936638753], [0.0, 2.0, 0.0], [0.88737932936638753, 0.0, 2.4610400479494858]],[1.0, 2.0, 3.0]) \
1013     , ("case23",[[3.0, 0.0, 0.0], [0.0, 1.153410887977139, -0.36038311201720391], [0.0, -0.36038311201720391, 1.8465891120228608]],[1.0, 2.0, 3.0]) \
1014     , ("case24",[[1.5928567395431172, 0.67348185484323142, -0.51356980156651744], [0.67348185484323153, 2.6000847801882254, -0.033486506584313548], [-0.51356980156651744, -0.033486506584313541, 1.8070584802686565]],[1.0, 2.0, 3.0]) \
1015     , ("case25",[[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 30000.0]],[1.0, 2.0, 30000.0]) \
1016     , ("case26",[[1.0, 0.0, 0.0], [0.0, 21183.286995177881, -13665.625800132779], [0.0, -13665.625800132779, 8818.7130048221279]],[1.0, 2.0, 30000.0]) \
1017     , ("case27",[[8085.1298007817086, 0.0, 13310.246250831115], [0.0, 2.0, 0.0], [13310.246250831115, 0.0, 21915.870199218316]],[1.0, 2.0, 30000.0]) \
1018     , ("case28",[[30000.0, 0.0, 0.0], [0.0, 1.153410887977139, -0.36038311201720391], [0.0, -0.36038311201720391, 1.8465891120228608]],[1.0, 2.0, 30000.0]) \
1019     , ("case29",[[7140.1907849945546, 12308.774438213351, -3419.2256841313947], [12308.774438213351, 21223.762934183575, -5894.4478052274408], [-3419.2256841313947, -5894.4478052274408, 1639.0462808218595]],[1.0, 2.0, 30000.0]) \
1020     ]
1021    
1022     dim=3
1023    
1024     alltests= \
1025     [ ("case0",[[0.0]],[0.0]) \
1026     , ("case1",[[1.0]],[1.0]) \
1027     ]
1028     dim=1
1029    
1030 gross 587 for case in ["constData","taggedData","expandedData"]:
1031     n=0
1032     while n<len(alltests):
1033     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1034     tname="test_eigenvalues_and_eigenvectors_%s_dim%s_%s"%(case,dim,alltests[n][0])
1035     text+=" def %s(self):\n"%tname
1036     a_0=numarray.array(alltests[n][1],numarray.Float64)
1037     ev_0=numarray.array(alltests[n][2],numarray.Float64)
1038     if case in ["taggedData", "expandedData"]:
1039 gross 588 if n+1<len(alltests):
1040     a1_0=numarray.array(alltests[n+1][1],numarray.Float64)
1041     ev1_0=numarray.array(alltests[n+1][2],numarray.Float64)
1042     else:
1043     a1_0=numarray.array(alltests[0][1],numarray.Float64)
1044     ev1_0=numarray.array(alltests[0][2],numarray.Float64)
1045 gross 587 n+=2
1046     else:
1047     a1_0=a_0
1048     ev1_0=ev_0
1049     n+=1
1050     text+=mkText(case,"arg",a_0,a1_0)
1051     text+=" res=eigenvalues_and_eigenvectors(arg)\n"
1052     text+=mkText(case,"ref_ev",ev_0,ev1_0)
1053     text+=mkTypeAndShapeTest(case,(dim,),"res[0]"," for eigenvalues")
1054     text+=mkTypeAndShapeTest(case,(dim,dim),"res[1]"," for eigenvectors")
1055     text+=" self.failUnless(Lsup(res[0]-ref_ev)<=self.RES_TOL*Lsup(ref_ev),\"wrong eigenvalues\")\n"
1056     for i in range(dim):
1057     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)
1058     text+=" self.failUnless(Lsup(matrixmult(transpose(res[1]),res[1])-kronecker(%s))<=self.RES_TOL,\"eigenvectors are not orthonormal\")\n"%dim
1059     if case == "taggedData" :
1060     t_prog_with_tags+=text
1061     else:
1062     t_prog+=text
1063     print test_header
1064     print t_prog
1065     print t_prog_with_tags
1066     print test_tail
1067     1/0
1068     #=======================================================================================================
1069 gross 550 # get slices
1070 gross 536 #=======================================================================================================
1071     from esys.escript import *
1072 gross 550 for case0 in ["constData","taggedData","expandedData"]:
1073     for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,3,3)]:
1074     # get perm:
1075     if len(sh0)==2:
1076     check=[[1,0]]
1077     elif len(sh0)==3:
1078     check=[[1,0,2],
1079     [1,2,0],
1080     [2,1,0],
1081     [2,0,2],
1082     [0,2,1]]
1083     elif len(sh0)==4:
1084     check=[[0,1,3,2],
1085     [0,2,1,3],
1086     [0,2,3,1],
1087     [0,3,2,1],
1088     [0,3,1,2] ,
1089     [1,0,2,3],
1090     [1,0,3,2],
1091     [1,2,0,3],
1092     [1,2,3,0],
1093     [1,3,2,0],
1094     [1,3,0,2],
1095     [2,0,1,3],
1096     [2,0,3,1],
1097     [2,1,0,3],
1098     [2,1,3,0],
1099     [2,3,1,0],
1100     [2,3,0,1],
1101     [3,0,1,2],
1102     [3,0,2,1],
1103     [3,1,0,2],
1104     [3,1,2,0],
1105     [3,2,1,0],
1106     [3,2,0,1]]
1107     else:
1108     check=[]
1109    
1110     # create the test cases:
1111     processed=[]
1112     l=["R","U","L","P","C","N"]
1113     c=[""]
1114     for i in range(len(sh0)):
1115     tmp=[]
1116     for ci in c:
1117     tmp+=[ci+li for li in l]
1118     c=tmp
1119     # SHUFFLE
1120     c2=[]
1121     while len(c)>0:
1122     i=int(random.random()*len(c))
1123     c2.append(c[i])
1124     del c[i]
1125     c=c2
1126     for ci in c:
1127     t=""
1128     sh=()
1129     sl=()
1130     for i in range(len(ci)):
1131     if ci[i]=="R":
1132     s="%s:%s"%(1,sh0[i]-1)
1133     sl=sl+(slice(1,sh0[i]-1),)
1134     sh=sh+(sh0[i]-2,)
1135     if ci[i]=="U":
1136     s=":%s"%(sh0[i]-1)
1137     sh=sh+(sh0[i]-1,)
1138     sl=sl+(slice(0,sh0[i]-1),)
1139     if ci[i]=="L":
1140     s="2:"
1141     sh=sh+(sh0[i]-2,)
1142     sl=sl+(slice(2,sh0[i]),)
1143     if ci[i]=="P":
1144     s="%s"%(int(sh0[i]/2))
1145     sl=sl+(int(sh0[i]/2),)
1146     if ci[i]=="C":
1147     s=":"
1148     sh=sh+(sh0[i],)
1149     sl=sl+(slice(0,sh0[i]),)
1150     if ci[i]=="N":
1151     s=""
1152     sh=sh+(sh0[i],)
1153     if len(s)>0:
1154     if not t=="": t+=","
1155     t+=s
1156     if len(sl)==1: sl=sl[0]
1157     N_found=False
1158     noN_found=False
1159     process=len(t)>0
1160     for i in ci:
1161     if i=="N":
1162     if not noN_found and N_found: process=False
1163     N_found=True
1164     else:
1165     if N_found: process=False
1166     noNfound=True
1167     # is there a similar one processed allready
1168     if process and ci.find("N")==-1:
1169     for ci2 in processed:
1170     for chi in check:
1171     is_perm=True
1172     for i in range(len(chi)):
1173     if not ci[i]==ci2[chi[i]]: is_perm=False
1174     if is_perm: process=False
1175     # if not process: print ci," rejected"
1176     if process:
1177     processed.append(ci)
1178     for case1 in ["array","constData","taggedData","expandedData"]:
1179     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1180     tname="test_setslice_%s_rank%s_%s_%s"%(case0,len(sh0),case1,ci)
1181     text+=" def %s(self):\n"%tname
1182     a_0=makeNumberedArray(sh0)
1183     if case0 in ["taggedData", "expandedData"]:
1184     a1_0=makeNumberedArray(sh0)
1185     else:
1186     a1_0=a_0*1.
1187    
1188     a_1=makeNumberedArray(sh)
1189     if case1 in ["taggedData", "expandedData"]:
1190     a1_1=makeNumberedArray(sh)
1191     else:
1192     a1_1=a_1*1.
1193    
1194     text+=mkText(case0,"arg",a_0,a1_0)
1195     text+=mkText(case1,"val",a_1,a1_1)
1196     text+=" arg[%s]=val\n"%t
1197     a_0.__setitem__(sl,a_1)
1198     a1_0.__setitem__(sl,a1_1)
1199     if Lsup(a_0-a1_0)==0.:
1200     text+=mkText("constData","ref",a_0,a1_0)
1201     else:
1202     text+=mkText("expandedData","ref",a_0,a1_0)
1203     text+=" self.failUnless(Lsup(arg-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"
1204    
1205     if case0 == "taggedData" or case1 == "taggedData":
1206     t_prog_with_tags+=text
1207     else:
1208     t_prog+=text
1209    
1210     print test_header
1211     # print t_prog
1212     print t_prog_with_tags
1213     print test_tail
1214     1/0
1215    
1216     #=======================================================================================================
1217     # (non)symmetric part
1218     #=======================================================================================================
1219     from esys.escript import *
1220 gross 536 for name in ["symmetric", "nonsymmetric"]:
1221     f=1.
1222     if name=="nonsymmetric": f=-1
1223     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1224     for sh0 in [ (3,3), (2,3,2,3)]:
1225     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1226     tname="test_%s_%s_rank%s"%(name,case0,len(sh0))
1227     text+=" def %s(self):\n"%tname
1228     a_0=makeNumberedArray(sh0,s=1.)
1229     r_0=(a_0+f*transpose(a_0))/2.
1230     if case0 in ["taggedData", "expandedData"]:
1231     a1_0=makeNumberedArray(sh0,s=-1.)
1232     r1_0=(a1_0+f*transpose(a1_0))/2.
1233     else:
1234     a1_0=a_0
1235     r1_0=r_0
1236     text+=mkText(case0,"arg",a_0,a1_0)
1237     text+=" res=%s(arg)\n"%name
1238     if case0=="Symbol":
1239     text+=mkText("array","s",a_0,a1_0)
1240     text+=" sub=res.substitute({arg:s})\n"
1241     res="sub"
1242     text+=mkText("array","ref",r_0,r1_0)
1243     else:
1244     res="res"
1245     text+=mkText(case0,"ref",r_0,r1_0)
1246     text+=mkTypeAndShapeTest(case0,sh0,"res")
1247     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1248    
1249     if case0 == "taggedData" :
1250     t_prog_with_tags+=text
1251     else:
1252     t_prog+=text
1253     print test_header
1254     print t_prog
1255     # print t_prog_with_tags
1256     print test_tail
1257     1/0
1258    
1259     #=======================================================================================================
1260 gross 530 # eigenvalues
1261     #=======================================================================================================
1262     import numarray.linear_algebra
1263     name="eigenvalues"
1264     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1265     for sh0 in [ (1,1), (2,2), (3,3)]:
1266     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1267     tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1268     text+=" def %s(self):\n"%tname
1269     a_0=makeArray(sh0,[-1.,1])
1270     a_0=(a_0+numarray.transpose(a_0))/2.
1271     ev=numarray.linear_algebra.eigenvalues(a_0)
1272     ev.sort()
1273     if case0 in ["taggedData", "expandedData"]:
1274     a1_0=makeArray(sh0,[-1.,1])
1275     a1_0=(a1_0+numarray.transpose(a1_0))/2.
1276     ev1=numarray.linear_algebra.eigenvalues(a1_0)
1277     ev1.sort()
1278     else:
1279     a1_0=a_0
1280     ev1=ev
1281     text+=mkText(case0,"arg",a_0,a1_0)
1282     text+=" res=%s(arg)\n"%name
1283     if case0=="Symbol":
1284     text+=mkText("array","s",a_0,a1_0)
1285     text+=" sub=res.substitute({arg:s})\n"
1286     res="sub"
1287     text+=mkText("array","ref",ev,ev1)
1288     else:
1289     res="res"
1290     text+=mkText(case0,"ref",ev,ev1)
1291     text+=mkTypeAndShapeTest(case0,(sh0[0],),"res")
1292     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1293    
1294     if case0 == "taggedData" :
1295     t_prog_with_tags+=text
1296     else:
1297     t_prog+=text
1298     print test_header
1299     # print t_prog
1300     print t_prog_with_tags
1301     print test_tail
1302     1/0
1303 jgs 154
1304 gross 517 #=======================================================================================================
1305 gross 550 # get slices
1306 gross 517 #=======================================================================================================
1307     for case0 in ["constData","taggedData","expandedData","Symbol"]:
1308     for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,5,3)]:
1309     # get perm:
1310     if len(sh0)==2:
1311     check=[[1,0]]
1312     elif len(sh0)==3:
1313     check=[[1,0,2],
1314     [1,2,0],
1315     [2,1,0],
1316     [2,0,2],
1317     [0,2,1]]
1318     elif len(sh0)==4:
1319     check=[[0,1,3,2],
1320     [0,2,1,3],
1321     [0,2,3,1],
1322     [0,3,2,1],
1323     [0,3,1,2] ,
1324     [1,0,2,3],
1325     [1,0,3,2],
1326     [1,2,0,3],
1327     [1,2,3,0],
1328     [1,3,2,0],
1329     [1,3,0,2],
1330     [2,0,1,3],
1331     [2,0,3,1],
1332     [2,1,0,3],
1333     [2,1,3,0],
1334     [2,3,1,0],
1335     [2,3,0,1],
1336     [3,0,1,2],
1337     [3,0,2,1],
1338     [3,1,0,2],
1339     [3,1,2,0],
1340     [3,2,1,0],
1341     [3,2,0,1]]
1342     else:
1343     check=[]
1344    
1345     # create the test cases:
1346     processed=[]
1347     l=["R","U","L","P","C","N"]
1348     c=[""]
1349     for i in range(len(sh0)):
1350     tmp=[]
1351     for ci in c:
1352     tmp+=[ci+li for li in l]
1353     c=tmp
1354     # SHUFFLE
1355     c2=[]
1356     while len(c)>0:
1357     i=int(random.random()*len(c))
1358     c2.append(c[i])
1359     del c[i]
1360     c=c2
1361     for ci in c:
1362     t=""
1363     sh=()
1364     for i in range(len(ci)):
1365     if ci[i]=="R":
1366     s="%s:%s"%(1,sh0[i]-1)
1367     sh=sh+(sh0[i]-2,)
1368     if ci[i]=="U":
1369     s=":%s"%(sh0[i]-1)
1370     sh=sh+(sh0[i]-1,)
1371     if ci[i]=="L":
1372     s="2:"
1373     sh=sh+(sh0[i]-2,)
1374     if ci[i]=="P":
1375     s="%s"%(int(sh0[i]/2))
1376     if ci[i]=="C":
1377     s=":"
1378     sh=sh+(sh0[i],)
1379     if ci[i]=="N":
1380     s=""
1381     sh=sh+(sh0[i],)
1382     if len(s)>0:
1383     if not t=="": t+=","
1384     t+=s
1385     N_found=False
1386     noN_found=False
1387     process=len(t)>0
1388     for i in ci:
1389     if i=="N":
1390     if not noN_found and N_found: process=False
1391     N_found=True
1392     else:
1393     if N_found: process=False
1394     noNfound=True
1395     # is there a similar one processed allready
1396     if process and ci.find("N")==-1:
1397     for ci2 in processed:
1398     for chi in check:
1399     is_perm=True
1400     for i in range(len(chi)):
1401     if not ci[i]==ci2[chi[i]]: is_perm=False
1402     if is_perm: process=False
1403     # if not process: print ci," rejected"
1404     if process:
1405     processed.append(ci)
1406     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1407     tname="test_getslice_%s_rank%s_%s"%(case0,len(sh0),ci)
1408     text+=" def %s(self):\n"%tname
1409     a_0=makeNumberedArray(sh0,s=1)
1410     if case0 in ["taggedData", "expandedData"]:
1411     a1_0=makeNumberedArray(sh0,s=-1.)
1412     else:
1413     a1_0=a_0
1414     r=eval("a_0[%s]"%t)
1415     r1=eval("a1_0[%s]"%t)
1416     text+=mkText(case0,"arg",a_0,a1_0)
1417     text+=" res=arg[%s]\n"%t
1418     if case0=="Symbol":
1419     text+=mkText("array","s",a_0,a1_0)
1420     text+=" sub=res.substitute({arg:s})\n"
1421     res="sub"
1422     text+=mkText("array","ref",r,r1)
1423     else:
1424     res="res"
1425     text+=mkText(case0,"ref",r,r1)
1426     text+=mkTypeAndShapeTest(case0,sh,"res")
1427     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1428    
1429     if case0 == "taggedData" :
1430     t_prog_with_tags+=text
1431     else:
1432     t_prog+=text
1433    
1434     print test_header
1435     # print t_prog
1436     print t_prog_with_tags
1437     print test_tail
1438     1/0
1439     #============================================================================================
1440 gross 291 def innerTEST(arg0,arg1):
1441     if isinstance(arg0,float):
1442     out=numarray.array(arg0*arg1)
1443     else:
1444     out=(arg0*arg1).sum()
1445     return out
1446    
1447     def outerTEST(arg0,arg1):
1448     if isinstance(arg0,float):
1449     out=numarray.array(arg0*arg1)
1450     elif isinstance(arg1,float):
1451     out=numarray.array(arg0*arg1)
1452     else:
1453     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1454     return out
1455    
1456     def tensorProductTest(arg0,arg1,sh_s):
1457     if isinstance(arg0,float):
1458     out=numarray.array(arg0*arg1)
1459     elif isinstance(arg1,float):
1460     out=numarray.array(arg0*arg1)
1461     elif len(sh_s)==0:
1462     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1463     else:
1464     l=len(sh_s)
1465     sh0=arg0.shape[:arg0.rank-l]
1466     sh1=arg1.shape[l:]
1467     ls,l0,l1=1,1,1
1468     for i in sh_s: ls*=i
1469     for i in sh0: l0*=i
1470     for i in sh1: l1*=i
1471     out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
1472     out2=numarray.zeros((l0,l1),numarray.Float)
1473     for i0 in range(l0):
1474     for i1 in range(l1):
1475     for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
1476     out=out2.resize(sh0+sh1)
1477     return out
1478    
1479     def testMatrixMult(arg0,arg1,sh_s):
1480     return numarray.matrixmultiply(arg0,arg1)
1481    
1482    
1483     def testTensorMult(arg0,arg1,sh_s):
1484     if len(arg0)==2:
1485     return numarray.matrixmultiply(arg0,arg1)
1486     else:
1487     if arg1.rank==4:
1488     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
1489     for i0 in range(arg0.shape[0]):
1490     for i1 in range(arg0.shape[1]):
1491     for i2 in range(arg0.shape[2]):
1492     for i3 in range(arg0.shape[3]):
1493     for j2 in range(arg1.shape[2]):
1494     for j3 in range(arg1.shape[3]):
1495     out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
1496     elif arg1.rank==3:
1497     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
1498     for i0 in range(arg0.shape[0]):
1499     for i1 in range(arg0.shape[1]):
1500     for i2 in range(arg0.shape[2]):
1501     for i3 in range(arg0.shape[3]):
1502     for j2 in range(arg1.shape[2]):
1503     out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
1504     elif arg1.rank==2:
1505     out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
1506     for i0 in range(arg0.shape[0]):
1507     for i1 in range(arg0.shape[1]):
1508     for i2 in range(arg0.shape[2]):
1509     for i3 in range(arg0.shape[3]):
1510     out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
1511     return out
1512 gross 313
1513     def testReduce(arg0,init_val,test_expr,post_expr):
1514     out=init_val
1515     if isinstance(arg0,float):
1516     out=eval(test_expr.replace("%a1%","arg0"))
1517     elif arg0.rank==0:
1518     out=eval(test_expr.replace("%a1%","arg0"))
1519     elif arg0.rank==1:
1520     for i0 in range(arg0.shape[0]):
1521     out=eval(test_expr.replace("%a1%","arg0[i0]"))
1522     elif arg0.rank==2:
1523     for i0 in range(arg0.shape[0]):
1524     for i1 in range(arg0.shape[1]):
1525     out=eval(test_expr.replace("%a1%","arg0[i0,i1]"))
1526     elif arg0.rank==3:
1527     for i0 in range(arg0.shape[0]):
1528     for i1 in range(arg0.shape[1]):
1529     for i2 in range(arg0.shape[2]):
1530     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2]"))
1531     elif arg0.rank==4:
1532     for i0 in range(arg0.shape[0]):
1533     for i1 in range(arg0.shape[1]):
1534     for i2 in range(arg0.shape[2]):
1535     for i3 in range(arg0.shape[3]):
1536     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))
1537     return eval(post_expr)
1538 gross 396
1539     def clipTEST(arg0,mn,mx):
1540     if isinstance(arg0,float):
1541     return max(min(arg0,mx),mn)
1542     out=numarray.zeros(arg0.shape,numarray.Float64)
1543     if arg0.rank==1:
1544     for i0 in range(arg0.shape[0]):
1545     out[i0]=max(min(arg0[i0],mx),mn)
1546     elif arg0.rank==2:
1547     for i0 in range(arg0.shape[0]):
1548     for i1 in range(arg0.shape[1]):
1549     out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
1550     elif arg0.rank==3:
1551     for i0 in range(arg0.shape[0]):
1552     for i1 in range(arg0.shape[1]):
1553     for i2 in range(arg0.shape[2]):
1554     out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
1555     elif arg0.rank==4:
1556     for i0 in range(arg0.shape[0]):
1557     for i1 in range(arg0.shape[1]):
1558     for i2 in range(arg0.shape[2]):
1559     for i3 in range(arg0.shape[3]):
1560     out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
1561     return out
1562     def minimumTEST(arg0,arg1):
1563     if isinstance(arg0,float):
1564     if isinstance(arg1,float):
1565     if arg0>arg1:
1566     return arg1
1567     else:
1568     return arg0
1569     else:
1570     arg0=numarray.ones(arg1.shape)*arg0
1571     else:
1572     if isinstance(arg1,float):
1573     arg1=numarray.ones(arg0.shape)*arg1
1574     out=numarray.zeros(arg0.shape,numarray.Float64)
1575     if arg0.rank==0:
1576     if arg0>arg1:
1577     out=arg1
1578     else:
1579     out=arg0
1580     elif arg0.rank==1:
1581     for i0 in range(arg0.shape[0]):
1582     if arg0[i0]>arg1[i0]:
1583     out[i0]=arg1[i0]
1584     else:
1585     out[i0]=arg0[i0]
1586     elif arg0.rank==2:
1587     for i0 in range(arg0.shape[0]):
1588     for i1 in range(arg0.shape[1]):
1589     if arg0[i0,i1]>arg1[i0,i1]:
1590     out[i0,i1]=arg1[i0,i1]
1591     else:
1592     out[i0,i1]=arg0[i0,i1]
1593     elif arg0.rank==3:
1594     for i0 in range(arg0.shape[0]):
1595     for i1 in range(arg0.shape[1]):
1596     for i2 in range(arg0.shape[2]):
1597     if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
1598     out[i0,i1,i2]=arg1[i0,i1,i2]
1599     else:
1600     out[i0,i1,i2]=arg0[i0,i1,i2]
1601     elif arg0.rank==4:
1602     for i0 in range(arg0.shape[0]):
1603     for i1 in range(arg0.shape[1]):
1604     for i2 in range(arg0.shape[2]):
1605     for i3 in range(arg0.shape[3]):
1606     if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
1607     out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
1608     else:
1609     out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
1610     return out
1611 gross 443
1612 gross 439 def unrollLoops(a,b,o,arg,tap="",x="x"):
1613 gross 437 out=""
1614     if a.rank==1:
1615     z=""
1616     for i99 in range(a.shape[0]):
1617     if not z=="": z+="+"
1618     if o=="1":
1619 gross 439 z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
1620 gross 437 else:
1621 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
1622 gross 437
1623     out+=tap+"%s=%s\n"%(arg,z)
1624    
1625     elif a.rank==2:
1626     for i0 in range(a.shape[0]):
1627     z=""
1628     for i99 in range(a.shape[1]):
1629     if not z=="": z+="+"
1630     if o=="1":
1631     z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
1632     else:
1633 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
1634 gross 437
1635     out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
1636     elif a.rank==3:
1637     for i0 in range(a.shape[0]):
1638     for i1 in range(a.shape[1]):
1639     z=""
1640     for i99 in range(a.shape[2]):
1641     if not z=="": z+="+"
1642     if o=="1":
1643 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
1644 gross 437 else:
1645 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
1646 gross 437
1647     out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
1648     elif a.rank==4:
1649     for i0 in range(a.shape[0]):
1650     for i1 in range(a.shape[1]):
1651     for i2 in range(a.shape[2]):
1652     z=""
1653     for i99 in range(a.shape[3]):
1654     if not z=="": z+="+"
1655     if o=="1":
1656 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
1657 gross 437 else:
1658 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
1659 gross 437
1660     out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
1661     elif a.rank==5:
1662     for i0 in range(a.shape[0]):
1663     for i1 in range(a.shape[1]):
1664     for i2 in range(a.shape[2]):
1665     for i3 in range(a.shape[3]):
1666     z=""
1667     for i99 in range(a.shape[4]):
1668     if not z=="": z+="+"
1669     if o=="1":
1670 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
1671 gross 437 else:
1672 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)
1673 gross 437
1674     out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
1675     return out
1676    
1677     def unrollLoopsOfGrad(a,b,o,arg,tap=""):
1678     out=""
1679     if a.rank==1:
1680     for i99 in range(a.shape[0]):
1681     if o=="1":
1682     out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
1683     else:
1684     out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
1685    
1686     elif a.rank==2:
1687     for i0 in range(a.shape[0]):
1688     for i99 in range(a.shape[1]):
1689     if o=="1":
1690     out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
1691     else:
1692     out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
1693    
1694     elif a.rank==3:
1695     for i0 in range(a.shape[0]):
1696     for i1 in range(a.shape[1]):
1697     for i99 in range(a.shape[2]):
1698     if o=="1":
1699     out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
1700     else:
1701     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])
1702    
1703     elif a.rank==4:
1704     for i0 in range(a.shape[0]):
1705     for i1 in range(a.shape[1]):
1706     for i2 in range(a.shape[2]):
1707     for i99 in range(a.shape[3]):
1708     if o=="1":
1709     out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
1710     else:
1711     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])
1712 gross 438 return out
1713 gross 441 def unrollLoopsOfDiv(a,b,o,arg,tap=""):
1714     out=tap+arg+"="
1715     if o=="1":
1716     z=0.
1717     for i99 in range(a.shape[0]):
1718     z+=b[i99,i99]+a[i99,i99]
1719     out+="(%s)"%z
1720     else:
1721     z=0.
1722     for i99 in range(a.shape[0]):
1723     z+=b[i99,i99]
1724     if i99>0: out+="+"
1725     out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
1726     out+="+(%s)"%z
1727     return out
1728    
1729 gross 437 def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
1730     if where=="Function":
1731     xfac_o=1.
1732     xfac_op=0.
1733     z_fac=1./2.
1734     z_fac_s=""
1735     zo_fac_s="/(o+1.)"
1736     elif where=="FunctionOnBoundary":
1737     xfac_o=1.
1738     xfac_op=0.
1739     z_fac=1.
1740     z_fac_s="*dim"
1741     zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
1742     elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1743     xfac_o=0.
1744     xfac_op=1.
1745     z_fac=1./2.
1746     z_fac_s=""
1747     zo_fac_s="/(o+1.)"
1748     out=""
1749     if a.rank==1:
1750     zo=0.
1751     zop=0.
1752     z=0.
1753     for i99 in range(a.shape[0]):
1754     if i99==0:
1755     zo+= xfac_o*a[i99]
1756     zop+= xfac_op*a[i99]
1757     else:
1758     zo+=a[i99]
1759     z+=b[i99]
1760    
1761     out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
1762     if zop==0.:
1763     out+="\n"
1764     else:
1765     out+="+(%s)*0.5**o\n"%zop
1766     elif a.rank==2:
1767     for i0 in range(a.shape[0]):
1768     zo=0.
1769     zop=0.
1770     z=0.
1771     for i99 in range(a.shape[1]):
1772     if i99==0:
1773     zo+= xfac_o*a[i0,i99]
1774     zop+= xfac_op*a[i0,i99]
1775     else:
1776     zo+=a[i0,i99]
1777     z+=b[i0,i99]
1778    
1779     out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
1780     if zop==0.:
1781     out+="\n"
1782     else:
1783     out+="+(%s)*0.5**o\n"%zop
1784     elif a.rank==3:
1785