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

Annotation of /trunk/escript/py_src/generatediff

Parent Directory Parent Directory | Revision Log Revision Log


Revision 443 - (hide annotations)
Fri Jan 20 06:22:38 2006 UTC (13 years, 9 months ago) by gross
Original Path: trunk/escript/py_src/generateutil
File size: 128530 byte(s)
L2 function added
1 jgs 154 #!/usr/bin/python
2 gross 283 # $Id$
3 jgs 154
4     """
5     program generates parts of the util.py and the test_util.py script
6     """
7 gross 157 test_header=""
8     test_header+="import unittest\n"
9     test_header+="import numarray\n"
10     test_header+="from esys.escript import *\n"
11     test_header+="from esys.finley import Rectangle\n"
12     test_header+="class Test_util2(unittest.TestCase):\n"
13     test_header+=" RES_TOL=1.e-7\n"
14     test_header+=" def setUp(self):\n"
15     test_header+=" self.__dom =Rectangle(11,11,2)\n"
16     test_header+=" self.functionspace = FunctionOnBoundary(self.__dom)\n"
17     test_tail=""
18     test_tail+="suite = unittest.TestSuite()\n"
19     test_tail+="suite.addTest(unittest.makeSuite(Test_util2))\n"
20     test_tail+="unittest.TextTestRunner(verbosity=2).run(suite)\n"
21    
22     case_set=["float","array","constData","taggedData","expandedData","Symbol"]
23     shape_set=[ (),(2,), (4,5), (6,2,2),(3,2,3,4)]
24    
25 jgs 154 t_prog=""
26 gross 157 t_prog_with_tags=""
27     t_prog_failing=""
28     u_prog=""
29 jgs 154
30 gross 157 def wherepos(arg):
31     if arg>0.:
32     return 1.
33     else:
34     return 0.
35    
36 gross 313
37 gross 157 class OPERATOR:
38     def __init__(self,nickname,rng=[-1000.,1000],test_expr="",math_expr=None,
39     numarray_expr="",symbol_expr=None,diff=None,name=""):
40     self.nickname=nickname
41     self.rng=rng
42     self.test_expr=test_expr
43     if math_expr==None:
44     self.math_expr=test_expr.replace("%a1%","arg")
45     else:
46     self.math_expr=math_expr
47     self.numarray_expr=numarray_expr
48     self.symbol_expr=symbol_expr
49     self.diff=diff
50     self.name=name
51    
52 jgs 154 import random
53     import numarray
54     import math
55     finc=1.e-6
56    
57 gross 157 def getResultCaseForBin(case0,case1):
58     if case0=="float":
59     if case1=="float":
60     case="float"
61     elif case1=="array":
62     case="array"
63     elif case1=="constData":
64     case="constData"
65     elif case1=="taggedData":
66     case="taggedData"
67     elif case1=="expandedData":
68     case="expandedData"
69     elif case1=="Symbol":
70     case="Symbol"
71     else:
72     raise ValueError,"unknown case1=%s"%case1
73     elif case0=="array":
74     if case1=="float":
75     case="array"
76     elif case1=="array":
77     case="array"
78     elif case1=="constData":
79     case="constData"
80     elif case1=="taggedData":
81     case="taggedData"
82     elif case1=="expandedData":
83     case="expandedData"
84     elif case1=="Symbol":
85     case="Symbol"
86     else:
87     raise ValueError,"unknown case1=%s"%case1
88     elif case0=="constData":
89     if case1=="float":
90     case="constData"
91     elif case1=="array":
92     case="constData"
93     elif case1=="constData":
94     case="constData"
95     elif case1=="taggedData":
96     case="taggedData"
97     elif case1=="expandedData":
98     case="expandedData"
99     elif case1=="Symbol":
100     case="Symbol"
101     else:
102     raise ValueError,"unknown case1=%s"%case1
103     elif case0=="taggedData":
104     if case1=="float":
105     case="taggedData"
106     elif case1=="array":
107     case="taggedData"
108     elif case1=="constData":
109     case="taggedData"
110     elif case1=="taggedData":
111     case="taggedData"
112     elif case1=="expandedData":
113     case="expandedData"
114     elif case1=="Symbol":
115     case="Symbol"
116     else:
117     raise ValueError,"unknown case1=%s"%case1
118     elif case0=="expandedData":
119     if case1=="float":
120     case="expandedData"
121     elif case1=="array":
122     case="expandedData"
123     elif case1=="constData":
124     case="expandedData"
125     elif case1=="taggedData":
126     case="expandedData"
127     elif case1=="expandedData":
128     case="expandedData"
129     elif case1=="Symbol":
130     case="Symbol"
131     else:
132     raise ValueError,"unknown case1=%s"%case1
133     elif case0=="Symbol":
134     if case1=="float":
135     case="Symbol"
136     elif case1=="array":
137     case="Symbol"
138     elif case1=="constData":
139     case="Symbol"
140     elif case1=="taggedData":
141     case="Symbol"
142     elif case1=="expandedData":
143     case="Symbol"
144     elif case1=="Symbol":
145     case="Symbol"
146     else:
147     raise ValueError,"unknown case1=%s"%case1
148     else:
149     raise ValueError,"unknown case0=%s"%case0
150     return case
151 jgs 154
152 gross 157
153 jgs 154 def makeArray(shape,rng):
154     l=rng[1]-rng[0]
155     out=numarray.zeros(shape,numarray.Float64)
156     if len(shape)==0:
157 gross 157 out=l*random.random()+rng[0]
158 jgs 154 elif len(shape)==1:
159     for i0 in range(shape[0]):
160     out[i0]=l*random.random()+rng[0]
161     elif len(shape)==2:
162     for i0 in range(shape[0]):
163     for i1 in range(shape[1]):
164     out[i0,i1]=l*random.random()+rng[0]
165     elif len(shape)==3:
166     for i0 in range(shape[0]):
167     for i1 in range(shape[1]):
168     for i2 in range(shape[2]):
169     out[i0,i1,i2]=l*random.random()+rng[0]
170     elif len(shape)==4:
171     for i0 in range(shape[0]):
172     for i1 in range(shape[1]):
173     for i2 in range(shape[2]):
174     for i3 in range(shape[3]):
175     out[i0,i1,i2,i3]=l*random.random()+rng[0]
176 gross 437 elif len(shape)==5:
177     for i0 in range(shape[0]):
178     for i1 in range(shape[1]):
179     for i2 in range(shape[2]):
180     for i3 in range(shape[3]):
181     for i4 in range(shape[4]):
182     out[i0,i1,i2,i3,i4]=l*random.random()+rng[0]
183 jgs 154 else:
184 gross 437 raise SystemError,"rank is restricted to 5"
185 jgs 154 return out
186    
187    
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    
1024 gross 441 #=======================================================================================================
1025 gross 443 # L2
1026     #=======================================================================================================
1027     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1028     for data in ["Data","Symbol"]:
1029     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1030     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1031     tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
1032     text+=" def %s(self):\n"%tname
1033     text+=" \"\"\"\n"
1034     text+=" tests L2-norm of %s on the %s\n\n"%(data,where)
1035     text+=" assumptions: self.domain supports integration on %s\n"%where
1036     text+=" \"\"\"\n"
1037     text+=" dim=self.domain.getDim()\n"
1038     text+=" w=%s(self.domain)\n"%where
1039     text+=" x=w.getX()\n"
1040     o="1"
1041     if len(sh)>0:
1042     sh_2=sh[:len(sh)-1]+(2,)
1043     sh_3=sh[:len(sh)-1]+(3,)
1044     b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
1045     b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
1046     else:
1047     sh_2=()
1048     sh_3=()
1049     b_2=makeArray(sh,[-1.,1])
1050     b_3=makeArray(sh,[-1.,1])
1051    
1052     if data=="Symbol":
1053     val="s"
1054     res="sub"
1055     else:
1056     val="arg"
1057     res="res"
1058     text+=" if dim==2:\n"
1059     if data=="Symbol":
1060     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
1061    
1062     text+=" %s=Data(0,%s,w)\n"%(val,sh_2)
1063     text+=unrollLoopsSimplified(b_2,val,tap=" ")
1064     text+=unrollLoopsOfL2(b_2,where,"ref",tap=" ")
1065     text+="\n else:\n"
1066     if data=="Symbol":
1067     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
1068     text+=" %s=Data(0,%s,w)\n"%(val,sh_3)
1069     text+=unrollLoopsSimplified(b_3,val,tap=" ")
1070     text+=unrollLoopsOfL2(b_3,where,"ref",tap=" ")
1071     text+="\n res=L2(arg)\n"
1072     if data=="Symbol":
1073     text+=" sub=res.substitute({arg:s})\n"
1074     text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1075     text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1076     else:
1077     text+=" self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1078     text+=" self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
1079     t_prog+=text
1080     print t_prog
1081     1/0
1082    
1083     #=======================================================================================================
1084 gross 441 # div
1085     #=======================================================================================================
1086     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1087     for data in ["Data","Symbol"]:
1088     for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1089     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1090     tname="test_div_on%s_from%s_%s"%(where,data,case)
1091     text+=" def %s(self):\n"%tname
1092     text+=" \"\"\"\n"
1093     text+=" tests divergence of %s on the %s\n\n"%(data,where)
1094     text+=" assumptions: %s(self.domain) exists\n"%case
1095     text+=" self.domain supports div on %s\n"%where
1096     text+=" \"\"\"\n"
1097     if case=="ReducedSolution":
1098     text+=" o=1\n"
1099     o="1"
1100     else:
1101     text+=" o=self.order\n"
1102     o="o"
1103     text+=" dim=self.domain.getDim()\n"
1104     text+=" w_ref=%s(self.domain)\n"%where
1105     text+=" x_ref=w_ref.getX()\n"
1106     text+=" w=%s(self.domain)\n"%case
1107     text+=" x=w.getX()\n"
1108     a_2=makeArray((2,2),[-1.,1])
1109     b_2=makeArray((2,2),[-1.,1])
1110     a_3=makeArray((3,3),[-1.,1])
1111     b_3=makeArray((3,3),[-1.,1])
1112     if data=="Symbol":
1113     text+=" arg=Symbol(shape=(dim,),dim=dim)\n"
1114     val="s"
1115     res="sub"
1116     else:
1117     val="arg"
1118     res="res"
1119     text+=" %s=Vector(0,w)\n"%val
1120     text+=" if dim==2:\n"
1121     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1122     text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap=" ")
1123     text+="\n else:\n"
1124    
1125     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1126     text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap=" ")
1127     text+="\n res=div(arg,where=w_ref)\n"
1128     if data=="Symbol":
1129     text+=" sub=res.substitute({arg:s})\n"
1130     text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1131     text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1132     text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
1133     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1134 gross 437
1135 gross 441
1136     t_prog+=text
1137     print t_prog
1138     1/0
1139    
1140 gross 439 #=======================================================================================================
1141     # interpolation
1142     #=======================================================================================================
1143     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
1144     for data in ["Data","Symbol"]:
1145     for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1146     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1147     if where==case or \
1148     ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
1149     ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
1150     (case=="ContinuousFunction" and where in ["Solution","ReducedSolution"]) or \
1151     (case=="Solution" and where=="ReducedSolution") :
1152    
1153    
1154     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1155     tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1156     text+=" def %s(self):\n"%tname
1157     text+=" \"\"\"\n"
1158     text+=" tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
1159     text+=" assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
1160     text+=" \"\"\"\n"
1161     if case=="ReducedSolution" or where=="ReducedSolution":
1162     text+=" o=1\n"
1163     o="1"
1164     else:
1165     text+=" o=self.order\n"
1166     o="o"
1167     text+=" dim=self.domain.getDim()\n"
1168     text+=" w_ref=%s(self.domain)\n"%where
1169     text+=" x_ref=w_ref.getX()\n"
1170     text+=" w=%s(self.domain)\n"%case
1171     text+=" x=w.getX()\n"
1172     a_2=makeArray(sh+(2,),[-1.,1])
1173     b_2=makeArray(sh+(2,),[-1.,1])
1174     a_3=makeArray(sh+(3,),[-1.,1])
1175     b_3=makeArray(sh+(3,),[-1.,1])
1176     if data=="Symbol":
1177     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1178     val="s"
1179     res="sub"
1180     else:
1181     val="arg"
1182     res="res"
1183     text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1184     text+=" ref=Data(0,%s,w_ref)\n"%str(sh)
1185     text+=" if dim==2:\n"
1186     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1187     text+=unrollLoops(a_2,b_2,o,"ref",tap=" ",x="x_ref")
1188     text+=" else:\n"
1189    
1190     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1191     text+=unrollLoops(a_3,b_3,o,"ref",tap=" ",x="x_ref")
1192     text+=" res=interpolate(arg,where=w_ref)\n"
1193     if data=="Symbol":
1194     text+=" sub=res.substitute({arg:s})\n"
1195     text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1196     text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
1197     text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1198     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1199     t_prog+=text
1200     print test_header
1201     print t_prog
1202     print test_tail
1203     1/0
1204 gross 437
1205 gross 429 #=======================================================================================================
1206 gross 438 # grad
1207 gross 437 #=======================================================================================================
1208     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1209     for data in ["Data","Symbol"]:
1210 gross 438 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1211     for sh in [ (),(2,), (4,5), (6,2,2)]:
1212 gross 437 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1213 gross 438 tname="test_grad_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1214 gross 437 text+=" def %s(self):\n"%tname
1215     text+=" \"\"\"\n"
1216 gross 438 text+=" tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
1217 gross 437 text+=" assumptions: %s(self.domain) exists\n"%case
1218 gross 438 text+=" self.domain supports gardient on %s\n"%where
1219 gross 437 text+=" \"\"\"\n"
1220     if case=="ReducedSolution":
1221     text+=" o=1\n"
1222     o="1"
1223     else:
1224     text+=" o=self.order\n"
1225     o="o"
1226     text+=" dim=self.domain.getDim()\n"
1227     text+=" w_ref=%s(self.domain)\n"%where
1228 gross 438 text+=" x_ref=w_ref.getX()\n"
1229 gross 437 text+=" w=%s(self.domain)\n"%case
1230     text+=" x=w.getX()\n"
1231     a_2=makeArray(sh+(2,),[-1.,1])
1232     b_2=makeArray(sh+(2,),[-1.,1])
1233     a_3=makeArray(sh+(3,),[-1.,1])
1234     b_3=makeArray(sh+(3,),[-1.,1])
1235     if data=="Symbol":
1236 gross 438 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1237 gross 437 val="s"
1238     res="sub"
1239     else:
1240     val="arg"
1241     res="res"
1242     text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1243 gross 438 text+=" ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1244 gross 437 text+=" if dim==2:\n"
1245     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1246 gross 438 text+=unrollLoopsOfGrad(a_2,b_2,o,"ref",tap=" ")
1247 gross 437 text+=" else:\n"
1248    
1249     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1250 gross 438 text+=unrollLoopsOfGrad(a_3,b_3,o,"ref",tap=" ")
1251     text+=" res=grad(arg,where=w_ref)\n"
1252 gross 437 if data=="Symbol":
1253     text+=" sub=res.substitute({arg:s})\n"
1254 gross 438 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1255     text+=" self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1256 gross 437 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1257    
1258    
1259     t_prog+=text
1260     print test_header
1261     print t_prog
1262     print test_tail
1263     1/0
1264 gross 438
1265    
1266 gross 437 #=======================================================================================================
1267 gross 438 # integrate
1268 gross 437 #=======================================================================================================
1269 gross 438 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1270 gross 437 for data in ["Data","Symbol"]:
1271 gross 438 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1272     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1273     if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:
1274 gross 437 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1275 gross 438 tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1276 gross 437 text+=" def %s(self):\n"%tname
1277     text+=" \"\"\"\n"
1278 gross 438 text+=" tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
1279     text+=" assumptions: %s(self.domain) exists\n"%case
1280     text+=" self.domain supports integral on %s\n"%where
1281    
1282 gross 437 text+=" \"\"\"\n"
1283     if case=="ReducedSolution":
1284     text+=" o=1\n"
1285     o="1"
1286     else:
1287     text+=" o=self.order\n"
1288     o="o"
1289     text+=" dim=self.domain.getDim()\n"
1290     text+=" w_ref=%s(self.domain)\n"%where
1291     text+=" w=%s(self.domain)\n"%case
1292     text+=" x=w.getX()\n"
1293     a_2=makeArray(sh+(2,),[-1.,1])
1294     b_2=makeArray(sh+(2,),[-1.,1])
1295     a_3=makeArray(sh+(3,),[-1.,1])
1296     b_3=makeArray(sh+(3,),[-1.,1])
1297     if data=="Symbol":
1298 gross 438 text+=" arg=Symbol(shape=%s)\n"%str(sh)
1299 gross 437 val="s"
1300     res="sub"
1301     else:
1302     val="arg"
1303     res="res"
1304 gross 438
1305 gross 437 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1306 gross 438 if not len(sh)==0:
1307     text+=" ref=numarray.zeros(%s,numarray.Float)\n"%str(sh)
1308 gross 437 text+=" if dim==2:\n"
1309     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1310 gross 438 text+=unrollLoopsOfInteriorIntegral(a_2,b_2,where,"ref",tap=" ")
1311 gross 437 text+=" else:\n"
1312    
1313     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1314 gross 438 text+=unrollLoopsOfInteriorIntegral(a_3,b_3,where,"ref",tap=" ")
1315     if case in ["ContinuousFunction","Solution","ReducedSolution"]:
1316     text+=" res=integrate(arg,where=w_ref)\n"
1317     else:
1318     text+=" res=integrate(arg)\n"
1319    
1320 gross 437 if data=="Symbol":
1321     text+=" sub=res.substitute({arg:s})\n"
1322 gross 438 if len(sh)==0 and data=="Data":
1323     text+=" self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%res
1324     else:
1325     if data=="Symbol":
1326     text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1327     text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1328     else:
1329     text+=" self.failUnless(isinstance(res,numarray.NumArray),\"wrong type of result.\")\n"
1330     text+=" self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
1331 gross 437 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1332    
1333    
1334     t_prog+=text
1335     print test_header
1336     print t_prog
1337     print test_tail
1338     1/0
1339     #=======================================================================================================
1340 gross 433 # inverse
1341     #=======================================================================================================
1342     name="inverse"
1343     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1344     for sh0 in [ (1,1), (2,2), (3,3)]:
1345     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1346     tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1347     text+=" def %s(self):\n"%tname
1348     a_0=makeArray(sh0,[-1.,1])
1349     for i in range(sh0[0]): a_0[i,i]+=2
1350     if case0 in ["taggedData", "expandedData"]:
1351     a1_0=makeArray(sh0,[-1.,1])
1352     for i in range(sh0[0]): a1_0[i,i]+=3
1353     else:
1354     a1_0=a_0
1355    
1356     text+=mkText(case0,"arg",a_0,a1_0)
1357     text+=" res=%s(arg)\n"%name
1358     if case0=="Symbol":
1359     text+=mkText("array","s",a_0,a1_0)
1360     text+=" sub=res.substitute({arg:s})\n"
1361     res="sub"
1362     ref="s"
1363     else:
1364     ref="arg"
1365     res="res"
1366     text+=mkTypeAndShapeTest(case0,sh0,"res")
1367     text+=" self.failUnless(Lsup(matrixmult(%s,%s)-kronecker(%s))<=self.RES_TOL,\"wrong result\")\n"%(res,ref,sh0[0])
1368    
1369     if case0 == "taggedData" :
1370     t_prog_with_tags+=text
1371     else:
1372     t_prog+=text
1373    
1374     print test_header
1375     # print t_prog
1376     print t_prog_with_tags
1377     print test_tail
1378     1/0
1379    
1380     #=======================================================================================================
1381 gross 429 # trace
1382     #=======================================================================================================
1383     def traceTest(r,offset):
1384     sh=r.shape
1385     r1=1
1386     for i in range(offset): r1*=sh[i]
1387     r2=1
1388     for i in range(offset+2,len(sh)): r2*=sh[i]
1389     r_s=numarray.reshape(r,(r1,sh[offset],sh[offset],r2))
1390     s=numarray.zeros([r1,r2],numarray.Float)
1391     for i1 in range(r1):
1392     for i2 in range(r2):
1393     for j in range(sh[offset]): s[i1,i2]+=r_s[i1,j,j,i2]
1394     return s.resize(sh[:offset]+sh[offset+2:])
1395     name,tt="trace",traceTest
1396     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1397     for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
1398     for offset in range(len(sh0)-1):
1399     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1400     tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1401     text+=" def %s(self):\n"%tname
1402     sh_t=list(sh0)
1403     sh_t[offset+1]=sh_t[offset]
1404     sh_t=tuple(sh_t)
1405     sh_r=[]
1406     for i in range(offset): sh_r.append(sh0[i])
1407     for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1408     sh_r=tuple(sh_r)
1409     a_0=makeArray(sh_t,[-1.,1])
1410     if case0 in ["taggedData", "expandedData"]:
1411     a1_0=makeArray(sh_t,[-1.,1])
1412     else:
1413     a1_0=a_0
1414     r=tt(a_0,offset)
1415     r1=tt(a1_0,offset)
1416     text+=mkText(case0,"arg",a_0,a1_0)
1417     text+=" res=%s(arg,%s)\n"%(name,offset)
1418     if case0=="Symbol":
1419     text+=mkText("array","s",a_0,a1_0)
1420     text+=" sub=res.substitute({arg:s})\n"
1421     res="sub"
1422     text+=mkText("array","ref",r,r1)
1423     else:
1424     res="res"
1425     text+=mkText(case0,"ref",r,r1)
1426     text+=mkTypeAndShapeTest(case0,sh_r,"res")
1427     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1428    
1429     if case0 == "taggedData" :
1430     t_prog_with_tags+=text
1431     else:
1432     t_prog+=text
1433 gross 396
1434 gross 429 print test_header
1435     # print t_prog
1436     print t_prog_with_tags
1437     print test_tail
1438     1/0
1439 gross 396
1440 gross 157 #=======================================================================================================
1441 gross 396 # clip
1442     #=======================================================================================================
1443     oper_L=[["clip",clipTEST]]
1444     for oper in oper_L:
1445     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1446     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1447     if len(sh0)==0 or not case0=="float":
1448     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1449     tname="test_%s_%s_rank%s"%(oper[0],case0,len(sh0))
1450     text+=" def %s(self):\n"%tname
1451     a_0=makeArray(sh0,[-1.,1])
1452     if case0 in ["taggedData", "expandedData"]:
1453     a1_0=makeArray(sh0,[-1.,1])
1454     else:
1455     a1_0=a_0
1456    
1457     r=oper[1](a_0,-0.3,0.5)
1458     r1=oper[1](a1_0,-0.3,0.5)
1459     text+=mkText(case0,"arg",a_0,a1_0)
1460     text+=" res=%s(arg,-0.3,0.5)\n"%oper[0]
1461     if case0=="Symbol":
1462     text+=mkText("array","s",a_0,a1_0)
1463     text+=" sub=res.substitute({arg:s})\n"
1464     res="sub"
1465     text+=mkText("array","ref",r,r1)
1466     else:
1467     res="res"
1468     text+=mkText(case0,"ref",r,r1)
1469     text+=mkTypeAndShapeTest(case0,sh0,"res")
1470     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1471    
1472     if case0 == "taggedData" :
1473     t_prog_with_tags+=text
1474     else:
1475     t_prog+=text
1476    
1477     print test_header
1478     # print t_prog
1479     print t_prog_with_tags
1480     print test_tail
1481     1/0
1482 gross 429
1483 gross 396 #=======================================================================================================
1484     # maximum, minimum, clipping
1485     #=======================================================================================================
1486     oper_L=[ ["maximum",maximumTEST],
1487     ["minimum",minimumTEST]]
1488     for oper in oper_L:
1489     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1490     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1491     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1492     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1493     if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1494     and (sh0==sh1 or len(sh0)==0 or len(sh1)==0) :
1495     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1496    
1497     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1498     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1499     text+=" def %s(self):\n"%tname
1500     a_0=makeArray(sh0,[-1.,1])
1501     if case0 in ["taggedData", "expandedData"]:
1502     a1_0=makeArray(sh0,[-1.,1])
1503     else:
1504     a1_0=a_0
1505    
1506     a_1=makeArray(sh1,[-1.,1])
1507     if case1 in ["taggedData", "expandedData"]:
1508     a1_1=makeArray(sh1,[-1.,1])
1509     else:
1510     a1_1=a_1
1511     r=oper[1](a_0,a_1)
1512     r1=oper[1](a1_0,a1_1)
1513     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1514     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1515     text+=" res=%s(arg0,arg1)\n"%oper[0]
1516     case=getResultCaseForBin(case0,case1)
1517     if case=="Symbol":
1518     c0_res,c1_res=case0,case1
1519     subs="{"
1520     if case0=="Symbol":
1521     text+=mkText("array","s0",a_0,a1_0)
1522     subs+="arg0:s0"
1523     c0_res="array"
1524     if case1=="Symbol":
1525     text+=mkText("array","s1",a_1,a1_1)
1526     if not subs.endswith("{"): subs+=","
1527     subs+="arg1:s1"
1528     c1_res="array"
1529     subs+="}"
1530     text+=" sub=res.substitute(%s)\n"%subs
1531     res="sub"
1532     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1533     else:
1534     res="res"
1535     text+=mkText(case,"ref",r,r1)
1536     if len(sh0)>len(sh1):
1537     text+=mkTypeAndShapeTest(case,sh0,"res")
1538     else:
1539     text+=mkTypeAndShapeTest(case,sh1,"res")
1540     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1541    
1542     if case0 == "taggedData" or case1 == "taggedData":
1543     t_prog_with_tags+=text
1544     else:
1545     t_prog+=text
1546    
1547     print test_header
1548     # print t_prog
1549     print t_prog_with_tags
1550     print test_tail
1551     1/0
1552    
1553    
1554     #=======================================================================================================
1555     # outer inner
1556     #=======================================================================================================
1557     oper=["outer",outerTEST]
1558     # oper=["inner",innerTEST]
1559     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1560     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1561     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1562     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1563     if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1564     and len(sh0+sh1)<5:
1565     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1566    
1567     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1568     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1569     text+=" def %s(self):\n"%tname
1570     a_0=makeArray(sh0,[-1.,1])
1571     if case0 in ["taggedData", "expandedData"]:
1572     a1_0=makeArray(sh0,[-1.,1])
1573     else:
1574     a1_0=a_0
1575    
1576     a_1=makeArray(sh1,[-1.,1])
1577     if case1 in ["taggedData", "expandedData"]:
1578     a1_1=makeArray(sh1,[-1.,1])
1579     else:
1580     a1_1=a_1
1581     r=oper[1](a_0,a_1)
1582     r1=oper[1](a1_0,a1_1)
1583     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1584     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1585     text+=" res=%s(arg0,arg1)\n"%oper[0]
1586     case=getResultCaseForBin(case0,case1)
1587     if case=="Symbol":
1588     c0_res,c1_res=case0,case1
1589     subs="{"
1590     if case0=="Symbol":
1591     text+=mkText("array","s0",a_0,a1_0)
1592     subs+="arg0:s0"
1593     c0_res="array"
1594     if case1=="Symbol":
1595     text+=mkText("array","s1",a_1,a1_1)
1596     if not subs.endswith("{"): subs+=","
1597     subs+="arg1:s1"
1598     c1_res="array"
1599     subs+="}"
1600     text+=" sub=res.substitute(%s)\n"%subs
1601     res="sub"
1602     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1603     else:
1604     res="res"
1605     text+=mkText(case,"ref",r,r1)
1606     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1607     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1608    
1609     if case0 == "taggedData" or case1 == "taggedData":
1610     t_prog_with_tags+=text
1611     else:
1612     t_prog+=text
1613    
1614     print test_header
1615     # print t_prog
1616     print t_prog_with_tags
1617     print test_tail
1618     1/0
1619    
1620     #=======================================================================================================
1621 gross 313 # local reduction
1622     #=======================================================================================================
1623     for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],
1624     ["maxval",-1.e99,"max(out,%a1%)","out"],
1625     ["minval",1.e99,"min(out,%a1%)","out"] ]:
1626     for case in case_set:
1627     for sh in shape_set:
1628     if not case=="float" or len(sh)==0:
1629     text=""
1630     text+=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1631     tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))
1632     text+=" %s(self):\n"%tname
1633     a=makeArray(sh,[-1.,1.])
1634     a1=makeArray(sh,[-1.,1.])
1635     r1=testReduce(a1,oper[1],oper[2],oper[3])
1636     r=testReduce(a,oper[1],oper[2],oper[3])
1637    
1638     text+=mkText(case,"arg",a,a1)
1639     text+=" res=%s(arg)\n"%oper[0]
1640     if case=="Symbol":
1641     text+=mkText("array","s",a,a1)
1642     text+=" sub=res.substitute({arg:s})\n"
1643     text+=mkText("array","ref",r,r1)
1644     res="sub"
1645     else:
1646     text+=mkText(case,"ref",r,r1)
1647     res="res"
1648     if oper[0]=="length":
1649     text+=mkTypeAndShapeTest(case,(),"res")
1650     else:
1651     if case=="float" or case=="array":
1652     text+=mkTypeAndShapeTest("float",(),"res")
1653     else:
1654     text+=mkTypeAndShapeTest(case,(),"res")
1655     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1656     if case == "taggedData":
1657     t_prog_with_tags+=text
1658     else:
1659     t_prog+=text
1660     print test_header
1661     # print t_prog
1662     print t_prog_with_tags
1663     print test_tail
1664     1/0
1665    
1666     #=======================================================================================================
1667 gross 291 # tensor multiply
1668 gross 157 #=======================================================================================================
1669 gross 291 # oper=["generalTensorProduct",tensorProductTest]
1670     # oper=["matrixmult",testMatrixMult]
1671     oper=["tensormult",testTensorMult]
1672    
1673     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1674     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1675     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1676     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1677     for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
1678     if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
1679     and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
1680     # 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
1681     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
1682     case=getResultCaseForBin(case0,case1)
1683 gross 157 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1684     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1685 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))
1686     #tname="test_matrixmult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1687     tname="test_tensormult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1688     # if tname=="test_generalTensorProduct_array_rank1_array_rank2_offset1":
1689     # print tnametest_generalTensorProduct_Symbol_rank1_Symbol_rank3_offset1
1690     text+=" def %s(self):\n"%tname
1691     a_0=makeArray(sh0+sh_s,[-1.,1])
1692     if case0 in ["taggedData", "expandedData"]:
1693     a1_0=makeArray(sh0+sh_s,[-1.,1])
1694     else:
1695     a1_0=a_0
1696    
1697     a_1=makeArray(sh_s+sh1,[-1.,1])
1698     if case1 in ["taggedData", "expandedData"]:
1699     a1_1=makeArray(sh_s+sh1,[-1.,1])
1700     else:
1701     a1_1=a_1
1702     r=oper[1](a_0,a_1,sh_s)
1703     r1=oper[1](a1_0,a1_1,sh_s)
1704     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1705     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1706     #text+=" res=matrixmult(arg0,arg1)\n"
1707     text+=" res=tensormult(arg0,arg1)\n"
1708     #text+=" res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))
1709     if case=="Symbol":
1710     c0_res,c1_res=case0,case1
1711     subs="{"
1712     if case0=="Symbol":
1713     text+=mkText("array","s0",a_0,a1_0)
1714     subs+="arg0:s0"
1715     c0_res="array"
1716     if case1=="Symbol":
1717     text+=mkText("array","s1",a_1,a1_1)
1718     if not subs.endswith("{"): subs+=","
1719     subs+="arg1:s1"
1720     c1_res="array"
1721     subs+="}"
1722     text+=" sub=res.substitute(%s)\n"%subs
1723     res="sub"
1724     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1725     else:
1726     res="res"
1727     text+=mkText(case,"ref",r,r1)
1728     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1729     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1730     if case0 == "taggedData" or case1 == "taggedData":
1731     t_prog_with_tags+=text
1732     else:
1733     t_prog+=text
1734     print test_header
1735     # print t_prog
1736     print t_prog_with_tags
1737     print test_tail
1738     1/0
1739     #=======================================================================================================
1740 gross 157 # basic binary operation overloading (tests only!)
1741     #=======================================================================================================
1742     oper_range=[-5.,5.]
1743     for oper in [["add" ,"+",[-5.,5.]],
1744     ["sub" ,"-",[-5.,5.]],
1745     ["mult","*",[-5.,5.]],
1746     ["div" ,"/",[-5.,5.]],
1747     ["pow" ,"**",[0.01,5.]]]:
1748     for case0 in case_set:
1749     for sh0 in shape_set:
1750     for case1 in case_set:
1751     for sh1 in shape_set:
1752 gross 291 if not case0=="array" and \
1753     (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1754 gross 157 (sh0==() or sh1==() or sh1==sh0) and \
1755     not (case0 in ["float","array"] and case1 in ["float","array"]):
1756 gross 291 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1757 gross 157 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1758     tname="test_%s_overloaded_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1759     text+=" def %s(self):\n"%tname
1760     a_0=makeArray(sh0,oper[2])
1761     if case0 in ["taggedData", "expandedData"]:
1762     a1_0=makeArray(sh0,oper[2])
1763     else:
1764     a1_0=a_0
1765 jgs 154
1766 gross 157 a_1=makeArray(sh1,oper[2])
1767     if case1 in ["taggedData", "expandedData"]:
1768     a1_1=makeArray(sh1,oper[2])
1769 jgs 154 else:
1770 gross 157 a1_1=a_1
1771     r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1772     r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1773     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1774     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1775     text+=" res=arg0%sarg1\n"%oper[1]
1776    
1777     case=getResultCaseForBin(case0,case1)
1778     if case=="Symbol":
1779     c0_res,c1_res=case0,case1
1780     subs="{"
1781     if case0=="Symbol":
1782     text+=mkText("array","s0",a_0,a1_0)
1783     subs+="arg0:s0"
1784     c0_res="array"
1785     if case1=="Symbol":
1786     text+=mkText("array","s1",a_1,a1_1)
1787     if not subs.endswith("{"): subs+=","
1788     subs+="arg1:s1"
1789     c1_res="array"
1790     subs+="}"
1791     text+=" sub=res.substitute(%s)\n"%subs
1792     res="sub"
1793     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1794 jgs 154 else:
1795 gross 157 res="res"
1796     text+=mkText(case,"ref",r,r1)
1797     if isinstance(r,float):
1798     text+=mkTypeAndShapeTest(case,(),"res")
1799     else:
1800     text+=mkTypeAndShapeTest(case,r.shape,"res")
1801     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1802    
1803     if case0 in [ "constData","taggedData","expandedData"] and case1 == "Symbol":
1804     t_prog_failing+=text
1805     else:
1806     if case0 == "taggedData" or case1 == "taggedData":
1807     t_prog_with_tags+=text
1808     else:
1809     t_prog+=text
1810 jgs 154
1811 gross 157
1812     print test_header
1813 gross 291 # print t_prog
1814     # print t_prog_with_tags
1815     print t_prog_failing
1816     print test_tail
1817 jgs 154 1/0
1818 gross 291 #=======================================================================================================
1819     # basic binary operations (tests only!)
1820     #=======================================================================================================
1821     oper_range=[-5.,5.]
1822     for oper in [["add" ,"+",[-5.,5.]],
1823     ["mult","*",[-5.,5.]],
1824     ["quotient" ,"/",[-5.,5.]],
1825     ["power" ,"**",[0.01,5.]]]:
1826     for case0 in case_set:
1827     for case1 in case_set:
1828     for sh in shape_set:
1829     for sh_p in shape_set:
1830     if len(sh_p)>0:
1831     resource=[-1,1]
1832     else:
1833     resource=[1]
1834     for sh_d in resource:
1835     if sh_d>0:
1836     sh0=sh
1837     sh1=sh+sh_p
1838     else:
1839     sh1=sh
1840     sh0=sh+sh_p
1841    
1842     if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1843     len(sh0)<5 and len(sh1)<5:
1844     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1845     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1846     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1847     text+=" def %s(self):\n"%tname
1848     a_0=makeArray(sh0,oper[2])
1849     if case0 in ["taggedData", "expandedData"]:
1850     a1_0=makeArray(sh0,oper[2])
1851     else:
1852     a1_0=a_0
1853    
1854     a_1=makeArray