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

Legend:
Removed from v.3874  
changed lines
  Added in v.4252

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