Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

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

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