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

Annotation of /trunk/escript/py_src/generatediff

Parent Directory Parent Directory | Revision Log Revision Log


Revision 441 - (hide annotations)
Fri Jan 20 03:40:39 2006 UTC (13 years, 9 months ago) by gross
Original Path: trunk/escript/py_src/generateutil
File size: 123349 byte(s)
and finally test for the divergence operator added. div() has been modified to take Symbol arguments
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 439 def unrollLoops(a,b,o,arg,tap="",x="x"):
700 gross 437 out=""
701     if a.rank==1:
702     z=""
703     for i99 in range(a.shape[0]):
704     if not z=="": z+="+"
705     if o=="1":
706 gross 439 z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
707 gross 437 else:
708 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
709 gross 437
710     out+=tap+"%s=%s\n"%(arg,z)
711    
712     elif a.rank==2:
713     for i0 in range(a.shape[0]):
714     z=""
715     for i99 in range(a.shape[1]):
716     if not z=="": z+="+"
717     if o=="1":
718     z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
719     else:
720 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
721 gross 437
722     out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
723     elif a.rank==3:
724     for i0 in range(a.shape[0]):
725     for i1 in range(a.shape[1]):
726     z=""
727     for i99 in range(a.shape[2]):
728     if not z=="": z+="+"
729     if o=="1":
730 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
731 gross 437 else:
732 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
733 gross 437
734     out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
735     elif a.rank==4:
736     for i0 in range(a.shape[0]):
737     for i1 in range(a.shape[1]):
738     for i2 in range(a.shape[2]):
739     z=""
740     for i99 in range(a.shape[3]):
741     if not z=="": z+="+"
742     if o=="1":
743 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
744 gross 437 else:
745 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
746 gross 437
747     out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
748     elif a.rank==5:
749     for i0 in range(a.shape[0]):
750     for i1 in range(a.shape[1]):
751     for i2 in range(a.shape[2]):
752     for i3 in range(a.shape[3]):
753     z=""
754     for i99 in range(a.shape[4]):
755     if not z=="": z+="+"
756     if o=="1":
757 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
758 gross 437 else:
759 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)
760 gross 437
761     out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
762     return out
763    
764     def unrollLoopsOfGrad(a,b,o,arg,tap=""):
765     out=""
766     if a.rank==1:
767     for i99 in range(a.shape[0]):
768     if o=="1":
769     out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
770     else:
771     out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
772    
773     elif a.rank==2:
774     for i0 in range(a.shape[0]):
775     for i99 in range(a.shape[1]):
776     if o=="1":
777     out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
778     else:
779     out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
780    
781     elif a.rank==3:
782     for i0 in range(a.shape[0]):
783     for i1 in range(a.shape[1]):
784     for i99 in range(a.shape[2]):
785     if o=="1":
786     out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
787     else:
788     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])
789    
790     elif a.rank==4:
791     for i0 in range(a.shape[0]):
792     for i1 in range(a.shape[1]):
793     for i2 in range(a.shape[2]):
794     for i99 in range(a.shape[3]):
795     if o=="1":
796     out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
797     else:
798     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])
799 gross 438 return out
800 gross 441 def unrollLoopsOfDiv(a,b,o,arg,tap=""):
801    
802    
803    
804     out=tap+arg+"="
805     if o=="1":
806     z=0.
807     for i99 in range(a.shape[0]):
808     z+=b[i99,i99]+a[i99,i99]
809     out+="(%s)"%z
810     else:
811     z=0.
812     for i99 in range(a.shape[0]):
813     z+=b[i99,i99]
814     if i99>0: out+="+"
815     out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
816     out+="+(%s)"%z
817     return out
818    
819 gross 437 def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
820     if where=="Function":
821     xfac_o=1.
822     xfac_op=0.
823     z_fac=1./2.
824     z_fac_s=""
825     zo_fac_s="/(o+1.)"
826     elif where=="FunctionOnBoundary":
827     xfac_o=1.
828     xfac_op=0.
829     z_fac=1.
830     z_fac_s="*dim"
831     zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
832     elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
833     xfac_o=0.
834     xfac_op=1.
835     z_fac=1./2.
836     z_fac_s=""
837     zo_fac_s="/(o+1.)"
838     out=""
839     if a.rank==1:
840     zo=0.
841     zop=0.
842     z=0.
843     for i99 in range(a.shape[0]):
844     if i99==0:
845     zo+= xfac_o*a[i99]
846     zop+= xfac_op*a[i99]
847     else:
848     zo+=a[i99]
849     z+=b[i99]
850    
851     out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
852     if zop==0.:
853     out+="\n"
854     else:
855     out+="+(%s)*0.5**o\n"%zop
856     elif a.rank==2:
857     for i0 in range(a.shape[0]):
858     zo=0.
859     zop=0.
860     z=0.
861     for i99 in range(a.shape[1]):
862     if i99==0:
863     zo+= xfac_o*a[i0,i99]
864     zop+= xfac_op*a[i0,i99]
865     else:
866     zo+=a[i0,i99]
867     z+=b[i0,i99]
868    
869     out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
870     if zop==0.:
871     out+="\n"
872     else:
873     out+="+(%s)*0.5**o\n"%zop
874     elif a.rank==3:
875     for i0 in range(a.shape[0]):
876     for i1 in range(a.shape[1]):
877     zo=0.
878     zop=0.
879     z=0.
880     for i99 in range(a.shape[2]):
881     if i99==0:
882     zo+= xfac_o*a[i0,i1,i99]
883     zop+= xfac_op*a[i0,i1,i99]
884     else:
885     zo+=a[i0,i1,i99]
886     z+=b[i0,i1,i99]
887    
888     out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
889     if zop==0.:
890     out+="\n"
891     else:
892     out+="+(%s)*0.5**o\n"%zop
893     elif a.rank==4:
894     for i0 in range(a.shape[0]):
895     for i1 in range(a.shape[1]):
896     for i2 in range(a.shape[2]):
897     zo=0.
898     zop=0.
899     z=0.
900     for i99 in range(a.shape[3]):
901     if i99==0:
902     zo+= xfac_o*a[i0,i1,i2,i99]
903     zop+= xfac_op*a[i0,i1,i2,i99]
904    
905     else:
906     zo+=a[i0,i1,i2,i99]
907     z+=b[i0,i1,i2,i99]
908    
909     out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
910     if zop==0.:
911     out+="\n"
912     else:
913     out+="+(%s)*0.5**o\n"%zop
914    
915     elif a.rank==5:
916     for i0 in range(a.shape[0]):
917     for i1 in range(a.shape[1]):
918     for i2 in range(a.shape[2]):
919     for i3 in range(a.shape[3]):
920     zo=0.
921     zop=0.
922     z=0.
923     for i99 in range(a.shape[4]):
924     if i99==0:
925     zo+= xfac_o*a[i0,i1,i2,i3,i99]
926     zop+= xfac_op*a[i0,i1,i2,i3,i99]
927    
928     else:
929     zo+=a[i0,i1,i2,i3,i99]
930     z+=b[i0,i1,i2,i3,i99]
931     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)
932     if zop==0.:
933     out+="\n"
934     else:
935     out+="+(%s)*0.5**o\n"%zop
936    
937     return out
938 gross 441 #=======================================================================================================
939     # div
940     #=======================================================================================================
941     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
942     for data in ["Data","Symbol"]:
943     for case in ["ContinuousFunction","Solution","ReducedSolution"]:
944     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
945     tname="test_div_on%s_from%s_%s"%(where,data,case)
946     text+=" def %s(self):\n"%tname
947     text+=" \"\"\"\n"
948     text+=" tests divergence of %s on the %s\n\n"%(data,where)
949     text+=" assumptions: %s(self.domain) exists\n"%case
950     text+=" self.domain supports div on %s\n"%where
951     text+=" \"\"\"\n"
952     if case=="ReducedSolution":
953     text+=" o=1\n"
954     o="1"
955     else:
956     text+=" o=self.order\n"
957     o="o"
958     text+=" dim=self.domain.getDim()\n"
959     text+=" w_ref=%s(self.domain)\n"%where
960     text+=" x_ref=w_ref.getX()\n"
961     text+=" w=%s(self.domain)\n"%case
962     text+=" x=w.getX()\n"
963     a_2=makeArray((2,2),[-1.,1])
964     b_2=makeArray((2,2),[-1.,1])
965     a_3=makeArray((3,3),[-1.,1])
966     b_3=makeArray((3,3),[-1.,1])
967     if data=="Symbol":
968     text+=" arg=Symbol(shape=(dim,),dim=dim)\n"
969     val="s"
970     res="sub"
971     else:
972     val="arg"
973     res="res"
974     text+=" %s=Vector(0,w)\n"%val
975     text+=" if dim==2:\n"
976     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
977     text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap=" ")
978     text+="\n else:\n"
979    
980     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
981     text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap=" ")
982     text+="\n res=div(arg,where=w_ref)\n"
983     if data=="Symbol":
984     text+=" sub=res.substitute({arg:s})\n"
985     text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
986     text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
987     text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
988     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
989 gross 437
990 gross 441
991     t_prog+=text
992     print t_prog
993     1/0
994    
995 gross 439 #=======================================================================================================
996     # interpolation
997     #=======================================================================================================
998     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
999     for data in ["Data","Symbol"]:
1000     for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1001     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1002     if where==case or \
1003     ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
1004     ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
1005     (case=="ContinuousFunction" and where in ["Solution","ReducedSolution"]) or \
1006     (case=="Solution" and where=="ReducedSolution") :
1007    
1008    
1009     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1010     tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1011     text+=" def %s(self):\n"%tname
1012     text+=" \"\"\"\n"
1013     text+=" tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
1014     text+=" assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
1015     text+=" \"\"\"\n"
1016     if case=="ReducedSolution" or where=="ReducedSolution":
1017     text+=" o=1\n"
1018     o="1"
1019     else:
1020     text+=" o=self.order\n"
1021     o="o"
1022     text+=" dim=self.domain.getDim()\n"
1023     text+=" w_ref=%s(self.domain)\n"%where
1024     text+=" x_ref=w_ref.getX()\n"
1025     text+=" w=%s(self.domain)\n"%case
1026     text+=" x=w.getX()\n"
1027     a_2=makeArray(sh+(2,),[-1.,1])
1028     b_2=makeArray(sh+(2,),[-1.,1])
1029     a_3=makeArray(sh+(3,),[-1.,1])
1030     b_3=makeArray(sh+(3,),[-1.,1])
1031     if data=="Symbol":
1032     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1033     val="s"
1034     res="sub"
1035     else:
1036     val="arg"
1037     res="res"
1038     text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1039     text+=" ref=Data(0,%s,w_ref)\n"%str(sh)
1040     text+=" if dim==2:\n"
1041     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1042     text+=unrollLoops(a_2,b_2,o,"ref",tap=" ",x="x_ref")
1043     text+=" else:\n"
1044    
1045     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1046     text+=unrollLoops(a_3,b_3,o,"ref",tap=" ",x="x_ref")
1047     text+=" res=interpolate(arg,where=w_ref)\n"
1048     if data=="Symbol":
1049     text+=" sub=res.substitute({arg:s})\n"
1050     text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1051     text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
1052     text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1053     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1054     t_prog+=text
1055     print test_header
1056     print t_prog
1057     print test_tail
1058     1/0
1059 gross 437
1060 gross 429 #=======================================================================================================
1061 gross 438 # grad
1062 gross 437 #=======================================================================================================
1063     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1064     for data in ["Data","Symbol"]:
1065 gross 438 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1066     for sh in [ (),(2,), (4,5), (6,2,2)]:
1067 gross 437 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1068 gross 438 tname="test_grad_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1069 gross 437 text+=" def %s(self):\n"%tname
1070     text+=" \"\"\"\n"
1071 gross 438 text+=" tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
1072 gross 437 text+=" assumptions: %s(self.domain) exists\n"%case
1073 gross 438 text+=" self.domain supports gardient on %s\n"%where
1074 gross 437 text+=" \"\"\"\n"
1075     if case=="ReducedSolution":
1076     text+=" o=1\n"
1077     o="1"
1078     else:
1079     text+=" o=self.order\n"
1080     o="o"
1081     text+=" dim=self.domain.getDim()\n"
1082     text+=" w_ref=%s(self.domain)\n"%where
1083 gross 438 text+=" x_ref=w_ref.getX()\n"
1084 gross 437 text+=" w=%s(self.domain)\n"%case
1085     text+=" x=w.getX()\n"
1086     a_2=makeArray(sh+(2,),[-1.,1])
1087     b_2=makeArray(sh+(2,),[-1.,1])
1088     a_3=makeArray(sh+(3,),[-1.,1])
1089     b_3=makeArray(sh+(3,),[-1.,1])
1090     if data=="Symbol":
1091 gross 438 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1092 gross 437 val="s"
1093     res="sub"
1094     else:
1095     val="arg"
1096     res="res"
1097     text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1098 gross 438 text+=" ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1099 gross 437 text+=" if dim==2:\n"
1100     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1101 gross 438 text+=unrollLoopsOfGrad(a_2,b_2,o,"ref",tap=" ")
1102 gross 437 text+=" else:\n"
1103    
1104     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1105 gross 438 text+=unrollLoopsOfGrad(a_3,b_3,o,"ref",tap=" ")
1106     text+=" res=grad(arg,where=w_ref)\n"
1107 gross 437 if data=="Symbol":
1108     text+=" sub=res.substitute({arg:s})\n"
1109 gross 438 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1110     text+=" self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1111 gross 437 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1112    
1113    
1114     t_prog+=text
1115     print test_header
1116     print t_prog
1117     print test_tail
1118     1/0
1119 gross 438
1120    
1121 gross 437 #=======================================================================================================
1122 gross 438 # integrate
1123 gross 437 #=======================================================================================================
1124 gross 438 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1125 gross 437 for data in ["Data","Symbol"]:
1126 gross 438 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1127     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1128     if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:
1129 gross 437 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1130 gross 438 tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1131 gross 437 text+=" def %s(self):\n"%tname
1132     text+=" \"\"\"\n"
1133 gross 438 text+=" tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
1134     text+=" assumptions: %s(self.domain) exists\n"%case
1135     text+=" self.domain supports integral on %s\n"%where
1136    
1137 gross 437 text+=" \"\"\"\n"
1138     if case=="ReducedSolution":
1139     text+=" o=1\n"
1140     o="1"
1141     else:
1142     text+=" o=self.order\n"
1143     o="o"
1144     text+=" dim=self.domain.getDim()\n"
1145     text+=" w_ref=%s(self.domain)\n"%where
1146     text+=" w=%s(self.domain)\n"%case
1147     text+=" x=w.getX()\n"
1148     a_2=makeArray(sh+(2,),[-1.,1])
1149     b_2=makeArray(sh+(2,),[-1.,1])
1150     a_3=makeArray(sh+(3,),[-1.,1])
1151     b_3=makeArray(sh+(3,),[-1.,1])
1152     if data=="Symbol":
1153 gross 438 text+=" arg=Symbol(shape=%s)\n"%str(sh)
1154 gross 437 val="s"
1155     res="sub"
1156     else:
1157     val="arg"
1158     res="res"
1159 gross 438
1160 gross 437 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1161 gross 438 if not len(sh)==0:
1162     text+=" ref=numarray.zeros(%s,numarray.Float)\n"%str(sh)
1163 gross 437 text+=" if dim==2:\n"
1164     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1165 gross 438 text+=unrollLoopsOfInteriorIntegral(a_2,b_2,where,"ref",tap=" ")
1166 gross 437 text+=" else:\n"
1167    
1168     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1169 gross 438 text+=unrollLoopsOfInteriorIntegral(a_3,b_3,where,"ref",tap=" ")
1170     if case in ["ContinuousFunction","Solution","ReducedSolution"]:
1171     text+=" res=integrate(arg,where=w_ref)\n"
1172     else:
1173     text+=" res=integrate(arg)\n"
1174    
1175 gross 437 if data=="Symbol":
1176     text+=" sub=res.substitute({arg:s})\n"
1177 gross 438 if len(sh)==0 and data=="Data":
1178     text+=" self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%res
1179     else:
1180     if data=="Symbol":
1181     text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1182     text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1183     else:
1184     text+=" self.failUnless(isinstance(res,numarray.NumArray),\"wrong type of result.\")\n"
1185     text+=" self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
1186 gross 437 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1187    
1188    
1189     t_prog+=text
1190     print test_header
1191     print t_prog
1192     print test_tail
1193     1/0
1194     #=======================================================================================================
1195 gross 433 # inverse
1196     #=======================================================================================================
1197     name="inverse"
1198     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1199     for sh0 in [ (1,1), (2,2), (3,3)]:
1200     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1201     tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1202     text+=" def %s(self):\n"%tname
1203     a_0=makeArray(sh0,[-1.,1])
1204     for i in range(sh0[0]): a_0[i,i]+=2
1205     if case0 in ["taggedData", "expandedData"]:
1206     a1_0=makeArray(sh0,[-1.,1])
1207     for i in range(sh0[0]): a1_0[i,i]+=3
1208     else:
1209     a1_0=a_0
1210    
1211     text+=mkText(case0,"arg",a_0,a1_0)
1212     text+=" res=%s(arg)\n"%name
1213     if case0=="Symbol":
1214     text+=mkText("array","s",a_0,a1_0)
1215     text+=" sub=res.substitute({arg:s})\n"
1216     res="sub"
1217     ref="s"
1218     else:
1219     ref="arg"
1220     res="res"
1221     text+=mkTypeAndShapeTest(case0,sh0,"res")
1222     text+=" self.failUnless(Lsup(matrixmult(%s,%s)-kronecker(%s))<=self.RES_TOL,\"wrong result\")\n"%(res,ref,sh0[0])
1223    
1224     if case0 == "taggedData" :
1225     t_prog_with_tags+=text
1226     else:
1227     t_prog+=text
1228    
1229     print test_header
1230     # print t_prog
1231     print t_prog_with_tags
1232     print test_tail
1233     1/0
1234    
1235     #=======================================================================================================
1236 gross 429 # trace
1237     #=======================================================================================================
1238     def traceTest(r,offset):
1239     sh=r.shape
1240     r1=1
1241     for i in range(offset): r1*=sh[i]
1242     r2=1
1243     for i in range(offset+2,len(sh)): r2*=sh[i]
1244     r_s=numarray.reshape(r,(r1,sh[offset],sh[offset],r2))
1245     s=numarray.zeros([r1,r2],numarray.Float)
1246     for i1 in range(r1):
1247     for i2 in range(r2):
1248     for j in range(sh[offset]): s[i1,i2]+=r_s[i1,j,j,i2]
1249     return s.resize(sh[:offset]+sh[offset+2:])
1250     name,tt="trace",traceTest
1251     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1252     for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
1253     for offset in range(len(sh0)-1):
1254     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1255     tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1256     text+=" def %s(self):\n"%tname
1257     sh_t=list(sh0)
1258     sh_t[offset+1]=sh_t[offset]
1259     sh_t=tuple(sh_t)
1260     sh_r=[]
1261     for i in range(offset): sh_r.append(sh0[i])
1262     for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1263     sh_r=tuple(sh_r)
1264     a_0=makeArray(sh_t,[-1.,1])
1265     if case0 in ["taggedData", "expandedData"]:
1266     a1_0=makeArray(sh_t,[-1.,1])
1267     else:
1268     a1_0=a_0
1269     r=tt(a_0,offset)
1270     r1=tt(a1_0,offset)
1271     text+=mkText(case0,"arg",a_0,a1_0)
1272     text+=" res=%s(arg,%s)\n"%(name,offset)
1273     if case0=="Symbol":
1274     text+=mkText("array","s",a_0,a1_0)
1275     text+=" sub=res.substitute({arg:s})\n"
1276     res="sub"
1277     text+=mkText("array","ref",r,r1)
1278     else:
1279     res="res"
1280     text+=mkText(case0,"ref",r,r1)
1281     text+=mkTypeAndShapeTest(case0,sh_r,"res")
1282     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1283    
1284     if case0 == "taggedData" :
1285     t_prog_with_tags+=text
1286     else:
1287     t_prog+=text
1288 gross 396
1289 gross 429 print test_header
1290     # print t_prog
1291     print t_prog_with_tags
1292     print test_tail
1293     1/0
1294 gross 396
1295 gross 157 #=======================================================================================================
1296 gross 396 # clip
1297     #=======================================================================================================
1298     oper_L=[["clip",clipTEST]]
1299     for oper in oper_L:
1300     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1301     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1302     if len(sh0)==0 or not case0=="float":
1303     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1304     tname="test_%s_%s_rank%s"%(oper[0],case0,len(sh0))
1305     text+=" def %s(self):\n"%tname
1306     a_0=makeArray(sh0,[-1.,1])
1307     if case0 in ["taggedData", "expandedData"]:
1308     a1_0=makeArray(sh0,[-1.,1])
1309     else:
1310     a1_0=a_0
1311    
1312     r=oper[1](a_0,-0.3,0.5)
1313     r1=oper[1](a1_0,-0.3,0.5)
1314     text+=mkText(case0,"arg",a_0,a1_0)
1315     text+=" res=%s(arg,-0.3,0.5)\n"%oper[0]
1316     if case0=="Symbol":
1317     text+=mkText("array","s",a_0,a1_0)
1318     text+=" sub=res.substitute({arg:s})\n"
1319     res="sub"
1320     text+=mkText("array","ref",r,r1)
1321     else:
1322     res="res"
1323     text+=mkText(case0,"ref",r,r1)
1324     text+=mkTypeAndShapeTest(case0,sh0,"res")
1325     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1326    
1327     if case0 == "taggedData" :
1328     t_prog_with_tags+=text
1329     else:
1330     t_prog+=text
1331    
1332     print test_header
1333     # print t_prog
1334     print t_prog_with_tags
1335     print test_tail
1336     1/0
1337 gross 429
1338 gross 396 #=======================================================================================================
1339     # maximum, minimum, clipping
1340     #=======================================================================================================
1341     oper_L=[ ["maximum",maximumTEST],
1342     ["minimum",minimumTEST]]
1343     for oper in oper_L:
1344     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1345     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1346     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1347     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1348     if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1349     and (sh0==sh1 or len(sh0)==0 or len(sh1)==0) :
1350     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1351    
1352     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1353     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1354     text+=" def %s(self):\n"%tname
1355     a_0=makeArray(sh0,[-1.,1])
1356     if case0 in ["taggedData", "expandedData"]:
1357     a1_0=makeArray(sh0,[-1.,1])
1358     else:
1359     a1_0=a_0
1360    
1361     a_1=makeArray(sh1,[-1.,1])
1362     if case1 in ["taggedData", "expandedData"]:
1363     a1_1=makeArray(sh1,[-1.,1])
1364     else:
1365     a1_1=a_1
1366     r=oper[1](a_0,a_1)
1367     r1=oper[1](a1_0,a1_1)
1368     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1369     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1370     text+=" res=%s(arg0,arg1)\n"%oper[0]
1371     case=getResultCaseForBin(case0,case1)
1372     if case=="Symbol":
1373     c0_res,c1_res=case0,case1
1374     subs="{"
1375     if case0=="Symbol":
1376     text+=mkText("array","s0",a_0,a1_0)
1377     subs+="arg0:s0"
1378     c0_res="array"
1379     if case1=="Symbol":
1380     text+=mkText("array","s1",a_1,a1_1)
1381     if not subs.endswith("{"): subs+=","
1382     subs+="arg1:s1"
1383     c1_res="array"
1384     subs+="}"
1385     text+=" sub=res.substitute(%s)\n"%subs
1386     res="sub"
1387     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1388     else:
1389     res="res"
1390     text+=mkText(case,"ref",r,r1)
1391     if len(sh0)>len(sh1):
1392     text+=mkTypeAndShapeTest(case,sh0,"res")
1393     else:
1394     text+=mkTypeAndShapeTest(case,sh1,"res")
1395     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1396    
1397     if case0 == "taggedData" or case1 == "taggedData":
1398     t_prog_with_tags+=text
1399     else:
1400     t_prog+=text
1401    
1402     print test_header
1403     # print t_prog
1404     print t_prog_with_tags
1405     print test_tail
1406     1/0
1407    
1408    
1409     #=======================================================================================================
1410     # outer inner
1411     #=======================================================================================================
1412     oper=["outer",outerTEST]
1413     # oper=["inner",innerTEST]
1414     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1415     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1416     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1417     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1418     if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1419     and len(sh0+sh1)<5:
1420     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1421    
1422     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1423     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1424     text+=" def %s(self):\n"%tname
1425     a_0=makeArray(sh0,[-1.,1])
1426     if case0 in ["taggedData", "expandedData"]:
1427     a1_0=makeArray(sh0,[-1.,1])
1428     else:
1429     a1_0=a_0
1430    
1431     a_1=makeArray(sh1,[-1.,1])
1432     if case1 in ["taggedData", "expandedData"]:
1433     a1_1=makeArray(sh1,[-1.,1])
1434     else:
1435     a1_1=a_1
1436     r=oper[1](a_0,a_1)
1437     r1=oper[1](a1_0,a1_1)
1438     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1439     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1440     text+=" res=%s(arg0,arg1)\n"%oper[0]
1441     case=getResultCaseForBin(case0,case1)
1442     if case=="Symbol":
1443     c0_res,c1_res=case0,case1
1444     subs="{"
1445     if case0=="Symbol":
1446     text+=mkText("array","s0",a_0,a1_0)
1447     subs+="arg0:s0"
1448     c0_res="array"
1449     if case1=="Symbol":
1450     text+=mkText("array","s1",a_1,a1_1)
1451     if not subs.endswith("{"): subs+=","
1452     subs+="arg1:s1"
1453     c1_res="array"
1454     subs+="}"
1455     text+=" sub=res.substitute(%s)\n"%subs
1456     res="sub"
1457     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1458     else:
1459     res="res"
1460     text+=mkText(case,"ref",r,r1)
1461     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1462     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1463    
1464     if case0 == "taggedData" or case1 == "taggedData":
1465     t_prog_with_tags+=text
1466     else:
1467     t_prog+=text
1468    
1469     print test_header
1470     # print t_prog
1471     print t_prog_with_tags
1472     print test_tail
1473     1/0
1474    
1475     #=======================================================================================================
1476 gross 313 # local reduction
1477     #=======================================================================================================
1478     for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],
1479     ["maxval",-1.e99,"max(out,%a1%)","out"],
1480     ["minval",1.e99,"min(out,%a1%)","out"] ]:
1481     for case in case_set:
1482     for sh in shape_set:
1483     if not case=="float" or len(sh)==0:
1484     text=""
1485     text+=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1486     tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))
1487     text+=" %s(self):\n"%tname
1488     a=makeArray(sh,[-1.,1.])
1489     a1=makeArray(sh,[-1.,1.])
1490     r1=testReduce(a1,oper[1],oper[2],oper[3])
1491     r=testReduce(a,oper[1],oper[2],oper[3])
1492    
1493     text+=mkText(case,"arg",a,a1)
1494     text+=" res=%s(arg)\n"%oper[0]
1495     if case=="Symbol":
1496     text+=mkText("array","s",a,a1)
1497     text+=" sub=res.substitute({arg:s})\n"
1498     text+=mkText("array","ref",r,r1)
1499     res="sub"
1500     else:
1501     text+=mkText(case,"ref",r,r1)
1502     res="res"
1503     if oper[0]=="length":
1504     text+=mkTypeAndShapeTest(case,(),"res")
1505     else:
1506     if case=="float" or case=="array":
1507     text+=mkTypeAndShapeTest("float",(),"res")
1508     else:
1509     text+=mkTypeAndShapeTest(case,(),"res")
1510     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1511     if case == "taggedData":
1512     t_prog_with_tags+=text
1513     else:
1514     t_prog+=text
1515     print test_header
1516     # print t_prog
1517     print t_prog_with_tags
1518     print test_tail
1519     1/0
1520    
1521     #=======================================================================================================
1522 gross 291 # tensor multiply
1523 gross 157 #=======================================================================================================
1524 gross 291 # oper=["generalTensorProduct",tensorProductTest]
1525     # oper=["matrixmult",testMatrixMult]
1526     oper=["tensormult",testTensorMult]
1527    
1528     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1529     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1530     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1531     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1532     for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
1533     if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
1534     and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
1535     # 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
1536     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
1537     case=getResultCaseForBin(case0,case1)
1538 gross 157 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1539     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1540 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))
1541     #tname="test_matrixmult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1542     tname="test_tensormult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1543     # if tname=="test_generalTensorProduct_array_rank1_array_rank2_offset1":
1544     # print tnametest_generalTensorProduct_Symbol_rank1_Symbol_rank3_offset1
1545     text+=" def %s(self):\n"%tname
1546     a_0=makeArray(sh0+sh_s,[-1.,1])
1547     if case0 in ["taggedData", "expandedData"]:
1548     a1_0=makeArray(sh0+sh_s,[-1.,1])
1549     else:
1550     a1_0=a_0
1551    
1552     a_1=makeArray(sh_s+sh1,[-1.,1])
1553     if case1 in ["taggedData", "expandedData"]:
1554     a1_1=makeArray(sh_s+sh1,[-1.,1])
1555     else:
1556     a1_1=a_1
1557     r=oper[1](a_0,a_1,sh_s)
1558     r1=oper[1](a1_0,a1_1,sh_s)
1559     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1560     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1561     #text+=" res=matrixmult(arg0,arg1)\n"
1562     text+=" res=tensormult(arg0,arg1)\n"
1563     #text+=" res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))
1564     if case=="Symbol":
1565     c0_res,c1_res=case0,case1
1566     subs="{"
1567     if case0=="Symbol":
1568     text+=mkText("array","s0",a_0,a1_0)
1569     subs+="arg0:s0"
1570     c0_res="array"
1571     if case1=="Symbol":
1572     text+=mkText("array","s1",a_1,a1_1)
1573     if not subs.endswith("{"): subs+=","
1574     subs+="arg1:s1"
1575     c1_res="array"
1576     subs+="}"
1577     text+=" sub=res.substitute(%s)\n"%subs
1578     res="sub"
1579     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1580     else:
1581     res="res"
1582     text+=mkText(case,"ref",r,r1)
1583     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1584     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1585     if case0 == "taggedData" or case1 == "taggedData":
1586     t_prog_with_tags+=text
1587     else:
1588     t_prog+=text
1589     print test_header
1590     # print t_prog
1591     print t_prog_with_tags
1592     print test_tail
1593     1/0
1594     #=======================================================================================================
1595 gross 157 # basic binary operation overloading (tests only!)
1596     #=======================================================================================================
1597     oper_range=[-5.,5.]
1598     for oper in [["add" ,"+",[-5.,5.]],
1599     ["sub" ,"-",[-5.,5.]],
1600     ["mult","*",[-5.,5.]],
1601     ["div" ,"/",[-5.,5.]],
1602     ["pow" ,"**",[0.01,5.]]]:
1603     for case0 in case_set:
1604     for sh0 in shape_set:
1605     for case1 in case_set:
1606     for sh1 in shape_set:
1607 gross 291 if not case0=="array" and \
1608     (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1609 gross 157 (sh0==() or sh1==() or sh1==sh0) and \
1610     not (case0 in ["float","array"] and case1 in ["float","array"]):
1611 gross 291 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1612 gross 157 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1613     tname="test_%s_overloaded_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1614     text+=" def %s(self):\n"%tname
1615     a_0=makeArray(sh0,oper[2])
1616     if case0 in ["taggedData", "expandedData"]:
1617     a1_0=makeArray(sh0,oper[2])
1618     else:
1619     a1_0=a_0
1620 jgs 154
1621 gross 157 a_1=makeArray(sh1,oper[2])
1622     if case1 in ["taggedData", "expandedData"]:
1623     a1_1=makeArray(sh1,oper[2])
1624 jgs 154 else:
1625 gross 157 a1_1=a_1
1626     r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1627     r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1628     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1629     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1630     text+=" res=arg0%sarg1\n"%oper[1]
1631    
1632     case=getResultCaseForBin(case0,case1)
1633     if case=="Symbol":
1634     c0_res,c1_res=case0,case1
1635     subs="{"
1636     if case0=="Symbol":
1637     text+=mkText("array","s0",a_0,a1_0)
1638     subs+="arg0:s0"
1639     c0_res="array"
1640     if case1=="Symbol":
1641     text+=mkText("array","s1",a_1,a1_1)
1642     if not subs.endswith("{"): subs+=","
1643     subs+="arg1:s1"
1644     c1_res="array"
1645     subs+="}"
1646     text+=" sub=res.substitute(%s)\n"%subs
1647     res="sub"
1648     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1649 jgs 154 else:
1650 gross 157 res="res"
1651     text+=mkText(case,"ref",r,r1)
1652     if isinstance(r,float):
1653     text+=mkTypeAndShapeTest(case,(),"res")
1654     else:
1655     text+=mkTypeAndShapeTest(case,r.shape,"res")
1656     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1657    
1658     if case0 in [ "constData","taggedData","expandedData"] and case1 == "Symbol":
1659     t_prog_failing+=text
1660     else:
1661     if case0 == "taggedData" or case1 == "taggedData":
1662     t_prog_with_tags+=text
1663     else:
1664     t_prog+=text
1665 jgs 154
1666 gross 157
1667     print test_header
1668 gross 291 # print t_prog
1669     # print t_prog_with_tags
1670     print t_prog_failing
1671     print test_tail
1672 jgs 154 1/0
1673 gross 291 #=======================================================================================================
1674     # basic binary operations (tests only!)
1675     #=======================================================================================================
1676     oper_range=[-5.,5.]
1677     for oper in [["add" ,"+",[-5.,5.]],
1678     ["mult","*",[-5.,5.]],
1679     ["quotient" ,"/",[-5.,5.]],
1680     ["power" ,"**",[0.01,5.]]]:
1681     for case0 in case_set:
1682     for case1 in case_set:
1683     for sh in shape_set:
1684     for sh_p in shape_set:
1685     if len(sh_p)>0:
1686     resource=[-1,1]
1687     else:
1688     resource=[1]
1689     for sh_d in resource:
1690     if sh_d>0:
1691     sh0=sh
1692     sh1=sh+sh_p
1693     else:
1694     sh1=sh
1695     sh0=sh+sh_p
1696    
1697     if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1698     len(sh0)<5 and len(sh1)<5:
1699     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1700     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1701     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1702     text+=" def %s(self):\n"%tname
1703     a_0=makeArray(sh0,oper[2])
1704     if case0 in ["taggedData", "expandedData"]:
1705     a1_0=makeArray(sh0,oper[2])
1706     else:
1707     a1_0=a_0
1708    
1709     a_1=makeArray(sh1,oper[2])
1710     if case1 in ["taggedData", "expandedData"]:
1711     a1_1=makeArray(sh1,oper[2])
1712     else:
1713     a1_1=a_1
1714     r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1715     r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1716     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1717     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1718     text+=" res=%s(arg0,arg1)\n"%oper[0]
1719    
1720     case=getResultCaseForBin(case0,case1)
1721     if case=="Symbol":
1722     c0_res,c1_res=case0,case1
1723     subs="{"
1724     if case0=="Symbol":
1725     text+=mkText("array","s0",a_0,a1_0)
1726     subs+="arg0:s0"
1727     c0_res="array"
1728     if case1=="Symbol":
1729     text+=mkText("array","s1",a_1,a1_1)
1730     if not subs.endswith("{"): subs+=","
1731     subs+="arg1:s1"
1732     c1_res="array"
1733     subs+="}"
1734     text+=" sub=res.substitute(%s)\n"%subs
1735     res="sub"
1736     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1737     else:
1738     res="res"
1739     text+=mkText(case,"ref",r,r1)
1740     if isinstance(r,float):
1741     text+=mkTypeAndShapeTest(case,(),"res")
1742     else:
1743     text+=mkTypeAndShapeTest(case,r.shape,"res")
1744     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1745    
1746     if case0 == "taggedData" or case1 == "taggedData":
1747     t_prog_with_tags+=text
1748     else:
1749     t_prog+=text
1750     print test_header
1751     # print t_prog
1752     print t_prog_with_tags
1753     print test_tail
1754     1/0
1755    
1756 gross 157 # print t_prog_with_tagsoper_range=[-5.,5.]
1757     for oper in [["add" ,"+",[-5.,5.]],
1758     ["sub" ,"-",[-5.,5.]],
1759     ["mult","*",[-5.,5.]],
1760     ["div" ,"/",[-5.,5.]],
1761     ["pow" ,"**",[0.01,5.]]]:
1762     for case0 in case_set:
1763     for sh0 in shape_set:
1764     for case1 in case_set:
1765     for sh1 in shape_set:
1766     if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1767     (sh0==() or sh1==() or sh1==sh0) and \
1768     not (case0 in ["float","array"] and case1 in ["float","array"]):
1769     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1770     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1771     text+=" def %s(self):\n"%tname
1772     a_0=makeArray(sh0,oper[2])
1773     if case0 in ["taggedData", "expandedData"]:
1774     a1_0=makeArray(sh0,oper[2])
1775     else:
1776     a1_0=a_0
1777 jgs 154
1778 gross 157 a_1=makeArray(sh1,oper[2])
1779     if case1 in ["taggedData", "expandedData"]:
1780     a1_1=makeArray(sh1,oper[2])
1781 jgs 154 else:
1782 gross 157 a1_1=a_1
1783     r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1784     r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1785     text+=mkText(case0,"arg0",a_0,a1_0)
1786     text+=mkText(case1,"arg1",a_1,a1_1)
1787     text+=" res=arg0%sarg1\n"%oper[1]
1788    
1789     case=getResultCaseForBin(case0,case1)
1790     if case=="Symbol":
1791     c0_res,c1_res=case0,case1
1792     subs="{"
1793     if case0=="Symbol":
1794     text+=mkText("array","s0",a_0,a1_0)
1795     subs+="arg0:s0"
1796     c0_res="array"
1797     if case1=="Symbol":
1798     text+=mkText("array","s1",a_1,a1_1)
1799     if not subs.endswith("{"): subs+=","
1800     subs+="arg1:s1"
1801     c1_res="array"
1802     subs+="}"
1803     text+=" sub=res.substitute(%s)\n"%subs
1804     res="sub"
1805     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1806 jgs 154 else:
1807 gross 157 res="res"
1808     text+=mkText(case,"ref",r,r1)
1809     if isinstance(r,float):
1810     text+=mkTypeAndShapeTest(case,(),"res")
1811     else:
1812     text+=mkTypeAndShapeTest(case,r.shape,"res")
1813     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1814    
1815     if case0 in [ "constData","taggedData","expandedData"] and case1 == "Symbol":
1816     t_prog_failing+=text
1817     else:
1818     if case0 == "taggedData" or case1 == "taggedData":
1819     t_prog_with_tags+=text
1820     else:
1821     t_prog+=text
1822 jgs 154
1823 gross 157
1824     # print u_prog
1825     # 1/0
1826     print test_header
1827 jgs 154 print t_prog
1828 gross 157 # print t_prog_with_tags
1829     # print t_prog_failing
1830     print test_tail
1831     # print t_prog_failing
1832     print test_tail
1833 jgs 154
1834 gross 157 #=======================================================================================================
1835     # unary operations:
1836     #=======================================================================================================
1837     func= [
1838     OPERATOR(nickname="log10",\
1839     rng=[1.e-3,100.],\
1840     test_expr="math.log10(%a1%)",\
1841     math_expr="math.log10(%a1%)",\
1842     numarray_expr="numarray.log10(%a1%)",\
1843     symbol_expr="log(%a1%)/log(10.)",\
1844     name="base-10 logarithm"),
1845     OPERATOR(nickname="wherePositive",\
1846     rng=[-100.,100.],\
1847     test_expr="wherepos(%a1%)",\
1848     math_expr="if arg>0:\n return 1.\nelse:\n return 0.",
1849     numarray_expr="numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))",\
1850     name="mask of positive values"),
1851     OPERATOR(nickname="whereNegative",\
1852     rng=[-100.,100.],\
1853     test_expr="wherepos(-%a1%)",\
1854     math_expr="if arg<0:\n return 1.\nelse:\n return 0.",
1855     numarray_expr="numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))",\
1856     name="mask of positive values"),
1857     OPERATOR(nickname="whereNonNegative",\
1858     rng=[-100.,100.],\
1859     test_expr="1-wherepos(-%a1%)", \
1860     math_expr="if arg<0:\n return 0.\nelse:\n return 1.",
1861     numarray_expr="numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))",\
1862     symbol_expr="1-wherePositive(%a1%)",\
1863     name="mask of non-negative values"),
1864     OPERATOR(nickname="whereNonPositive",\
1865     rng=[-100.,100.],\
1866     test_expr="1-wherepos(%a1%)",\
1867     math_expr="if arg>0:\n return 0.\nelse:\n return 1.",
1868     numarray_expr="numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))",\
1869     symbol_expr="1-whereNegative(%a1%)",\
1870