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

Contents of /trunk/escript/py_src/generateutil

Parent Directory Parent Directory | Revision Log Revision Log


Revision 493 - (show annotations)
Fri Feb 3 02:18:45 2006 UTC (13 years, 9 months ago) by gross
File size: 130886 byte(s)
transpose added (just for me)
1 #!/usr/bin/python
2 # $Id$
3
4 """
5 program generates parts of the util.py and the test_util.py script
6 """
7 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 t_prog=""
26 t_prog_with_tags=""
27 t_prog_failing=""
28 u_prog=""
29
30 def wherepos(arg):
31 if arg>0.:
32 return 1.
33 else:
34 return 0.
35
36
37 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 import random
53 import numarray
54 import math
55 finc=1.e-6
56
57 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
152
153 def makeArray(shape,rng):
154 l=rng[1]-rng[0]
155 out=numarray.zeros(shape,numarray.Float64)
156 if len(shape)==0:
157 out=l*random.random()+rng[0]
158 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 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 else:
184 raise SystemError,"rank is restricted to 5"
185 return out
186
187
188 def makeResult(val,test_expr):
189 if isinstance(val,float):
190 out=eval(test_expr.replace("%a1%","val"))
191 elif len(val.shape)==0:
192 out=eval(test_expr.replace("%a1%","val"))
193 elif len(val.shape)==1:
194 out=numarray.zeros(val.shape,numarray.Float64)
195 for i0 in range(val.shape[0]):
196 out[i0]=eval(test_expr.replace("%a1%","val[i0]"))
197 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 out[i0,i1]=eval(test_expr.replace("%a1%","val[i0,i1]"))
202 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 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val[i0,i1,i2]"))
208 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 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val[i0,i1,i2,i3]"))
215 else:
216 raise SystemError,"rank is restricted to 4"
217 return out
218
219 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
436
437 def mkText(case,name,a,a1=None,use_tagging_for_expanded_data=False):
438 t_out=""
439 if case=="float":
440 if isinstance(a,float):
441 t_out+=" %s=%s\n"%(name,a)
442 elif a.rank==0:
443 t_out+=" %s=%s\n"%(name,a)
444 else:
445 t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
446 elif case=="array":
447 if isinstance(a,float):
448 t_out+=" %s=numarray.array(%s)\n"%(name,a)
449 elif a.rank==0:
450 t_out+=" %s=numarray.array(%s)\n"%(name,a)
451 else:
452 t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
453 elif case=="constData":
454 if isinstance(a,float):
455 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
456 elif a.rank==0:
457 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 elif a.rank==0:
465 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 t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
470 elif case=="expandedData":
471 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 elif a.rank==0:
476 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 elif a.rank==0:
487 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 elif case=="Symbol":
491 if isinstance(a,float):
492 t_out+=" %s=Symbol(shape=())\n"%(name)
493 elif a.rank==0:
494 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 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
515 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
528 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
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
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
700 def unrollLoops(a,b,o,arg,tap="",x="x"):
701 out=""
702 if a.rank==1:
703 z=""
704 for i99 in range(a.shape[0]):
705 if not z=="": z+="+"
706 if o=="1":
707 z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
708 else:
709 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
710
711 out+=tap+"%s=%s\n"%(arg,z)
712
713 elif a.rank==2:
714 for i0 in range(a.shape[0]):
715 z=""
716 for i99 in range(a.shape[1]):
717 if not z=="": z+="+"
718 if o=="1":
719 z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
720 else:
721 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
722
723 out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
724 elif a.rank==3:
725 for i0 in range(a.shape[0]):
726 for i1 in range(a.shape[1]):
727 z=""
728 for i99 in range(a.shape[2]):
729 if not z=="": z+="+"
730 if o=="1":
731 z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
732 else:
733 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
734
735 out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
736 elif a.rank==4:
737 for i0 in range(a.shape[0]):
738 for i1 in range(a.shape[1]):
739 for i2 in range(a.shape[2]):
740 z=""
741 for i99 in range(a.shape[3]):
742 if not z=="": z+="+"
743 if o=="1":
744 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
745 else:
746 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
747
748 out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
749 elif a.rank==5:
750 for i0 in range(a.shape[0]):
751 for i1 in range(a.shape[1]):
752 for i2 in range(a.shape[2]):
753 for i3 in range(a.shape[3]):
754 z=""
755 for i99 in range(a.shape[4]):
756 if not z=="": z+="+"
757 if o=="1":
758 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
759 else:
760 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
761
762 out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
763 return out
764
765 def unrollLoopsOfGrad(a,b,o,arg,tap=""):
766 out=""
767 if a.rank==1:
768 for i99 in range(a.shape[0]):
769 if o=="1":
770 out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
771 else:
772 out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
773
774 elif a.rank==2:
775 for i0 in range(a.shape[0]):
776 for i99 in range(a.shape[1]):
777 if o=="1":
778 out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
779 else:
780 out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
781
782 elif a.rank==3:
783 for i0 in range(a.shape[0]):
784 for i1 in range(a.shape[1]):
785 for i99 in range(a.shape[2]):
786 if o=="1":
787 out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
788 else:
789 out+=tap+"%s[%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99],i99,b[i0,i1,i99])
790
791 elif a.rank==4:
792 for i0 in range(a.shape[0]):
793 for i1 in range(a.shape[1]):
794 for i2 in range(a.shape[2]):
795 for i99 in range(a.shape[3]):
796 if o=="1":
797 out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
798 else:
799 out+=tap+"%s[%s,%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99],i99,b[i0,i1,i2,i99])
800 return out
801 def unrollLoopsOfDiv(a,b,o,arg,tap=""):
802 out=tap+arg+"="
803 if o=="1":
804 z=0.
805 for i99 in range(a.shape[0]):
806 z+=b[i99,i99]+a[i99,i99]
807 out+="(%s)"%z
808 else:
809 z=0.
810 for i99 in range(a.shape[0]):
811 z+=b[i99,i99]
812 if i99>0: out+="+"
813 out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
814 out+="+(%s)"%z
815 return out
816
817 def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
818 if where=="Function":
819 xfac_o=1.
820 xfac_op=0.
821 z_fac=1./2.
822 z_fac_s=""
823 zo_fac_s="/(o+1.)"
824 elif where=="FunctionOnBoundary":
825 xfac_o=1.
826 xfac_op=0.
827 z_fac=1.
828 z_fac_s="*dim"
829 zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
830 elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
831 xfac_o=0.
832 xfac_op=1.
833 z_fac=1./2.
834 z_fac_s=""
835 zo_fac_s="/(o+1.)"
836 out=""
837 if a.rank==1:
838 zo=0.
839 zop=0.
840 z=0.
841 for i99 in range(a.shape[0]):
842 if i99==0:
843 zo+= xfac_o*a[i99]
844 zop+= xfac_op*a[i99]
845 else:
846 zo+=a[i99]
847 z+=b[i99]
848
849 out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
850 if zop==0.:
851 out+="\n"
852 else:
853 out+="+(%s)*0.5**o\n"%zop
854 elif a.rank==2:
855 for i0 in range(a.shape[0]):
856 zo=0.
857 zop=0.
858 z=0.
859 for i99 in range(a.shape[1]):
860 if i99==0:
861 zo+= xfac_o*a[i0,i99]
862 zop+= xfac_op*a[i0,i99]
863 else:
864 zo+=a[i0,i99]
865 z+=b[i0,i99]
866
867 out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
868 if zop==0.:
869 out+="\n"
870 else:
871 out+="+(%s)*0.5**o\n"%zop
872 elif a.rank==3:
873 for i0 in range(a.shape[0]):
874 for i1 in range(a.shape[1]):
875 zo=0.
876 zop=0.
877 z=0.
878 for i99 in range(a.shape[2]):
879 if i99==0:
880 zo+= xfac_o*a[i0,i1,i99]
881 zop+= xfac_op*a[i0,i1,i99]
882 else:
883 zo+=a[i0,i1,i99]
884 z+=b[i0,i1,i99]
885
886 out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
887 if zop==0.:
888 out+="\n"
889 else:
890 out+="+(%s)*0.5**o\n"%zop
891 elif a.rank==4:
892 for i0 in range(a.shape[0]):
893 for i1 in range(a.shape[1]):
894 for i2 in range(a.shape[2]):
895 zo=0.
896 zop=0.
897 z=0.
898 for i99 in range(a.shape[3]):
899 if i99==0:
900 zo+= xfac_o*a[i0,i1,i2,i99]
901 zop+= xfac_op*a[i0,i1,i2,i99]
902
903 else:
904 zo+=a[i0,i1,i2,i99]
905 z+=b[i0,i1,i2,i99]
906
907 out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
908 if zop==0.:
909 out+="\n"
910 else:
911 out+="+(%s)*0.5**o\n"%zop
912
913 elif a.rank==5:
914 for i0 in range(a.shape[0]):
915 for i1 in range(a.shape[1]):
916 for i2 in range(a.shape[2]):
917 for i3 in range(a.shape[3]):
918 zo=0.
919 zop=0.
920 z=0.
921 for i99 in range(a.shape[4]):
922 if i99==0:
923 zo+= xfac_o*a[i0,i1,i2,i3,i99]
924 zop+= xfac_op*a[i0,i1,i2,i3,i99]
925
926 else:
927 zo+=a[i0,i1,i2,i3,i99]
928 z+=b[i0,i1,i2,i3,i99]
929 out+=tap+"%s[%s,%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,i3,zo,zo_fac_s,z*z_fac,z_fac_s)
930 if zop==0.:
931 out+="\n"
932 else:
933 out+="+(%s)*0.5**o\n"%zop
934
935 return out
936 def unrollLoopsSimplified(b,arg,tap=""):
937 out=""
938 if isinstance(b,float) or b.rank==0:
939 out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
940
941 elif b.rank==1:
942 for i0 in range(b.shape[0]):
943 out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
944 elif b.rank==2:
945 for i0 in range(b.shape[0]):
946 for i1 in range(b.shape[1]):
947 out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
948 elif b.rank==3:
949 for i0 in range(b.shape[0]):
950 for i1 in range(b.shape[1]):
951 for i2 in range(b.shape[2]):
952 out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
953 elif b.rank==4:
954 for i0 in range(b.shape[0]):
955 for i1 in range(b.shape[1]):
956 for i2 in range(b.shape[2]):
957 for i3 in range(b.shape[3]):
958 out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
959 return out
960
961 def unrollLoopsOfL2(b,where,arg,tap=""):
962 out=""
963 z=[]
964 if isinstance(b,float) or b.rank==0:
965 z.append(b**2)
966 elif b.rank==1:
967 for i0 in range(b.shape[0]):
968 z.append(b[i0]**2)
969 elif b.rank==2:
970 for i1 in range(b.shape[1]):
971 s=0
972 for i0 in range(b.shape[0]):
973 s+=b[i0,i1]**2
974 z.append(s)
975 elif b.rank==3:
976 for i2 in range(b.shape[2]):
977 s=0
978 for i0 in range(b.shape[0]):
979 for i1 in range(b.shape[1]):
980 s+=b[i0,i1,i2]**2
981 z.append(s)
982
983 elif b.rank==4:
984 for i3 in range(b.shape[3]):
985 s=0
986 for i0 in range(b.shape[0]):
987 for i1 in range(b.shape[1]):
988 for i2 in range(b.shape[2]):
989 s+=b[i0,i1,i2,i3]**2
990 z.append(s)
991 if where=="Function":
992 xfac_o=1.
993 xfac_op=0.
994 z_fac_s=""
995 zo_fac_s=""
996 zo_fac=1./3.
997 elif where=="FunctionOnBoundary":
998 xfac_o=1.
999 xfac_op=0.
1000 z_fac_s="*dim"
1001 zo_fac_s="*(2.*dim+1.)/3."
1002 zo_fac=1.
1003 elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1004 xfac_o=0.
1005 xfac_op=1.
1006 z_fac_s=""
1007 zo_fac_s=""
1008 zo_fac=1./3.
1009 zo=0.
1010 zop=0.
1011 for i99 in range(len(z)):
1012 if i99==0:
1013 zo+=xfac_o*z[i99]
1014 zop+=xfac_op*z[i99]
1015 else:
1016 zo+=z[i99]
1017 out+=tap+"%s=sqrt((%s)%s"%(arg,zo*zo_fac,zo_fac_s)
1018 if zop==0.:
1019 out+=")\n"
1020 else:
1021 out+="+(%s))\n"%(zop*0.5**2)
1022 return out
1023 #=======================================================================================================
1024 # transpose
1025 #=======================================================================================================
1026 def transposeTest(r,offset):
1027 if isinstance(r,float): return r
1028 s=r.shape
1029 s1=1
1030 for i in s[:offset]: s1*=i
1031 s2=1
1032 for i in s[offset:]: s2*=i
1033 out=numarray.reshape(r,(s1,s2))
1034 out.transpose()
1035 return numarray.resize(out,s[offset:]+s[:offset])
1036
1037 name,tt="transpose",transposeTest
1038 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1039 for sh0 in [ (), (3,), (4,5), (6,2,2),(3,2,3,4)]:
1040 for offset in range(len(sh0)+1):
1041 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1042 tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1043 text+=" def %s(self):\n"%tname
1044 sh_t=sh0[offset:]+sh0[:offset]
1045
1046 # sh_t=list(sh0)
1047 # sh_t[offset+1]=sh_t[offset]
1048 # sh_t=tuple(sh_t)
1049 # sh_r=[]
1050 # for i in range(offset): sh_r.append(sh0[i])
1051 # for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1052 # sh_r=tuple(sh_r)
1053
1054 a_0=makeArray(sh0,[-1.,1])
1055 if case0 in ["taggedData", "expandedData"]:
1056 a1_0=makeArray(sh0,[-1.,1])
1057 else:
1058 a1_0=a_0
1059 r=tt(a_0,offset)
1060 r1=tt(a1_0,offset)
1061 text+=mkText(case0,"arg",a_0,a1_0)
1062 text+=" res=%s(arg,%s)\n"%(name,offset)
1063 if case0=="Symbol":
1064 text+=mkText("array","s",a_0,a1_0)
1065 text+=" sub=res.substitute({arg:s})\n"
1066 res="sub"
1067 text+=mkText("array","ref",r,r1)
1068 else:
1069 res="res"
1070 text+=mkText(case0,"ref",r,r1)
1071 text+=mkTypeAndShapeTest(case0,sh_t,"res")
1072 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1073
1074 if case0 == "taggedData" :
1075 t_prog_with_tags+=text
1076 else:
1077 t_prog+=text
1078
1079 print test_header
1080 # print t_prog
1081 print t_prog_with_tags
1082 print test_tail
1083 1/0
1084 #=======================================================================================================
1085 # L2
1086 #=======================================================================================================
1087 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1088 for data in ["Data","Symbol"]:
1089 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1090 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1091 tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
1092 text+=" def %s(self):\n"%tname
1093 text+=" \"\"\"\n"
1094 text+=" tests L2-norm of %s on the %s\n\n"%(data,where)
1095 text+=" assumptions: self.domain supports integration on %s\n"%where
1096 text+=" \"\"\"\n"
1097 text+=" dim=self.domain.getDim()\n"
1098 text+=" w=%s(self.domain)\n"%where
1099 text+=" x=w.getX()\n"
1100 o="1"
1101 if len(sh)>0:
1102 sh_2=sh[:len(sh)-1]+(2,)
1103 sh_3=sh[:len(sh)-1]+(3,)
1104 b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
1105 b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
1106 else:
1107 sh_2=()
1108 sh_3=()
1109 b_2=makeArray(sh,[-1.,1])
1110 b_3=makeArray(sh,[-1.,1])
1111
1112 if data=="Symbol":
1113 val="s"
1114 res="sub"
1115 else:
1116 val="arg"
1117 res="res"
1118 text+=" if dim==2:\n"
1119 if data=="Symbol":
1120 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
1121
1122 text+=" %s=Data(0,%s,w)\n"%(val,sh_2)
1123 text+=unrollLoopsSimplified(b_2,val,tap=" ")
1124 text+=unrollLoopsOfL2(b_2,where,"ref",tap=" ")
1125 text+="\n else:\n"
1126 if data=="Symbol":
1127 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
1128 text+=" %s=Data(0,%s,w)\n"%(val,sh_3)
1129 text+=unrollLoopsSimplified(b_3,val,tap=" ")
1130 text+=unrollLoopsOfL2(b_3,where,"ref",tap=" ")
1131 text+="\n res=L2(arg)\n"
1132 if data=="Symbol":
1133 text+=" sub=res.substitute({arg:s})\n"
1134 text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1135 text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1136 else:
1137 text+=" self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1138 text+=" self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
1139 t_prog+=text
1140 print t_prog
1141 1/0
1142
1143 #=======================================================================================================
1144 # div
1145 #=======================================================================================================
1146 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1147 for data in ["Data","Symbol"]:
1148 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1149 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1150 tname="test_div_on%s_from%s_%s"%(where,data,case)
1151 text+=" def %s(self):\n"%tname
1152 text+=" \"\"\"\n"
1153 text+=" tests divergence of %s on the %s\n\n"%(data,where)
1154 text+=" assumptions: %s(self.domain) exists\n"%case
1155 text+=" self.domain supports div on %s\n"%where
1156 text+=" \"\"\"\n"
1157 if case=="ReducedSolution":
1158 text+=" o=1\n"
1159 o="1"
1160 else:
1161 text+=" o=self.order\n"
1162 o="o"
1163 text+=" dim=self.domain.getDim()\n"
1164 text+=" w_ref=%s(self.domain)\n"%where
1165 text+=" x_ref=w_ref.getX()\n"
1166 text+=" w=%s(self.domain)\n"%case
1167 text+=" x=w.getX()\n"
1168 a_2=makeArray((2,2),[-1.,1])
1169 b_2=makeArray((2,2),[-1.,1])
1170 a_3=makeArray((3,3),[-1.,1])
1171 b_3=makeArray((3,3),[-1.,1])
1172 if data=="Symbol":
1173 text+=" arg=Symbol(shape=(dim,),dim=dim)\n"
1174 val="s"
1175 res="sub"
1176 else:
1177 val="arg"
1178 res="res"
1179 text+=" %s=Vector(0,w)\n"%val
1180 text+=" if dim==2:\n"
1181 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1182 text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap=" ")
1183 text+="\n else:\n"
1184
1185 text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1186 text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap=" ")
1187 text+="\n res=div(arg,where=w_ref)\n"
1188 if data=="Symbol":
1189 text+=" sub=res.substitute({arg:s})\n"
1190 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1191 text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1192 text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
1193 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1194
1195
1196 t_prog+=text
1197 print t_prog
1198 1/0
1199
1200 #=======================================================================================================
1201 # interpolation
1202 #=======================================================================================================
1203 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
1204 for data in ["Data","Symbol"]:
1205 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1206 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1207 if where==case or \
1208 ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
1209 ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
1210 (case=="ContinuousFunction" and where in ["Solution","ReducedSolution"]) or \
1211 (case=="Solution" and where=="ReducedSolution") :
1212
1213
1214 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1215 tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1216 text+=" def %s(self):\n"%tname
1217 text+=" \"\"\"\n"
1218 text+=" tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
1219 text+=" assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
1220 text+=" \"\"\"\n"
1221 if case=="ReducedSolution" or where=="ReducedSolution":
1222 text+=" o=1\n"
1223 o="1"
1224 else:
1225 text+=" o=self.order\n"
1226 o="o"
1227 text+=" dim=self.domain.getDim()\n"
1228 text+=" w_ref=%s(self.domain)\n"%where
1229 text+=" x_ref=w_ref.getX()\n"
1230 text+=" w=%s(self.domain)\n"%case
1231 text+=" x=w.getX()\n"
1232 a_2=makeArray(sh+(2,),[-1.,1])
1233 b_2=makeArray(sh+(2,),[-1.,1])
1234 a_3=makeArray(sh+(3,),[-1.,1])
1235 b_3=makeArray(sh+(3,),[-1.,1])
1236 if data=="Symbol":
1237 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1238 val="s"
1239 res="sub"
1240 else:
1241 val="arg"
1242 res="res"
1243 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1244 text+=" ref=Data(0,%s,w_ref)\n"%str(sh)
1245 text+=" if dim==2:\n"
1246 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1247 text+=unrollLoops(a_2,b_2,o,"ref",tap=" ",x="x_ref")
1248 text+=" else:\n"
1249
1250 text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1251 text+=unrollLoops(a_3,b_3,o,"ref",tap=" ",x="x_ref")
1252 text+=" res=interpolate(arg,where=w_ref)\n"
1253 if data=="Symbol":
1254 text+=" sub=res.substitute({arg:s})\n"
1255 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1256 text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
1257 text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1258 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1259 t_prog+=text
1260 print test_header
1261 print t_prog
1262 print test_tail
1263 1/0
1264
1265 #=======================================================================================================
1266 # grad
1267 #=======================================================================================================
1268 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1269 for data in ["Data","Symbol"]:
1270 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1271 for sh in [ (),(2,), (4,5), (6,2,2)]:
1272 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1273 tname="test_grad_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1274 text+=" def %s(self):\n"%tname
1275 text+=" \"\"\"\n"
1276 text+=" tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
1277 text+=" assumptions: %s(self.domain) exists\n"%case
1278 text+=" self.domain supports gardient on %s\n"%where
1279 text+=" \"\"\"\n"
1280 if case=="ReducedSolution":
1281 text+=" o=1\n"
1282 o="1"
1283 else:
1284 text+=" o=self.order\n"
1285 o="o"
1286 text+=" dim=self.domain.getDim()\n"
1287 text+=" w_ref=%s(self.domain)\n"%where
1288 text+=" x_ref=w_ref.getX()\n"
1289 text+=" w=%s(self.domain)\n"%case
1290 text+=" x=w.getX()\n"
1291 a_2=makeArray(sh+(2,),[-1.,1])
1292 b_2=makeArray(sh+(2,),[-1.,1])
1293 a_3=makeArray(sh+(3,),[-1.,1])
1294 b_3=makeArray(sh+(3,),[-1.,1])
1295 if data=="Symbol":
1296 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1297 val="s"
1298 res="sub"
1299 else:
1300 val="arg"
1301 res="res"
1302 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1303 text+=" ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1304 text+=" if dim==2:\n"
1305 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1306 text+=unrollLoopsOfGrad(a_2,b_2,o,"ref",tap=" ")
1307 text+=" else:\n"
1308
1309 text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1310 text+=unrollLoopsOfGrad(a_3,b_3,o,"ref",tap=" ")
1311 text+=" res=grad(arg,where=w_ref)\n"
1312 if data=="Symbol":
1313 text+=" sub=res.substitute({arg:s})\n"
1314 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1315 text+=" self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1316 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1317
1318
1319 t_prog+=text
1320 print test_header
1321 print t_prog
1322 print test_tail
1323 1/0
1324
1325
1326 #=======================================================================================================
1327 # integrate
1328 #=======================================================================================================
1329 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1330 for data in ["Data","Symbol"]:
1331 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1332 for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1333 if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:
1334 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1335 tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1336 text+=" def %s(self):\n"%tname
1337 text+=" \"\"\"\n"
1338 text+=" tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
1339 text+=" assumptions: %s(self.domain) exists\n"%case
1340 text+=" self.domain supports integral on %s\n"%where
1341
1342 text+=" \"\"\"\n"
1343 if case=="ReducedSolution":
1344 text+=" o=1\n"
1345 o="1"
1346 else:
1347 text+=" o=self.order\n"
1348 o="o"
1349 text+=" dim=self.domain.getDim()\n"
1350 text+=" w_ref=%s(self.domain)\n"%where
1351 text+=" w=%s(self.domain)\n"%case
1352 text+=" x=w.getX()\n"
1353 a_2=makeArray(sh+(2,),[-1.,1])
1354 b_2=makeArray(sh+(2,),[-1.,1])
1355 a_3=makeArray(sh+(3,),[-1.,1])
1356 b_3=makeArray(sh+(3,),[-1.,1])
1357 if data=="Symbol":
1358 text+=" arg=Symbol(shape=%s)\n"%str(sh)
1359 val="s"
1360 res="sub"
1361 else:
1362 val="arg"
1363 res="res"
1364
1365 text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1366 if not len(sh)==0:
1367 text+=" ref=numarray.zeros(%s,numarray.Float)\n"%str(sh)
1368 text+=" if dim==2:\n"
1369 text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1370 text+=unrollLoopsOfInteriorIntegral(a_2,b_2,where,"ref",tap=" ")
1371 text+=" else:\n"
1372
1373 text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1374 text+=unrollLoopsOfInteriorIntegral(a_3,b_3,where,"ref",tap=" ")
1375 if case in ["ContinuousFunction","Solution","ReducedSolution"]:
1376 text+=" res=integrate(arg,where=w_ref)\n"
1377 else:
1378 text+=" res=integrate(arg)\n"
1379
1380 if data=="Symbol":
1381 text+=" sub=res.substitute({arg:s})\n"
1382 if len(sh)==0 and data=="Data":
1383 text+=" self.failUnless(isinstance(%s,float),\"wrong type of result.\")\n"%res
1384 else:
1385 if data=="Symbol":
1386 text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1387 text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1388 else:
1389 text+=" self.failUnless(isinstance(res,numarray.NumArray),\"wrong type of result.\")\n"
1390 text+=" self.failUnlessEqual(res.shape,%s,\"wrong shape of result.\")\n"%str(sh)
1391 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1392
1393
1394 t_prog+=text
1395 print test_header
1396 print t_prog
1397 print test_tail
1398 1/0
1399 #=======================================================================================================
1400 # inverse
1401 #=======================================================================================================
1402 name="inverse"
1403 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1404 for sh0 in [ (1,1), (2,2), (3,3)]:
1405 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1406 tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
1407 text+=" def %s(self):\n"%tname
1408 a_0=makeArray(sh0,[-1.,1])
1409 for i in range(sh0[0]): a_0[i,i]+=2
1410 if case0 in ["taggedData", "expandedData"]:
1411 a1_0=makeArray(sh0,[-1.,1])
1412 for i in range(sh0[0]): a1_0[i,i]+=3
1413 else:
1414 a1_0=a_0
1415
1416 text+=mkText(case0,"arg",a_0,a1_0)
1417 text+=" res=%s(arg)\n"%name
1418 if case0=="Symbol":
1419 text+=mkText("array","s",a_0,a1_0)
1420 text+=" sub=res.substitute({arg:s})\n"
1421 res="sub"
1422 ref="s"
1423 else:
1424 ref="arg"
1425 res="res"
1426 text+=mkTypeAndShapeTest(case0,sh0,"res")
1427 text+=" self.failUnless(Lsup(matrixmult(%s,%s)-kronecker(%s))<=self.RES_TOL,\"wrong result\")\n"%(res,ref,sh0[0])
1428
1429 if case0 == "taggedData" :
1430 t_prog_with_tags+=text
1431 else:
1432 t_prog+=text
1433
1434 print test_header
1435 # print t_prog
1436 print t_prog_with_tags
1437 print test_tail
1438 1/0
1439
1440 #=======================================================================================================
1441 # trace
1442 #=======================================================================================================
1443 def traceTest(r,offset):
1444 sh=r.shape
1445 r1=1
1446 for i in range(offset): r1*=sh[i]
1447 r2=1
1448 for i in range(offset+2,len(sh)): r2*=sh[i]
1449 r_s=numarray.reshape(r,(r1,sh[offset],sh[offset],r2))
1450 s=numarray.zeros([r1,r2],numarray.Float)
1451 for i1 in range(r1):
1452 for i2 in range(r2):
1453 for j in range(sh[offset]): s[i1,i2]+=r_s[i1,j,j,i2]
1454 return s.resize(sh[:offset]+sh[offset+2:])
1455 name,tt="trace",traceTest
1456 for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1457 for sh0 in [ (4,5), (6,2,2),(3,2,3,4)]:
1458 for offset in range(len(sh0)-1):
1459 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1460 tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1461 text+=" def %s(self):\n"%tname
1462 sh_t=list(sh0)
1463 sh_t[offset+1]=sh_t[offset]
1464 sh_t=tuple(sh_t)
1465 sh_r=[]
1466 for i in range(offset): sh_r.append(sh0[i])
1467 for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1468 sh_r=tuple(sh_r)
1469 a_0=makeArray(sh_t,[-1.,1])
1470 if case0 in ["taggedData", "expandedData"]:
1471 a1_0=makeArray(sh_t,[-1.,1])
1472 else:
1473 a1_0=a_0
1474 r=tt(a_0,offset)
1475 r1=tt(a1_0,offset)
1476 text+=mkText(case0,"arg",a_0,a1_0)
1477 text+=" res=%s(arg,%s)\n"%(name,offset)
1478 if case0=="Symbol":
1479 text+=mkText("array","s",a_0,a1_0)
1480 text+=" sub=res.substitute({arg:s})\n"
1481 res="sub"
1482 text+=mkText("array","ref",r,r1)
1483 else:
1484 res="res"
1485 text+=mkText(case0,"ref",r,r1)
1486 text+=mkTypeAndShapeTest(case0,sh_r,"res")
1487 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1488
1489 if case0 == "taggedData" :
1490 t_prog_with_tags+=text
1491 else:
1492 t_prog+=text
1493
1494 print test_header
1495 # print t_prog
1496 print t_prog_with_tags
1497 print test_tail
1498 1/0
1499
1500 #=======================================================================================================
1501 # clip
1502 #=======================================================================================================
1503 oper_L=[["clip",clipTEST]]
1504 for oper in oper_L:
1505 for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1506 for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1507 if len(sh0)==0 or not case0=="float":
1508 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1509 tname="test_%s_%s_rank%s"%(oper[0],case0,len(sh0))
1510 text+=" def %s(self):\n"%tname
1511 a_0=makeArray(sh0,[-1.,1])
1512 if case0 in ["taggedData", "expandedData"]:
1513 a1_0=makeArray(sh0,[-1.,1])
1514 else:
1515 a1_0=a_0
1516
1517 r=oper[1](a_0,-0.3,0.5)
1518 r1=oper[1](a1_0,-0.3,0.5)
1519 text+=mkText(case0,"arg",a_0,a1_0)
1520 text+=" res=%s(arg,-0.3,0.5)\n"%oper[0]
1521 if case0=="Symbol":
1522 text+=mkText("array","s",a_0,a1_0)
1523 text+=" sub=res.substitute({arg:s})\n"
1524 res="sub"
1525 text+=mkText("array","ref",r,r1)
1526 else:
1527 res="res"
1528 text+=mkText(case0,"ref",r,r1)
1529 text+=mkTypeAndShapeTest(case0,sh0,"res")
1530 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1531
1532 if case0 == "taggedData" :
1533 t_prog_with_tags+=text
1534 else:
1535 t_prog+=text
1536
1537 print test_header
1538 # print t_prog
1539 print t_prog_with_tags
1540 print test_tail
1541 1/0
1542
1543 #=======================================================================================================
1544 # maximum, minimum, clipping
1545 #=======================================================================================================
1546 oper_L=[ ["maximum",maximumTEST],
1547 ["minimum",minimumTEST]]
1548 for oper in oper_L:
1549 for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1550 for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1551 for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1552 for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1553 if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1554 and (sh0==sh1 or len(sh0)==0 or len(sh1)==0) :
1555 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1556
1557 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1558 tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1559 text+=" def %s(self):\n"%tname
1560 a_0=makeArray(sh0,[-1.,1])
1561 if case0 in ["taggedData", "expandedData"]:
1562 a1_0=makeArray(sh0,[-1.,1])
1563 else:
1564 a1_0=a_0
1565
1566 a_1=makeArray(sh1,[-1.,1])
1567 if case1 in ["taggedData", "expandedData"]:
1568 a1_1=makeArray(sh1,[-1.,1])
1569 else:
1570 a1_1=a_1
1571 r=oper[1](a_0,a_1)
1572 r1=oper[1](a1_0,a1_1)
1573 text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1574 text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1575 text+=" res=%s(arg0,arg1)\n"%oper[0]
1576 case=getResultCaseForBin(case0,case1)
1577 if case=="Symbol":
1578 c0_res,c1_res=case0,case1
1579 subs="{"
1580 if case0=="Symbol":
1581 text+=mkText("array","s0",a_0,a1_0)
1582 subs+="arg0:s0"
1583 c0_res="array"
1584 if case1=="Symbol":
1585 text+=mkText("array","s1",a_1,a1_1)
1586 if not subs.endswith("{"): subs+=","
1587 subs+="arg1:s1"
1588 c1_res="array"
1589 subs+="}"
1590 text+=" sub=res.substitute(%s)\n"%subs
1591 res="sub"
1592 text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1593 else:
1594 res="res"
1595 text+=mkText(case,"ref",r,r1)
1596 if len(sh0)>len(sh1):
1597 text+=mkTypeAndShapeTest(case,sh0,"res")
1598 else:
1599 text+=mkTypeAndShapeTest(case,sh1,"res")
1600 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1601
1602 if case0 == "taggedData" or case1 == "taggedData":
1603 t_prog_with_tags+=text
1604 else:
1605 t_prog+=text
1606
1607 print test_header
1608 # print t_prog
1609 print t_prog_with_tags
1610 print test_tail
1611 1/0
1612
1613
1614 #=======================================================================================================
1615 # outer inner
1616 #=======================================================================================================
1617 oper=["outer",outerTEST]
1618 # oper=["inner",innerTEST]
1619 for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1620 for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1621 for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1622 for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1623 if (len(sh0)==0 or not case0=="float") and (len(sh1)==0 or not case1=="float") \
1624 and len(sh0+sh1)<5:
1625 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1626
1627 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1628 tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1629 text+=" def %s(self):\n"%tname
1630 a_0=makeArray(sh0,[-1.,1])
1631 if case0 in ["taggedData", "expandedData"]:
1632 a1_0=makeArray(sh0,[-1.,1])
1633 else:
1634 a1_0=a_0
1635
1636 a_1=makeArray(sh1,[-1.,1])
1637 if case1 in ["taggedData", "expandedData"]:
1638 a1_1=makeArray(sh1,[-1.,1])
1639 else:
1640 a1_1=a_1
1641 r=oper[1](a_0,a_1)
1642 r1=oper[1](a1_0,a1_1)
1643 text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1644 text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1645 text+=" res=%s(arg0,arg1)\n"%oper[0]
1646 case=getResultCaseForBin(case0,case1)
1647 if case=="Symbol":
1648 c0_res,c1_res=case0,case1
1649 subs="{"
1650 if case0=="Symbol":
1651 text+=mkText("array","s0",a_0,a1_0)
1652 subs+="arg0:s0"
1653 c0_res="array"
1654 if case1=="Symbol":
1655 text+=mkText("array","s1",a_1,a1_1)
1656 if not subs.endswith("{"): subs+=","
1657 subs+="arg1:s1"
1658 c1_res="array"
1659 subs+="}"
1660 text+=" sub=res.substitute(%s)\n"%subs
1661 res="sub"
1662 text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1663 else:
1664 res="res"
1665 text+=mkText(case,"ref",r,r1)
1666 text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1667 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1668
1669 if case0 == "taggedData" or case1 == "taggedData":
1670 t_prog_with_tags+=text
1671 else:
1672 t_prog+=text
1673
1674 print test_header
1675 # print t_prog
1676 print t_prog_with_tags
1677 print test_tail
1678 1/0
1679
1680 #=======================================================================================================
1681 # local reduction
1682 #=======================================================================================================
1683 for oper in [["length",0.,"out+%a1%**2","math.sqrt(out)"],
1684 ["maxval",-1.e99,"max(out,%a1%)","out"],
1685 ["minval",1.e99,"min(out,%a1%)","out"] ]:
1686 for case in case_set:
1687 for sh in shape_set:
1688 if not case=="float" or len(sh)==0:
1689 text=""
1690 text+=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1691 tname="def test_%s_%s_rank%s"%(oper[0],case,len(sh))
1692 text+=" %s(self):\n"%tname
1693 a=makeArray(sh,[-1.,1.])
1694 a1=makeArray(sh,[-1.,1.])
1695 r1=testReduce(a1,oper[1],oper[2],oper[3])
1696 r=testReduce(a,oper[1],oper[2],oper[3])
1697
1698 text+=mkText(case,"arg",a,a1)
1699 text+=" res=%s(arg)\n"%oper[0]
1700 if case=="Symbol":
1701 text+=mkText("array","s",a,a1)
1702 text+=" sub=res.substitute({arg:s})\n"
1703 text+=mkText("array","ref",r,r1)
1704 res="sub"
1705 else:
1706 text+=mkText(case,"ref",r,r1)
1707 res="res"
1708 if oper[0]=="length":
1709 text+=mkTypeAndShapeTest(case,(),"res")
1710 else:
1711 if case=="float" or case=="array":
1712 text+=mkTypeAndShapeTest("float",(),"res")
1713 else:
1714 text+=mkTypeAndShapeTest(case,(),"res")
1715 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1716 if case == "taggedData":
1717 t_prog_with_tags+=text
1718 else:
1719 t_prog+=text
1720 print test_header
1721 # print t_prog
1722 print t_prog_with_tags
1723 print test_tail
1724 1/0
1725
1726 #=======================================================================================================
1727 # tensor multiply
1728 #=======================================================================================================
1729 # oper=["generalTensorProduct",tensorProductTest]
1730 # oper=["matrixmult",testMatrixMult]
1731 oper=["tensormult",testTensorMult]
1732
1733 for case0 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1734 for sh0 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1735 for case1 in ["float","array","Symbol","constData","taggedData","expandedData"]:
1736 for sh1 in [ (),(2,), (4,5), (6,2,2),(3,2,3,4)]:
1737 for sh_s in [ (),(3,), (2,3), (2,4,3),(4,2,3,2)]:
1738 if (len(sh0+sh_s)==0 or not case0=="float") and (len(sh1+sh_s)==0 or not case1=="float") \
1739 and len(sh0+sh1)<5 and len(sh0+sh_s)<5 and len(sh1+sh_s)<5:
1740 # if len(sh_s)==1 and len(sh0+sh_s)==2 and (len(sh_s+sh1)==1 or len(sh_s+sh1)==2)): # test for matrixmult
1741 if ( len(sh_s)==1 and len(sh0+sh_s)==2 and ( len(sh1+sh_s)==2 or len(sh1+sh_s)==1 )) or (len(sh_s)==2 and len(sh0+sh_s)==4 and (len(sh1+sh_s)==2 or len(sh1+sh_s)==3 or len(sh1+sh_s)==4)): # test for tensormult
1742 case=getResultCaseForBin(case0,case1)
1743 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1744 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1745 # tname="test_generalTensorProduct_%s_rank%s_%s_rank%s_offset%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1),len(sh_s))
1746 #tname="test_matrixmult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1747 tname="test_tensormult_%s_rank%s_%s_rank%s"%(case0,len(sh0+sh_s),case1,len(sh_s+sh1))
1748 # if tname=="test_generalTensorProduct_array_rank1_array_rank2_offset1":
1749 # print tnametest_generalTensorProduct_Symbol_rank1_Symbol_rank3_offset1
1750 text+=" def %s(self):\n"%tname
1751 a_0=makeArray(sh0+sh_s,[-1.,1])
1752 if case0 in ["taggedData", "expandedData"]:
1753 a1_0=makeArray(sh0+sh_s,[-1.,1])
1754 else:
1755 a1_0=a_0
1756
1757 a_1=makeArray(sh_s+sh1,[-1.,1])
1758 if case1 in ["taggedData", "expandedData"]:
1759 a1_1=makeArray(sh_s+sh1,[-1.,1])
1760 else:
1761 a1_1=a_1
1762 r=oper[1](a_0,a_1,sh_s)
1763 r1=oper[1](a1_0,a1_1,sh_s)
1764 text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1765 text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1766 #text+=" res=matrixmult(arg0,arg1)\n"
1767 text+=" res=tensormult(arg0,arg1)\n"
1768 #text+=" res=generalTensorProduct(arg0,arg1,offset=%s)\n"%(len(sh_s))
1769 if case=="Symbol":
1770 c0_res,c1_res=case0,case1
1771 subs="{"
1772 if case0=="Symbol":
1773 text+=mkText("array","s0",a_0,a1_0)
1774 subs+="arg0:s0"
1775 c0_res="array"
1776 if case1=="Symbol":
1777 text+=mkText("array","s1",a_1,a1_1)
1778 if not subs.endswith("{"): subs+=","
1779 subs+="arg1:s1"
1780 c1_res="array"
1781 subs+="}"
1782 text+=" sub=res.substitute(%s)\n"%subs
1783 res="sub"
1784 text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1785 else:
1786 res="res"
1787 text+=mkText(case,"ref",r,r1)
1788 text+=mkTypeAndShapeTest(case,sh0+sh1,"res")
1789 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1790 if case0 == "taggedData" or case1 == "taggedData":
1791 t_prog_with_tags+=text
1792 else:
1793 t_prog+=text
1794 print test_header
1795 # print t_prog
1796 print t_prog_with_tags
1797 print test_tail
1798 1/0
1799 #=======================================================================================================
1800 # basic binary operation overloading (tests only!)
1801 #=======================================================================================================
1802 oper_range=[-5.,5.]
1803 for oper in [["add" ,"+",[-5.,5.]],
1804 ["sub" ,"-",[-5.,5.]],
1805 ["mult","*",[-5.,5.]],
1806 ["div" ,"/",[-5.,5.]],
1807 ["pow" ,"**",[0.01,5.]]]:
1808 for case0 in case_set:
1809 for sh0 in shape_set:
1810 for case1 in case_set:
1811 for sh1 in shape_set:
1812 if not case0=="array" and \
1813 (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1814 (sh0==() or sh1==() or sh1==sh0) and \
1815 not (case0 in ["float","array"] and case1 in ["float","array"]):
1816 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1817 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1818 tname="test_%s_overloaded_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1819 text+=" def %s(self):\n"%tname
1820 a_0=makeArray(sh0,oper[2])
1821 if case0 in ["taggedData", "expandedData"]:
1822 a1_0=makeArray(sh0,oper[2])
1823 else:
1824 a1_0=a_0
1825
1826 a_1=makeArray(sh1,oper[2])
1827 if case1 in ["taggedData", "expandedData"]:
1828 a1_1=makeArray(sh1,oper[2])
1829 else:
1830 a1_1=a_1
1831 r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1832 r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1833 text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1834 text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1835 text+=" res=arg0%sarg1\n"%oper[1]
1836
1837 case=getResultCaseForBin(case0,case1)
1838 if case=="Symbol":
1839 c0_res,c1_res=case0,case1
1840 subs="{"
1841 if case0=="Symbol":
1842 text+=mkText("array","s0",a_0,a1_0)
1843 subs+="arg0:s0"
1844 c0_res="array"
1845 if case1=="Symbol":
1846 text+=mkText("array","s1",a_1,a1_1)
1847 if not subs.endswith("{"): subs+=","
1848 subs+="arg1:s1"
1849 c1_res="array"
1850 subs+="}"
1851 text+=" sub=res.substitute(%s)\n"%subs
1852 res="sub"
1853 text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1854 else:
1855 res="res"
1856 text+=mkText(case,"ref",r,r1)
1857 if isinstance(r,float):
1858 text+=mkTypeAndShapeTest(case,(),"res")
1859 else:
1860 text+=mkTypeAndShapeTest(case,r.shape,"res")
1861 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1862
1863 if case0 in [ "constData","taggedData","expandedData"] and case1 == "Symbol":
1864 t_prog_failing+=text
1865 else:
1866 if case0 == "taggedData" or case1 == "taggedData":
1867 t_prog_with_tags+=text
1868 else:
1869 t_prog+=text
1870
1871
1872 print test_header
1873 # print t_prog
1874 # print t_prog_with_tags
1875 print t_prog_failing
1876 print test_tail
1877 1/0
1878 #=======================================================================================================
1879 # basic binary operations (tests only!)
1880 #=======================================================================================================
1881 oper_range=[-5.,5.]
1882 for oper in [["add" ,"+",[-5.,5.]],
1883 ["mult","*",[-5.,5.]],
1884 ["quotient" ,"/",[-5.,5.]],
1885 ["power" ,"**",[0.01,5.]]]:
1886 for case0 in case_set:
1887 for case1 in case_set:
1888 for sh in shape_set:
1889 for sh_p in shape_set:
1890 if len(sh_p)>0:
1891 resource=[-1,1]
1892 else:
1893 resource=[1]
1894 for sh_d in resource:
1895 if sh_d>0:
1896 sh0=sh
1897 sh1=sh+sh_p
1898 else:
1899 sh1=sh
1900 sh0=sh+sh_p
1901
1902 if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1903 len(sh0)<5 and len(sh1)<5:
1904 use_tagging_for_expanded_data= case0=="taggedData" or case1=="taggedData"
1905 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1906 tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1907 text+=" def %s(self):\n"%tname
1908 a_0=makeArray(sh0,oper[2])
1909 if case0 in ["taggedData", "expandedData"]:
1910 a1_0=makeArray(sh0,oper[2])
1911 else:
1912 a1_0=a_0
1913
1914 a_1=makeArray(sh1,oper[2])
1915 if case1 in ["taggedData", "expandedData"]:
1916 a1_1=makeArray(sh1,oper[2])
1917 else:
1918 a1_1=a_1
1919 r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1920 r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1921 text+=mkText(case0,"arg0",a_0,a1_0,use_tagging_for_expanded_data)
1922 text+=mkText(case1,"arg1",a_1,a1_1,use_tagging_for_expanded_data)
1923 text+=" res=%s(arg0,arg1)\n"%oper[0]
1924
1925 case=getResultCaseForBin(case0,case1)
1926 if case=="Symbol":
1927 c0_res,c1_res=case0,case1
1928 subs="{"
1929 if case0=="Symbol":
1930 text+=mkText("array","s0",a_0,a1_0)
1931 subs+="arg0:s0"
1932 c0_res="array"
1933 if case1=="Symbol":
1934 text+=mkText("array","s1",a_1,a1_1)
1935 if not subs.endswith("{"): subs+=","
1936 subs+="arg1:s1"
1937 c1_res="array"
1938 subs+="}"
1939 text+=" sub=res.substitute(%s)\n"%subs
1940 res="sub"
1941 text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
1942 else:
1943 res="res"
1944 text+=mkText(case,"ref",r,r1)
1945 if isinstance(r,float):
1946 text+=mkTypeAndShapeTest(case,(),"res")
1947 else:
1948 text+=mkTypeAndShapeTest(case,r.shape,"res")
1949 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1950
1951 if case0 == "taggedData" or case1 == "taggedData":
1952 t_prog_with_tags+=text
1953 else:
1954 t_prog+=text
1955 print test_header
1956 # print t_prog
1957 print t_prog_with_tags
1958 print test_tail
1959 1/0
1960
1961 # print t_prog_with_tagsoper_range=[-5.,5.]
1962 for oper in [["add" ,"+",[-5.,5.]],
1963 ["sub" ,"-",[-5.,5.]],
1964 ["mult","*",[-5.,5.]],
1965 ["div" ,"/",[-5.,5.]],
1966 ["pow" ,"**",[0.01,5.]]]:
1967 for case0 in case_set:
1968 for sh0 in shape_set:
1969 for case1 in case_set:
1970 for sh1 in shape_set:
1971 if (not case0=="float" or len(sh0)==0) and (not case1=="float" or len(sh1)==0) and \
1972 (sh0==() or sh1==() or sh1==sh0) and \
1973 not (case0 in ["float","array"] and case1 in ["float","array"]):
1974 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1975 tname="test_%s_%s_rank%s_%s_rank%s"%(oper[0],case0,len(sh0),case1,len(sh1))
1976 text+=" def %s(self):\n"%tname
1977 a_0=makeArray(sh0,oper[2])
1978 if case0 in ["taggedData", "expandedData"]:
1979 a1_0=makeArray(sh0,oper[2])
1980 else:
1981 a1_0=a_0
1982
1983 a_1=makeArray(sh1,oper[2])
1984 if case1 in ["taggedData", "expandedData"]:
1985 a1_1=makeArray(sh1,oper[2])
1986 else:
1987 a1_1=a_1
1988 r1=makeResult2(a1_0,a1_1,"%a1%"+oper[1]+"%a2%")
1989 r=makeResult2(a_0,a_1,"%a1%"+oper[1]+"%a2%")
1990 text+=mkText(case0,"arg0",a_0,a1_0)
1991 text+=mkText(case1,"arg1",a_1,a1_1)
1992 text+=" res=arg0%sarg1\n"%oper[1]
1993
1994 case=getResultCaseForBin(case0,case1)
1995 if case=="Symbol":
1996 c0_res,c1_res=case0,case1
1997 subs="{"
1998 if case0=="Symbol":
1999 text+=mkText("array","s0",a_0,a1_0)
2000 subs+="arg0:s0"
2001 c0_res="array"
2002 if case1=="Symbol":
2003 text+=mkText("array","s1",a_1,a1_1)
2004 if not subs.endswith("{"): subs+=","
2005 subs+="arg1:s1"
2006 c1_res="array"
2007 subs+="}"
2008 text+=" sub=res.substitute(%s)\n"%subs
2009 res="sub"
2010 text+=mkText(getResultCaseForBin(c0_res,c1_res),"ref",r,r1)
2011 else:
2012 res="res"
2013 text+=mkText(case,"ref",r,r1)
2014 if isinstance(r,float):
2015 text+=mkTypeAndShapeTest(case,(),"res")
2016 else:
2017 text+=mkTypeAndShapeTest(case,r.shape,"res")
2018 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
2019
2020 if case0 in [ "constData","taggedData","expandedData"] and case1 == "Symbol":
2021 t_prog_failing+=text
2022 else:
2023 if case0 == "taggedData" or case1 == "taggedData":
2024 t_prog_with_tags+=text
2025 else:
2026 t_prog+=text
2027
2028
2029 # print u_prog
2030 # 1/0
2031 print test_header
2032 print t_prog
2033 # print t_prog_with_tags
2034 # print t_prog_failing
2035 print test_tail
2036 # print t_prog_failing
2037 print test_tail
2038
2039 #=======================================================================================================
2040 # unary operations:
2041 #=======================================================================================================
2042 func= [
2043 OPERATOR(nickname="log10",\
2044 rng=[1.e-3,100.],\
2045 test_expr="math.log10(%a1%)",\
2046 math_expr="math.log10(%a1%)",\
2047 numarray_expr="numarray.log10(%a1%)",\
2048 symbol_expr="log(%a1%)/log(10.)",\
2049 name="base-10 logarithm"),
2050 OPERATOR(nickname="wherePositive",\
2051 rng=[-100.,100.],\
2052 test_expr="wherepos(%a1%)",\
2053 math_expr="if arg>0:\n return 1.\nelse:\n return 0.",
2054 numarray_expr="numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))",\
2055 name="mask of positive values"),
2056 OPERATOR(nickname="whereNegative",\
2057 rng=[-100.,100.],\
2058 test_expr="wherepos(-%a1%)",\
2059 math_expr="if arg<0:\n return 1.\nelse:\n return 0.",
2060 numarray_expr="numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))",\
2061 name="mask of positive values"),
2062 OPERATOR(nickname="whereNonNegative",\
2063 rng=[-100.,100.],\
2064 test_expr="1-wherepos(-%a1%)", \
2065 math_expr="if arg<0:\n return 0.\nelse:\n return 1.",
2066 numarray_expr="numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))",\
2067 symbol_expr="1-wherePositive(%a1%)",\
2068 name="mask of non-negative values"),
2069 OPERATOR(nickname="whereNonPositive",\
2070 rng=[-100.,100.],\
2071 test_expr="1-wherepos(%a1%)",\
2072 math_expr="if arg>0:\n return 0.\nelse:\n return 1.",
2073 numarray_expr="numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))",\
2074 symbol_expr="1-whereNegative(%a1%)",\
2075 name="mask of non-positive values"),
2076 OPERATOR(nickname="whereZero",\
2077 rng=[-100.,100.],\
2078 test_expr="1-wherepos(%a1%)-wherepos(-%a1%)",\
2079 math_expr="if abs(%a1%)<=tol:\n return 1.\nelse:\n return 0.",
2080 numarray_expr="numarray.less_equal(abs(%a1%)-tol,numarray.zeros(arg.shape,numarray.Float))",\
2081 name="mask of zero entries"),
2082 OPERATOR(nickname="whereNonZero",\
2083 rng=[-100.,100.],\
2084 test_expr="wherepos(%a1%)+wherepos(-%a1%)",\
2085 math_expr="if abs(%a1%)>tol:\n return 1.\nelse:\n return 0.",\
2086 numarray_expr="numarray.greater(abs(%a1%)-tol,numarray.zeros(arg.shape,numarray.Float))",\
2087 symbol_expr="1-whereZero(arg,tol)",\
2088 name="mask of values different from zero"),
2089 OPERATOR(nickname="sin",\
2090 rng=[-100.,100.],\
2091 test_expr="math.sin(%a1%)",
2092 numarray_expr="numarray.sin(%a1%)",\
2093 diff="cos(%a1%)",\
2094 name="sine"),
2095 OPERATOR(nickname="cos",\
2096 rng=[-100.,100.],\
2097 test_expr="math.cos(%a1%)",
2098 numarray_expr="numarray.cos(%a1%)",\
2099 diff="-sin(%a1%)",
2100 name="cosine"),
2101 OPERATOR(nickname="tan",\
2102 rng=[-100.,100.],\
2103 test_expr="math.tan(%a1%)",
2104 numarray_expr="numarray.tan(%a1%)",\
2105 diff="1./cos(%a1%)**2",
2106 name="tangent"),
2107 OPERATOR(nickname="asin",\
2108 rng=[-0.99,0.99],\
2109 test_expr="math.asin(%a1%)",
2110 numarray_expr="numarray.arcsin(%a1%)",
2111 diff="1./sqrt(1.-%a1%**2)",
2112 name="inverse sine"),
2113 OPERATOR(nickname="acos",\
2114 rng=[-0.99,0.99],\
2115 test_expr="math.acos(%a1%)",
2116 numarray_expr="numarray.arccos(%a1%)",
2117 diff="-1./sqrt(1.-%a1%**2)",
2118 name="inverse cosine"),
2119 OPERATOR(nickname="atan",\
2120 rng=[-100.,100.],\
2121 test_expr="math.atan(%a1%)",
2122 numarray_expr="numarray.arctan(%a1%)",
2123 diff="1./(1+%a1%**2)",
2124 name="inverse tangent"),
2125 OPERATOR(nickname="sinh",\
2126 rng=[-5,5],\
2127 test_expr="math.sinh(%a1%)",
2128 numarray_expr="numarray.sinh(%a1%)",
2129 diff="cosh(%a1%)",
2130 name="hyperbolic sine"),
2131 OPERATOR(nickname="cosh",\
2132 rng=[-5.,5.],
2133 test_expr="math.cosh(%a1%)",
2134 numarray_expr="numarray.cosh(%a1%)",
2135 diff="sinh(%a1%)",
2136 name="hyperbolic cosine"),
2137 OPERATOR(nickname="tanh",\
2138 rng=[-5.,5.],
2139 test_expr="math.tanh(%a1%)",
2140 numarray_expr="numarray.tanh(%a1%)",
2141 diff="1./cosh(%a1%)**2",
2142 name="hyperbolic tangent"),
2143 OPERATOR(nickname="asinh",\
2144 rng=[-100.,100.], \
2145 test_expr="numarray.arcsinh(%a1%)",
2146 math_expr="numarray.arcsinh(%a1%)",
2147 numarray_expr="numarray.arcsinh(%a1%)",
2148 diff="1./sqrt(%a1%**2+1)",
2149 name="inverse hyperbolic sine"),
2150 OPERATOR(nickname="acosh",\
2151 rng=[1.001,100.],\
2152 test_expr="numarray.arccosh(%a1%)",
2153 math_expr="numarray.arccosh(%a1%)",
2154 numarray_expr="numarray.arccosh(%a1%)",
2155 diff="1./sqrt(%a1%**2-1)",
2156 name="inverse hyperolic cosine"),
2157 OPERATOR(nickname="atanh",\
2158 rng=[-0.99,0.99], \
2159 test_expr="numarray.arctanh(%a1%)",
2160 math_expr="numarray.arctanh(%a1%)",
2161 numarray_expr="numarray.arctanh(%a1%)",
2162 diff="1./(1.-%a1%**2)",
2163 name="inverse hyperbolic tangent"),
2164 OPERATOR(nickname="exp",\
2165 rng=[-5.,5.],
2166 test_expr="math.exp(%a1%)",
2167 numarray_expr="numarray.exp(%a1%)",
2168 diff="self",
2169 name="exponential"),
2170 OPERATOR(nickname="sqrt",\
2171 rng=[1.e-3,100.],\
2172 test_expr="math.sqrt(%a1%)",
2173 numarray_expr="numarray.sqrt(%a1%)",
2174 diff="0.5/self",
2175 name="square root"),
2176 OPERATOR(nickname="log", \
2177 rng=[1.e-3,100.],\
2178 test_expr="math.log(%a1%)",
2179 numarray_expr="numarray.log(%a1%)",
2180 diff="1./arg",
2181 name="natural logarithm"),
2182 OPERATOR(nickname="sign",\
2183 rng=[-100.,100.], \
2184 math_expr="if %a1%>0:\n return 1.\nelif %a1%<0:\n return -1.\nelse:\n return 0.",
2185 test_expr="wherepos(%a1%)-wherepos(-%a1%)",
2186 numarray_expr="numarray.sign(%a1%)",
2187 symbol_expr="wherePositive(%a1%)-whereNegative(%a1%)",\
2188 name="sign"),
2189 OPERATOR(nickname="abs",\
2190 rng=[-100.,100.], \
2191 math_expr="if %a1%>0:\n return %a1% \nelif %a1%<0:\n return -(%a1%)\nelse:\n return 0.",
2192 test_expr="wherepos(%a1%)*(%a1%)-wherepos(-%a1%)*(%a1%)",
2193 numarray_expr="abs(%a1%)",
2194 diff="sign(%a1%)",
2195 name="absolute value")
2196
2197 ]
2198 for f in func:
2199 symbol_name=f.nickname[0].upper()+f.nickname[1:]
2200 if f.nickname!="abs":
2201 u_prog+="def %s(arg):\n"%f.nickname
2202 u_prog+=" \"\"\"\n"
2203 u_prog+=" returns %s of argument arg\n\n"%f.name
2204 u_prog+=" @param arg: argument\n"
2205 u_prog+=" @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.\n"
2206 u_prog+=" @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.\n"
2207 u_prog+=" @raises TypeError: if the type of the argument is not expected.\n"
2208 u_prog+=" \"\"\"\n"
2209 u_prog+=" if isinstance(arg,numarray.NumArray):\n"
2210 u_prog+=mkCode(f.numarray_expr,["arg"],2*" ")
2211 u_prog+=" elif isinstance(arg,escript.Data):\n"
2212 u_prog+=mkCode("arg._%s()"%f.nickname,[],2*" ")
2213 u_prog+=" elif isinstance(arg,float):\n"
2214 u_prog+=mkCode(f.math_expr,["arg"],2*" ")
2215 u_prog+=" elif isinstance(arg,int):\n"
2216 u_prog+=mkCode(f.math_expr,["float(arg)"],2*" ")
2217 u_prog+=" elif isinstance(arg,Symbol):\n"
2218 if f.symbol_expr==None:
2219 u_prog+=mkCode("%s_Symbol(arg)"%symbol_name,[],2*" ")
2220 else:
2221 u_prog+=mkCode(f.symbol_expr,["arg"],2*" ")
2222 u_prog+=" else:\n"
2223 u_prog+=" raise TypeError,\"%s: Unknown argument type.\"\n\n"%f.nickname
2224 if f.symbol_expr==None:
2225 u_prog+="class %s_Symbol(DependendSymbol):\n"%symbol_name
2226 u_prog+=" \"\"\"\n"
2227 u_prog+=" L{Symbol} representing the result of the %s function\n"%f.name
2228 u_prog+=" \"\"\"\n"
2229 u_prog+=" def __init__(self,arg):\n"
2230 u_prog+=" \"\"\"\n"
2231 u_prog+=" initialization of %s L{Symbol} with argument arg\n"%f.nickname
2232 u_prog+=" @param arg: argument of function\n"
2233 u_prog+=" @type arg: typically L{Symbol}.\n"
2234 u_prog+=" \"\"\"\n"
2235 u_prog+=" DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())\n"
2236 u_prog+="\n"
2237
2238 u_prog+=" def getMyCode(self,argstrs,format=\"escript\"):\n"
2239 u_prog+=" \"\"\"\n"
2240 u_prog+=" returns a program code that can be used to evaluate the symbol.\n\n"
2241
2242 u_prog+=" @param argstrs: gives for each argument a string representing the argument for the evaluation.\n"
2243 u_prog+=" @type argstrs: C{str} or a C{list} of length 1 of C{str}.\n"
2244 u_prog+=" @param format: specifies the format to be used. At the moment only \"escript\" ,\"text\" and \"str\" are supported.\n"
2245 u_prog+=" @type format: C{str}\n"
2246 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"
2247 u_prog+=" @rtype: C{str}\n"
2248 u_prog+=" @raise: NotImplementedError: if the requested format is not available\n"
2249 u_prog+=" \"\"\"\n"
2250 u_prog+=" if isinstance(argstrs,list):\n"
2251 u_prog+=" argstrs=argstrs[0]\n"
2252 u_prog+=" if format==\"escript\" or format==\"str\" or format==\"text\":\n"
2253 u_prog+=" return \"%s(%%s)\"%%argstrs\n"%f.nickname
2254 u_prog+=" else:\n"
2255 u_prog+=" raise NotImplementedError,\"%s_Symbol does not provide program code for format %%s.\"%%format\n"%symbol_name
2256 u_prog+="\n"
2257
2258 u_prog+=" def substitute(self,argvals):\n"
2259 u_prog+=" \"\"\"\n"
2260 u_prog+=" assigns new values to symbols in the definition of the symbol.\n"
2261 u_prog+=" The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.\n"
2262 u_prog+="\n"
2263 u_prog+=" @param argvals: new values assigned to symbols\n"
2264 u_prog+=" @type argvals: C{dict} with keywords of type L{Symbol}.\n"
2265 u_prog+=" @return: result of the substitution process. Operations are executed as much as possible.\n"
2266 u_prog+=" @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution\n"
2267 u_prog+=" @raise TypeError: if a value for a L{Symbol} cannot be substituted.\n"
2268 u_prog+=" \"\"\"\n"
2269 u_prog+=" if argvals.has_key(self):\n"
2270 u_prog+=" arg=argvals[self]\n"
2271 u_prog+=" if self.isAppropriateValue(arg):\n"
2272 u_prog+=" return arg\n"
2273 u_prog+=" else:\n"
2274 u_prog+=" raise TypeError,\"%s: new value is not appropriate.\"%str(self)\n"
2275 u_prog+=" else:\n"
2276 u_prog+=" arg=self.getSubstitutedArguments(argvals)[0]\n"
2277 u_prog+=" return %s(arg)\n\n"%f.nickname
2278 if not f.diff==None:
2279 u_prog+=" def diff(self,arg):\n"
2280 u_prog+=" \"\"\"\n"
2281 u_prog+=" differential of this object\n"
2282 u_prog+="\n"
2283 u_prog+=" @param arg: the derivative is calculated with respect to arg\n"
2284 u_prog+=" @type arg: L{escript.Symbol}\n"
2285 u_prog+=" @return: derivative with respect to C{arg}\n"
2286 u_prog+=" @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray} are possible.\n"
2287 u_prog+=" \"\"\"\n"
2288 u_prog+=" if arg==self:\n"
2289 u_prog+=" return identity(self.getShape())\n"
2290 u_prog+=" else:\n"
2291 u_prog+=" myarg=self.getArgument()[0]\n"
2292 u_prog+=" val=matchShape(%s,self.getDifferentiatedArguments(arg)[0])\n"%f.diff.replace("%a1%","myarg")
2293 u_prog+=" return val[0]*val[1]\n\n"
2294
2295 for case in case_set:
2296 for sh in shape_set:
2297 if not case=="float" or len(sh)==0:
2298 text=""
2299 text+=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
2300 tname