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 4269, Tue Jul 26 14:28:03 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 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        [n, m] = fty.get_shape(ityp1)
369        rtn = 0
370        for i in range(n):
371                for j in range(m):
372                    rtn+=exp1[i][j]*exp2[i][j]
373        return rtn
374    
375    
376    
377    #evaluate det
378    def fn_det(fld):
379        exp = field.get_data(fld)
380        ityp = field.get_ty(fld)
381        rtn=[]
382        if(field.is_Matrix(fld)):
383            [n, m] = fty.get_shape(ityp)
384            if (n!=m):
385                raise Exception("matrix is not identitical")
386            a = exp[0][0]
387            d = exp[1][1]
388            c = exp[1][0]
389            b = exp[0][1]
390            if(n==2):
391                return a*d-b*c
392            elif(n==3):
393                a = exp[0][0]
394                b = exp[0][1]
395                c = exp[0][2]
396                d = exp[1][0]
397                e = exp[1][1]
398                f = exp[1][2]
399                g = exp[2][0]
400                h = exp[2][1]
401                i = exp[2][2]
402                return a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)
403            else:
404                raise "unsupported matrix field"
405        else:
406            raise "unsupported trace"
407    
408    
409    #evaluate outer product
410    def fn_outer(fld1, fld2):
411        exp1 = field.get_data(fld1)
412        ityp1 = field.get_ty(fld1)
413        exp2 = field.get_data(fld2)
414        ityp2 = field.get_ty(fld2)
415        ashape = fty.get_shape(ityp1)
416        bshape = fty.get_shape(ityp2)
417        x= "ashape", ashape, "bshape", bshape
418        #print "exp1",exp1,"ityp1",ityp1.name,"-length",len(exp1)
419        #print "exp2",exp2,"ityp2",ityp2.name,"-length",len(exp2)
420        rtn = []
421        if(fty.is_Vector(ityp1)):
422            [n1] = fty.get_shape(ityp1)
423            if(fty.is_Vector(ityp2)):
424                #both vectors
425                [n2] = fty.get_shape(ityp2)
426                for i in  range(n1):
427                    tmpI = []
428                    for j in range(n2):
429                        tmpI.append(exp1[i]*exp2[j])
430                    rtn.append(tmpI)
431                return rtn
432            elif(fty.is_Matrix(ityp2)):
433                [n2,n3] = fty.get_shape(ityp2)
434                for i in  range(n1):
435                    tmpI = []
436                    for j in range(n2):
437                        tmpJ = []
438                        for k in range(n3):
439                            tmpJ.append(exp1[i]*exp2[j][k])
440                        tmpI.append(tmpJ)
441                    rtn.append(tmpI)
442                return rtn
443        elif(fty.is_Matrix(ityp1)):
444            [n1,n2] = fty.get_shape(ityp1)
445            if(fty.is_Vector(ityp2)):
446                [n3] = fty.get_shape(ityp2)
447                for i in  range(n1):
448                    tmpI = []
449                    for j in range(n2):
450                        tmpJ = []
451                        for k in range(n3):
452                            tmpJ.append(exp1[i][j]*exp2[k])
453                        tmpI.append(tmpJ)
454                    rtn.append(tmpI)
455                return rtn
456            elif(fty.is_Matrix(ityp2)):
457                [n3, n4] = 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                            tmpK = []
464                            for l in range(n4):
465                                tmpK.append(exp1[i][j]*exp2[k][l])
466                            tmpJ.append(tmpK)
467                        tmpI.append(tmpJ)
468                    rtn.append(tmpI)
469                return rtn
470            else:
471                raise "outer product is not supported"
472        else:
473            raise "outer product is not supported"
474    #evaluate inner product
475    def fn_inner(fld1, fld2):
476        exp1 = field.get_data(fld1)
477        ityp1 = field.get_ty(fld1)
478        exp2 = field.get_data(fld2)
479        ityp2 = field.get_ty(fld2)
480        ashape = fty.get_shape(ityp1)
481        bshape = fty.get_shape(ityp2)
482        x= "ashape", ashape, "bshape", bshape
483        if(fty.is_Vector(ityp1)):
484            [a1] = fty.get_shape(ityp1)
485            if(fty.is_Vector(ityp2)):
486                #length of vetors
487                rtn=0
488                [b1] = fty.get_shape(ityp2)
489                if(a1!=b1):
490                    raise x
491                for s in  range(a1):
492                    rtn += exp1[s]*exp2[s]
493                return rtn
494            elif(fty.is_Matrix(ityp2)):
495                [b1,b2] = fty.get_shape(ityp2)
496                rtn=[]
497                if(a1!=b1):
498                    raise x
499                for i in  range(b2):
500                    sumrtn=0
501                    for s in  range(a1):
502                        sumrtn +=  exp1[s]*exp2[s][i]
503                    rtn.append(sumrtn)
504                return rtn
505            elif(fty.is_Ten3(ityp2)):
506                [b1, b2, b3] = fty.get_shape(ityp2)
507                rtn = []
508                if(a1!=b1):
509                    raise x
510                for i in  range(b2):
511                    tmpJ = []
512                    for j in  range(b3):
513                        sumrtn=0
514                        for s in  range(a1):
515                            sumrtn +=  exp1[s]*exp2[s][i][j]
516                        tmpJ.append(sumrtn)
517                    rtn.append(tmpJ)
518                return rtn
519            else:
520                raise "inner product is not supported"
521        elif(fty.is_Matrix(ityp1)):
522    
523            [a1,a2] = fty.get_shape(ityp1)
524    
525            if(fty.is_Vector(ityp2)):
526                [b1] = fty.get_shape(ityp2)
527                if(a2!=b1):
528                    raise x
529                rtn=[]
530                for i in  range(a1):
531                    sumrtn=0
532                    for s in  range(a2):
533                        sumrtn += exp1[i][s]*exp2[s]
534                    rtn.append(sumrtn)
535                return rtn
536            if(fty.is_Matrix(ityp2)):
537                [b1,b2] = fty.get_shape(ityp2)
538                rtn=[]
539                if(a2!=b1):
540                    raise x
541                for i in  range(a1):
542                    rtnj = []
543                    for j in  range(b2):
544                        sumrtn=0
545                        for s in  range(a2):
546                            sumrtn += exp1[i][s]*exp2[s][j]
547                        rtnj.append(sumrtn)
548                    rtn.append(rtnj)
549                return rtn
550            else:
551                raise "inner product is not supported"
552        elif(fty.is_Ten3(ityp1)):
553            [a1,a2, a3] = ashape
554            if(fty.is_Vector(ityp2)):
555                [b1] = bshape
556                if(a3!=b1):
557                    raise x
558                rtn=[]
559                for i in  range(a1):
560                    tmpI=[]
561                    for j in  range(a2):
562                        sumrtn=0
563                        for s in  range(a3):
564                            sumrtn += exp1[i][j][s]*exp2[s]
565                        tmpI.append(sumrtn)
566                    rtn.append(tmpI)
567                return rtn
568            if(fty.is_Matrix(ityp2)):
569                [b1,b2] = bshape
570                if(a3!=b1):
571                    raise x
572                rtn=[]
573                for i in  range(a1):
574                    tmpI=[]
575                    for j in  range(a2):
576                        tmpJ = []
577                        for k in range(b2):
578                            sumrtn=0
579                            for s in  range(a3):
580                                sumrtn += exp1[i][j][s]*exp2[s][k]
581                            tmpJ.append(sumrtn)
582                        tmpI.append(tmpJ)
583                    rtn.append(tmpI)
584                return rtn
585            else:
586                raise "inner product is not supported"
587        else:
588            raise "inner product is not supported"
589    
590    
591    # ***************************  generic apply operators ***************************
592    #unary operator on a vector
593    def applyToVector(vec, unary):
594        rtn = []
595        for v in vec:
596            rtn.append(unary(v))
597        return rtn
598    #binary operator on a vector
599    def applyToVectors(vecA, vecB,  binary):
600        rtn = []
601        for (a,b) in zip(vecA,vecB):
602            x= binary(a,b)
603            rtn.append(x)
604        return rtn
605    
606    def applyToM(vec, unary):
607        rtn = []
608        for i in vec:
609            tmpI = []
610            for v in i:
611                tmpI.append(unary(v))
612            rtn.append(tmpI)
613        return rtn
614    
615    def applyToMs(vecA,vecB, unary):
616        rtn = []
617        for (a,b) in zip(vecA,vecB):
618            tmpI = []
619            for (u,v) in zip(a,b):
620                tmpI.append(unary(u, v))
621            rtn.append(tmpI)
622        return rtn
623    
624    
625    def applyToT3(vec, unary):
626        rtn = []
627        for i in vec:
628            tmpI = []
629            for j in i:
630                tmpJ = []
631                for v in j:
632                    tmpJ.append(unary(v))
633                tmpI.append(tmpJ)
634            rtn.append(tmpI)
635        return rtn
636    
637    # ***************************  apply to scalars or vectors ***************************
638    #apply operators to expression
639    # return output types and expression
640    # unary operator
641    # exp: scalar types
642    
643    def applyToExp_U_S(fn_name, fld):
644        exp = field.get_data(fld)
645        dim = field.get_dim(fld)
646        #print fn_name
647        if(op_probe==fn_name): #probing
648            return  exp
649        elif(op_negation==fn_name): #negation
650            return fn_negation(exp)
651        elif(op_cosine==fn_name): #cosine
652            return cos(exp)
653        elif(op_sine==fn_name): #sine
654            return sin(exp)
655        elif(op_tangent==fn_name): # tangent
656            return tan(exp)
657        elif(op_asine==fn_name): #asine
658            return asin(exp)
659        elif(op_acosine==fn_name): #acos
660            return acos(exp)
661        elif(op_atangent==fn_name): #atan
662            return atan(exp)
663        elif(op_sqrt==fn_name): #sqrt
664            return sqrt(exp)
665        elif(op_gradient==fn_name):
666            return fn_grad(exp, dim)
667        else:
668            raise Exception("unsupported unary operator on scalar field:"+ fn_name.name)
669    
670    # unary operator
671    # exp: vector  types
672    def applyToExp_U_V(fn_name, fld):
673        exp = field.get_data(fld)
674        if(op_probe==fn_name): #probing
675            return exp
676        elif(op_negation==fn_name): #negation
677            return applyToVector(exp, fn_negation)
678        elif(op_divergence==fn_name):
679            return fn_divergence(fld)
680        elif(op_curl==fn_name):
681            return fn_curl(fld)
682        elif(op_jacob==fn_name): #jacob
683            return fn_jacob(fld)
684        elif(op_slicev0==fn_name) :
685            return fn_slicev0(fld)
686        elif(op_slicev1==fn_name):
687            return fn_slicev1(fld)
688        elif(op_grad==op_grad):
689            return fn_grad(fld)
690        else:
691            raise Exception("unsupported unary operator:"+ fn_name.name)
692    
693    def applyToExp_U_M(fn_name, fld):
694        exp = field.get_data(fld)
695        if(op_probe==fn_name): #probing
696            return exp
697        elif(op_negation==fn_name): #negation
698            return applyToM(exp, fn_negation)
699        elif(op_jacob==fn_name): #jacob
700            return fn_jacob(fld)
701        elif(op_slicem0==fn_name) :
702            return fn_slicem0(fld)
703        elif(op_slicem1==fn_name):
704            return fn_slicem1(fld)
705        elif(op_trace == fn_name):
706            return fn_trace(fld)
707        elif(op_transpose==fn_name):
708            return fn_transpose(fld)
709        elif(op_det==fn_name):
710            return fn_det(fld)
711        else:
712            raise Exception("unsupported unary operator:"+ fn_name.name)
713    
714    def applyToExp_U_T3(fn_name, fld):
715        exp = field.get_data(fld)
716        if(op_probe==fn_name): #probing
717            return exp
718        elif(op_negation==fn_name): #negation
719            return applyToT3(exp, fn_negation)
720        elif(op_jacob==fn_name): #jacob
721            return fn_jacob(fld)
722        elif(op_slicet0==fn_name) :
723            return fn_slicet0(fld)
724        elif(op_slicet1==fn_name):
725            return fn_slicet1(fld)
726        else:
727            raise Exception("unsupported unary operator:"+ fn_name.name)
728    
729    # binary operator
730    # exp: scalar types
731    def applyToExp_B_S(e):
732        fn_name=e.opr
733        print "b.s"
734        (fld1,fld2) =  apply.get_binary(e)
735        exp1 = field.get_data(fld1)
736        exp2 = field.get_data(fld2)
737        #print fn_name
738        if(op_add==fn_name):#addition
739            return fn_add(exp1,exp2)
740        elif(op_subtract==fn_name):#subtract
741            return fn_subtract(exp1,exp2)
742        elif(op_modulate==fn_name):#modulate
743            return fn_modulate(exp1,exp2)
744        elif(op_scale==fn_name): #scaling
745            return fn_scaling(fld1,fld2)
746        elif(op_division==fn_name): #division
747            return fn_division(fld1,fld2)
748        else:
749            raise Exception("unsupported binary operator on scalar fields:"+ fn_name.name)
750    
751    # binary, args do not need to have the same shape
752    def applyToExp_B_uneven(e):
753        fn_name=e.opr
754        (fld1,fld2) =  apply.get_binary(e)
755        exp1 = field.get_data(fld1)
756        exp2 = field.get_data(fld2)
757        if(op_outer==fn_name):
758            return fn_outer(fld1, fld2)
759        elif(op_inner==fn_name):
760            return fn_inner(fld1, fld2)
761        elif(op_scale==fn_name): #scaling
762            return fn_scaling(fld1,fld2)
763        else:
764            raise Exception("unsupported unary operator:"+fn_name.name)
765    
766    
767    # binary operator
768    # args have the same shape
769    def applyToExp_B_V(e):
770        print "b.v"
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_add==fn_name):#addition
776            return applyToVectors(exp1, exp2,  fn_add)
777        elif(op_subtract==fn_name):#subtract
778            return  applyToVectors(exp1, exp2, fn_subtract)
779        elif(op_modulate==fn_name):#modulate
780            return applyToVectors(exp1,exp2 ,fn_modulate)
781        elif(op_cross==fn_name):
782            return fn_cross(fld1, fld2)
783        else:
784           return applyToExp_B_uneven(e)
785    
786    def applyToExp_B_M(e):
787        print "b. m"
788        fn_name=e.opr
789        (fld1,fld2) =  apply.get_binary(e)
790        exp1 = field.get_data(fld1)
791        exp2 = field.get_data(fld2)
792        if(op_add==fn_name):#addition
793            return applyToMs(exp1, exp2,  fn_add)
794        elif(op_subtract==fn_name):#subtract
795            return  applyToMs(exp1, exp2, fn_subtract)
796        elif(op_modulate==fn_name):#modulate
797            return applyToMs(exp1,exp2,fn_modulate)
798        elif(op_doubledot==fn_name):#op_doubledot
799            return fn_doubledot (fld1,fld2)
800        else:
801            return applyToExp_B_uneven(e)
802    
803    
804    # ***************************  unary/binary operators ***************************
805    def unary(e):
806        #apply.toStr(e)
807        fld =apply.get_unary(e)
808        fn_name=e.opr
809        exp = field.get_data(fld)
810        dim = field.get_dim(fld)
811        if(op_norm==fn_name):#norm
812            return fn_norm(fld, dim)
813        if(op_normalize==fn_name):#normalize
814            x= fn_normalize(fld, dim)
815            return x
816        elif (field.is_Scalar(fld)): # input is a scalar field
817            return applyToExp_U_S(fn_name, fld)
818        elif(field.is_Vector(fld)): # input is a vector field
819            return applyToExp_U_V(fn_name, fld)
820        elif(field.is_Matrix(fld)): # input is a vector field
821            return applyToExp_U_M(fn_name, fld)
822        else:
823            return applyToExp_U_T3(fn_name, fld)
824    
825    def binary(e):
826        (f, g) =apply.get_binary(e)
827        fn_name = e.opr
828        # type is checked elsewhere or does not matter
829        if(op_division==fn_name): #division
830            return fn_division(f, g)
831        elif (field.is_Scalar(f) and field.is_Scalar(g)): # input is a scalar field
832            return applyToExp_B_S(e)
833        elif (field.is_Vector(f)):# input is a vector field
834            if(field.is_Vector(g)):
835                return applyToExp_B_V(e)
836            else: # input is a vector field, _
837                return applyToExp_B_V(e)
838        elif (field.is_Matrix(f)):# input is a matrix field
839            if(field.is_Matrix(g)):
840                return applyToExp_B_M(e)
841            else:
842                return  applyToExp_B_V(e)
843        else:
844             return applyToExp_B_V(e)
845    
846    def applyUnaryOnce(oexp_inner,app_inner,app_outer):
847        #print "applyUnaryOnce"
848        #apply.toStr(app_inner)
849        oty_inner = apply.get_oty(app_inner)
850        oty_outer = apply.get_oty(app_outer)
851        opr_outer = app_outer.opr
852        #print "oexp_inner",oexp_inner,"opr_outer",opr_outer.name
853        lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
854        app_tmp = apply("tmp", opr_outer, lhs_tmp, None, oty_outer, true, true)
855        oexp_tmp =unary(app_tmp)
856        #print " oexp_tmp", oexp_tmp
857        return (oty_outer, oexp_tmp)
858    
859    def applyBinaryOnce(oexp_inner,app_inner,app_outer,rhs):
860        oty_inner = apply.get_oty(app_inner)
861        oty_outer = apply.get_oty(app_outer)
862        opr_outer = app_outer.opr
863    
864        lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
865    
866        app_tmp = apply("tmp", opr_outer, lhs_tmp, rhs, oty_outer, true,true)
867        oexp_tmp =binary(app_tmp)
868        return (oty_outer, oexp_tmp)
869    
870    
871    # operators with scalar field and vector field
872    def sort(e):
873        #apply.toStr(e)
874        arity = apply.get_arity(e)
875        if(e.isrootlhs): # is root
876            #print "sort is a root"
877            oty = apply.get_oty(e)
878            if (arity ==1):
879                return (oty, unary(e))
880            elif (arity ==2): # binary operator
881                return (oty, binary(e))
882            else:
883                raise Exception ("arity is not supported: "+str(arity_outer))
884        else:
885            app_outer = e
886            arity_outer = arity
887            #print "app_outer",app_outer.opr.name
888            if (arity_outer==1):  #assuming both arity
889                app_inner = apply.get_unary(app_outer)
890                #print "outer(1) app_inner:",app_inner.opr.name
891                arity_inner=  app_inner.opr.arity
892                if (arity_inner==1):
893                    oexp_inner = unary(app_inner)
894    
895                    (oty_outer, oexp_tmp) =  applyUnaryOnce(oexp_inner ,app_inner, app_outer)
896    
897                    return (oty_outer, oexp_tmp)
898                elif(arity_inner==2):
899                    oexp_inner = binary(app_inner)
900                    (oty_outer, oexp_tmp) =  applyUnaryOnce(oexp_inner ,app_inner, app_outer)
901    
902                    return (oty_outer, oexp_tmp)
903                else:
904                    raise Exception ("arity is not supported: "+str(arity_outer))
905            elif (arity_outer==2):
906                (app_inner, G) = apply.get_binary(app_outer)
907                arity_inner=  app_inner.opr.arity
908                #print "outer(2) app_inner",app_inner.opr.name
909                if (arity_inner==1):
910                    oexp_inner = unary(app_inner)
911                    rhs = G
912                    (oty_outer, oexp_tmp) =  applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
913                    return (oty_outer, oexp_tmp)
914                elif(arity_inner==2):
915                    oexp_inner = binary(app_inner)
916                    #print "applying binary frst time", oexp_inner
917                    rhs = G
918                    (oty_outer, oexp_tmp) =  applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
919                    #print "applying binary second time",  oexp_tmp
920                    return (oty_outer, oexp_tmp)
921                else:
922                    raise Exception ("arity is not supported: "+str(arity_outer))
923            else:
924                raise Exception ("arity is not supported: "+str(arity_outer))
925    
926    # ***************************  evaluate at positions ***************************
927    #evaluate scalar field exp
928    def eval_d1(pos0, exp):
929        #print "eval vec d1"
930        #print "exp:",exp
931        #print "pos0",pos0
932  #evaluate field defined by coefficients at position  #evaluate field defined by coefficients at position
933  def eval(x,y,a,b,c,d,e,f,g,h,i):      exp = exp.subs(x,pos0)
934        #print "exp",exp
935        return exp
936    
937    def eval_d2(pos0, pos1, exp):
938        #print "exp:",exp
939      #evaluate field defined by coefficients at position      #evaluate field defined by coefficients at position
940      t0 = b*x + c*y      exp = exp.subs(x,pos0)
941      t1 = d*x*x + g*y*x*x      exp = exp.subs(y,pos1)
942      t2 = f*y*y + h*x*y*y      return exp
943      t3 = a + e*x*y + i*(y*y)*(x*x)  
944      t4= t0+t1+t2+t3  def eval_d3(pos0, pos1, pos2, exp):
945      return t4      #evaluate field defined by coefficients at position
946        exp = exp.subs(x,pos0)
947  def negation(x,y,a,b,c,d,e,f,g,h,i):      exp = exp.subs(y,pos1)
948      return (-1*eval(x,y,a,b,c,d,e,f,g,h,i))      exp = exp.subs(z,pos2)
949        return exp
950  #gradient of field defined by coefficients at position  
951  def grad(x,y,a,b,c,d,e,f,g,h,i):  #evaluate vector field [exp]
952      x0 = b  def eval_vec_d1(pos0, vec):
953      x1 = d*2*x + g*y*2*x      #print "eval vec d1"
954      x2 = h*y*y      rtn = []
955      x3 = e*y + i*(y*y)*(2*x)      for v in vec:
956      x4= x0+x1+x2+x3 #xaxis          rtn.append(eval_d1(pos0, v))
957      y0 = c      return rtn
958      y1 = g*x*x  
959      y2 = f*2*y + h*x*2*y  #evaluate vector field [exp,exp]
960      y3 = e*x + i*(y*2)*(x*x)  def eval_vec_d2(pos0, pos1, vec):
961      y4= y0+y1+y2+y3 #yaxis      #print "eval_vec_d2 vec:",vec
962      return [x4,y4]      rtn = []
963        for v in vec:
964            rtn.append(eval_d2(pos0, pos1, v))
965        return rtn
966    
967    def eval_ten3_d1(pos0,ten3):
968        rtn = []
969        for i in ten3:
970            for j in i:
971                for v in j:
972                    rtn.append(eval_d1(pos0, v))
973        return rtn
974    
975    
976    
977    #evaluate vector field [exp,exp]
978    def eval_mat_d1(pos0, mat):
979        #print "eval_vec_d2 vec:",vec
980        rtn = []
981        for i in mat:
982            for v in i:
983                rtn.append(eval_d1(pos0, v))
984        return rtn
985    
986    #evaluate vector field [exp,exp]
987    def eval_mat_d2(pos0, pos1, mat):
988        #print "eval_vec_d2 vec:",vec
989        rtn = []
990        #print "eval_mat_d2 mat",mat
991        for i in mat:
992            for v in i:
993                rtn.append(eval_d2(pos0, pos1, v))
994        return rtn
995    
996    def eval_ten3_d2(pos0, pos1, ten3):
997        #print "eval_vec_d2 vec:",vec
998        rtn = []
999        for i in ten3:
1000            for j in i:
1001                for v in j:
1002                    rtn.append(eval_d2(pos0, pos1, v))
1003        return rtn
1004    
1005    
1006    
1007    #evaluate vector field [exp,exp]
1008    def eval_vec_d3(pos0, pos1, pos2, vec):
1009        rtn = []
1010        for v in vec:
1011            rtn.append(eval_d3(pos0, pos1, pos2, v))
1012        return rtn
1013    
1014    
1015    #evaluate vector field [exp,exp]
1016    def eval_mat_d3(pos0, pos1, pos2, mat):
1017        #print "eval_vec_d2 vec:",vec
1018        rtn = []
1019        for i in mat:
1020            for v in i:
1021                rtn.append(eval_d3(pos0, pos1, pos2, v))
1022        return rtn
1023    
1024  def iter(k, pos, base, xsq, ysq, diag):  def eval_ten3_d3(pos0, pos1, pos2,ten3):
1025        rtn = []
1026        for i in ten3:
1027            for j in i:
1028                for v in j:
1029                    rtn.append(eval_d3(pos0, pos1, pos2, v))
1030        return rtn
1031    
1032    
1033    
1034    def iter_d1(k, pos, exp):
1035      corr = []      corr = []
1036      b = base[0]      for x in pos:
1037      c = base[1]          val = k(x, exp)
1038      d = xsq[0]          corr.append(val)
1039      g = xsq[1]      return corr
1040      f = ysq[0]  
1041      h = ysq[1]  def iter_d2(k, pos, exp):
1042      a = diag[0]      corr = []
1043      e = diag[1]      #print "iter expr:", exp
1044      i = diag[2]      #print "pos", pos
1045      for p in pos:      for p in pos:
1046            #print "p", p
1047          x=p[0]          x=p[0]
1048          y=p[1]          y=p[1]
1049          val = k(x,y,a,b,c,d,e,f,g,h,i)          val = k(x,y,exp)
1050          corr.append(val)          corr.append(val)
1051      return corr      return corr
1052    
1053    def iter_d3(k, pos, exp):
1054        corr = []
1055        #print "iter exp:", exp
1056        for p in pos:
1057            x=p[0]
1058            y=p[1]
1059            z=p[2]
1060            val = k(x,y,z, exp)
1061            #print "pos: ",x,y,z, " val:", val
1062            corr.append(val)
1063        return corr
1064    
1065  def unary(ex, pos,base,xsq,ysq,diag):  def probeField(otyp1,pos, ortn):
1066      if(ex==op_probe): #probing      dim = fty.get_dim(otyp1)
1067          return iter(eval,pos,base,xsq,ysq,diag)      #print "output type"+otyp1.name
1068      elif(ex==op_negation): #negation      if (dim==1):
1069          return iter(negation,pos,base,xsq,ysq,diag)          def get_k():
1070      elif(ex==op_gradient): #gradient              if (fty.is_ScalarField(otyp1)): # output is a scalar field
1071          return iter(grad,pos,base,xsq,ysq,diag)                  #print "s_d1"
1072                    return eval_d1
1073                elif (fty.is_VectorField(otyp1)):
1074                    #print "v_d1"
1075                    return eval_vec_d1
1076                elif (fty.is_Matrix(otyp1)):
1077                    return eval_mat_d1
1078                elif(fty.is_Ten3(otyp1)):
1079                    return eval_ten3_d1
1080      else:      else:
1081          raise ("no operator value for "+str(ex))                  raise "error"+otyp1.name
1082            return iter_d1(get_k(), pos, ortn)
1083        elif (dim==2):
1084            def get_k():
1085                if (fty.is_ScalarField(otyp1)): # output is a scalar field
1086                    return eval_d2
1087                elif (fty.is_VectorField(otyp1)):
1088                    return eval_vec_d2
1089                elif (fty.is_Matrix(otyp1)):
1090                    return eval_mat_d2
1091                elif(fty.is_Ten3(otyp1)):
1092                    return eval_ten3_d2
1093                else:
1094                    raise "error"+otyp1.name
1095            return iter_d2(get_k(), pos, ortn)
1096        elif (dim==3):
1097            def get_k():
1098                if (fty.is_ScalarField(otyp1)): # output is a scalar field
1099                    return eval_d3
1100                elif (fty.is_VectorField(otyp1)):
1101                    return eval_vec_d3
1102                elif (fty.is_Matrix(otyp1)):
1103                    return eval_mat_d3
1104                elif(fty.is_Ten3(otyp1)):
1105                    return eval_ten3_d3
1106                else:
1107                    raise "error"+otyp1.name
1108            return iter_d3(get_k(), pos, ortn)
1109        else:
1110            raise "unsupported dimension"
1111    
1112    # ***************************  main  ***************************
1113    
1114    # operators with scalar field and vector field
1115    def eval(app, pos):
1116        #print "evalname",app.name
1117        #apply.toStr(app)
1118        (otyp1, ortn) = sort(app) #apply operations to expressions
1119        #print "ortn",ortn
1120        rtn = probeField(otyp1, pos, ortn) #evaluate expression at positions
1121        #print "rtn", rtn
1122        return rtn

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

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