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

Annotation of /trunk/escript/py_src/generatediff

Parent Directory Parent Directory | Revision Log Revision Log


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