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

Annotation of /trunk/escript/py_src/generatediff

Parent Directory Parent Directory | Revision Log Revision Log


Revision 396 - (hide annotations)
Wed Dec 21 05:08:25 2005 UTC (13 years, 10 months ago) by gross
Original Path: trunk/escript/py_src/generateutil
File size: 98895 byte(s)
tests for clip, maximum, minimum 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     else:
177     raise SystemError,"rank is restricted to 4"
178     return out
179    
180    
181 gross 157 def makeResult(val,test_expr):
182 jgs 154 if isinstance(val,float):
183 gross 157 out=eval(test_expr.replace("%a1%","val"))
184 jgs 154 elif len(val.shape)==0:
185 gross 157 out=eval(test_expr.replace("%a1%","val"))
186 jgs 154 elif len(val.shape)==1:
187     out=numarray.zeros(val.shape,numarray.Float64)
188     for i0 in range(val.shape[0]):
189 gross 157 out[i0]=eval(test_expr.replace("%a1%","val[i0]"))
190 jgs 154 elif len(val.shape)==2:
191     out=numarray.zeros(val.shape,numarray.Float64)
192     for i0 in range(val.shape[0]):
193     for i1 in range(val.shape[1]):
194 gross 157 out[i0,i1]=eval(test_expr.replace("%a1%","val[i0,i1]"))
195 jgs 154 elif len(val.shape)==3:
196     out=numarray.zeros(val.shape,numarray.Float64)
197     for i0 in range(val.shape[0]):
198     for i1 in range(val.shape[1]):
199     for i2 in range(val.shape[2]):
200 gross 157 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val[i0,i1,i2]"))
201 jgs 154 elif len(val.shape)==4:
202     out=numarray.zeros(val.shape,numarray.Float64)
203     for i0 in range(val.shape[0]):
204     for i1 in range(val.shape[1]):
205     for i2 in range(val.shape[2]):
206     for i3 in range(val.shape[3]):
207 gross 157 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val[i0,i1,i2,i3]"))
208 jgs 154 else:
209     raise SystemError,"rank is restricted to 4"
210     return out
211    
212 gross 157 def makeResult2(val0,val1,test_expr):
213     if isinstance(val0,float):
214     if isinstance(val1,float):
215     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
216     elif len(val1.shape)==0:
217     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
218     elif len(val1.shape)==1:
219     out=numarray.zeros(val1.shape,numarray.Float64)
220     for i0 in range(val1.shape[0]):
221     out[i0]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0]"))
222     elif len(val1.shape)==2:
223     out=numarray.zeros(val1.shape,numarray.Float64)
224     for i0 in range(val1.shape[0]):
225     for i1 in range(val1.shape[1]):
226     out[i0,i1]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1]"))
227     elif len(val1.shape)==3:
228     out=numarray.zeros(val1.shape,numarray.Float64)
229     for i0 in range(val1.shape[0]):
230     for i1 in range(val1.shape[1]):
231     for i2 in range(val1.shape[2]):
232     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2]"))
233     elif len(val1.shape)==4:
234     out=numarray.zeros(val1.shape,numarray.Float64)
235     for i0 in range(val1.shape[0]):
236     for i1 in range(val1.shape[1]):
237     for i2 in range(val1.shape[2]):
238     for i3 in range(val1.shape[3]):
239     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2,i3]"))
240     else:
241     raise SystemError,"rank of val1 is restricted to 4"
242     elif len(val0.shape)==0:
243     if isinstance(val1,float):
244     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
245     elif len(val1.shape)==0:
246     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
247     elif len(val1.shape)==1:
248     out=numarray.zeros(val1.shape,numarray.Float64)
249     for i0 in range(val1.shape[0]):
250     out[i0]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0]"))
251     elif len(val1.shape)==2:
252     out=numarray.zeros(val1.shape,numarray.Float64)
253     for i0 in range(val1.shape[0]):
254     for i1 in range(val1.shape[1]):
255     out[i0,i1]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1]"))
256     elif len(val1.shape)==3:
257     out=numarray.zeros(val1.shape,numarray.Float64)
258     for i0 in range(val1.shape[0]):
259     for i1 in range(val1.shape[1]):
260     for i2 in range(val1.shape[2]):
261     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2]"))
262     elif len(val1.shape)==4:
263     out=numarray.zeros(val1.shape,numarray.Float64)
264     for i0 in range(val1.shape[0]):
265     for i1 in range(val1.shape[1]):
266     for i2 in range(val1.shape[2]):
267     for i3 in range(val1.shape[3]):
268     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2,i3]"))
269     else:
270     raise SystemError,"rank of val1 is restricted to 4"
271     elif len(val0.shape)==1:
272     if isinstance(val1,float):
273     out=numarray.zeros(val0.shape,numarray.Float64)
274     for i0 in range(val0.shape[0]):
275     out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1"))
276     elif len(val1.shape)==0:
277     out=numarray.zeros(val0.shape,numarray.Float64)
278     for i0 in range(val0.shape[0]):
279     out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1"))
280     elif len(val1.shape)==1:
281     out=numarray.zeros(val0.shape,numarray.Float64)
282     for i0 in range(val0.shape[0]):
283     out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0]"))
284     elif len(val1.shape)==2:
285     out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
286     for i0 in range(val0.shape[0]):
287     for j1 in range(val1.shape[1]):
288     out[i0,j1]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1]"))
289     elif len(val1.shape)==3:
290     out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
291     for i0 in range(val0.shape[0]):
292     for j1 in range(val1.shape[1]):
293     for j2 in range(val1.shape[2]):
294     out[i0,j1,j2]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1,j2]"))
295     elif len(val1.shape)==4:
296     out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
297     for i0 in range(val0.shape[0]):
298     for j1 in range(val1.shape[1]):
299     for j2 in range(val1.shape[2]):
300     for j3 in range(val1.shape[3]):
301     out[i0,j1,j2,j3]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1,j2,j3]"))
302     else:
303     raise SystemError,"rank of val1 is restricted to 4"
304     elif len(val0.shape)==2:
305     if isinstance(val1,float):
306     out=numarray.zeros(val0.shape,numarray.Float64)
307     for i0 in range(val0.shape[0]):
308     for i1 in range(val0.shape[1]):
309     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1"))
310     elif len(val1.shape)==0:
311     out=numarray.zeros(val0.shape,numarray.Float64)
312     for i0 in range(val0.shape[0]):
313     for i1 in range(val0.shape[1]):
314     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1"))
315     elif len(val1.shape)==1:
316     out=numarray.zeros(val0.shape,numarray.Float64)
317     for i0 in range(val0.shape[0]):
318     for i1 in range(val0.shape[1]):
319     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0]"))
320     elif len(val1.shape)==2:
321     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
322     for i0 in range(val0.shape[0]):
323     for i1 in range(val0.shape[1]):
324     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1]"))
325     elif len(val1.shape)==3:
326     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
327     for i0 in range(val0.shape[0]):
328     for i1 in range(val0.shape[1]):
329     for j2 in range(val1.shape[2]):
330     out[i0,i1,j2]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1,j2]"))
331     elif len(val1.shape)==4:
332     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
333     for i0 in range(val0.shape[0]):
334     for i1 in range(val0.shape[1]):
335     for j2 in range(val1.shape[2]):
336     for j3 in range(val1.shape[3]):
337     out[i0,i1,j2,j3]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1,j2,j3]"))
338     else:
339     raise SystemError,"rank of val1 is restricted to 4"
340     elif len(val0.shape)==3:
341     if isinstance(val1,float):
342     out=numarray.zeros(val0.shape,numarray.Float64)
343     for i0 in range(val0.shape[0]):
344     for i1 in range(val0.shape[1]):
345     for i2 in range(val0.shape[2]):
346     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1"))
347     elif len(val1.shape)==0:
348     out=numarray.zeros(val0.shape,numarray.Float64)
349     for i0 in range(val0.shape[0]):
350     for i1 in range(val0.shape[1]):
351     for i2 in range(val0.shape[2]):
352     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1"))
353     elif len(val1.shape)==1:
354     out=numarray.zeros(val0.shape,numarray.Float64)
355     for i0 in range(val0.shape[0]):
356     for i1 in range(val0.shape[1]):
357     for i2 in range(val0.shape[2]):
358     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0]"))
359     elif len(val1.shape)==2:
360     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
361     for i0 in range(val0.shape[0]):
362     for i1 in range(val0.shape[1]):
363     for i2 in range(val0.shape[2]):
364     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1]"))
365     elif len(val1.shape)==3:
366     out=numarray.zeros(val0.shape,numarray.Float64)
367     for i0 in range(val0.shape[0]):
368     for i1 in range(val0.shape[1]):
369     for i2 in range(val0.shape[2]):
370     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1,i2]"))
371     elif len(val1.shape)==4:
372     out=numarray.zeros(val0.shape+val1.shape[3:],numarray.Float64)
373     for i0 in range(val0.shape[0]):
374     for i1 in range(val0.shape[1]):
375     for i2 in range(val0.shape[2]):
376     for j3 in range(val1.shape[3]):
377     out[i0,i1,i2,j3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1,i2,j3]"))
378     else:
379     raise SystemError,"rank of val1 is rargs[1]estricted to 4"
380     elif len(val0.shape)==4:
381     if isinstance(val1,float):
382     out=numarray.zeros(val0.shape,numarray.Float64)
383     for i0 in range(val0.shape[0]):
384     for i1 in range(val0.shape[1]):
385     for i2 in range(val0.shape[2]):
386     for i3 in range(val0.shape[3]):
387     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1"))
388     elif len(val1.shape)==0:
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)==1:
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[i0]"))
402     elif len(val1.shape)==2:
403     out=numarray.zeros(val0.shape+val1.shape[2:],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,i1]"))
409     elif len(val1.shape)==3:
410     out=numarray.zeros(val0.shape,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,i2]"))
416     elif len(val1.shape)==4:
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,i3]"))
423     else:
424     raise SystemError,"rank of val1 is restricted to 4"
425     else:
426     raise SystemError,"rank is restricted to 4"
427     return out
428 jgs 154
429 gross 157
430     def mkText(case,name,a,a1=None,use_tagging_for_expanded_data=False):
431 jgs 154 t_out=""
432 gross 157 if case=="float":
433 jgs 154 if isinstance(a,float):
434     t_out+=" %s=%s\n"%(name,a)
435 gross 291 elif a.rank==0:
436 jgs 154 t_out+=" %s=%s\n"%(name,a)
437     else:
438     t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
439 gross 157 elif case=="array":
440     if isinstance(a,float):
441     t_out+=" %s=numarray.array(%s)\n"%(name,a)
442 gross 291 elif a.rank==0:
443 gross 157 t_out+=" %s=numarray.array(%s)\n"%(name,a)
444     else:
445     t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
446 jgs 154 elif case=="constData":
447     if isinstance(a,float):
448     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
449 gross 291 elif a.rank==0:
450 jgs 154 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
451     else:
452     t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
453     elif case=="taggedData":
454     if isinstance(a,float):
455     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
456     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
457 gross 291 elif a.rank==0:
458 jgs 154 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
459     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
460     else:
461     t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
462 gross 157 t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
463 jgs 154 elif case=="expandedData":
464 gross 157 if use_tagging_for_expanded_data:
465     if isinstance(a,float):
466     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
467     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
468 gross 291 elif a.rank==0:
469 gross 157 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
470     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
471     else:
472     t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
473     t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
474     t_out+=" %s.expand()\n"%name
475     else:
476     t_out+=" msk_%s=whereNegative(self.functionspace.getX()[0]-0.5)\n"%name
477     if isinstance(a,float):
478     t_out+=" %s=msk_%s*(%s)+(1.-msk_%s)*(%s)\n"%(name,name,a,name,a1)
479 gross 291 elif a.rank==0:
480 gross 157 t_out+=" %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a,name,a1)
481     else:
482     t_out+=" %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a.tolist(),name,a1.tolist())
483 jgs 154 elif case=="Symbol":
484     if isinstance(a,float):
485     t_out+=" %s=Symbol(shape=())\n"%(name)
486 gross 291 elif a.rank==0:
487 jgs 154 t_out+=" %s=Symbol(shape=())\n"%(name)
488     else:
489     t_out+=" %s=Symbol(shape=%s)\n"%(name,str(a.shape))
490    
491     return t_out
492    
493 gross 157 def mkTypeAndShapeTest(case,sh,argstr):
494     text=""
495     if case=="float":
496     text+=" self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%argstr
497     elif case=="array":
498     text+=" self.failUnless(isinstance(%s,numarray.NumArray),\"wrong type of result.\")\n"%argstr
499     text+=" self.failUnlessEqual(%s.shape,%s,\"wrong shape of result.\")\n"%(argstr,str(sh))
500     elif case in ["constData","taggedData","expandedData"]:
501     text+=" self.failUnless(isinstance(%s,Data),\"wrong type of result.\")\n"%argstr
502     text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result.\")\n"%(argstr,str(sh))
503     elif case=="Symbol":
504     text+=" self.failUnless(isinstance(%s,Symbol),\"wrong type of result.\")\n"%argstr
505     text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result.\")\n"%(argstr,str(sh))
506     return text
507 jgs 154
508 gross 157 def mkCode(txt,args=[],intend=""):
509     s=txt.split("\n")
510     if len(s)>1:
511     out=""
512     for l in s:
513     out+=intend+l+"\n"
514     else:
515     out="%sreturn %s\n"%(intend,txt)
516     c=1
517     for r in args:
518     out=out.replace("%%a%s%%"%c,r)
519     return out
520 jgs 154
521 gross 291 def innerTEST(arg0,arg1):
522     if isinstance(arg0,float):
523     out=numarray.array(arg0*arg1)
524     else:
525     out=(arg0*arg1).sum()
526     return out
527    
528     def outerTEST(arg0,arg1):
529     if isinstance(arg0,float):
530     out=numarray.array(arg0*arg1)
531     elif isinstance(arg1,float):
532     out=numarray.array(arg0*arg1)
533     else:
534     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
535     return out
536    
537     def tensorProductTest(arg0,arg1,sh_s):
538     if isinstance(arg0,float):
539     out=numarray.array(arg0*arg1)
540     elif isinstance(arg1,float):
541     out=numarray.array(arg0*arg1)
542     elif len(sh_s)==0:
543     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
544     else:
545     l=len(sh_s)
546     sh0=arg0.shape[:arg0.rank-l]
547     sh1=arg1.shape[l:]
548     ls,l0,l1=1,1,1
549     for i in sh_s: ls*=i
550     for i in sh0: l0*=i
551     for i in sh1: l1*=i
552     out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
553     out2=numarray.zeros((l0,l1),numarray.Float)
554     for i0 in range(l0):
555     for i1 in range(l1):
556     for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
557     out=out2.resize(sh0+sh1)
558     return out
559    
560     def testMatrixMult(arg0,arg1,sh_s):
561     return numarray.matrixmultiply(arg0,arg1)
562    
563    
564     def testTensorMult(arg0,arg1,sh_s):
565     if len(arg0)==2:
566     return numarray.matrixmultiply(arg0,arg1)
567     else:
568     if arg1.rank==4:
569     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
570     for i0 in range(arg0.shape[0]):
571     for i1 in range(arg0.shape[1]):
572     for i2 in range(arg0.shape[2]):
573     for i3 in range(arg0.shape[3]):
574     for j2 in range(arg1.shape[2]):
575     for j3 in range(arg1.shape[3]):
576     out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
577     elif arg1.rank==3:
578     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
579     for i0 in range(arg0.shape[0]):
580     for i1 in range(arg0.shape[1]):
581     for i2 in range(arg0.shape[2]):
582     for i3 in range(arg0.shape[3]):
583     for j2 in range(arg1.shape[2]):
584     out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
585     elif arg1.rank==2:
586     out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
587     for i0 in range(arg0.shape[0]):
588     for i1 in range(arg0.shape[1]):
589     for i2 in range(arg0.shape[2]):
590     for i3 in range(arg0.shape[3]):
591     out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
592     return out
593 gross 313
594     def testReduce(arg0,init_val,test_expr,post_expr):
595     out=init_val
596     if isinstance(arg0,float):
597     out=eval(test_expr.replace("%a1%","arg0"))
598     elif arg0.rank==0:
599     out=eval(test_expr.replace("%a1%","arg0"))
600     elif arg0.rank==1:
601     for i0 in range(arg0.shape[0]):
602     out=eval(test_expr.replace("%a1%","arg0[i0]"))
603     elif arg0.rank==2:
604     for i0 in range(arg0.shape[0]):
605     for i1 in range(arg0.shape[1]):
606     out=eval(test_expr.replace("%a1%","arg0[i0,i1]"))
607     elif arg0.rank==3:
608     for i0 in range(arg0.shape[0]):
609     for i1 in range(arg0.shape[1]):
610     for i2 in range(arg0.shape[2]):
611     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2]"))
612     elif arg0.rank==4:
613     for i0 in range(arg0.shape[0]):
614     for i1 in range(arg0.shape[1]):
615     for i2 in range(arg0.shape[2]):
616     for i3 in range(arg0.shape[3]):
617     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))
618     return eval(post_expr)
619 gross 396
620     def clipTEST(arg0,mn,mx):
621     if isinstance(arg0,float):
622     return max(min(arg0,mx),mn)
623     out=numarray.zeros(arg0.shape,numarray.Float64)
624     if arg0.rank==1:
625     for i0 in range(arg0.shape[0]):
626     out[i0]=max(min(arg0[i0],mx),mn)
627     elif arg0.rank==2:
628     for i0 in range(arg0.shape[0]):
629     for i1 in range(arg0.shape[1]):
630     out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
631     elif arg0.rank==3:
632     for i0 in range(arg0.shape[0]):
633     for i1 in range(arg0.shape[1]):
634     for i2 in range(arg0.shape[2]):
635     out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
636     elif arg0.rank==4:
637     for i0 in range(arg0.shape[0]):
638     for i1 in range(arg0.shape[1]):
639     for i2 in range(arg0.shape[2]):
640     for i3 in range(arg0.shape[3]):
641     out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
642     return out
643     def minimumTEST(arg0,arg1):
644     if isinstance(arg0,float):
645     if isinstance(arg1,float):
646     if arg0>arg1:
647     return arg1
648     else:
649     return arg0
650     else:
651     arg0=numarray.ones(arg1.shape)*arg0
652     else:
653     if isinstance(arg1,float):
654     arg1=numarray.ones(arg0.shape)*arg1
655     out=numarray.zeros(arg0.shape,numarray.Float64)
656     if arg0.rank==0:
657     if arg0>arg1:
658     out=arg1
659     else:
660     out=arg0
661     elif arg0.rank==1:
662     for i0 in range(arg0.shape[0]):
663     if arg0[i0]>arg1[i0]:
664     out[i0]=arg1[i0]
665     else:
666     out[i0]=arg0[i0]
667     elif arg0.rank==2:
668     for i0 in range(arg0.shape[0]):
669     for i1 in range(arg0.shape[1]):
670     if arg0[i0,i1]>arg1[i0,i1]:
671     out[i0,i1]=arg1[i0,i1]
672     else:
673     out[i0,i1]=arg0[i0,i1]
674     elif arg0.rank==3:
675     for i0 in range(arg0.shape[0]):
676     for i1 in range(arg0.shape[1]):
677     for i2 in range(arg0.shape[2]):
678     if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
679     out[i0,i1,i2]=arg1[i0,i1,i2]
680     else:
681     out[i0,i1,i2]=arg0[i0,i1,i2]
682     elif arg0.rank==4:
683     for i0 in range(arg0.shape[0]):
684     for i1 in range(arg0.shape[1]):
685     for i2 in range(arg0.shape[2]):
686     for i3 in range(arg0.shape[3]):
687     if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
688     out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
689     else:
690     out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
691     return out
692    
693    
694 gross 157 #=======================================================================================================
695 gross 396 # clip
696     #=======================================================================================================
697     oper_L=[["clip",clipTEST]]
698     for oper in oper_L:
699     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
700     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
701     if len(sh0)==0 or not case0=="float":
702     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
703     tname="test_%s_%s_rank%s"%(oper[0],case0,len(sh0))
704     text+=" def %s(self):\n"%tname
705     a_0=makeArray(sh0,[-1.,1])
706     if case0 in ["taggedData", "expandedData"]:
707     a1_0=makeArray(sh0,[-1.,1])
708     else:
709     a1_0=a_0
710    
711     r=oper[1](a_0,-0.3,0.5)
712     r1=oper[1](a1_0,-0.3,0.5)
713     text+=mkText(case0,"arg",a_0,a1_0)
714     text+=" res=%s(arg,-0.3,0.5)\n"%oper[0]
715     if case0=="Symbol":
716     text+=mkText("array","s",a_0,a1_0)
717     text+=" sub=res.substitute({arg:s})\n"
718     res="sub"
719     text+=mkText("array","ref",r,r1)
720     else:
721     res="res"
722     text+=mkText(case0,"ref",r,r1)
723     text+=mkTypeAndShapeTest(case0,sh0,"res")
724     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
725    
726     if case0 == "taggedData" :
727     t_prog_with_tags+=text
728     else:
729     t_prog+=text
730    
731     print test_header
732     # print t_prog
733     print t_prog_with_tags
734     print test_tail
735     1/0
736     #=======================================================================================================
737     # maximum, minimum, clipping
738     #=======================================================================================================
739     oper_L=[ ["maximum",maximumTEST],
740     ["minimum",minimumTEST]]
741     for oper in oper_L:
742     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
743     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
744     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
745     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
746     if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
747     and (sh0==sh1 or len(sh0)==0 or len(sh1)==0) :
748     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
749    
750     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
751     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
752     text+=" def %s(self):\n"%tname
753     a_0=makeArray(sh0,[-1.,1])
754     if case0 in ["taggedData", "expandedData"]:
755     a1_0=makeArray(sh0,[-1.,1])
756     else:
757     a1_0=a_0
758    
759     a_1=makeArray(sh1,[-1.,1])
760     if case1 in ["taggedData", "expandedData"]:
761     a1_1=makeArray(sh1,[-1.,1])
762     else:
763     a1_1=a_1
764     r=oper[1](a_0,a_1)
765     r1=oper[1](a1_0,a1_1)
766     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
767     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
768     text+=" res=%s(arg0,arg1)\n"%oper[0]
769     case=getResultCaseForBin(case0,case1)
770     if case=="Symbol":
771     c0_res,c1_res=case0,case1
772     subs="{"
773     if case0=="Symbol":
774     text+=mkText("array","s0",a_0,a1_0)
775     subs+="arg0:s0"
776     c0_res="array"
777     if case1=="Symbol":
778     text+=mkText("array","s1",a_1,a1_1)
779     if not subs.endswith("{"): subs+=","
780     subs+="arg1:s1"
781     c1_res="array"
782     subs+="}"
783     text+=" sub=res.substitute(%s)\n"%subs
784     res="sub"
785     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
786     else:
787     res="res"
788     text+=mkText(case,"ref",r,r1)
789     if len(sh0)>len(sh1):
790     text+=mkTypeAndShapeTest(case,sh0,"res")
791     else:
792     text+=mkTypeAndShapeTest(case,sh1,"res")
793     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
794    
795     if case0 == "taggedData" or case1 == "taggedData":
796     t_prog_with_tags+=text
797     else:
798     t_prog+=text
799    
800     print test_header
801     # print t_prog
802     print t_prog_with_tags
803     print test_tail
804     1/0
805    
806    
807     #=======================================================================================================
808     # outer inner
809     #=======================================================================================================
810     oper=["outer",outerTEST]
811     # oper=["inner",innerTEST]
812     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
813     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
814     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
815     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
816     if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
817     and len(sh0+sh1)<5:
818     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
819    
820     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
821     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
822     text+=" def %s(self):\n"%tname
823     a_0=makeArray(sh0,[-1.,1])
824     if case0 in ["taggedData", "expandedData"]:
825     a1_0=makeArray(sh0,[-1.,1])
826     else:
827     a1_0=a_0
828    
829     a_1=makeArray(sh1,[-1.,1])
830     if case1 in ["taggedData", "expandedData"]:
831     a1_1=makeArray(sh1,[-1.,1])
832     else:
833     a1_1=a_1
834     r=oper[1](a_0,a_1)
835     r1=oper[1](a1_0,a1_1)
836     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
837     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
838     text+=" res=%s(arg0,arg1)\n"%oper[0]
839     case=getResultCaseForBin(case0,case1)
840     if case=="Symbol":
841     c0_res,c1_res=case0,case1
842     subs="{"
843     if case0=="Symbol":
844     text+=mkText("array","s0",a_0,a1_0)
845     subs+="arg0:s0"
846     c0_res="array"
847     if case1=="Symbol":
848     text+=mkText("array","s1",a_1,a1_1)
849     if not subs.endswith("{"): subs+=","
850     subs+="arg1:s1"
851     c1_res="array"
852     subs+="}"
853     text+=" sub=res.substitute(%s)\n"%subs
854     res="sub"
855     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
856     else:
857     res="res"
858     text+=mkText(case,"ref",r,r1)
859     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
860     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
861    
862     if case0 == "taggedData" or case1 == "taggedData":
863     t_prog_with_tags+=text
864     else:
865     t_prog+=text
866    
867     print test_header
868     # print t_prog
869     print t_prog_with_tags
870     print test_tail
871     1/0
872    
873     #=======================================================================================================
874 gross 313 # local reduction
875     #=======================================================================================================
876     for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],
877     ["maxval",-1.e99,"max(out,%a1%)","out"],
878     ["minval",1.e99,"min(out,%a1%)","out"] ]:
879     for case in case_set:
880     for sh in shape_set:
881     if not case=="float" or len(sh)==0:
882     text=""
883     text+=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
884     tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))
885     text+=" %s(self):\n"%tname
886     a=makeArray(sh,[-1.,1.])
887     a1=makeArray(sh,[-1.,1.])
888     r1=testReduce(a1,oper[1],oper[2],oper[3])
889     r=testReduce(a,oper[1],oper[2],oper[3])
890    
891     text+=mkText(case,"arg",a,a1)
892     text+=" res=%s(arg)\n"%oper[0]
893     if case=="Symbol":
894     text+=mkText("array","s",a,a1)
895     text+=" sub=res.substitute({arg:s})\n"
896     text+=mkText("array","ref",r,r1)
897     res="sub"
898     else:
899     text+=mkText(case,"ref",r,r1)
900     res="res"
901     if oper[0]=="length":
902     text+=mkTypeAndShapeTest(case,(),"res")
903     else:
904     if case=="float" or case=="array":
905     text+=mkTypeAndShapeTest("float",(),"res")
906     else:
907     text+=mkTypeAndShapeTest(case,(),"res")
908     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
909     if case == "taggedData":
910     t_prog_with_tags+=text
911     else:
912     t_prog+=text
913     print test_header
914     # print t_prog
915     print t_prog_with_tags
916     print test_tail
917     1/0
918    
919     #=======================================================================================================
920 gross 291 # tensor multiply
921 gross 157 #=======================================================================================================
922 gross 291 # oper=["generalTensorProduct",tensorProductTest]
923     # oper=["matrixmult",testMatrixMult]
924     oper=["tensormult",testTensorMult]
925    
926     for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
927     for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
928     for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
929     for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
930     for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
931     if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
932     and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
933     # 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
934     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
935     case=getResultCaseForBin(case0,case1)
936 gross 157 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
937     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
938 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))
939     #tname="test_matrixmult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
940     tname="test_tensormult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
941     # if tname=="test_generalTensorProduct_array_rank1_array_rank2_offset1":
942     # print tnametest_generalTensorProduct_Symbol_rank1_Symbol_rank3_offset1
943     text+=" def %s(self):\n"%tname
944     a_0=makeArray(sh0+sh_s,[-1.,1])
945     if case0 in ["taggedData", "expandedData"]:
946     a1_0=makeArray(sh0+sh_s,[-1.,1])
947     else:
948     a1_0=a_0
949    
950     a_1=makeArray(sh_s+sh1,[-1.,1])
951     if case1 in ["taggedData", "expandedData"]:
952     a1_1=makeArray(sh_s+sh1,[-1.,1])
953     else:
954     a1_1=a_1
955     r=oper[1](a_0,a_1,sh_s)
956     r1=oper[1](a1_0,a1_1,sh_s)
957     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
958     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
959     #text+=" res=matrixmult(arg0,arg1)\n"
960     text+=" res=tensormult(arg0,arg1)\n"
961     #text+=" res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))
962     if case=="Symbol":
963     c0_res,c1_res=case0,case1
964     subs="{"
965     if case0=="Symbol":
966     text+=mkText("array","s0",a_0,a1_0)
967     subs+="arg0:s0"
968     c0_res="array"
969     if case1=="Symbol":
970     text+=mkText("array","s1",a_1,a1_1)
971     if not subs.endswith("{"): subs+=","
972     subs+="arg1:s1"
973     c1_res="array"
974     subs+="}"
975     text+=" sub=res.substitute(%s)\n"%subs
976     res="sub"
977     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
978     else:
979     res="res"
980     text+=mkText(case,"ref",r,r1)
981     text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
982     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
983     if case0 == "taggedData" or case1 == "taggedData":
984     t_prog_with_tags+=text
985     else:
986     t_prog+=text
987     print test_header
988     # print t_prog
989     print t_prog_with_tags
990     print test_tail
991     1/0
992     #=======================================================================================================
993 gross 157 # basic binary operation overloading (tests only!)
994     #=======================================================================================================
995     oper_range=[-5.,5.]
996     for oper in [["add" ,"+",[-5.,5.]],
997     ["sub" ,"-",[-5.,5.]],
998     ["mult","*",[-5.,5.]],
999     ["div" ,"/",[-5.,5.]],
1000     ["pow" ,"**",[0.01,5.]]]:
1001     for case0 in case_set:
1002     for sh0 in shape_set:
1003     for case1 in case_set:
1004     for sh1 in shape_set:
1005 gross 291 if not case0=="array" and \
1006     (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1007 gross 157 (sh0==() or sh1==() or sh1==sh0) and \
1008     not (case0 in ["float","array"] and case1 in ["float","array"]):
1009 gross 291 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1010 gross 157 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1011     tname="test_%s_overloaded_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1012     text+=" def %s(self):\n"%tname
1013     a_0=makeArray(sh0,oper[2])
1014     if case0 in ["taggedData", "expandedData"]:
1015     a1_0=makeArray(sh0,oper[2])
1016     else:
1017     a1_0=a_0
1018 jgs 154
1019 gross 157 a_1=makeArray(sh1,oper[2])
1020     if case1 in ["taggedData", "expandedData"]:
1021     a1_1=makeArray(sh1,oper[2])
1022 jgs 154 else:
1023 gross 157 a1_1=a_1
1024     r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1025     r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1026     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1027     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1028     text+=" res=arg0%sarg1\n"%oper[1]
1029    
1030     case=getResultCaseForBin(case0,case1)
1031     if case=="Symbol":
1032     c0_res,c1_res=case0,case1
1033     subs="{"
1034     if case0=="Symbol":
1035     text+=mkText("array","s0",a_0,a1_0)
1036     subs+="arg0:s0"
1037     c0_res="array"
1038     if case1=="Symbol":
1039     text+=mkText("array","s1",a_1,a1_1)
1040     if not subs.endswith("{"): subs+=","
1041     subs+="arg1:s1"
1042     c1_res="array"
1043     subs+="}"
1044     text+=" sub=res.substitute(%s)\n"%subs
1045     res="sub"
1046     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1047 jgs 154 else:
1048 gross 157 res="res"
1049     text+=mkText(case,"ref",r,r1)
1050     if isinstance(r,float):
1051     text+=mkTypeAndShapeTest(case,(),"res")
1052     else:
1053     text+=mkTypeAndShapeTest(case,r.shape,"res")
1054     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1055    
1056     if case0 in [ "constData","taggedData","expandedData"] and case1 == "Symbol":
1057     t_prog_failing+=text
1058     else:
1059     if case0 == "taggedData" or case1 == "taggedData":
1060     t_prog_with_tags+=text
1061     else:
1062     t_prog+=text
1063 jgs 154
1064 gross 157
1065     print test_header
1066 gross 291 # print t_prog
1067     # print t_prog_with_tags
1068     print t_prog_failing
1069     print test_tail
1070 jgs 154 1/0
1071 gross 291 #=======================================================================================================
1072     # basic binary operations (tests only!)
1073     #=======================================================================================================
1074     oper_range=[-5.,5.]
1075     for oper in [["add" ,"+",[-5.,5.]],
1076     ["mult","*",[-5.,5.]],
1077     ["quotient" ,"/",[-5.,5.]],
1078     ["power" ,"**",[0.01,5.]]]:
1079     for case0 in case_set:
1080     for case1 in case_set:
1081     for sh in shape_set:
1082     for sh_p in shape_set:
1083     if len(sh_p)>0:
1084     resource=[-1,1]
1085     else:
1086     resource=[1]
1087     for sh_d in resource:
1088     if sh_d>0:
1089     sh0=sh
1090     sh1=sh+sh_p
1091     else:
1092     sh1=sh
1093     sh0=sh+sh_p
1094    
1095     if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1096     len(sh0)<5 and len(sh1)<5:
1097     use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1098     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1099     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1100     text+=" def %s(self):\n"%tname
1101     a_0=makeArray(sh0,oper[2])
1102     if case0 in ["taggedData", "expandedData"]:
1103     a1_0=makeArray(sh0,oper[2])
1104     else:
1105     a1_0=a_0
1106    
1107     a_1=makeArray(sh1,oper[2])
1108     if case1 in ["taggedData", "expandedData"]:
1109     a1_1=makeArray(sh1,oper[2])
1110     else:
1111     a1_1=a_1
1112     r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1113     r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1114     text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1115     text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1116     text+=" res=%s(arg0,arg1)\n"%oper[0]
1117    
1118     case=getResultCaseForBin(case0,case1)
1119     if case=="Symbol":
1120     c0_res,c1_res=case0,case1
1121     subs="{"
1122     if case0=="Symbol":
1123     text+=mkText("array","s0",a_0,a1_0)
1124     subs+="arg0:s0"
1125     c0_res="array"
1126     if case1=="Symbol":
1127     text+=mkText("array","s1",a_1,a1_1)
1128     if not subs.endswith("{"): subs+=","
1129     subs+="arg1:s1"
1130     c1_res="array"
1131     subs+="}"
1132     text+=" sub=res.substitute(%s)\n"%subs
1133     res="sub"
1134     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1135     else:
1136     res="res"
1137     text+=mkText(case,"ref",r,r1)
1138     if isinstance(r,float):
1139     text+=mkTypeAndShapeTest(case,(),"res")
1140     else:
1141     text+=mkTypeAndShapeTest(case,r.shape,"res")
1142     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1143    
1144     if case0 == "taggedData" or case1 == "taggedData":
1145     t_prog_with_tags+=text
1146     else:
1147     t_prog+=text
1148     print test_header
1149     # print t_prog
1150     print t_prog_with_tags
1151     print test_tail
1152     1/0
1153    
1154 gross 157 # print t_prog_with_tagsoper_range=[-5.,5.]
1155     for oper in [["add" ,"+",[-5.,5.]],
1156     ["sub" ,"-",[-5.,5.]],
1157     ["mult","*",[-5.,5.]],
1158     ["div" ,"/",[-5.,5.]],
1159     ["pow" ,"**",[0.01,5.]]]:
1160     for case0 in case_set:
1161     for sh0 in shape_set:
1162     for case1 in case_set:
1163     for sh1 in shape_set:
1164     if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1165     (sh0==() or sh1==() or sh1==sh0) and \
1166     not (case0 in ["float","array"] and case1 in ["float","array"]):
1167     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1168     tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1169     text+=" def %s(self):\n"%tname
1170     a_0=makeArray(sh0,oper[2])
1171     if case0 in ["taggedData", "expandedData"]:
1172     a1_0=makeArray(sh0,oper[2])
1173     else:
1174     a1_0=a_0
1175 jgs 154
1176 gross 157 a_1=makeArray(sh1,oper[2])
1177     if case1 in ["taggedData", "expandedData"]:
1178     a1_1=makeArray(sh1,oper[2])
1179 jgs 154 else:
1180 gross 157 a1_1=a_1
1181     r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1182     r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1183     text+=mkText(case0,"arg0",a_0,a1_0)
1184     text+=mkText(case1,"arg1",a_1,a1_1)
1185     text+=" res=arg0%sarg1\n"%oper[1]
1186    
1187     case=getResultCaseForBin(case0,case1)
1188     if case=="Symbol":
1189     c0_res,c1_res=case0,case1
1190     subs="{"
1191     if case0=="Symbol":
1192     text+=mkText("array","s0",a_0,a1_0)
1193     subs+="arg0:s0"
1194     c0_res="array"
1195     if case1=="Symbol":
1196     text+=mkText("array","s1",a_1,a1_1)
1197     if not subs.endswith("{"): subs+=","
1198     subs+="arg1:s1"
1199     c1_res="array"
1200     subs+="}"
1201     text+=" sub=res.substitute(%s)\n"%subs
1202     res="sub"
1203     text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1204 jgs 154 else:
1205 gross 157 res="res"
1206     text+=mkText(case,"ref",r,r1)
1207     if isinstance(r,float):
1208     text+=mkTypeAndShapeTest(case,(),"res")
1209     else:
1210     text+=mkTypeAndShapeTest(case,r.shape,"res")
1211     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1212    
1213     if case0 in [ "constData","taggedData","expandedData"] and case1 == "Symbol":
1214     t_prog_failing+=text
1215     else:
1216     if case0 == "taggedData" or case1 == "taggedData":
1217     t_prog_with_tags+=text
1218     else:
1219     t_prog+=text
1220 jgs 154
1221 gross 157
1222     # print u_prog
1223     # 1/0
1224     print test_header
1225 jgs 154 print t_prog
1226 gross 157 # print t_prog_with_tags
1227     # print t_prog_failing
1228     print test_tail
1229     # print t_prog_failing
1230     print test_tail
1231 jgs 154
1232 gross 157 #=======================================================================================================
1233     # unary operations:
1234     #=======================================================================================================
1235     func= [
1236     OPERATOR(nickname="log10",\
1237     rng=[1.e-3,100.],\
1238     test_expr="math.log10(%a1%)",\
1239     math_expr="math.log10(%a1%)",\
1240     numarray_expr="numarray.log10(%a1%)",\
1241     symbol_expr="log(%a1%)/log(10.)",\
1242     name="base-10 logarithm"),
1243     OPERATOR(nickname="wherePositive",\
1244     rng=[-100.,100.],\
1245     test_expr="wherepos(%a1%)",\
1246     math_expr="if arg>0:\n return 1.\nelse:\n return 0.",
1247     numarray_expr="numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))",\
1248     name="mask of positive values"),
1249     OPERATOR(nickname="whereNegative",\
1250     rng=[-100.,100.],\
1251     test_expr="wherepos(-%a1%)",\
1252     math_expr="if arg<0:\n return 1.\nelse:\n return 0.",
1253     numarray_expr="numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))",\
1254     name="mask of positive values"),
1255     OPERATOR(nickname="whereNonNegative",\
1256     rng=[-100.,100.],\
1257     test_expr="1-wherepos(-%a1%)", \
1258     math_expr="if arg<0:\n return 0.\nelse:\n return 1.",
1259     numarray_expr="numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))",\
1260     symbol_expr="1-wherePositive(%a1%)",\
1261     name="mask of non-negative values"),
1262     OPERATOR(nickname="whereNonPositive",\
1263     rng=[-100.,100.],\
1264     test_expr="1-wherepos(%a1%)",\
1265     math_expr="if arg>0:\n return 0.\nelse:\n return 1.",
1266     numarray_expr="numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))",\
1267     symbol_expr="1-whereNegative(%a1%)",\
1268     name="mask of non-positive values"),
1269     OPERATOR(nickname="whereZero",\
1270     rng=[-100.,100.],\
1271     test_expr="1-wherepos(%a1%)-wherepos(-%a1%)",\
1272     math_expr="if abs(%a1%)<=tol:\n return 1.\nelse:\n return 0.",
1273     numarray_expr="numarray.less_equal(abs(%a1%)-tol,numarray.zeros(arg.shape,numarray.Float))",\
1274     name="mask of zero entries"),
1275     OPERATOR(nickname="whereNonZero",\
1276     rng=[-100.,100.],\
1277     test_expr="wherepos(%a1%)+wherepos(-%a1%)",\
1278     math_expr="if abs(%a1%)>tol:\n return 1.\nelse:\n return 0.",\
1279     numarray_expr="numarray.greater(abs(%a1%)-tol,numarray.zeros(arg.shape,numarray.Float))",\
1280     symbol_expr="1-whereZero(arg,tol)",\
1281     name="mask of values different from zero"),
1282     OPERATOR(nickname="sin",\
1283     rng=[-100.,100.],\
1284     test_expr="math.sin(%a1%)",
1285     numarray_expr="numarray.sin(%a1%)",\
1286     diff="cos(%a1%)",\
1287     name="sine"),
1288     OPERATOR(nickname="cos",\
1289     rng=[-100.,100.],\
1290     test_expr="math.cos(%a1%)",
1291     numarray_expr="numarray.cos(%a1%)",\
1292     diff="-sin(%a1%)",
1293     name="cosine"),
1294     OPERATOR(nickname="tan",\
1295     rng=[-100.,100.],\
1296     test_expr="math.tan(%a1%)",
1297     numarray_expr="numarray.tan(%a1%)",\
1298     diff="1./cos(%a1%)**2",
1299     name="tangent"),
1300     OPERATOR(nickname="asin",\
1301     rng=[-0.99,0.99],\
1302     test_expr="math.asin(%a1%)",
1303     numarray_expr="numarray.arcsin(%a1%)",
1304     diff="1./sqrt(1.-%a1%**2)",
1305     name="inverse sine"),
1306     OPERATOR(nickname="acos",\
1307     rng=[-0.99,0.99],\
1308     test_expr="math.acos(%a1%)",
1309     numarray_expr="numarray.arccos(%a1%)",
1310     diff="-1./sqrt(1.-%a1%**2)",
1311     name="inverse cosine"),
1312     OPERATOR(nickname="atan",\
1313     rng=[-100.,100.],\
1314     test_expr="math.atan(%a1%)",
1315     numarray_expr="numarray.arctan(%a1%)",
1316     diff="1./(1+%a1%**2)",
1317     name="inverse tangent"),
1318     OPERATOR(nickname="sinh",\
1319     rng=[-5,5],\
1320     test_expr="math.sinh(%a1%)",
1321     numarray_expr="numarray.sinh(%a1%)",
1322     diff="cosh(%a1%)",
1323     name="hyperbolic sine"),
1324     OPERATOR(nickname="cosh",\
1325     rng=[-5.,5.],
1326     test_expr="math.cosh(%a1%)",
1327     numarray_expr="numarray.cosh(%a1%)",
1328     diff="sinh(%a1%)",
1329     name="hyperbolic cosine"),
1330     OPERATOR(nickname="tanh",\
1331     rng=[-5.,5.],
1332     test_expr="math.tanh(%a1%)",
1333     numarray_expr="numarray.tanh(%a1%)",
1334     diff="1./cosh(%a1%)**2",
1335     name="hyperbolic tangent"),
1336     OPERATOR(nickname="asinh",\
1337     rng=[-100.,100.], \
1338     test_expr="numarray.arcsinh(%a1%)",
1339     math_expr="numarray.arcsinh(%a1%)",
1340     numarray_expr="numarray.arcsinh(%a1%)",
1341     diff="1./sqrt(%a1%**2+1)",
1342     name="inverse hyperbolic sine"),
1343     OPERATOR(nickname="acosh",\
1344     rng=[1.001,100.],\
1345     test_expr="numarray.arccosh(%a1%)",
1346     math_expr="numarray.arccosh(%a1%)",
1347     numarray_expr="numarray.arccosh(%a1%)",
1348     diff="1./sqrt(%a1%**2-1)",
1349     name="inverse hyperolic cosine"),
1350     OPERATOR(nickname="atanh",\
1351     rng=[-0.99,0.99], \
1352     test_expr="numarray.arctanh(%a1%)",
1353     math_expr="numarray.arctanh(%a1%)",
1354     numarray_expr="numarray.arctanh(%a1%)",
1355     diff="1./(1.-%a1%**2)",
1356     name="inverse hyperbolic tangent"),
1357     OPERATOR(nickname="exp",\
1358     rng=[-5.,5.],
1359     test_expr="math.exp(%a1%)",
1360     numarray_expr="numarray.exp(%a1%)",
1361     diff="self",
1362     name="exponential"),
1363     OPERATOR(nickname="sqrt",\
1364     rng=[1.e-3,100.],\
1365     test_expr="math.sqrt(%a1%)",
1366     numarray_expr="numarray.sqrt(%a1%)",
1367     diff="0.5/self",
1368     name="square root"),
1369     OPERATOR(nickname="log", \
1370     rng=[1.e-3,100.],\
1371     test_expr="math.log(%a1%)",
1372     numarray_expr="numarray.log(%a1%)",
1373     diff="1./arg",
1374     name="natural logarithm"),
1375     OPERATOR(nickname="sign",\
1376     rng=[-100.,100.], \
1377     math_expr="if %a1%>0:\n return 1.\nelif %a1%<0:\n return -1.\nelse:\n return 0.",
1378     test_expr="wherepos(%a1%)-wherepos(-%a1%)",
1379     numarray_expr="numarray.sign(%a1%)",
1380     symbol_expr="wherePositive(%a1%)-whereNegative(%a1%)",\
1381     name="sign"),
1382     OPERATOR(nickname="abs",\
1383     rng=[-100.,100.], \
1384     math_expr="if %a1%>0:\n return %a1% \nelif %a1%<0:\n return -(%a1%)\nelse:\n return 0.",
1385     test_expr="wherepos(%a1%)*(%a1%)-wherepos(-%a1%)*(%a1%)",
1386     numarray_expr="abs(%a1%)",
1387     diff="sign(%a1%)",
1388     name="absolute value")
1389    
1390     ]
1391     for f in func:
1392     symbol_name=f.nickname[0].upper()+f.nickname[1:]
1393     if f.nickname!="abs":
1394     u_prog+="def %s(arg):\n"%f.nickname
1395     u_prog+=" \"\"\"\n"
1396     u_prog+=" returns %s of argument arg\n\n"%f.name
1397     u_prog+=" @param arg: argument\n"
1398     u_prog+=" @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.\n"
1399     u_prog+=" @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.\n"
1400     u_prog+=" @raises TypeError: if the type of the argument is not expected.\n"
1401     u_prog+=" \"\"\"\n"
1402     u_prog+=" if isinstance(arg,numarray.NumArray):\n"
1403     u_prog+=mkCode(f.numarray_expr,["arg"],2*" ")
1404     u_prog+=" elif isinstance(arg,escript.Data):\n"
1405     u_prog+=mkCode("arg._%s()"%f.nickname,[],2*" ")
1406     u_prog+=" elif isinstance(arg,float):\n"
1407     u_prog+=mkCode(f.math_expr,["arg"],2*" ")
1408     u_prog+=" elif isinstance(arg,int):\n"
1409     u_prog+=mkCode(f.math_expr,["float(arg)"],2*" ")
1410     u_prog+=" elif isinstance(arg,Symbol):\n"
1411     if f.symbol_expr==None:
1412     u_prog+=mkCode("%s_Symbol(arg)"%symbol_name,[],2*" ")
1413     else:
1414     u_prog+=mkCode(f.symbol_expr,["arg"],2*" ")
1415     u_prog+=" else:\n"
1416     u_prog+=" raise TypeError,\"%s: Unknown argument type.\"\n\n"%f.nickname
1417     if f.symbol_expr==None:
1418     u_prog+="class %s_Symbol(DependendSymbol):\n"%symbol_name
1419     u_prog+=" \"\"\"\n"
1420     u_prog+=" L{Symbol} representing the result of the %s function\n"%f.name
1421     u_prog+=" \"\"\"\n"
1422     u_prog+=" def __init__(self,arg):\n"
1423     u_prog+=" \"\"\"\n"
1424     u_prog+=" initialization of %s L{Symbol} with argument arg\n"%f.nickname
1425     u_prog+=" @param arg: argument of function\n"
1426     u_prog+=" @type arg: typically L{Symbol}.\n"
1427     u_prog+=" \"\"\"\n"
1428     u_prog+=" DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())\n"
1429     u_prog+="\n"
1430    
1431     u_prog+=" def getMyCode(self,argstrs,format=\"escript\"):\n"
1432     u_prog+=" \"\"\"\n"
1433     u_prog+=" returns a program code that can be used to evaluate the symbol.\n\n"
1434 jgs 154
1435 gross 157 u_prog+=" @param argstrs: gives for each argument a string representing the argument for the evaluation.\n"
1436     u_prog+=" @type argstrs: C{str} or a C{list} of length 1 of C{str}.\n"
1437     u_prog+=" @param format: specifies the format to be used. At the moment only \"escript\" ,\"text\" and \"str\" are supported.\n"
1438     u_prog+=" @type format: C{str}\n"
1439     u_prog+=" @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.\n"
1440     u_prog+=" @rtype: C{str}\n"
1441     u_prog+=" @raise: NotImplementedError: if the requested format is not available\n"
1442     u_prog+=" \"\"\"\n"
1443     u_prog+=" if isinstance(argstrs,list):\n"
1444     u_prog+=" argstrs=argstrs[0]\n"
1445     u_prog+=" if format==\"escript\" or format==\"str\" or format==\"text\":\n"
1446     u_prog+=" return \"%s(%%s)\"%%argstrs\n"%f.nickname
1447     u_prog+=" else:\n"
1448     u_prog+=" raise NotImplementedError,\"%s_Symbol does not provide program code for format %%s.\"%%format\n"%symbol_name
1449     u_prog+="\n"
1450 jgs 154
1451 gross 157 u_prog+=" def substitute(self,argvals):\n"
1452     u_prog+=" \"\"\"\n"
1453     u_prog+=" assigns new values to symbols in the definition of the symbol.\n"
1454     u_prog+=" The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.\n"
1455     u_prog+="\n"
1456     u_prog+=" @param argvals: new values assigned to symbols\n"
1457     u_prog+=" @type argvals: C{dict} with keywords of type L{Symbol}.\n"
1458     u_prog+=" @return: result of the substitution process. Operations are executed as much as possible.\n"
1459     u_prog+=" @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution\n"
1460     u_prog+=" @raise TypeError: if a value for a L{Symbol} cannot be substituted.\n"
1461     u_prog+=" \"\"\"\n"
1462     u_prog+=" if argvals.has_key(self):\n"
1463     u_prog+=" arg=argvals[self]\n"
1464     u_prog+=" if self.isAppropriateValue(arg):\n"
1465     u_prog+=" return arg\n"
1466     u_prog+=" else:\n"
1467     u_prog+=" raise TypeError,\"%s: new value is not appropriate.\"%str(self)\n"
1468     u_prog+=" else:\n"
1469     u_prog+=" arg=self.getSubstitutedArguments(argvals)[0]\n"
1470     u_prog+=" return %s(arg)\n\n"%f.nickname
1471     if not f.diff==None:
1472     u_prog+=" def diff(self,arg):\n"
1473     u_prog+=" \"\"\"\n"
1474     u_prog+=" differential of this object\n"
1475     u_prog+="\n"
1476     u_prog+=" @param arg: the derivative is calculated with respect to arg\n"
1477     u_prog+=" @type arg: L{escript.Symbol}\n"
1478     u_prog+=" @return: derivative with respect to C{arg}\n"
1479     u_prog+=" @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray} are possible.\n"
1480     u_prog+=" \"\"\"\n"
1481     u_prog+=" if arg==self:\n"
1482     u_prog+=" return identity(self.getShape())\n"
1483     u_prog+=" else:\n"
1484     u_prog+=" myarg=self.getArgument()[0]\n"
1485     u_prog+=" val=matchShape(%s,self.getDifferentiatedArguments(arg)[0])\n"%f.diff.replace("%a1%","myarg")
1486     u_prog+=" return val[0]*val[1]\n\n"
1487 jgs 154
1488 gross 157 for case in case_set:
1489     for sh in shape_set:
1490     if not case=="float" or len(sh)==0:
1491     text=""
1492     text+=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1493     tname="def test_%s_%s_rank%s"%(f.nickname,case,len(sh))
1494     text+=" %s(self):\n"%tname
1495     a=makeArray(sh,f.rng)
1496     a1=makeArray(sh,f.rng)
1497     r1=makeResult(a1,f.test_expr)
1498     r=makeResult(a,f.test_expr)
1499    
1500     text+=mkText(case,"arg",a,a1)
1501     text+=" res=%s(arg)\n"%f.nickname
1502     if case=="Symbol":
1503     text+=mkText("array","s",a,a1)
1504     text+=" sub=res.substitute({arg:s})\n"
1505     text+=mkText("array","ref",r,r1)
1506     res="sub"
1507     else:
1508     text+=mkText(case,"ref",r,r1)
1509     res="res"
1510     text+=mkTypeAndShapeTest(case,sh,"res")
1511     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1512     if case == "taggedData":
1513     t_prog_with_tags+=text
1514     else:
1515     t_prog+=text
1516    
1517     #=========== END OF GOOD CODE +++++++++++++++++++++++++++
1518 jgs 154
1519 gross 157 1/0
1520 jgs 154
1521 gross 157 def X():
1522     if args=="float":
1523     a=makeArray(sh,f[RANGE])
1524 jgs 154 r=makeResult(a,f)
1525 gross 157 t_prog+=" arg=%s\n"%a[0]
1526     t_prog+=" ref=%s\n"%r[0]
1527     t_prog+=" res=%s(%a1%)\n"%f.nickname
1528     t_prog+=" self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1529     t_prog+=" self.failUnless(Lsup(res-ref)<=self.tol*Lsup(ref),\"wrong result\")\n"
1530     elif args == "array":
1531     a=makeArray(sh,f[RANGE])
1532     r=makeResult(a,f)
1533 jgs 154 if len(sh)==0:
1534 gross 157 t_prog+=" arg=numarray.array(%s)\n"%a[0]
1535     t_prog+=" ref=numarray.array(%s)\n"%r[0]
1536 jgs 154 else:
1537     t_prog+=" arg=numarray.array(%s)\n"%a.tolist()
1538     t_prog+=" ref=numarray.array(%s)\n"%r.tolist()
1539 gross 157 t_prog+=" res=%s(%a1%)\n"%f.nickname
1540     t_prog+=" self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
1541     t_prog+=" self.failUnless(Lsup(res-ref)<=self.tol*Lsup(ref),\"wrong result\")\n"
1542 jgs 154 elif args== "constData":
1543 gross 157 a=makeArray(sh,f[RANGE])
1544 jgs 154 r=makeResult(a,f)
1545     if len(sh)==0:
1546     t_prog+=" arg=Data(%s,self.functionspace)\n"%(a)
1547     t_prog+=" ref=%s\n"%r
1548     else:
1549     t_prog+=" arg=Data(numarray.array(%s),self.functionspace)\n"%(a.tolist())
1550     t_prog+=" ref=numarray.array(%s)\n"%r.tolist()
1551 gross 157 t_prog+=" res=%s(%a1%)\n"%f.nickname
1552 jgs 154 t_prog+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1553     t_prog+=" self.failUnless(Lsup(res-ref)<=self.tol*Lsup(ref),\"wrong result\")\n"
1554     elif args in [ "taggedData","expandedData"]:
1555 gross 157 a=makeArray(sh,f[RANGE])
1556 jgs 154 r=makeResult(a,f)
1557 gross 157 a1=makeArray(sh,f[RANGE])
1558 jgs 154 r1=makeResult(a1,f)
1559     if len(sh)==0:
1560     if args=="expandedData":
1561     t_prog+=" arg=Data(%s,self.functionspace,True)\n"%(a)
1562     t_prog+=" ref=Data(%s,self.functionspace,True)\n"%(r)
1563     else:
1564     t_prog+=" arg=Data(%s,self.functionspace)\n"%(a)
1565     t_prog+=" ref=Data(%s,self.functionspace)\n"%(r)
1566     t_prog+=" arg.setTaggedValue(1,%s)\n"%a
1567     t_prog+=" ref.setTaggedValue(1,%s)\n"%r1
1568     else:
1569     if args=="expandedData":
1570     t_prog+=" arg=Data(numarray.array(%s),self.functionspace,True)\n"%(a.tolist())
1571     t_prog+=" ref=Data(numarray.array(%s),self.functionspace,True)\n"%(r.tolist())
1572     else:
1573     t_prog+=" arg=Data(numarray.array(%s),self.functionspace)\n"%(a.tolist())
1574     t_prog+=" ref=Data(numarray.array(%s),self.functionspace)\n"%(r.tolist())
1575     t_prog+=" arg.setTaggedValue(1,%s)\n"%a1.tolist()
1576     t_prog+=" ref.setTaggedValue(1,%s)\n"%r1.tolist()
1577 gross 157 t_prog+=" res=%s(%a1%)\n"%f.nickname
1578 jgs 154 t_prog+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1579     t_prog+=" self.failUnless(Lsup(res-ref)<=self.tol*Lsup(ref),\"wrong result\")\n"
1580     elif args=="Symbol":
1581     t_prog+=" arg=Symbol(shape=%s)\n"%str(sh)
1582 gross 157 t_prog+=" v=%s(%a1%)\n"%f.nickname
1583 jgs 154 t_prog+=" self.failUnlessRaises(ValueError,v.substitute,Symbol(shape=(1,1)),\"illegal shape of substitute not identified.\")\n"
1584 gross 157 a=makeArray(sh,f[RANGE])
1585 jgs 154 r=makeResult(a,f)
1586     if len(sh)==0:
1587     t_prog+=" res=v.substitute({arg : %s})\n"%a
1588     t_prog+=" ref=%s\n"%r
1589     t_prog+=" self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1590     else:
1591     t_prog+=" res=v.substitute({arg : numarray.array(%s)})\n"%a.tolist()
1592     t_prog+=" ref=numarray.array(%s)\n"%r.tolist()
1593     t_prog+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of substitution result.\")\n"%str(sh)
1594     t_prog+=" self.failUnless(Lsup(res-ref)<=self.tol*Lsup(ref),\"wrong result\")\n"
1595    
1596     if len(sh)==0:
1597     t_prog+=" # test derivative with respect to itself:\n"
1598     t_prog+=" dvdv=v.diff(v)\n"
1599     t_prog+=" self.failUnlessEqual(dvdv,1.,\"derivative with respect to self is not 1.\")\n"
1600     elif len(sh)==1:
1601     t_prog+=" # test derivative with respect to itself:\n"
1602     t_prog+=" dvdv=v.diff(v)\n"
1603     t_prog+=" self.failUnlessEqual(dvdv.shape,%s,\"shape of derivative with respect is wrong\")\n"%str(sh+sh)
1604     for i0_l in range(sh[0]):
1605     for i0_r in range(sh[0]):
1606     if i0_l == i0_r:
1607     v=1.
1608     else:
1609     v=0.
1610     t_prog+=" self.failUnlessEqual(dvdv[%s,%s],%s,\"derivative with respect to self: [%s,%s] is not %s\")\n"%(i0_l,i0_r,v,i0_l,i0_r,v)
1611     elif len(sh)==2:
1612     t_prog+=" # test derivative with respect to itself:\n"
1613     t_prog+=" dvdv=v.diff(v)\n"
1614     t_prog+=" self.failUnlessEqual(dvdv.shape,%s,\"shape of derivative with respect is wrong\")\n"%str(sh+sh)
1615     for i0_l in range(sh[0]):
1616     for i0_r in range(sh[0]):
1617     for i1_l in range(sh[1]):
1618     for i1_r in range(sh[1]):
1619     if i0_l == i0_r and i1_l == i1_r:
1620     v=1.
1621     else:
1622     v=0.
1623     t_prog+=" self.failUnlessEqual(dvdv[%s,%s,%s,%s],%s,\"derivative with respect to self: [%s,%s,%s,%s] is not %s\")\n"%(i0_l,i1_l,i0_r,i1_r,v,i0_l,i1_l,i0_r,i1_r,v)
1624    
1625     for sh_in in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1626     if len(sh_in)+len(sh)<=4:
1627    
1628     t_prog+=" # test derivative with shape %s as argument\n"%str(sh_in)
1629     trafo=makeArray(sh+sh_in,[0,1.])
1630 gross 157 a_in=makeArray(sh_in,f[RANGE])
1631 jgs 154 t_prog+=" arg_in=Symbol(shape=%s)\n"%str(sh_in)
1632     t_prog+=" arg2=Symbol(shape=%s)\n"%str(sh)
1633    
1634     if len(sh)==0:
1635     t_prog+=" arg2="
1636     if len(sh_in)==0:
1637     t_prog+="%s*arg_in\n"%trafo
1638     elif len(sh_in)==1:
1639     for i0 in range(sh_in[0]):
1640     if i0>0: t_prog+="+"
1641     t_prog+="%s*arg_in[%s]"%(trafo[i0],i0)
1642     t_prog+="\n"
1643     elif len(sh_in)==2:
1644     for i0 in range(sh_in[0]):
1645     for i1 in range(sh_in[1]):
1646     if i0+i1>0: t_prog+="+"
1647     t_prog+="%s*arg_in[%s,%s]"%(trafo[i0,i1],i0,i1)
1648    
1649     elif len(sh_in)==3:
1650     for i0 in range(sh_in[0]):
1651     for i1 in range(sh_in[1]):
1652     for i2 in range(sh_in[2]):
1653     if i0+i1+i2>0: t_prog+="+"
1654     t_prog+="%s*arg_in[%s,%s,%s]"%(trafo[i0,i1,i2],i0,i1,i2)
1655     elif len(sh_in)==4:
1656     for i0 in range(sh_in[0]):
1657     for i1 in range(sh_in[1]):
1658     for i2 in range(sh_in[2]):
1659     for i3 in range(sh_in[3]):
1660     if i0+i1+i2+i3>0: t_prog+="+"
1661     t_prog+="%s*arg_in[%s,%s,%s,%s]"%(trafo[i0,i1,i2,i3],i0,i1,i2,i3)
1662     t_prog+="\n"
1663     elif len(sh)==1:
1664     for j0 in range(sh[0]):
1665     t_prog+=" arg2[%s]="%j0
1666     if len(sh_in)==0:
1667     t_prog+="%s*arg_in"%trafo[j0]
1668     elif len(sh_in)==1:
1669     for i0 in range(sh_in[0]):
1670     if i0>0: t_prog+="+"
1671     t_prog+="%s*arg_in[%s]"%(trafo[j0,i0],i0)
1672     elif len(sh_in)==2:
1673     for i0 in range(sh_in[0]):
1674     for i1 in range(sh_in[1]):
1675     if i0+i1>0: t_prog+="+"
1676     t_prog+="%s*arg_in[%s,%s]"%(trafo[j0,i0,i1],i0,i1)
1677     elif len(sh_in)==3:
1678     for i0 in range(sh_in[0]):
1679     for i1 in range(sh_in[1]):
1680     for i2 in range(sh_in[2]):
1681     if i0+i1+i2>0: t_prog+="+"
1682     t_prog+="%s*arg_in[%s,%s,%s]"%(trafo[j0,i0,i1,i2],i0,i1,i2)
1683     t_prog+="\n"
1684     elif len(sh)==2:
1685     for j0 in range(sh[0]):
1686     for j1 in range(sh[1]):
1687     t_prog+=" arg2[%s,%s]="%(j0,j1)
1688     if len(sh_in)==0:
1689     t_prog+="%s*arg_in"%trafo[j0,j1]
1690     elif len(sh_in)==1:
1691     for i0 in range(sh_in[0]):
1692     if i0>0: t_prog+="+"
1693     t_prog+="%s*arg_in[%s]"%(trafo[j0,j1,i0],i0)
1694     elif len(sh_in)==2:
1695     for i0 in range(sh_in[0]):
1696     for i1 in range(sh_in[1]):
1697     if i0+i1>0: t_prog+="+"
1698     t_prog+="%s*arg_in[%s,%s]"%(trafo[j0,j1,i0,i1],i0,i1)
1699     t_prog+="\n"
1700     elif len(sh)==3:
1701     for j0 in range(sh[0]):
1702     for j1 in range(sh[1]):
1703     for j2 in range(sh[2]):
1704     t_prog+=" arg2[%s,%s,%s]="%(j0,j1,j2)
1705     if len(sh_in)==0:
1706     t_prog+="%s*arg_in"%trafo[j0,j1,j2]
1707     elif len(sh_in)==1:
1708     for i0 in range(sh_in[0]):
1709     if i0>0: t_prog+="+"
1710     t_prog+="%s*arg_in[%s]"%(trafo[j0,j1,j2,i0],i0)
1711     t_prog+="\n"
1712     elif len(sh)==4:
1713     for j0 in range(sh[0]):
1714     for j1 in range(sh[1]):
1715     for j2 in range(sh[2]):
1716     for j3 in range(sh[3]):
1717     t_prog+=" arg2[%s,%s,%s,%s]="%(j0,j1,j2,j3)
1718     if len(sh_in)==0:
1719     t_prog+="%s*arg_in"%trafo[j0,j1,j2,j3]
1720     t_prog+="\n"
1721     t_prog+=" dvdin=v.substitute({arg : arg2}).diff(arg_in)\n"
1722     if len(sh_in)==0:
1723     t_prog+=" res_in=dvdin.substitute({arg_in : %s})\n"%a_in
1724     else:
1725     t_prog+=" res_in=dvdin.substitute({arg : numarray.array(%s)})\n"%a_in.tolist()
1726    
1727     if len(sh)==0:
1728     if len(sh_in)==0:
1729     ref_diff=(makeResult(trafo*a_in+finc,f)-makeResult(trafo*a_in,f))/finc
1730     t_prog+=" self.failUnlessAlmostEqual(dvdin,%s,self.places,\"%s-derivative: wrong derivative\")\n"%(ref_diff,str(sh_in))
1731     elif len(sh_in)==1:
1732     s=0
1733     for k0 in range(sh_in[0]):
1734     s+=trafo[k0]*a_in[k0]
1735     for i0 in range(sh_in[0]):
1736     ref_diff=(makeResult(s+trafo[i0]*finc,f)-makeResult(s,f))/finc
1737     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s],%s,self.places,\"%s-derivative: wrong derivative with respect of %s\")\n"%(i0,ref_diff,str(sh_in),str(i0))
1738     elif len(sh_in)==2:
1739     s=0
1740     for k0 in range(sh_in[0]):
1741     for k1 in range(sh_in[1]):
1742     s+=trafo[k0,k1]*a_in[k0,k1]
1743     for i0 in range(sh_in[0]):
1744     for i1 in range(sh_in[1]):
1745     ref_diff=(makeResult(s+trafo[i0,i1]*finc,f)-makeResult(s,f))/finc
1746     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s],%s,self.places,\"%s-derivative: wrong derivative with respect of %s\")\n"%(i0,i1,ref_diff,str(sh_in),str((i0,i1)))
1747    
1748     elif len(sh_in)==3:
1749     s=0
1750     for k0 in range(sh_in[0]):
1751     for k1 in range(sh_in[1]):
1752     for k2 in range(sh_in[2]):
1753     s+=trafo[k0,k1,k2]*a_in[k0,k1,k2]
1754     for i0 in range(sh_in[0]):
1755     for i1 in range(sh_in[1]):
1756     for i2 in range(sh_in[2]):
1757     ref_diff=(makeResult(s+trafo[i0,i1,i2]*finc,f)-makeResult(s,f))/finc
1758     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative with respect of %s\")\n"%(i0,i1,i2,ref_diff,str(sh_in),str((i0,i1,i2)))
1759     elif len(sh_in)==4:
1760     s=0
1761     for k0 in range(sh_in[0]):
1762     for k1 in range(sh_in[1]):
1763     for k2 in range(sh_in[2]):
1764     for k3 in range(sh_in[3]):
1765     s+=trafo[k0,k1,k2,k3]*a_in[k0,k1,k2,k3]
1766     for i0 in range(sh_in[0]):
1767     for i1 in range(sh_in[1]):
1768     for i2 in range(sh_in[2]):
1769     for i3 in range(sh_in[3]):
1770     ref_diff=(makeResult(s+trafo[i0,i1,i2,i3]*finc,f)-makeResult(s,f))/finc
1771     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative with respect of %s\")\n"%(i0,i1,i2,i3,ref_diff,str(sh_in),str((i0,i1,i2,i3)))
1772     elif len(sh)==1:
1773     for j0 in range(sh[0]):
1774     if len(sh_in)==0:
1775     ref_diff=(makeResult(trafo[j0]*a_in+finc,f)-makeResult(trafo[j0]*a_in,f))/finc
1776     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s],%s,self.places,\"%s-derivative: wrong derivative of %s\")\n"%(j0,ref_diff,str(sh_in),j0)
1777     elif len(sh_in)==1:
1778     s=0
1779     for k0 in range(sh_in[0]):
1780     s+=trafo[j0,k0]*a_in[k0]
1781     for i0 in range(sh_in[0]):
1782     ref_diff=(makeResult(s+trafo[j0,i0]*finc,f)-makeResult(s,f))/finc
1783     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s with respect of %s\")\n"%(j0,i0,ref_diff,str(sh_in),str(j0),str(i0))
1784     elif len(sh_in)==2:
1785     s=0
1786     for k0 in range(sh_in[0]):
1787     for k1 in range(sh_in[1]):
1788     s+=trafo[j0,k0,k1]*a_in[k0,k1]
1789     for i0 in range(sh_in[0]):
1790     for i1 in range(sh_in[1]):
1791     ref_diff=(makeResult(s+trafo[j0,i0,i1]*finc,f)-makeResult(s,f))/finc
1792     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s with respect of %s\")\n"%(j0,i0,i1,ref_diff,str(sh_in),str(j0),str((i0,i1)))
1793    
1794     elif len(sh_in)==3:
1795    
1796     s=0
1797     for k0 in range(sh_in[0]):
1798     for k1 in range(sh_in[1]):
1799     for k2 in range(sh_in[2]):
1800     s+=trafo[j0,k0,k1,k2]*a_in[k0,k1,k2]
1801    
1802     for i0 in range(sh_in[0]):
1803     for i1 in range(sh_in[1]):
1804     for i2 in range(sh_in[2]):
1805     ref_diff=(makeResult(s+trafo[j0,i0,i1,i2]*finc,f)-makeResult(s,f))/finc
1806     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s with respect of %s\")\n"%(j0,i0,i1,i2,ref_diff,str(sh_in),str(j0),str((i0,i1,i2)))
1807     elif len(sh)==2:
1808     for j0 in range(sh[0]):
1809     for j1 in range(sh[1]):
1810     if len(sh_in)==0:
1811     ref_diff=(makeResult(trafo[j0,j1]*a_in+finc,f)-makeResult(trafo[j0,j1]*a_in,f))/finc
1812     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s\")\n"%(j0,j1,ref_diff,str(sh_in),str((j0,j1)))
1813     elif len(sh_in)==1:
1814     s=0
1815     for k0 in range(sh_in[0]):
1816     s+=trafo[j0,j1,k0]*a_in[k0]
1817     for i0 in range(sh_in[0]):
1818     ref_diff=(makeResult(s+trafo[j0,j1,i0]*finc,f)-makeResult(s,f))/finc
1819     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s with respect of %s\")\n"%(j0,j1,i0,ref_diff,str(sh_in),str((j0,j1)),str(i0))
1820    
1821     elif len(sh_in)==2:
1822     s=0
1823     for k0 in range(sh_in[0]):
1824     for k1 in range(sh_in[1]):
1825     s+=trafo[j0,j1,k0,k1]*a_in[k0,k1]
1826     for i0 in range(sh_in[0]):
1827     for i1 in range(sh_in[1]):
1828     ref_diff=(makeResult(s+trafo[j0,j1,i0,i1]*finc,f)-makeResult(s,f))/finc
1829     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s with respect of %s\")\n"%(j0,j1,i0,i1,ref_diff,str(sh_in),str((j0,j1)),str((i0,i1)))
1830     elif len(sh)==3:
1831     for j0 in range(sh[0]):
1832     for j1 in range(sh[1]):
1833     for j2 in range(sh[2]):
1834     if len(sh_in)==0:
1835     ref_diff=(makeResult(trafo[j0,j1,j2]*a_in+finc,f)-makeResult(trafo[j0,j1,j2]*a_in,f))/finc
1836     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s\")\n"%(j0,j1,j2,ref_diff,str(sh_in),str((j0,j1,j2)))
1837     elif len(sh_in)==1:
1838     s=0
1839     for k0 in range(sh_in[0]):
1840     s+=trafo[j0,j1,j2,k0]*a_in[k0]
1841     for i0 in range(sh_in[0]):
1842     ref_diff=(makeResult(s+trafo[j0,j1,j2,i0]*finc,f)-makeResult(s,f))/finc
1843     t_prog+=" self.failUnlessAlmostEqual(dvdin[%s,%s,%s,%s],%s,self.places,\"%s-derivative: wrong derivative of %s with respect of %s\")\n"%(j0,j1,j2,i0,ref_diff,str(sh_in),str((j0,j1,j2)),i0)
1844     elif len(sh)==4:
1845     for j0 in range(sh