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

Annotation of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


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