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

Annotation of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


Revision 785 - (hide annotations)
Tue Jul 25 03:48:10 2006 UTC (16 years, 2 months ago) by gross
File size: 164408 byte(s)
some new functions around tensor product added.
1 jgs 154 #!/usr/bin/python
2 gross 283 # $Id$
3 jgs 154
4     """
5     program generates parts of the util.py and the test_util.py script
6     """
7 gross 157 test_header=""
8     test_header+="import unittest\n"
9     test_header+="import numarray\n"
10     test_header+="from esys.escript import *\n"
11     test_header+="from esys.finley import Rectangle\n"
12     test_header+="class Test_util2(unittest.TestCase):\n"
13     test_header+=" RES_TOL=1.e-7\n"
14     test_header+=" def setUp(self):\n"
15     test_header+=" self.__dom =Rectangle(11,11,2)\n"
16     test_header+=" self.functionspace = FunctionOnBoundary(self.__dom)\n"
17     test_tail=""
18     test_tail+="suite = unittest.TestSuite()\n"
19     test_tail+="suite.addTest(unittest.makeSuite(Test_util2))\n"
20     test_tail+="unittest.TextTestRunner(verbosity=2).run(suite)\n"
21    
22     case_set=["float","array","constData","taggedData","expandedData","Symbol"]
23     shape_set=[ (),(2,), (4,5), (6,2,2),(3,2,3,4)]
24    
25 jgs 154 t_prog=""
26 gross 157 t_prog_with_tags=""
27     t_prog_failing=""
28     u_prog=""
29 jgs 154
30 gross 157 def wherepos(arg):
31     if arg>0.:
32     return 1.
33     else:
34     return 0.
35    
36 gross 313
37 gross 157 class OPERATOR:
38     def __init__(self,nickname,rng=[-1000.,1000],test_expr="",math_expr=None,
39     numarray_expr="",symbol_expr=None,diff=None,name=""):
40     self.nickname=nickname
41     self.rng=rng
42     self.test_expr=test_expr
43     if math_expr==None:
44     self.math_expr=test_expr.replace("%a1%","arg")
45     else:
46     self.math_expr=math_expr
47     self.numarray_expr=numarray_expr
48     self.symbol_expr=symbol_expr
49     self.diff=diff
50     self.name=name
51    
52 jgs 154 import random
53     import numarray
54     import math
55     finc=1.e-6
56    
57 gross 157 def getResultCaseForBin(case0,case1):
58     if case0=="float":
59     if case1=="float":
60     case="float"
61     elif case1=="array":
62     case="array"
63     elif case1=="constData":
64     case="constData"
65     elif case1=="taggedData":
66     case="taggedData"
67     elif case1=="expandedData":
68     case="expandedData"
69     elif case1=="Symbol":
70     case="Symbol"
71     else:
72     raise ValueError,"unknown case1=%s"%case1
73     elif case0=="array":
74     if case1=="float":
75     case="array"
76     elif case1=="array":
77     case="array"
78     elif case1=="constData":
79     case="constData"
80     elif case1=="taggedData":
81     case="taggedData"
82     elif case1=="expandedData":
83     case="expandedData"
84     elif case1=="Symbol":
85     case="Symbol"
86     else:
87     raise ValueError,"unknown case1=%s"%case1
88     elif case0=="constData":
89     if case1=="float":
90     case="constData"
91     elif case1=="array":
92     case="constData"
93     elif case1=="constData":
94     case="constData"
95     elif case1=="taggedData":
96     case="taggedData"
97     elif case1=="expandedData":
98     case="expandedData"
99     elif case1=="Symbol":
100     case="Symbol"
101     else:
102     raise ValueError,"unknown case1=%s"%case1
103     elif case0=="taggedData":
104     if case1=="float":
105     case="taggedData"
106     elif case1=="array":
107     case="taggedData"
108     elif case1=="constData":
109     case="taggedData"
110     elif case1=="taggedData":
111     case="taggedData"
112     elif case1=="expandedData":
113     case="expandedData"
114     elif case1=="Symbol":
115     case="Symbol"
116     else:
117     raise ValueError,"unknown case1=%s"%case1
118     elif case0=="expandedData":
119     if case1=="float":
120     case="expandedData"
121     elif case1=="array":
122     case="expandedData"
123     elif case1=="constData":
124     case="expandedData"
125     elif case1=="taggedData":
126     case="expandedData"
127     elif case1=="expandedData":
128     case="expandedData"
129     elif case1=="Symbol":
130     case="Symbol"
131     else:
132     raise ValueError,"unknown case1=%s"%case1
133     elif case0=="Symbol":
134     if case1=="float":
135     case="Symbol"
136     elif case1=="array":
137     case="Symbol"
138     elif case1=="constData":
139     case="Symbol"
140     elif case1=="taggedData":
141     case="Symbol"
142     elif case1=="expandedData":
143     case="Symbol"
144     elif case1=="Symbol":
145     case="Symbol"
146     else:
147     raise ValueError,"unknown case1=%s"%case1
148     else:
149     raise ValueError,"unknown case0=%s"%case0
150     return case
151 jgs 154
152 gross 157
153 jgs 154 def makeArray(shape,rng):
154     l=rng[1]-rng[0]
155     out=numarray.zeros(shape,numarray.Float64)
156 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 785 #=================================================
584     def testMatrixMult(arg0,arg1,sh_s):
585     return numarray.matrixmultiply(arg0,arg1)
586 gross 587
587 gross 785 def shapeTestMatrixMult(sh0,sh1,sh_s):
588     if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2):
589     return sh0+sh_s,sh_s+sh1
590 gross 587
591 gross 785 def testTensorMult(arg0,arg1,sh_s):
592     if arg0.rank==2:
593     return numarray.matrixmultiply(arg0,arg1)
594     else:
595     if arg1.rank==4:
596     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
597     for i0 in range(arg0.shape[0]):
598     for i1 in range(arg0.shape[1]):
599     for i2 in range(arg0.shape[2]):
600     for i3 in range(arg0.shape[3]):
601     for j2 in range(arg1.shape[2]):
602     for j3 in range(arg1.shape[3]):
603     out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
604     elif arg1.rank==3:
605     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
606     for i0 in range(arg0.shape[0]):
607     for i1 in range(arg0.shape[1]):
608     for i2 in range(arg0.shape[2]):
609     for i3 in range(arg0.shape[3]):
610     for j2 in range(arg1.shape[2]):
611     out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
612     elif arg1.rank==2:
613     out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
614     for i0 in range(arg0.shape[0]):
615     for i1 in range(arg0.shape[1]):
616     for i2 in range(arg0.shape[2]):
617     for i3 in range(arg0.shape[3]):
618     out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
619     return out
620    
621     def shapeTestTensorMult(sh0,sh1,sh_s):
622     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)):
623     return sh0+sh_s,sh_s+sh1
624    
625     def tensorProductTest(arg0,arg1,sh_s):
626     if isinstance(arg0,float):
627     out=numarray.array(arg0*arg1)
628     elif isinstance(arg1,float):
629     out=numarray.array(arg0*arg1)
630     elif len(sh_s)==0:
631     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
632     else:
633     l=len(sh_s)
634     sh0=arg0.shape[:arg0.rank-l]
635     sh1=arg1.shape[l:]
636     ls,l0,l1=1,1,1
637     for i in sh_s: ls*=i
638     for i in sh0: l0*=i
639     for i in sh1: l1*=i
640     out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
641     out2=numarray.zeros((l0,l1),numarray.Float)
642     for i0 in range(l0):
643     for i1 in range(l1):
644     for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
645     out=out2.resize(sh0+sh1)
646     return out
647     def tensorProductShape(sh0,sh1,sh_s):
648     return sh0+sh_s,sh_s+sh1
649    
650     def transpose(r,offset):
651     if isinstance(r,float): return r
652     s=r.shape
653     s1=1
654     for i in s[:offset]: s1*=i
655     s2=1
656     for i in s[offset:]: s2*=i
657     out=numarray.reshape(r,(s1,s2))
658     out.transpose()
659     return numarray.resize(out,s[offset:]+s[:offset])
660    
661     def shapeTestTMatrixMult(sh0,sh1,sh_s):
662     if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2):
663     return sh_s+sh0,sh_s+sh1
664    
665     def testTMatrixMult(arg0,arg1,sh_s):
666     return numarray.matrixmultiply(transpose(arg0,1),arg1)
667    
668     def shapeTestTTensorMult(sh0,sh1,sh_s):
669     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)):
670     return sh_s+sh0,sh_s+sh1
671    
672     def testTTensorMult(arg0,arg1,sh_s):
673     if arg0.rank==2:
674     return numarray.matrixmultiply(transpose(arg0,1),arg1)
675     else:
676     if arg1.rank==4:
677     out=numarray.zeros((arg0.shape[2],arg0.shape[3],arg1.shape[2],arg1.shape[3]),numarray.Float)
678     for i0 in range(arg0.shape[2]):
679     for i1 in range(arg0.shape[3]):
680     for i2 in range(arg0.shape[0]):
681     for i3 in range(arg0.shape[1]):
682     for j2 in range(arg1.shape[2]):
683     for j3 in range(arg1.shape[3]):
684     out[i0,i1,j2,j3]+=arg0[i2,i3,i0,i1]*arg1[i2,i3,j2,j3]
685     elif arg1.rank==3:
686     out=numarray.zeros((arg0.shape[2],arg0.shape[3],arg1.shape[2]),numarray.Float)
687     for i0 in range(arg0.shape[2]):
688     for i1 in range(arg0.shape[3]):
689     for i2 in range(arg0.shape[0]):
690     for i3 in range(arg0.shape[1]):
691     for j2 in range(arg1.shape[2]):
692     out[i0,i1,j2]+=arg0[i2,i3,i0,i1]*arg1[i2,i3,j2]
693     elif arg1.rank==2:
694     out=numarray.zeros((arg0.shape[2],arg0.shape[3]),numarray.Float)
695     for i0 in range(arg0.shape[2]):
696     for i1 in range(arg0.shape[3]):
697     for i2 in range(arg0.shape[0]):
698     for i3 in range(arg0.shape[1]):
699     out[i0,i1]+=arg0[i2,i3,i0,i1]*arg1[i2,i3]
700     return out
701    
702     def tensorTProductShape(sh0,sh1,sh_s):
703     return sh_s+sh0,sh_s+sh1
704    
705     def tensorTProductTest(arg0,arg1,sh_s):
706     if isinstance(arg0,float):
707     out=numarray.array(arg0*arg1)
708     elif isinstance(arg1,float):
709     out=numarray.array(arg0*arg1)
710     elif len(sh_s)==0:
711     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
712     else:
713     l=len(sh_s)
714     sh0=arg0.shape[l:]
715     sh1=arg1.shape[l:]
716     ls,l0,l1=1,1,1
717     for i in sh_s: ls*=i
718     for i in sh0: l0*=i
719     for i in sh1: l1*=i
720     out1=numarray.outerproduct(arg0,arg1).resize((ls,l0,ls,l1))
721     out2=numarray.zeros((l0,l1),numarray.Float)
722     for i0 in range(l0):
723     for i1 in range(l1):
724     for i in range(ls): out2[i0,i1]+=out1[i,i0,i,i1]
725     out=out2.resize(sh0+sh1)
726     return out
727    
728     def shapeTestMatrixTMult(sh0,sh1,sh_s):
729     if len(sh_s)==1 and len(sh0+sh_s)==2 and len(sh_s+sh1)==2:
730     return sh0+sh_s,sh1+sh_s
731    
732     def testMatrixTMult(arg0,arg1,sh_s):
733     return numarray.matrixmultiply(arg0,transpose(arg1,1))
734    
735     def shapeTestTensorTMult(sh0,sh1,sh_s):
736     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)):
737     return sh0+sh_s,sh1+sh_s
738    
739     def testTensorTMult(arg0,arg1,sh_s):
740     if arg0.rank==2:
741     return numarray.matrixmultiply(arg0,transpose(arg1,1))
742     else:
743     if arg1.rank==4:
744     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[0],arg1.shape[1]),numarray.Float)
745     for i0 in range(arg0.shape[0]):
746     for i1 in range(arg0.shape[1]):
747     for i2 in range(arg0.shape[2]):
748     for i3 in range(arg0.shape[3]):
749     for j0 in range(arg1.shape[0]):
750     for j1 in range(arg1.shape[1]):
751     out[i0,i1,j0,j1]+=arg0[i0,i1,i2,i3]*arg1[j0,j1,i2,i3]
752     elif arg1.rank==3:
753     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[0]),numarray.Float)
754     for i0 in range(arg0.shape[0]):
755     for i1 in range(arg0.shape[1]):
756     for i2 in range(arg0.shape[2]):
757     for i3 in range(arg0.shape[3]):
758     for j0 in range(arg1.shape[0]):
759     out[i0,i1,j0]+=arg0[i0,i1,i2,i3]*arg1[j0,i2,i3]
760     elif arg1.rank==2:
761     out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
762     for i0 in range(arg0.shape[0]):
763     for i1 in range(arg0.shape[1]):
764     for i2 in range(arg0.shape[2]):
765     for i3 in range(arg0.shape[3]):
766     out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
767     return out
768    
769     def tensorProductTTest(arg0,arg1,sh_s):
770     if isinstance(arg0,float):
771     out=numarray.array(arg0*arg1)
772     elif isinstance(arg1,float):
773     out=numarray.array(arg0*arg1)
774     elif len(sh_s)==0:
775     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
776     else:
777     l=len(sh_s)
778     sh0=arg0.shape[:arg0.rank-l]
779     sh1=arg1.shape[:arg1.rank-l]
780     ls,l0,l1=1,1,1
781     for i in sh_s: ls*=i
782     for i in sh0: l0*=i
783     for i in sh1: l1*=i
784     out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,l1,ls))
785     out2=numarray.zeros((l0,l1),numarray.Float)
786     for i0 in range(l0):
787     for i1 in range(l1):
788     for i in range(ls): out2[i0,i1]+=out1[i0,i,i1,i]
789     out=out2.resize(sh0+sh1)
790     return out
791    
792     def tensorProductTShape(sh0,sh1,sh_s):
793     return sh0+sh_s,sh1+sh_s
794    
795 gross 530 #=======================================================================================================
796 gross 785 # tensor multiply
797     #=======================================================================================================
798     # oper=
799     # oper=
800     for oper in [ ["tensor_mult",testTensorMult,shapeTestTensorMult], \
801     ["generalTensorProduct",tensorProductTest,tensorProductShape], \
802     ["matrix_mult",testMatrixMult,shapeTestMatrixMult], \
803     ["transposed_matrix_mult",testTMatrixMult,shapeTestTMatrixMult], \
804     ["transposed_tensor_mult",testTTensorMult,shapeTestTTensorMult], \
805     ["generalTransposedTensorProduct",tensorTProductTest,tensorTProductShape], \
806     ["matrix_transposed_mult",testMatrixTMult,shapeTestMatrixTMult], \
807     ["tensor_transposed_mult",testTensorTMult,shapeTestTensorTMult], \
808     ["generalTensorTransposedProduct",tensorProductTTest,tensorProductTShape] ]:
809    
810     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
811     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
812     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
813     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
814     for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
815     if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
816     and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
817     arg_shape=oper[2](sh0,sh1,sh_s)
818     if not arg_shape==None:
819     case=getResultCaseForBin(case0,case1)
820     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
821     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
822     if oper[0] in [ "generalTensorProduct", "generalTransposedTensorProduct", "generalTensorTransposedProduct"]:
823     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))
824    
825     res_text=" res=%s(arg0,arg1,axis_offset=%s)\n"%(oper[0],len(sh_s))
826     else:
827     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0+sh_s),case1,len(sh_s+sh1))
828    
829     res_text=" res=%s(arg0,arg1)\n"%oper[0]
830    
831     text+=" def %s(self):\n"%tname
832     a_0=makeArray(arg_shape[0],[-8,8])
833     if case0 in ["taggedData", "expandedData"]:
834     a1_0=makeArray(arg_shape[0],[-8,8])
835     else:
836     a1_0=a_0
837    
838     a_1=makeArray(arg_shape[1],[-8,8])
839     if case1 in ["taggedData", "expandedData"]:
840     a1_1=makeArray(arg_shape[1],[-8,8])
841     else:
842     a1_1=a_1
843     r=oper[1](a_0,a_1,sh_s)
844     r1=oper[1](a1_0,a1_1,sh_s)
845     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
846     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
847     text+=res_text
848     if case=="Symbol":
849     c0_res,c1_res=case0,case1
850     subs="{"
851     if case0=="Symbol":
852     text+=mkText("array","s0",a_0,a1_0)
853     subs+="arg0:s0"
854     c0_res="array"
855     if case1=="Symbol":
856     text+=mkText("array","s1",a_1,a1_1)
857     if not subs.endswith("{"): subs+=","
858     subs+="arg1:s1"
859     c1_res="array"
860     subs+="}"
861     text+=" sub=res.substitute(%s)\n"%subs
862     res="sub"
863     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
864     else:
865     res="res"
866     text+=mkText(case,"ref",r,r1)
867     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
868     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
869     if case0 == "taggedData" or case1 == "taggedData":
870     t_prog_with_tags+=text
871     else:
872     t_prog+=text
873     print test_header
874     print t_prog
875     # print t_prog_with_tags
876     print test_tail
877     1/0
878    
879    
880    
881     #=======================================================================================================
882 gross 587 # eigenvalues and eigen vectors 2D:
883     #=======================================================================================================
884     alltests= \
885     [ ("case0",[[0.0, 0.0], [0.0, 0.0]],[0.0, 0.0]) \
886     , ("case3",[[-1.0, 0.0], [0.0, -1.0]],[-1.0, -1.0]) \
887     , ("case5",[[-0.99999999999999967, -6.4606252205695602e-16], [-6.4606252205695602e-16, -0.99999999999999967]],[-1.0, -1.0]) \
888     , ("case6",[[0.0, 0.0], [0.0, 0.0001]],[0.0, 0.0001]) \
889     , ("case7",[[0.0001, 0.0], [0.0, 0.0]],[0.0, 0.0001]) \
890     , ("case8",[[6.0598371831785722e-06, 2.3859213977648625e-05], [2.3859213977648629e-05, 9.3940162816821425e-05]],[0.0, 0.0001]) \
891     , ("case9",[[1.0, 0.0], [0.0, 2.0]],[1.0, 2.0]) \
892     , ("case10",[[2.0, 0.0], [0.0, 1.0]],[1.0, 2.0]) \
893     , ("case11",[[1.0605983718317855, 0.23859213977648688], [0.23859213977648688, 1.9394016281682138]],[1.0, 2.0]) \
894     , ("case12",[[1.0, 0.0], [0.0, 1000000.0]],[1.0, 1000000.0]) \
895     , ("case13",[[1000000.0, 0.0], [0.0, 1.0]],[1.0, 1000000.0]) \
896     , ("case14",[[60599.311233413886, 238591.90118434647], [238591.90118434647, 939401.68876658613]],[1.0, 1000000.0]) \
897     ]
898     dim=2
899 gross 588
900    
901     alltests= \
902     [ ("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]) \
903     , ("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]) \
904     , ("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]) \
905     , ("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]) \
906     , ("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]) \
907     , ("case13",[[1.0, 0.0, 0.0], [0.0, 0.9, 0.], [0.0, 0., 0.9]],[0.9, 0.9, 1.0]) \
908     , ("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]) \
909     , ("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]) \
910     , ("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]) \
911     , ("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]) \
912     , ("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]) \
913     , ("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]) \
914     , ("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]) \
915     , ("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]) \
916     , ("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]) \
917     , ("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]) \
918     , ("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]) \
919     , ("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]) \
920     , ("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]) \
921     , ("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]) \
922     , ("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]) \
923     ]
924    
925     dim=3
926    
927     alltests= \
928     [ ("case0",[[0.0]],[0.0]) \
929     , ("case1",[[1.0]],[1.0]) \
930     ]
931     dim=1
932    
933 gross 587 for case in ["constData","taggedData","expandedData"]:
934     n=0
935     while n<len(alltests):
936     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
937     tname="test_eigenvalues_and_eigenvectors_%s_dim%s_%s"%(case,dim,alltests[n][0])
938     text+=" def %s(self):\n"%tname
939     a_0=numarray.array(alltests[n][1],numarray.Float64)
940     ev_0=numarray.array(alltests[n][2],numarray.Float64)
941     if case in ["taggedData", "expandedData"]:
942 gross 588 if n+1<len(alltests):
943     a1_0=numarray.array(alltests[n+1][1],numarray.Float64)
944     ev1_0=numarray.array(alltests[n+1][2],numarray.Float64)
945     else:
946     a1_0=numarray.array(alltests[0][1],numarray.Float64)
947     ev1_0=numarray.array(alltests[0][2],numarray.Float64)
948 gross 587 n+=2
949     else:
950     a1_0=a_0
951     ev1_0=ev_0
952     n+=1
953     text+=mkText(case,"arg",a_0,a1_0)
954     text+=" res=eigenvalues_and_eigenvectors(arg)\n"
955     text+=mkText(case,"ref_ev",ev_0,ev1_0)
956     text+=mkTypeAndShapeTest(case,(dim,),"res[0]"," for eigenvalues")
957     text+=mkTypeAndShapeTest(case,(dim,dim),"res[1]"," for eigenvectors")
958     text+=" self.failUnless(Lsup(res[0]-ref_ev)<=self.RES_TOL*Lsup(ref_ev),\"wrong eigenvalues\")\n"
959     for i in range(dim):
960     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)
961     text+=" self.failUnless(Lsup(matrixmult(transpose(res[1]),res[1])-kronecker(%s))<=self.RES_TOL,\"eigenvectors are not orthonormal\")\n"%dim
962     if case == "taggedData" :
963     t_prog_with_tags+=text
964     else:
965     t_prog+=text
966     print test_header
967     print t_prog
968     print t_prog_with_tags
969     print test_tail
970     1/0
971     #=======================================================================================================
972 gross 550 # get slices
973 gross 536 #=======================================================================================================
974     from esys.escript import *
975 gross 550 for case0 in ["constData","taggedData","expandedData"]:
976     for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,3,3)]:
977     # get perm:
978     if len(sh0)==2:
979     check=[[1,0]]
980     elif len(sh0)==3:
981     check=[[1,0,2],
982     [1,2,0],
983     [2,1,0],
984     [2,0,2],
985     [0,2,1]]
986     elif len(sh0)==4:
987     check=[[0,1,3,2],
988     [0,2,1,3],
989     [0,2,3,1],
990     [0,3,2,1],
991     [0,3,1,2] ,
992     [1,0,2,3],
993     [1,0,3,2],
994     [1,2,0,3],
995     [1,2,3,0],
996     [1,3,2,0],
997     [1,3,0,2],
998     [2,0,1,3],
999     [2,0,3,1],
1000     [2,1,0,3],
1001     [2,1,3,0],
1002     [2,3,1,0],
1003     [2,3,0,1],
1004     [3,0,1,2],
1005     [3,0,2,1],
1006     [3,1,0,2],
1007     [3,1,2,0],
1008     [3,2,1,0],
1009     [3,2,0,1]]
1010     else:
1011     check=[]
1012    
1013     # create the test cases:
1014     processed=[]
1015     l=["R","U","L","P","C","N"]
1016     c=[""]
1017     for i in range(len(sh0)):
1018     tmp=[]
1019     for ci in c:
1020     tmp+=[ci+li for li in l]
1021     c=tmp
1022     # SHUFFLE
1023     c2=[]
1024     while len(c)>0:
1025     i=int(random.random()*len(c))
1026     c2.append(c[i])
1027     del c[i]
1028     c=c2
1029     for ci in c:
1030     t=""
1031     sh=()
1032     sl=()
1033     for i in range(len(ci)):
1034     if ci[i]=="R":
1035     s="%s:%s"%(1,sh0[i]-1)
1036     sl=sl+(slice(1,sh0[i]-1),)
1037     sh=sh+(sh0[i]-2,)
1038     if ci[i]=="U":
1039     s=":%s"%(sh0[i]-1)
1040     sh=sh+(sh0[i]-1,)
1041     sl=sl+(slice(0,sh0[i]-1),)
1042     if ci[i]=="L":
1043     s="2:"
1044     sh=sh+(sh0[i]-2,)
1045     sl=sl+(slice(2,sh0[i]),)
1046     if ci[i]=="P":
1047     s="%s"%(int(sh0[i]/2))
1048     sl=sl+(int(sh0[i]/2),)
1049     if ci[i]=="C":
1050     s=":"
1051     sh=sh+(sh0[i],)
1052     sl=sl+(slice(0,sh0[i]),)
1053     if ci[i]=="N":
1054     s=""
1055     sh=sh+(sh0[i],)
1056     if len(s)>0:
1057     if not t=="": t+=","
1058     t+=s
1059     if len(sl)==1: sl=sl[0]
1060     N_found=False
1061     noN_found=False
1062     process=len(t)>0
1063     for i in ci:
1064     if i=="N":
1065     if not noN_found and N_found: process=False
1066     N_found=True
1067     else:
1068     if N_found: process=False
1069     noNfound=True
1070     # is there a similar one processed allready
1071     if process and ci.find("N")==-1:
1072     for ci2 in processed:
1073     for chi in check:
1074     is_perm=True
1075     for i in range(len(chi)):
1076     if not ci[i]==ci2[chi[i]]: is_perm=False
1077     if is_perm: process=False
1078     # if not process: print ci," rejected"
1079     if process:
1080     processed.append(ci)
1081     for case1 in ["array","constData","taggedData","expandedData"]:
1082     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1083     tname="test_setslice_%s_rank%s_%s_%s"%(case0,len(sh0),case1,ci)
1084     text+=" def %s(self):\n"%tname
1085     a_0=makeNumberedArray(sh0)
1086     if case0 in ["taggedData", "expandedData"]:
1087     a1_0=makeNumberedArray(sh0)
1088     else:
1089     a1_0=a_0*1.
1090    
1091     a_1=makeNumberedArray(sh)
1092     if case1 in ["taggedData", "expandedData"]:
1093     a1_1=makeNumberedArray(sh)
1094     else:
1095     a1_1=a_1*1.
1096    
1097     text+=mkText(case0,"arg",a_0,a1_0)
1098     text+=mkText(case1,"val",a_1,a1_1)
1099     text+=" arg[%s]=val\n"%t
1100     a_0.__setitem__(sl,a_1)
1101     a1_0.__setitem__(sl,a1_1)
1102     if Lsup(a_0-a1_0)==0.:
1103     text+=mkText("constData","ref",a_0,a1_0)
1104     else:
1105     text+=mkText("expandedData","ref",a_0,a1_0)
1106     text+=" self.failUnless(Lsup(arg-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"
1107    
1108     if case0 == "taggedData" or case1 == "taggedData":
1109     t_prog_with_tags+=text
1110     else:
1111     t_prog+=text
1112    
1113     print test_header
1114     # print t_prog
1115     print t_prog_with_tags
1116     print test_tail
1117     1/0
1118    
1119     #=======================================================================================================
1120     # (non)symmetric part
1121     #=======================================================================================================
1122     from esys.escript import *
1123 gross 536 for name in ["symmetric", "nonsymmetric"]:
1124     f=1.
1125     if name=="nonsymmetric": f=-1
1126     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1127     for sh0 in [ (3,3), (2,3,2,3)]:
1128     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1129     tname="test_%s_%s_rank%s"%(name,case0,len(sh0))
1130     text+=" def %s(self):\n"%tname
1131     a_0=makeNumberedArray(sh0,s=1.)
1132     r_0=(a_0+f*transpose(a_0))/2.
1133     if case0 in ["taggedData", "expandedData"]:
1134     a1_0=makeNumberedArray(sh0,s=-1.)
1135     r1_0=(a1_0+f*transpose(a1_0))/2.
1136     else:
1137     a1_0=a_0
1138     r1_0=r_0
1139     text+=mkText(case0,"arg",a_0,a1_0)
1140     text+=" res=%s(arg)\n"%name
1141     if case0=="Symbol":
1142     text+=mkText("array","s",a_0,a1_0)
1143     text+=" sub=res.substitute({arg:s})\n"
1144     res="sub"
1145     text+=mkText("array","ref",r_0,r1_0)
1146     else:
1147     res="res"
1148     text+=mkText(case0,"ref",r_0,r1_0)
1149     text+=mkTypeAndShapeTest(case0,sh0,"res")
1150     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1151    
1152     if case0 == "taggedData" :
1153     t_prog_with_tags+=text
1154     else:
1155     t_prog+=text
1156     print test_header
1157     print t_prog
1158     # print t_prog_with_tags
1159     print test_tail
1160     1/0
1161    
1162     #=======================================================================================================
1163 gross 530 # eigenvalues
1164     #=======================================================================================================
1165     import numarray.linear_algebra
1166     name="eigenvalues"
1167     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1168     for sh0 in [ (1,1), (2,2), (3,3)]:
1169     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1170     tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1171     text+=" def %s(self):\n"%tname
1172     a_0=makeArray(sh0,[-1.,1])
1173     a_0=(a_0+numarray.transpose(a_0))/2.
1174     ev=numarray.linear_algebra.eigenvalues(a_0)
1175     ev.sort()
1176     if case0 in ["taggedData", "expandedData"]:
1177     a1_0=makeArray(sh0,[-1.,1])
1178     a1_0=(a1_0+numarray.transpose(a1_0))/2.
1179     ev1=numarray.linear_algebra.eigenvalues(a1_0)
1180     ev1.sort()
1181     else:
1182     a1_0=a_0
1183     ev1=ev
1184     text+=mkText(case0,"arg",a_0,a1_0)
1185     text+=" res=%s(arg)\n"%name
1186     if case0=="Symbol":
1187     text+=mkText("array","s",a_0,a1_0)
1188     text+=" sub=res.substitute({arg:s})\n"
1189     res="sub"
1190     text+=mkText("array","ref",ev,ev1)
1191     else:
1192     res="res"
1193     text+=mkText(case0,"ref",ev,ev1)
1194     text+=mkTypeAndShapeTest(case0,(sh0[0],),"res")
1195     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1196    
1197     if case0 == "taggedData" :
1198     t_prog_with_tags+=text
1199     else:
1200     t_prog+=text
1201     print test_header
1202     # print t_prog
1203     print t_prog_with_tags
1204     print test_tail
1205     1/0
1206 jgs 154
1207 gross 517 #=======================================================================================================
1208 gross 550 # get slices
1209 gross 517 #=======================================================================================================
1210     for case0 in ["constData","taggedData","expandedData","Symbol"]:
1211     for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,5,3)]:
1212     # get perm:
1213     if len(sh0)==2:
1214     check=[[1,0]]
1215     elif len(sh0)==3:
1216     check=[[1,0,2],
1217     [1,2,0],
1218     [2,1,0],
1219     [2,0,2],
1220     [0,2,1]]
1221     elif len(sh0)==4:
1222     check=[[0,1,3,2],
1223     [0,2,1,3],
1224     [0,2,3,1],
1225     [0,3,2,1],
1226     [0,3,1,2] ,
1227     [1,0,2,3],
1228     [1,0,3,2],
1229     [1,2,0,3],
1230     [1,2,3,0],
1231     [1,3,2,0],
1232     [1,3,0,2],
1233     [2,0,1,3],
1234     [2,0,3,1],
1235     [2,1,0,3],
1236     [2,1,3,0],
1237     [2,3,1,0],
1238     [2,3,0,1],
1239     [3,0,1,2],
1240     [3,0,2,1],
1241     [3,1,0,2],
1242     [3,1,2,0],
1243     [3,2,1,0],
1244     [3,2,0,1]]
1245     else:
1246     check=[]
1247    
1248     # create the test cases:
1249     processed=[]
1250     l=["R","U","L","P","C","N"]
1251     c=[""]
1252     for i in range(len(sh0)):
1253     tmp=[]
1254     for ci in c:
1255     tmp+=[ci+li for li in l]
1256     c=tmp
1257     # SHUFFLE
1258     c2=[]
1259     while len(c)>0:
1260     i=int(random.random()*len(c))
1261     c2.append(c[i])
1262     del c[i]
1263     c=c2
1264     for ci in c:
1265     t=""
1266     sh=()
1267     for i in range(len(ci)):
1268     if ci[i]=="R":
1269     s="%s:%s"%(1,sh0[i]-1)
1270     sh=sh+(sh0[i]-2,)
1271     if ci[i]=="U":
1272     s=":%s"%(sh0[i]-1)
1273     sh=sh+(sh0[i]-1,)
1274     if ci[i]=="L":
1275     s="2:"
1276     sh=sh+(sh0[i]-2,)
1277     if ci[i]=="P":
1278     s="%s"%(int(sh0[i]/2))
1279     if ci[i]=="C":
1280     s=":"
1281     sh=sh+(sh0[i],)
1282     if ci[i]=="N":
1283     s=""
1284     sh=sh+(sh0[i],)
1285     if len(s)>0:
1286     if not t=="": t+=","
1287     t+=s
1288     N_found=False
1289     noN_found=False
1290     process=len(t)>0
1291     for i in ci:
1292     if i=="N":
1293     if not noN_found and N_found: process=False
1294     N_found=True
1295     else:
1296     if N_found: process=False
1297     noNfound=True
1298     # is there a similar one processed allready
1299     if process and ci.find("N")==-1:
1300     for ci2 in processed:
1301     for chi in check:
1302     is_perm=True
1303     for i in range(len(chi)):
1304     if not ci[i]==ci2[chi[i]]: is_perm=False
1305     if is_perm: process=False
1306     # if not process: print ci," rejected"
1307     if process:
1308     processed.append(ci)
1309     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1310     tname="test_getslice_%s_rank%s_%s"%(case0,len(sh0),ci)
1311     text+=" def %s(self):\n"%tname
1312     a_0=makeNumberedArray(sh0,s=1)
1313     if case0 in ["taggedData", "expandedData"]:
1314     a1_0=makeNumberedArray(sh0,s=-1.)
1315     else:
1316     a1_0=a_0
1317     r=eval("a_0[%s]"%t)
1318     r1=eval("a1_0[%s]"%t)
1319     text+=mkText(case0,"arg",a_0,a1_0)
1320     text+=" res=arg[%s]\n"%t
1321     if case0=="Symbol":
1322     text+=mkText("array","s",a_0,a1_0)
1323     text+=" sub=res.substitute({arg:s})\n"
1324     res="sub"
1325     text+=mkText("array","ref",r,r1)
1326     else:
1327     res="res"
1328     text+=mkText(case0,"ref",r,r1)
1329     text+=mkTypeAndShapeTest(case0,sh,"res")
1330     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1331    
1332     if case0 == "taggedData" :
1333     t_prog_with_tags+=text
1334     else:
1335     t_prog+=text
1336    
1337     print test_header
1338     # print t_prog
1339     print t_prog_with_tags
1340     print test_tail
1341     1/0
1342     #============================================================================================
1343 gross 291 def innerTEST(arg0,arg1):
1344     if isinstance(arg0,float):
1345     out=numarray.array(arg0*arg1)
1346     else:
1347     out=(arg0*arg1).sum()
1348     return out
1349    
1350     def outerTEST(arg0,arg1):
1351     if isinstance(arg0,float):
1352     out=numarray.array(arg0*arg1)
1353     elif isinstance(arg1,float):
1354     out=numarray.array(arg0*arg1)
1355     else:
1356     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1357     return out
1358    
1359     def tensorProductTest(arg0,arg1,sh_s):
1360     if isinstance(arg0,float):
1361     out=numarray.array(arg0*arg1)
1362     elif isinstance(arg1,float):
1363     out=numarray.array(arg0*arg1)
1364     elif len(sh_s)==0:
1365     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1366     else:
1367     l=len(sh_s)
1368     sh0=arg0.shape[:arg0.rank-l]
1369     sh1=arg1.shape[l:]
1370     ls,l0,l1=1,1,1
1371     for i in sh_s: ls*=i
1372     for i in sh0: l0*=i
1373     for i in sh1: l1*=i
1374     out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
1375     out2=numarray.zeros((l0,l1),numarray.Float)
1376     for i0 in range(l0):
1377     for i1 in range(l1):
1378     for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
1379     out=out2.resize(sh0+sh1)
1380     return out
1381    
1382     def testMatrixMult(arg0,arg1,sh_s):
1383     return numarray.matrixmultiply(arg0,arg1)
1384    
1385    
1386     def testTensorMult(arg0,arg1,sh_s):
1387     if len(arg0)==2:
1388     return numarray.matrixmultiply(arg0,arg1)
1389     else:
1390     if arg1.rank==4:
1391     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
1392     for i0 in range(arg0.shape[0]):
1393     for i1 in range(arg0.shape[1]):
1394     for i2 in range(arg0.shape[2]):
1395     for i3 in range(arg0.shape[3]):
1396     for j2 in range(arg1.shape[2]):
1397     for j3 in range(arg1.shape[3]):
1398     out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
1399     elif arg1.rank==3:
1400     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
1401     for i0 in range(arg0.shape[0]):
1402     for i1 in range(arg0.shape[1]):
1403     for i2 in range(arg0.shape[2]):
1404     for i3 in range(arg0.shape[3]):
1405     for j2 in range(arg1.shape[2]):
1406     out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
1407     elif arg1.rank==2:
1408     out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
1409     for i0 in range(arg0.shape[0]):
1410     for i1 in range(arg0.shape[1]):
1411     for i2 in range(arg0.shape[2]):
1412     for i3 in range(arg0.shape[3]):
1413     out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
1414     return out
1415 gross 313
1416     def testReduce(arg0,init_val,test_expr,post_expr):
1417     out=init_val
1418     if isinstance(arg0,float):
1419     out=eval(test_expr.replace("%a1%","arg0"))
1420     elif arg0.rank==0:
1421     out=eval(test_expr.replace("%a1%","arg0"))
1422     elif arg0.rank==1:
1423     for i0 in range(arg0.shape[0]):
1424     out=eval(test_expr.replace("%a1%","arg0[i0]"))
1425     elif arg0.rank==2:
1426     for i0 in range(arg0.shape[0]):
1427     for i1 in range(arg0.shape[1]):
1428     out=eval(test_expr.replace("%a1%","arg0[i0,i1]"))
1429     elif arg0.rank==3:
1430     for i0 in range(arg0.shape[0]):
1431     for i1 in range(arg0.shape[1]):
1432     for i2 in range(arg0.shape[2]):
1433     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2]"))
1434     elif arg0.rank==4:
1435     for i0 in range(arg0.shape[0]):
1436     for i1 in range(arg0.shape[1]):
1437     for i2 in range(arg0.shape[2]):
1438     for i3 in range(arg0.shape[3]):
1439     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))
1440     return eval(post_expr)
1441 gross 396
1442     def clipTEST(arg0,mn,mx):
1443     if isinstance(arg0,float):
1444     return max(min(arg0,mx),mn)
1445     out=numarray.zeros(arg0.shape,numarray.Float64)
1446     if arg0.rank==1:
1447     for i0 in range(arg0.shape[0]):
1448     out[i0]=max(min(arg0[i0],mx),mn)
1449     elif arg0.rank==2:
1450     for i0 in range(arg0.shape[0]):
1451     for i1 in range(arg0.shape[1]):
1452     out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
1453     elif arg0.rank==3:
1454     for i0 in range(arg0.shape[0]):
1455     for i1 in range(arg0.shape[1]):
1456     for i2 in range(arg0.shape[2]):
1457     out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
1458     elif arg0.rank==4:
1459     for i0 in range(arg0.shape[0]):
1460     for i1 in range(arg0.shape[1]):
1461     for i2 in range(arg0.shape[2]):
1462     for i3 in range(arg0.shape[3]):
1463     out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
1464     return out
1465     def minimumTEST(arg0,arg1):
1466     if isinstance(arg0,float):
1467     if isinstance(arg1,float):
1468     if arg0>arg1:
1469     return arg1
1470     else:
1471     return arg0
1472     else:
1473     arg0=numarray.ones(arg1.shape)*arg0
1474     else:
1475     if isinstance(arg1,float):
1476     arg1=numarray.ones(arg0.shape)*arg1
1477     out=numarray.zeros(arg0.shape,numarray.Float64)
1478     if arg0.rank==0:
1479     if arg0>arg1:
1480     out=arg1
1481     else:
1482     out=arg0
1483     elif arg0.rank==1:
1484     for i0 in range(arg0.shape[0]):
1485     if arg0[i0]>arg1[i0]:
1486     out[i0]=arg1[i0]
1487     else:
1488     out[i0]=arg0[i0]
1489     elif arg0.rank==2:
1490     for i0 in range(arg0.shape[0]):
1491     for i1 in range(arg0.shape[1]):
1492     if arg0[i0,i1]>arg1[i0,i1]:
1493     out[i0,i1]=arg1[i0,i1]
1494     else:
1495     out[i0,i1]=arg0[i0,i1]
1496     elif arg0.rank==3:
1497     for i0 in range(arg0.shape[0]):
1498     for i1 in range(arg0.shape[1]):
1499     for i2 in range(arg0.shape[2]):
1500     if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
1501     out[i0,i1,i2]=arg1[i0,i1,i2]
1502     else:
1503     out[i0,i1,i2]=arg0[i0,i1,i2]
1504     elif arg0.rank==4:
1505     for i0 in range(arg0.shape[0]):
1506     for i1 in range(arg0.shape[1]):
1507     for i2 in range(arg0.shape[2]):
1508     for i3 in range(arg0.shape[3]):
1509     if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
1510     out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
1511     else:
1512     out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
1513     return out
1514 gross 443
1515 gross 439 def unrollLoops(a,b,o,arg,tap="",x="x"):
1516 gross 437 out=""
1517     if a.rank==1:
1518     z=""
1519     for i99 in range(a.shape[0]):
1520     if not z=="": z+="+"
1521     if o=="1":
1522 gross 439 z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
1523 gross 437 else:
1524 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
1525 gross 437
1526     out+=tap+"%s=%s\n"%(arg,z)
1527    
1528     elif a.rank==2:
1529     for i0 in range(a.shape[0]):
1530     z=""
1531     for i99 in range(a.shape[1]):
1532     if not z=="": z+="+"
1533     if o=="1":
1534     z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
1535     else:
1536 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
1537 gross 437
1538     out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
1539     elif a.rank==3:
1540     for i0 in range(a.shape[0]):
1541     for i1 in range(a.shape[1]):
1542     z=""
1543     for i99 in range(a.shape[2]):
1544     if not z=="": z+="+"
1545     if o=="1":
1546 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
1547 gross 437 else:
1548 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
1549 gross 437
1550     out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
1551     elif a.rank==4:
1552     for i0 in range(a.shape[0]):
1553     for i1 in range(a.shape[1]):
1554     for i2 in range(a.shape[2]):
1555     z=""
1556     for i99 in range(a.shape[3]):
1557     if not z=="": z+="+"
1558     if o=="1":
1559 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
1560 gross 437 else:
1561 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
1562 gross 437
1563     out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
1564     elif a.rank==5:
1565     for i0 in range(a.shape[0]):
1566     for i1 in range(a.shape[1]):
1567     for i2 in range(a.shape[2]):
1568     for i3 in range(a.shape[3]):
1569     z=""
1570     for i99 in range(a.shape[4]):
1571     if not z=="": z+="+"
1572     if o=="1":
1573 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
1574 gross 437 else:
1575 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)
1576 gross 437
1577     out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
1578     return out
1579    
1580     def unrollLoopsOfGrad(a,b,o,arg,tap=""):
1581     out=""
1582     if a.rank==1:
1583     for i99 in range(a.shape[0]):
1584     if o=="1":
1585     out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
1586     else:
1587     out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
1588    
1589     elif a.rank==2:
1590     for i0 in range(a.shape[0]):
1591     for i99 in range(a.shape[1]):
1592     if o=="1":
1593     out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
1594     else:
1595     out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
1596    
1597     elif a.rank==3:
1598     for i0 in range(a.shape[0]):
1599     for i1 in range(a.shape[1]):
1600     for i99 in range(a.shape[2]):
1601     if o=="1":
1602     out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
1603     else:
1604     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])
1605    
1606     elif a.rank==4:
1607     for i0 in range(a.shape[0]):
1608     for i1 in range(a.shape[1]):
1609     for i2 in range(a.shape[2]):
1610     for i99 in range(a.shape[3]):
1611     if o=="1":
1612     out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
1613     else:
1614     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])
1615 gross 438 return out
1616 gross 441 def unrollLoopsOfDiv(a,b,o,arg,tap=""):
1617     out=tap+arg+"="
1618     if o=="1":
1619     z=0.
1620     for i99 in range(a.shape[0]):
1621     z+=b[i99,i99]+a[i99,i99]
1622     out+="(%s)"%z
1623     else:
1624     z=0.
1625     for i99 in range(a.shape[0]):
1626     z+=b[i99,i99]
1627     if i99>0: out+="+"
1628     out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
1629     out+="+(%s)"%z
1630     return out
1631    
1632 gross 437 def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
1633     if where=="Function":
1634     xfac_o=1.
1635     xfac_op=0.
1636     z_fac=1./2.
1637     z_fac_s=""
1638     zo_fac_s="/(o+1.)"
1639     elif where=="FunctionOnBoundary":
1640     xfac_o=1.
1641     xfac_op=0.
1642     z_fac=1.
1643     z_fac_s="*dim"
1644     zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
1645     elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1646     xfac_o=0.
1647     xfac_op=1.
1648     z_fac=1./2.
1649     z_fac_s=""
1650     zo_fac_s="/(o+1.)"
1651     out=""
1652     if a.rank==1:
1653     zo=0.
1654     zop=0.
1655     z=0.
1656     for i99 in range(a.shape[0]):
1657     if i99==0:
1658     zo+= xfac_o*a[i99]
1659     zop+= xfac_op*a[i99]
1660     else:
1661     zo+=a[i99]
1662     z+=b[i99]
1663    
1664     out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
1665     if zop==0.:
1666     out+="\n"
1667     else:
1668     out+="+(%s)*0.5**o\n"%zop
1669     elif a.rank==2:
1670     for i0 in range(a.shape[0]):
1671     zo=0.
1672     zop=0.
1673     z=0.
1674     for i99 in range(a.shape[1]):
1675     if i99==0:
1676     zo+= xfac_o*a[i0,i99]
1677     zop+= xfac_op*a[i0,i99]
1678     else:
1679     zo+=a[i0,i99]
1680     z+=b[i0,i99]
1681    
1682     out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
1683     if zop==0.:
1684     out+="\n"
1685     else:
1686     out+="+(%s)*0.5**o\n"%zop
1687     elif a.rank==3:
1688     for i0 in range(a.shape[0]):
1689     for i1 in range(a.shape[1]):
1690     zo=0.
1691     zop=0.
1692     z=0.
1693     for i99 in range(a.shape[2]):
1694     if i99==0:
1695     zo+= xfac_o*a[i0,i1,i99]
1696     zop+= xfac_op*a[i0,i1,i99]
1697     else:
1698     zo+=a[i0,i1,i99]
1699     z+=b[i0,i1,i99]
1700    
1701     out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
1702     if zop==0.:
1703     out+="\n"
1704     else:
1705     out+="+(%s)*0.5**o\n"%zop
1706     elif a.rank==4:
1707     for i0 in range(a.shape[0]):
1708     for i1 in range(a.shape[1]):
1709     for i2 in range(a.shape[2]):
1710     zo=0.
1711     zop=0.
1712     z=0.
1713     for i99 in range(a.shape[3]):
1714     if i99==0:
1715     zo+= xfac_o*a[i0,i1,i2,i99]
1716     zop+= xfac_op*a[i0,i1,i2,i99]
1717    
1718     else:
1719     zo+=a[i0,i1,i2,i99]
1720     z+=b[i0,i1,i2,i99]
1721    
1722     out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
1723     if zop==0.:
1724     out+="\n"
1725     else:
1726     out+="+(%s)*0.5**o\n"%zop
1727    
1728     elif a.rank==5:
1729     for i0 in range(a.shape[0]):
1730     for i1 in range(a.shape[1]):
1731     for i2 in range(a.shape[2]):
1732     for i3 in range(a.shape[3]):
1733     zo=0.
1734     zop=0.
1735     z=0.
1736     for i99 in range(a.shape[4]):
1737     if i99==0:
1738     zo+= xfac_o*a[i0,i1,i2,i3,i99]
1739     zop+= xfac_op*a[i0,i1,i2,i3,i99]
1740    
1741     else:
1742     zo+=a[i0,i1,i2,i3,i99]
1743     z+=b[i0,i1,i2,i3,i99]
1744     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)
1745     if zop==0.:
1746     out+="\n"
1747     else:
1748     out+="+(%s)*0.5**o\n"%zop
1749    
1750     return out
1751 gross 443 def unrollLoopsSimplified(b,arg,tap=""):
1752     out=""
1753     if isinstance(b,float) or b.rank==0:
1754     out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
1755    
1756     elif b.rank==1:
1757     for i0 in range(b.shape[0]):
1758     out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
1759     elif b.rank==2:
1760     for i0 in range(b.shape[0]):
1761     for i1 in range(b.shape[1]):
1762     out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
1763     elif b.rank==3:
1764     for i0 in range(b.shape[0]):
1765     for i1 in range(b.shape[1]):
1766     for i2 in range(b.shape[2]):
1767     out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
1768     elif b.rank==4:
1769     for i0 in range(b.shape[0]):
1770     for i1 in range(b.shape[1]):
1771     for i2 in range(b.shape[2]):
1772     for i3 in range(b.shape[3]):
1773     out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
1774     return out
1775    
1776     def unrollLoopsOfL2(b,where,arg,tap=""):
1777     out=""
1778     z=[]
1779     if isinstance(<