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

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

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