Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0