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

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

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