Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] Diff of /branches/ein16/synth/d2/test_eval.py
 [diderot] / branches / ein16 / synth / d2 / test_eval.py

# Diff of /branches/ein16/synth/d2/test_eval.py

revision 3865, Mon May 16 16:32:19 2016 UTC revision 4436, Fri Aug 19 19:35:37 2016 UTC
# Line 1  Line 1
1    import sympy
2    from sympy import *
3    #symbols
4    x,y,z =symbols('x y z')
5  import sys  import sys
6  import re  import re
7  from test_examples import value_probe  import math
8  from test_examples import value_negation  from obj_ty import *
9  from test_examples import value_gradient  from obj_apply import *
10  from test_examples import value_scalarFd2  from obj_operator import *
11  from test_examples import value_vec2Fd2  from obj_field import *
from test_examples import value_scalarT
from test_examples import value_vec2T

#fixed variables
op_probe = value_probe()
op_negation = value_negation()
ty_scalarF_d2 = value_scalarFd2()
ty_vec2F_d2 = value_vec2Fd2()
ty_scalarT = value_scalarT()
ty_vec2T = value_vec2T()

def printFunc(a,b,c,d,e,f,g,h,i):
if(b!=0):
print b, "*x"
if(c!=0):
print " +", c,"y"
if(d!=0):
print " +", d, "x^2"
if(g!=0):
print " + ", g, "y*x^2"
if(f!=0):
print " + ", f, "y^2"
if(h!=0):
print " + ", h, "x*y^2"
if(a!=0):
print " + ", a
if(e!=0):
print " + ", e, "x*y"
if(i!=0):
print " + ", i, "x^2*y^2"
12
13
14
15    # ***************************  unary operators ***************************
16    # binary operators
18        return exp1+exp2
19    def fn_subtract(exp1,exp2):
20        return exp1-exp2
21    def fn_modulate(exp1,exp2):
22        return exp1*exp2
23    # scaling operator
24    def fn_multiplication(exp_s, t):
25        exp_t = field.get_data(t)
26        ityp1 = field.get_ty(t)
27        shape1 = fty.get_shape(ityp1)
28        #print "exp_s",exp_s,"exp_t",exp_t
29        if(field.is_Scalar(t)):
30            return  exp_s*  exp_t
31        elif(field.is_Vector(t)):
32            [n1] =  shape1 #vector
33            rtn = []
34            for i in range(n1):
35                rtn.append(exp_s*exp_t[i])
36            return rtn
37        elif(field.is_Matrix(t)):
38            [n1,n2] =  shape1
39            rtn = []
40            for i in range(n1):
41                tmp = []
42                for j in range(n2):
43                    tmp.append(exp_s*exp_t[i][j])
44                rtn.append(tmp)
45            return rtn
46        elif(field.is_Ten3(t)):
47            [n1,n2,n3] =  shape1
48            rtn = []
49            for i in range(n1):
50                tmpI = []
51                for j in range(n2):
52                    tmpJ = []
53                    for k in range(n3):
54                        tmpJ.append(exp_s*exp_t[i][j][k])
55                    tmpI.append(tmpJ)
56                rtn.append(tmpI)
57            return rtn
58        else:
59            raise "unsupported scaling"
60
61    #scaling of a field
62    def fn_scaling(fld1, fld2):
63        def get_sca():
64            if(field.is_Scalar(fld1)):
65                return (fld1, fld2)
66            else:
67                return (fld2, fld1)
68        (s, t) = get_sca()
69        exp_s = field.get_data(s)
70        return fn_multiplication(exp_s, t)
71
72    #division of a field
73    def fn_division(t, s):
74        if(field.is_Scalar(s)):
75            exp_s = (1/field.get_data(s))
76            return fn_multiplication(exp_s, t)
77        else:
78            raise Exception ("err second arg in division should be a scalar")
79
80
81    # sine  of field
82    def fn_negation(exp):
83        return -1*exp
84
85    def fn_cross(fld1, fld2):
86        exp1 = field.get_data(fld1)
87        ityp1 = field.get_ty(fld1)
88        exp2 = field.get_data(fld2)
89        ityp2 = field.get_ty(fld2)
90        #print " exp1: ",exp1," exp2: ",exp2
91        # vectors
92        n1 = fty.get_vecLength(ityp1) #vector
93        n2 = fty.get_vecLength(ityp2)
94        if(n1==2):
95            return (exp1[0]*exp2[1]) -(exp1[1]*exp2[0])
96        elif(n1==3):
97            x0= (exp1[1]*exp2[2]) -(exp1[2]*exp2[1])
98            x1= (exp1[2]*exp2[0]) -(exp1[0]*exp2[2])
99            x2= (exp1[0]*exp2[1]) -(exp1[1]*exp2[0])
100            return [x0,x1,x2]
101        else:
102            raise "unsupported type for cross product"
103
106        if (dim==1):
107            return diff(exp,x)
108        elif (dim==2):
109            return [diff(exp,x), diff(exp,y)]
110        elif (dim==3):
111            return [diff(exp,x), diff(exp,y), diff(exp,z)]
112        else:
113            raise "dimension not supported"
114
115    #evaluate divergence
116    def fn_divergence(fld):
117        exp = field.get_data(fld)
118        ityp = field.get_ty(fld)
119        #print " exp1: ",exp1," exp2: ",exp2
120        # vectors
121        n1 = fty.get_vecLength(ityp) #vector
122        if(n1==2):
123            return diff(exp[0],x)+diff(exp[1],y)
124        elif(n1==3):
125
126            return diff(exp[0],x)+diff(exp[1],y)+diff(exp[2],z)
127        else:
128            raise "unsupported type for divergence"
129
130    #evaluate cross product
131    def fn_curl(fld):
132        exp = field.get_data(fld)
133        ityp = field.get_ty(fld)
134        dim = field.get_dim(fld)
135        n = fty.get_vecLength(ityp) #vector
136        if(n!=dim):
137            raise (" type not supported for curl")
138        if(n==2):
139           return diff(exp[1], x) - diff(exp[0], y)
140        elif(n==3):
141            x0= diff(exp[2],y) - diff(exp[1],z)
142            x1= diff(exp[0],z) - diff(exp[2],x)
143            x2= diff(exp[1],x) - diff(exp[0],y)
144            return [x0,x1,x2]
145        else:
146            raise "unsupported type for cross product"
147
148    #evaluate jacob
149    def fn_jacob(fld):
150        #print "inside jacob"
151        exp = field.get_data(fld)
152        ityp = field.get_ty(fld)
153        dim = field.get_dim(fld)
154        #print " exp1: ",exp1," exp2: ",exp2
155        shape = fty.get_shape(ityp)
156        # vectors
157        if(field.is_VectorField (fld)):
158            print "found a vector"
159            [n] = shape #vector
160            if(n!=dim):
161                raise (" type not supported for jacob")
162            if(dim==2):
163                return [[diff(exp[0],x), diff(exp[0],y)],
164                        [diff(exp[1],x), diff(exp[1],y)]]
165            elif(dim==3):
166                return  [[diff(exp[0],x), diff(exp[0],y), diff(exp[0],z)],
167                         [diff(exp[1],x), diff(exp[1],y), diff(exp[1],z)],
168                         [diff(exp[2],x), diff(exp[2],y), diff(exp[2],z)]]
169            else:
170                raise (" type not supported for jacob")
171        elif(field.is_Matrix (fld)):
172            #print "inside matrix"
173            [n,m] = shape #vector
174            if(n!=dim and n!=m):
175                raise (" type not supported for jacob")
176            elif(dim==2):
177                e_00 = exp[0][0]
178                e_10 = exp[1][0]
179                e_01 = exp[0][1]
180                e_11 = exp[1][1]
181                f_00 = [diff(e_00,x), diff(e_00,y)]
182                f_01 = [diff(e_01,x), diff(e_01,y)]
183                f_10 = [diff(e_10,x), diff(e_10,y)]
184                f_11 = [diff(e_11,x), diff(e_11,y)]
185                return [[f_00,f_01],[f_10,f_11]]
186            elif(dim==3):
187                e_00 = exp[0][0]
188                e_10 = exp[1][0]
189                e_20 = exp[2][0]
190                e_01 = exp[0][1]
191                e_11 = exp[1][1]
192                e_21 = exp[2][1]
193                e_02 = exp[0][2]
194                e_12 = exp[1][2]
195                e_22 = exp[2][2]
196                def diff_d3(e):
197                    return [diff(e,x), diff(e,y), diff(e,z)]
198                f_00 = diff_d3(e_00)
199                f_01 = diff_d3(e_01)
200                f_02 = diff_d3(e_02)
201                f_10 = diff_d3(e_10)
202                f_11 = diff_d3(e_11)
203                f_12 = diff_d3(e_12)
204                f_20 = diff_d3(e_20)
205                f_21 = diff_d3(e_21)
206                f_22 = diff_d3(e_22)
207                return [[f_00,f_01,f_02],[f_10,f_11,f_12],[f_20,f_21,f_22]]
208            else:
209                raise (" type not supported for jacob")
210        else:
211            raise "unsupported type for jacob"
212
213    #evaluate norm
214    def fn_norm(fld, dim):
215        exp = field.get_data(fld)
216        ityp = field.get_ty(fld)
217        dim = field.get_dim(fld)
218        #print " exp1: ",exp1," exp2: ",exp2
219        # vectors
220        def iter (es):
221            sum = 0
222            for i in es:
223                t=i*i
224                #print "t",t
225                sum+=t
226            #print "\nsum",sum
227            rtn  = sqrt(sum)
228            #print "\nrtn",rtn
229            return rtn
230        if(field.is_Scalar(fld)):
231            [] = fty.get_shape(ityp)
232            #print "scalar exp:",exp
233            t =(sqrt(exp*exp))
234            #print "t",t
235            return t
236        elif(field.is_Vector(fld)):
237            [n] = fty.get_shape(ityp)
238            rtn = []
239            for i in range(n):
240                rtn.append(exp[i])
241            return iter(rtn)
242        elif(field.is_Matrix(fld)):
243            [n, m] = fty.get_shape(ityp)
244            rtn = []
245            for i in range(n):
246                for j in range(m):
247                    rtn.append(exp[i][j])
248            return iter(rtn)
249        elif(field.is_Ten3(fld)):
250            [n, m, p] = fty.get_shape(ityp)
251            rtn = []
252            for i in range(n):
253                for j in range(m):
254                    for k in range(p):
255                        rtn.append(exp[i][j][k])
256            return iter(rtn)
257        else:
258            raise "unsupported type for norm"
259
260    #evaluate norm
261    def fn_normalize(fld, dim):
262        exp = field.get_data(fld)
263        ityp = field.get_ty(fld)
264        dim = field.get_dim(fld)
265        #print " exp1: ",exp1," exp2: ",exp2
266        norm = fn_norm(fld, dim)
267        if(field.is_Scalar(fld)):
268            #print "scal",exp
269            return exp
270        elif(field.is_Vector(fld)):
271            [n] = fty.get_shape(ityp)
272            rtn = []
273            for i in range(n):
274                rtn.append(exp[i]/norm)
275            #print "vec",rtn
276            return rtn
277        elif(field.is_Matrix(fld)):
278            [n, m] = fty.get_shape(ityp)
279            rtn = []
280            for i in range(n):
281                rtni = []
282                for j in range(m):
283                    rtni.append(exp[i][j]/norm)
284                rtn.append(rtni)
285                #print "matrix:",rtn
286            return rtn
287        elif(field.is_Ten3(fld)):
288            [n, m, p] = fty.get_shape(ityp)
289            rtn = []
290            for i in range(n):
291                rtni = []
292                for j in range(m):
293                    rtnj = []
294                    for k in range(p):
295                        rtnj.append(exp[i][j][k]/norm)
296                    rtni.append( rtnj)
297                rtn.append( rtni)
298    #print "ten3",rtn
299            return rtn
300        else:
301            raise "unsupported type for norm"
302
303    #[0]
304    def fn_slicev0(fld1):
305        exp1 = field.get_data(fld1)
306        ityp1 = field.get_ty(fld1)
307        return exp1[0]
308
309
310    #evaluate slice
311    #[1]
312    def fn_slicev1(fld1):
313        exp1 = field.get_data(fld1)
314        ityp1 = field.get_ty(fld1)
315        return exp1[1]
316
317    #evaluate slice
318    #[1,:]
319    def fn_slicem0(fld1):
320        exp1 = field.get_data(fld1)
321        ityp1 = field.get_ty(fld1)
322        rtn=[]
323        if(fty.is_Matrix(ityp1)):
324            [n2,n3] = fty.get_shape(ityp1)
325            for i in range(n3):
326                rtn.append(exp1[1][i])
327            return rtn
328        else:
329            raise "unsupported type for slice"
330
331    #evaluate slice
332    #[:,0]
333    def fn_slicem1(fld1):
334        exp1 = field.get_data(fld1)
335        ityp1 = field.get_ty(fld1)
336        rtn=[]
337        if(fty.is_Matrix(ityp1)):
338            [n2,n3] = fty.get_shape(ityp1)
339            for i in range(n2):
340                rtn.append(exp1[i][0])
341            return rtn
342        else:
343            raise "unsupported type for slice"
344
345    #evaluate slice
346    #[:,1,:]
347    def fn_slicet0(fld1):
348        exp1 = field.get_data(fld1)
349        ityp1 = field.get_ty(fld1)
350        rtn=[]
351        if(fty.is_Ten3(ityp1)):
352            [n1, n2,n3] = fty.get_shape(ityp1)
353            for i in range(n1):
354                rtnj=[]
355                for j in range(n3):
356                    rtnj.append(exp1[i][1][j])
357                rtn.append(rtnj)
358            return rtn
359        else:
360            raise "unsupported type for slice"
361
362    #evaluate slice
363    #[1,0,:]
364    def fn_slicet1(fld1):
365        exp1 = field.get_data(fld1)
366        ityp1 = field.get_ty(fld1)
367        rtn=[]
368        if(fty.is_Ten3(ityp1)):
369            [n1, n2, n3] = fty.get_shape(ityp1)
370            for i in range(n3):
371                rtn.append(exp1[1][0][i])
372            return rtn
373        else:
374            raise "unsupported type for slice"
375
376
377
378
379    #evaluate trace
380    def fn_trace(fld):
381        exp = field.get_data(fld)
382        ityp = field.get_ty(fld)
383        rtn=[]
384        if(field.is_Matrix(fld)):
385            [n, m] = fty.get_shape(ityp)
386            if (n!=m):
387                raise Exception("matrix is not identitical")
388            rtn = exp[0][0]+exp[1][1]
389            if(n==2):
390                return rtn
391            elif(n==3):
392                return rtn+exp[2][2]
393            else:
394                raise "unsupported matrix field"
395        else:
396            raise "unsupported trace"
397
398    #evaluate transpose
399    def fn_transpose(fld):
400        exp = field.get_data(fld)
401        ityp = field.get_ty(fld)
402        if(field.is_Matrix(fld)):
403            [n, m] = fty.get_shape(ityp)
404            rtn = []
405            for i in range(n):
406                rtni = []
407                for j in range(m):
408                    rtni.append(exp[j][i])
409                rtn.append(rtni)
410            return rtn
411        else:
412            raise "unsupported transpose"
413
414    def fn_doubledot(fld1, fld2):
415        exp1 = field.get_data(fld1)
416        exp2 = field.get_data(fld2)
417        ityp1 = field.get_ty(fld1)
418        if(field.is_Matrix(fld1)):
419            rtn = 0
420            [n, m] = fty.get_shape(ityp1)
421            for i in range(n):
422                for j in range(m):
423                    rtn+=exp1[i][j]*exp2[i][j]
424            return rtn
425        elif(field.is_Ten3(fld1)):
426            rtn = []
427            [n, m, l] = fty.get_shape(ityp1)
428            for i in range(n):
429                for j in range(m):
430                    for k in range(l):
431                        rtn.append(exp1[i][j][k]*exp2[i][j][k])
432            return rtn
433        else:
434            raise "unsupported double dot"
435
436
437
438    #evaluate det
439    def fn_det(fld):
440        exp = field.get_data(fld)
441        ityp = field.get_ty(fld)
442        rtn=[]
443        if(field.is_Matrix(fld)):
444            [n, m] = fty.get_shape(ityp)
445            if (n!=m):
446                raise Exception("matrix is not identitical")
447            a = exp[0][0]
448            d = exp[1][1]
449            c = exp[1][0]
450            b = exp[0][1]
451            if(n==2):
452                return a*d-b*c
453            elif(n==3):
454                a = exp[0][0]
455                b = exp[0][1]
456                c = exp[0][2]
457                d = exp[1][0]
458                e = exp[1][1]
459                f = exp[1][2]
460                g = exp[2][0]
461                h = exp[2][1]
462                i = exp[2][2]
463                return a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)
464            else:
465                raise "unsupported matrix field"
466        else:
467            raise "unsupported trace"
468
469    #evaluate det
470    def fn_inverse(fld):
471        exp = field.get_data(fld)
472        ityp = field.get_ty(fld)
473        rtn=[]
474        if(field.is_Matrix(fld)):
475            [n, m] = fty.get_shape(ityp)
476            if (n!=m):
477                raise Exception("matrix is not identitical")
478            a = exp[0][0]
479            d = exp[1][1]
480            c = exp[1][0]
481            b = exp[0][1]
482            if(n==2):
483                denom = a*d-b*c
484                return [[d/denom,-b/denom],[-c/denom, a/denom]]
485            elif(n==3):
486                a = exp[0][0]
487                b = exp[0][1]
488                c = exp[0][2]
489                d = exp[1][0]
490                e = exp[1][1]
491                f = exp[1][2]
492                g = exp[2][0]
493                h = exp[2][1]
494                i = exp[2][2]
495                denom  =  a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)
496                num1 = [(e*i - f*h), -(b*i - c*h), (b*f - c*e)]
497                num2 = [-(d*i - f*g), (a*i -c*g), -(a*f - c*d)]
498                num3 = [(d*h - e*g), -(a*h - b*g),(a*e - b*d)]
499                num =[num1,num2,num3]
500
501                rtn = []
502                for i in range(n):
503                    tmp = []
504                    for j in range(n):
505                        tmp.append(num[i][j]/denom)
506                    rtn.append(tmp)
507                return rtn
508            else:
509                raise "unsupported matrix field"
510        else:
511            raise "unsupported trace"
512
513
514
515
516
517    #evaluate outer product
518    def fn_outer(fld1, fld2):
519        exp1 = field.get_data(fld1)
520        ityp1 = field.get_ty(fld1)
521        exp2 = field.get_data(fld2)
522        ityp2 = field.get_ty(fld2)
523        ashape = fty.get_shape(ityp1)
524        bshape = fty.get_shape(ityp2)
525        x= "ashape", ashape, "bshape", bshape
526        #print "exp1",exp1,"ityp1",ityp1.name,"-length",len(exp1)
527        #print "exp2",exp2,"ityp2",ityp2.name,"-length",len(exp2)
528        rtn = []
529        if(fty.is_Vector(ityp1)):
530            [n1] = fty.get_shape(ityp1)
531            if(fty.is_Vector(ityp2)):
532                #both vectors
533                [n2] = fty.get_shape(ityp2)
534                for i in  range(n1):
535                    tmpI = []
536                    for j in range(n2):
537                        tmpI.append(exp1[i]*exp2[j])
538                    rtn.append(tmpI)
539                return rtn
540            elif(fty.is_Matrix(ityp2)):
541                [n2,n3] = fty.get_shape(ityp2)
542                for i in  range(n1):
543                    tmpI = []
544                    for j in range(n2):
545                        tmpJ = []
546                        for k in range(n3):
547                            tmpJ.append(exp1[i]*exp2[j][k])
548                        tmpI.append(tmpJ)
549                    rtn.append(tmpI)
550                return rtn
551        elif(fty.is_Matrix(ityp1)):
552            [n1,n2] = fty.get_shape(ityp1)
553            if(fty.is_Vector(ityp2)):
554                [n3] = fty.get_shape(ityp2)
555                for i in  range(n1):
556                    tmpI = []
557                    for j in range(n2):
558                        tmpJ = []
559                        for k in range(n3):
560                            tmpJ.append(exp1[i][j]*exp2[k])
561                        tmpI.append(tmpJ)
562                    rtn.append(tmpI)
563                return rtn
564            elif(fty.is_Matrix(ityp2)):
565                [n3, n4] = fty.get_shape(ityp2)
566                for i in  range(n1):
567                    tmpI = []
568                    for j in range(n2):
569                        tmpJ = []
570                        for k in range(n3):
571                            tmpK = []
572                            for l in range(n4):
573                                tmpK.append(exp1[i][j]*exp2[k][l])
574                            tmpJ.append(tmpK)
575                        tmpI.append(tmpJ)
576                    rtn.append(tmpI)
577                return rtn
578            else:
579                raise "outer product is not supported"
580        else:
581            raise "outer product is not supported"
582    #evaluate inner product
583    def fn_inner(fld1, fld2):
584        exp1 = field.get_data(fld1)
585        ityp1 = field.get_ty(fld1)
586        exp2 = field.get_data(fld2)
587        ityp2 = field.get_ty(fld2)
588        ashape = fty.get_shape(ityp1)
589        bshape = fty.get_shape(ityp2)
590        x= "ashape", ashape, "bshape", bshape
591        if(fty.is_Vector(ityp1)):
592            [a1] = fty.get_shape(ityp1)
593            if(fty.is_Vector(ityp2)):
594                #length of vetors
595                rtn=0
596                [b1] = fty.get_shape(ityp2)
597                if(a1!=b1):
598                    raise x
599                for s in  range(a1):
600                    rtn += exp1[s]*exp2[s]
601                return rtn
602            elif(fty.is_Matrix(ityp2)):
603                [b1,b2] = fty.get_shape(ityp2)
604                rtn=[]
605                if(a1!=b1):
606                    raise x
607                for i in  range(b2):
608                    sumrtn=0
609                    for s in  range(a1):
610                        sumrtn +=  exp1[s]*exp2[s][i]
611                    rtn.append(sumrtn)
612                return rtn
613            elif(fty.is_Ten3(ityp2)):
614                [b1, b2, b3] = fty.get_shape(ityp2)
615                rtn = []
616                if(a1!=b1):
617                    raise x
618                for i in  range(b2):
619                    tmpJ = []
620                    for j in  range(b3):
621                        sumrtn=0
622                        for s in  range(a1):
623                            sumrtn +=  exp1[s]*exp2[s][i][j]
624                        tmpJ.append(sumrtn)
625                    rtn.append(tmpJ)
626                return rtn
627            else:
628                raise "inner product is not supported"
629        elif(fty.is_Matrix(ityp1)):
630            [a1,a2] = fty.get_shape(ityp1)
631            if(fty.is_Vector(ityp2)):
632                [b1] = fty.get_shape(ityp2)
633                if(a2!=b1):
634                    raise x
635                rtn=[]
636                for i in  range(a1):
637                    sumrtn=0
638                    for s in  range(a2):
639                        sumrtn += exp1[i][s]*exp2[s]
640                    rtn.append(sumrtn)
641                return rtn
642            elif(fty.is_Matrix(ityp2)):
643                [b1,b2] = fty.get_shape(ityp2)
644                rtn=[]
645                if(a2!=b1):
646                    raise x
647                for i in  range(a1):
648                    rtnj = []
649                    for j in  range(b2):
650                        sumrtn=0
651                        for s in  range(a2):
652                            sumrtn += exp1[i][s]*exp2[s][j]
653                        rtnj.append(sumrtn)
654                    rtn.append(rtnj)
655                return rtn
656            elif(fty.is_Ten3(ityp2)):
657                [b1,b2, b3] = fty.get_shape(ityp2)
658                rtn=[]
659                if(a2!=b1):
660                    raise x
661                for i in  range(a1):
662                    rtnj = []
663                    for j in  range(b2):
664                        rtnk = []
665                        for k in range(b3):
666                            sumrtn=0
667                            for s in  range(a2):
668                                sumrtn += exp1[i][s]*exp2[s][j][k]
669                            rtnk.append(sumrtn)
670                        rtnj.append(rtnk)
671                    rtn.append(rtnj)
672                return rtn
673
674            else:
675                raise "inner product is not supported"
676        elif(fty.is_Ten3(ityp1)):
677            [a1,a2, a3] = ashape
678            if(fty.is_Vector(ityp2)):
679                [b1] = bshape
680                if(a3!=b1):
681                    raise x
682                rtn=[]
683                for i in  range(a1):
684                    tmpI=[]
685                    for j in  range(a2):
686                        sumrtn=0
687                        for s in  range(a3):
688                            sumrtn += exp1[i][j][s]*exp2[s]
689                        tmpI.append(sumrtn)
690                    rtn.append(tmpI)
691                return rtn
692            if(fty.is_Matrix(ityp2)):
693                [b1,b2] = bshape
694                if(a3!=b1):
695                    raise x
696                rtn=[]
697                for i in  range(a1):
698                    tmpI=[]
699                    for j in  range(a2):
700                        tmpJ = []
701                        for k in range(b2):
702                            sumrtn=0
703                            for s in  range(a3):
704                                sumrtn += exp1[i][j][s]*exp2[s][k]
705                            tmpJ.append(sumrtn)
706                        tmpI.append(tmpJ)
707                    rtn.append(tmpI)
708                return rtn
709            else:
710                raise "inner product is not supported"
711        else:
712            raise "inner product is not supported"
713
714
715    # ***************************  generic apply operators ***************************
716    #unary operator on a vector
717    def applyToVector(vec, unary):
718        rtn = []
719        for v in vec:
720            rtn.append(unary(v))
721        return rtn
722    #binary operator on a vector
723    def applyToVectors(vecA, vecB,  binary):
724        rtn = []
725        for (a,b) in zip(vecA,vecB):
726            x= binary(a,b)
727            rtn.append(x)
728        return rtn
729
730    def applyToM(vec, unary):
731        rtn = []
732        for i in vec:
733            tmpI = []
734            for v in i:
735                tmpI.append(unary(v))
736            rtn.append(tmpI)
737        return rtn
738
739    def applyToMs(vecA,vecB, unary):
740        rtn = []
741        for (a,b) in zip(vecA,vecB):
742            tmpI = []
743            for (u,v) in zip(a,b):
744                tmpI.append(unary(u, v))
745            rtn.append(tmpI)
746        return rtn
747
748
749    def applyToT3(vec, unary):
750        rtn = []
751        for i in vec:
752            tmpI = []
753            for j in i:
754                tmpJ = []
755                for v in j:
756                    tmpJ.append(unary(v))
757                tmpI.append(tmpJ)
758            rtn.append(tmpI)
759        return rtn
760
761    def applyToT3s(vecA, vecB, unary):
762        rtn = []
763        for (a,b) in zip(vecA, vecB):
764            tmpI = []
765            for (i,j) in zip(a,b):
766                tmpJ = []
767                for (u,v) in zip(i,j):
768                    tmpJ.append(unary(u, v))
769                tmpI.append(tmpJ)
770            rtn.append(tmpI)
771        return rtn
772
773
774    # ***************************  apply to scalars or vectors ***************************
775    #apply operators to expression
776    # return output types and expression
777    # unary operator
778    # exp: scalar types
779
780    def applyToExp_U_S(fn_name, fld):
781        exp = field.get_data(fld)
782        dim = field.get_dim(fld)
783        #print fn_name
784        if(op_copy==fn_name): #probing
785            return  exp
786        elif(op_negation==fn_name): #negation
787            return fn_negation(exp)
788        elif(op_cosine==fn_name): #cosine
789            return cos(exp)
790        elif(op_sine==fn_name): #sine
791            return sin(exp)
792        elif(op_tangent==fn_name): # tangent
793            return tan(exp)
794        elif(op_asine==fn_name): #asine mag(x)<=1
795            return asin(0.1*exp)
796        elif(op_acosine==fn_name): #acos  mag(x)<=1
797            return acos(0.1*exp)
798        elif(op_atangent==fn_name): #atan
799            return atan(exp)
800        elif(op_sqrt==fn_name): #sqrt
801            return sqrt(exp)
804        else:
805            raise Exception("unsupported unary operator on scalar field:"+ fn_name.name)
806
807    # unary operator
808    # exp: vector  types
809    def applyToExp_U_V(fn_name, fld):
810        exp = field.get_data(fld)
811        if(op_copy==fn_name): #probing
812            return exp
813        elif(op_negation==fn_name): #negation
814            return applyToVector(exp, fn_negation)
815        elif(op_divergence==fn_name):
816            return fn_divergence(fld)
817        elif(op_curl==fn_name):
818            return fn_curl(fld)
819        elif(op_jacob==fn_name): #jacob
820            return fn_jacob(fld)
821        elif(op_slicev0==fn_name) :
822            return fn_slicev0(fld)
823        elif(op_slicev1==fn_name):
824            return fn_slicev1(fld)
827        else:
828            raise Exception("unsupported unary operator:"+ fn_name.name)
829
830    def applyToExp_U_M(fn_name, fld):
831        exp = field.get_data(fld)
832        if(op_copy==fn_name): #probing
833            return exp
834        elif(op_negation==fn_name): #negation
835            return applyToM(exp, fn_negation)
836        elif(op_jacob==fn_name): #jacob
837            #print "app u-m"
838            x = fn_jacob(fld)
839            #print "x", x
840            return x
841        elif(op_slicem0==fn_name) :
842            return fn_slicem0(fld)
843        elif(op_slicem1==fn_name):
844            return fn_slicem1(fld)
845        elif(op_trace == fn_name):
846            return fn_trace(fld)
847        elif(op_transpose==fn_name):
848            return fn_transpose(fld)
849        elif(op_det==fn_name):
850            return fn_det(fld)
851        elif(op_inverse == fn_name):
852            return fn_inverse(fld)
853        else:
854            raise Exception("unsupported unary operator:"+ fn_name.name)
855
856    def applyToExp_U_T3(fn_name, fld):
857        exp = field.get_data(fld)
858        if(op_copy==fn_name): #probing
859            return exp
860        elif(op_negation==fn_name): #negation
861            return applyToT3(exp, fn_negation)
862        elif(op_jacob==fn_name): #jacob
863            return fn_jacob(fld)
864        elif(op_slicet0==fn_name) :
865            return fn_slicet0(fld)
866        elif(op_slicet1==fn_name):
867            return fn_slicet1(fld)
868        else:
869            raise Exception("unsupported unary operator:"+ fn_name.name)
870
871    # binary, args do not need to have the same shape
872    def applyToExp_B_rest(e):
873        fn_name=e.opr
874        (fld1,fld2) =  apply.get_binary(e)
875        exp1 = field.get_data(fld1)
876        exp2 = field.get_data(fld2)
877        if(op_outer==fn_name):
878            return fn_outer(fld1, fld2)
879        elif(op_inner==fn_name):
880            return fn_inner(fld1, fld2)
881        elif(op_scale==fn_name): #scaling
882            return fn_scaling(fld1,fld2)
883        else:
884            raise Exception("unsupported unary operator:"+fn_name.name)
885
886    # binary operator
887    # exp: scalar types
888    def applyToExp_B_S(e):
889        fn_name=e.opr
890
891        (fld1,fld2) =  apply.get_binary(e)
892        exp1 = field.get_data(fld1)
893        exp2 = field.get_data(fld2)
894        #print fn_name
897        elif(op_subtract==fn_name):#subtract
898            return fn_subtract(exp1,exp2)
899        elif(op_modulate==fn_name):#modulate
900            return fn_modulate(exp1,exp2)
901        elif(op_scale==fn_name): #scaling
902            return fn_scaling(fld1,fld2)
903        elif(op_division==fn_name): #division
904            return fn_division(fld1,fld2)
905        else:
906            raise Exception("unsupported binary operator on scalar fields:"+ fn_name.name)
907
908
909    # binary operator
910    # args have the same shape
911    def applyToExp_B_V(e):
912
913        fn_name=e.opr
914        (fld1,fld2) =  apply.get_binary(e)
915        exp1 = field.get_data(fld1)
916        exp2 = field.get_data(fld2)
919        elif(op_subtract==fn_name):#subtract
920            return  applyToVectors(exp1, exp2, fn_subtract)
921        elif(op_modulate==fn_name):#modulate
922            return applyToVectors(exp1,exp2 ,fn_modulate)
923        elif(op_cross==fn_name):
924            return fn_cross(fld1, fld2)
925        else:
926           return applyToExp_B_rest(e)
927
928    def applyToExp_B_M(e):
929        print "b. m"
930        fn_name=e.opr
931        (fld1,fld2) =  apply.get_binary(e)
932        exp1 = field.get_data(fld1)
933        exp2 = field.get_data(fld2)
936        elif(op_subtract==fn_name):#subtract
937            return  applyToMs(exp1, exp2, fn_subtract)
938        elif(op_modulate==fn_name):#modulate
939            return applyToMs(exp1,exp2,fn_modulate)
940        else:
941            return applyToExp_B_rest(e)
942
943    def applyToExp_B_T3(e):
944        print "b. m"
945        fn_name=e.opr
946        (fld1,fld2) =  apply.get_binary(e)
947        exp1 = field.get_data(fld1)
948        exp2 = field.get_data(fld2)
951        elif(op_subtract==fn_name):#subtract
952            return  applyToT3s(exp1, exp2, fn_subtract)
953        elif(op_modulate==fn_name):#modulate
954            return applyToT3s(exp1,exp2,fn_modulate)
955        else:
956            return applyToExp_B_rest(e)
957
958
959    # ***************************  unary/binary operators ***************************
960    def unary(e):
961        #apply.toStr(e)
962        fld =apply.get_unary(e)
963        fn_name=e.opr
964        exp = field.get_data(fld)
965        dim = field.get_dim(fld)
966        if(op_norm==fn_name):#norm
967            return fn_norm(fld, dim)
968        elif(op_normalize==fn_name):#normalize
969            x= fn_normalize(fld, dim)
970            return x
971        elif (field.is_Scalar(fld)): # input is a scalar field
972            return applyToExp_U_S(fn_name, fld)
973        elif(field.is_Vector(fld)): # input is a vector field
974            return applyToExp_U_V(fn_name, fld)
975        elif(field.is_Matrix(fld)): # input is a vector field
976            return applyToExp_U_M(fn_name, fld)
977        else:
978            return applyToExp_U_T3(fn_name, fld)
979
980    def binary(e):
981        (f, g) =apply.get_binary(e)
982        fn_name = e.opr
983        # type is checked elsewhere or does not matter
984        if(op_division==fn_name): #division
985            return fn_division(f, g)
986        elif(op_doubledot==fn_name):#op_doubledot
987            return fn_doubledot (f, g)
988        elif(op_outer==fn_name):
989            return fn_outer(f, g)
990        elif(op_inner==fn_name):
991            return fn_inner(f, g)
992        elif(op_scale==fn_name): #scaling
993            return fn_scaling(f,g)
994        elif (field.is_Scalar(f) and field.is_Scalar(g)): # input is a scalar field
995            return applyToExp_B_S(e)
996        elif (field.is_Vector(f) and field.is_Vector(g)):
997                return applyToExp_B_V(e)
998        elif (field.is_Matrix(f) and field.is_Matrix(g)):
999            return applyToExp_B_M(e)
1000        elif (field.is_Ten3(f) and field.is_Ten3(g)):
1001            return applyToExp_B_T3(e)
1002        else:
1003             return applyToExp_B_rest(e)
1004
1005    def applyUnaryOnce(oexp_inner,app_inner,app_outer):
1006        #print "applyUnaryOnce"
1007        #apply.toStr(app_inner)
1008        oty_inner = apply.get_oty(app_inner)
1009        oty_outer = apply.get_oty(app_outer)
1010        opr_outer = app_outer.opr
1011        #print "oexp_inner",oexp_inner,"opr_outer",opr_outer.name
1012        lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
1013        app_tmp = apply("tmp", opr_outer, lhs_tmp, None, oty_outer, true, true)
1014        oexp_tmp =unary(app_tmp)
1015        #print " oexp_tmp", oexp_tmp
1016        return (oty_outer, oexp_tmp)
1017
1018    def applyBinaryOnce(oexp_inner,app_inner,app_outer,rhs):
1019        oty_inner = apply.get_oty(app_inner)
1020        oty_outer = apply.get_oty(app_outer)
1021        opr_outer = app_outer.opr
1022
1023        lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
1024
1025        app_tmp = apply("tmp", opr_outer, lhs_tmp, rhs, oty_outer, true,true)
1026        oexp_tmp =binary(app_tmp)
1027        return (oty_outer, oexp_tmp)
1028
1029
1030    # operators with scalar field and vector field
1031    def sort(e):
1032        #apply.toStr(e)
1033        arity = apply.get_arity(e)
1034        if(e.isrootlhs): # is root
1035            #print "sort is a root"
1036            oty = apply.get_oty(e)
1037            if (arity ==1):
1038                return (oty, unary(e))
1039            elif (arity ==2): # binary operator
1040                return (oty, binary(e))
1041            else:
1042                raise Exception ("arity is not supported: "+str(arity_outer))
1043        else:
1044            app_outer = e
1045            arity_outer = arity
1046            #print "app_outer",app_outer.opr.name
1047            if (arity_outer==1):  #assuming both arity
1048                app_inner = apply.get_unary(app_outer)
1049                #print "outer(1) app_inner:",app_inner.opr.name
1050                arity_inner=  app_inner.opr.arity
1051                if (arity_inner==1):
1052                    oexp_inner = unary(app_inner)
1053
1054                    (oty_outer, oexp_tmp) =  applyUnaryOnce(oexp_inner ,app_inner, app_outer)
1055
1056                    return (oty_outer, oexp_tmp)
1057                elif(arity_inner==2):
1058                    oexp_inner = binary(app_inner)
1059                    (oty_outer, oexp_tmp) =  applyUnaryOnce(oexp_inner ,app_inner, app_outer)
1060
1061                    return (oty_outer, oexp_tmp)
1062                else:
1063                    raise Exception ("arity is not supported: "+str(arity_outer))
1064            elif (arity_outer==2):
1065                (app_inner, G) = apply.get_binary(app_outer)
1066                arity_inner=  app_inner.opr.arity
1067                #print "outer(2) app_inner",app_inner.opr.name
1068                if (arity_inner==1):
1069                    oexp_inner = unary(app_inner)
1070                    rhs = G
1071                    (oty_outer, oexp_tmp) =  applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
1072                    return (oty_outer, oexp_tmp)
1073                elif(arity_inner==2):
1074                    oexp_inner = binary(app_inner)
1075                    #print "applying binary frst time", oexp_inner
1076                    rhs = G
1077                    (oty_outer, oexp_tmp) =  applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
1078                    #print "applying binary second time",  oexp_tmp
1079                    return (oty_outer, oexp_tmp)
1080                else:
1081                    raise Exception ("arity is not supported: "+str(arity_outer))
1082            else:
1083                raise Exception ("arity is not supported: "+str(arity_outer))
1084
1085    # ***************************  evaluate at positions ***************************
1086    #evaluate scalar field exp
1087    def eval_d0(pos0, exp):
1088        return exp
1089
1090
1091    def eval_d1(pos0, exp):
1092        #print "eval vec d1"
1093        #print "exp:",exp
1094        #print "pos0",pos0
1095  #evaluate field defined by coefficients at position  #evaluate field defined by coefficients at position
1096  def eval(x,y,a,b,c,d,e,f,g,h,i):      exp = exp.subs(x,pos0)
1097        #print "exp",exp
1098        return exp
1099
1100    def eval_d2(pos0, pos1, exp):
1101        #print "exp:",exp
1102        #evaluate field defined by coefficients at position
1103        exp = exp.subs(x,pos0)
1104        exp = exp.subs(y,pos1)
1105        return exp
1106
1107    def eval_d3(pos0, pos1, pos2, exp):
1108      #evaluate field defined by coefficients at position      #evaluate field defined by coefficients at position
1109      t0 = b*x + c*y      exp = exp.subs(x,pos0)
1110      t1 = d*x*x + g*y*x*x      exp = exp.subs(y,pos1)
1111      t2 = f*y*y + h*x*y*y      exp = exp.subs(z,pos2)
1112      t3 = a + e*x*y + i*(y*y)*(x*x)      return exp
1113      t4= t0+t1+t2+t3
1114      return t4  #evaluate vector field [exp]
1115    def eval_vec_d1(pos0, vec):
1116  def negation(x,y,a,b,c,d,e,f,g,h,i):      #print "eval vec d1"
1117      return (-1*eval(x,y,a,b,c,d,e,f,g,h,i))      rtn = []
1118        for v in vec:
1119  #gradient of field defined by coefficients at position          rtn.append(eval_d1(pos0, v))
1121      x0 = b
1122      x1 = d*2*x + g*y*2*x  #evaluate vector field [exp,exp]
1123      x2 = h*y*y  def eval_vec_d2(pos0, pos1, vec):
1124      x3 = e*y + i*(y*y)*(2*x)      #print "eval_vec_d2 vec:",vec
1125      x4= x0+x1+x2+x3 #xaxis      rtn = []
1126      y0 = c      for v in vec:
1127      y1 = g*x*x          rtn.append(eval_d2(pos0, pos1, v))
1128      y2 = f*2*y + h*x*2*y      return rtn
1129      y3 = e*x + i*(y*2)*(x*x)
1130      y4= y0+y1+y2+y3 #yaxis  def eval_ten3_d1(pos0,ten3):
1131      return [x4,y4]      rtn = []
1132        for i in ten3:
1133            for j in i:
1134                for v in j:
1135                    rtn.append(eval_d1(pos0, v))
1136        return rtn
1137
1138
1139
1140    #evaluate vector field [exp,exp]
1141    def eval_mat_d1(pos0, mat):
1142        #print "eval_vec_d2 vec:",vec
1143        rtn = []
1144        for i in mat:
1145            for v in i:
1146                rtn.append(eval_d1(pos0, v))
1147        return rtn
1148
1149    #evaluate vector field [exp,exp]
1150    def eval_mat_d2(pos0, pos1, mat):
1151        #print "eval_vec_d2 vec:",vec
1152        rtn = []
1153        #print "eval_mat_d2 mat",mat
1154        for i in mat:
1155            for v in i:
1156                rtn.append(eval_d2(pos0, pos1, v))
1157        return rtn
1158
1159    def eval_ten3_d2(pos0, pos1, ten3):
1160        #print "eval_vec_d2 vec:",vec
1161        rtn = []
1162        for i in ten3:
1163            for j in i:
1164                for v in j:
1165                    rtn.append(eval_d2(pos0, pos1, v))
1166        return rtn
1167
1168
1169
1170    #evaluate vector field [exp,exp]
1171    def eval_vec_d3(pos0, pos1, pos2, vec):
1172        rtn = []
1173        for v in vec:
1174            rtn.append(eval_d3(pos0, pos1, pos2, v))
1175        return rtn
1176
1177
1178    #evaluate vector field [exp,exp]
1179    def eval_mat_d3(pos0, pos1, pos2, mat):
1180        #print "eval_vec_d2 vec:",vec
1181        rtn = []
1182        for i in mat:
1183            for v in i:
1184                rtn.append(eval_d3(pos0, pos1, pos2, v))
1185        return rtn
1186
1187    def eval_ten3_d3(pos0, pos1, pos2,ten3):
1188        rtn = []
1189        for i in ten3:
1190
1191            for j in i:
1192
1193  def iter(k, pos, base, xsq, ysq, diag):              for v in j:
1194
1195                    rtn.append(eval_d3(pos0, pos1, pos2, v))
1196        return rtn
1197
1198
1199
1200    def iter_d1(k, pos, exp):
1201        corr = []
1202        for x in pos:
1203            val = k(x, exp)
1204            corr.append(val)
1205        return corr
1206
1207    def iter_d2(k, pos, exp):
1208      corr = []      corr = []
1209      b = base[0]      #print "iter expr:", exp
1210      c = base[1]      #print "pos", pos
d = xsq[0]
g = xsq[1]
f = ysq[0]
h = ysq[1]
a = diag[0]
e = diag[1]
i = diag[2]
1211      for p in pos:      for p in pos:
1212            #print "p", p
1213          x=p[0]          x=p[0]
1214          y=p[1]          y=p[1]
1215          val = k(x,y,a,b,c,d,e,f,g,h,i)          val = k(x,y,exp)
1216          corr.append(val)          corr.append(val)
1217      return corr      return corr
1218
1219    def iter_d3(k, pos, exp):
1220        corr = []
1221        #print "iter exp:", exp
1222        for p in pos:
1223            x=p[0]
1224            y=p[1]
1225            z=p[2]
1226            val = k(x,y,z, exp)
1227            #print "pos: ",x,y,z, " val:", val
1228            corr.append(val)
1229        return corr
1230
1231    def probeField(otyp1,pos, ortn):
1232        dim = fty.get_dim(otyp1)
1233        #print "output type"+otyp1.name
1234        #print "inside probe field ortn", ortn
1235        if (dim==nonefield_dim):
1236            #still need to flatten
1237            rtn = []
1238            if (fty.is_Matrix(otyp1)):
1239                for i in ortn:
1240                    for j in i :
1241                        rtn.append(j)
1242                return [rtn]
1243            elif (fty.is_Ten3(otyp1)):
1244                for i in ortn:
1245                    for j in i :
1246                        for k in  j:
1247                            rtn.append(k)
1248                return [rtn]
1249
def unary(ex, pos,base,xsq,ysq,diag):
if(ex==op_probe): #probing
return iter(eval,pos,base,xsq,ysq,diag)
elif(ex==op_negation): #negation
return iter(negation,pos,base,xsq,ysq,diag)
1250      else:      else:
1251          raise ("no operator value for "+str(ex))              return [ortn]
1252
1253        elif (dim==1):
1254            def get_k():
1255                if (fty.is_ScalarField(otyp1)): # output is a scalar field
1256                    #print "s_d1"
1257                    return eval_d1
1258                elif (fty.is_VectorField(otyp1)):
1259                    #print "v_d1"
1260                    return eval_vec_d1
1261                elif (fty.is_Matrix(otyp1)):
1262                    return eval_mat_d1
1263                elif(fty.is_Ten3(otyp1)):
1264                    return eval_ten3_d1
1265                else:
1266                    raise "error"+otyp1.name
1267            return iter_d1(get_k(), pos, ortn)
1268        elif (dim==2):
1269            def get_k():
1270                if (fty.is_ScalarField(otyp1)): # output is a scalar field
1271                    return eval_d2
1272                elif (fty.is_VectorField(otyp1)):
1273                    return eval_vec_d2
1274                elif (fty.is_Matrix(otyp1)):
1275                    return eval_mat_d2
1276                elif(fty.is_Ten3(otyp1)):
1277                    return eval_ten3_d2
1278                else:
1279                    raise "error"+otyp1.name
1280            return iter_d2(get_k(), pos, ortn)
1281        elif (dim==3):
1282            def get_k():
1283                if (fty.is_ScalarField(otyp1)): # output is a scalar field
1284                    return eval_d3
1285                elif (fty.is_VectorField(otyp1)):
1286                    return eval_vec_d3
1287                elif (fty.is_Matrix(otyp1)):
1288                    return eval_mat_d3
1289                elif(fty.is_Ten3(otyp1)):
1290                    return eval_ten3_d3
1291                else:
1292                    raise "error"+otyp1.name
1293            return iter_d3(get_k(), pos, ortn)
1294        else:
1295            raise "unsupported dimension"
1296
1297    # ***************************  main  ***************************
1298
1299    # operators with scalar field and vector field
1300    def eval(app, pos):
1301        #print "evalname",app.name
1302        #apply.toStr(app)
1303        (otyp1, ortn) = sort(app) #apply operations to expressions
1304        #print "ortn",ortn
1305        rtn = probeField(otyp1, pos, ortn) #evaluate expression at positions
1306        #print "rtn", rtn
1307        return rtn

Legend:
 Removed from v.3865 changed lines Added in v.4436