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 4188, Mon Jul 11 23:21:11 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    # negation  of field
79    def fn_negation(exp):
80
81        return -1*exp
82
83
84    #evaluate cross product
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
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  = (sum)**0.5
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    #evaluate slice
254    def fn_slice(fld1):
255        exp1 = field.get_data(fld1)
256        ityp1 = field.get_ty(fld1)
257        rtn=[]
258        if(fty.is_Matrix(ityp1)):
259            [n2,n3] = fty.get_shape(ityp1)
260            for j in range(n3):
261                rtn.append(exp1[j][0])
262            return rtn
263        else:
264            raise "unsupported type for slice"
265
266    #evaluate trace
267    def fn_trace(fld):
268        exp = field.get_data(fld)
269        ityp = field.get_ty(fld)
270        rtn=[]
271        if(field.is_Matrix(fld)):
272            [n, m] = fty.get_shape(ityp)
273            if (n!=m):
274                raise Exception("matrix is not identitical")
275            rtn = exp[0][0]+exp[1][1]
276            if(n==2):
277                return rtn
278            elif(n==3):
279                return rtn+exp[2][2]
280            else:
281                raise "unsupported matrix field"
282        else:
283            raise "unsupported trace"
284
285    #evaluate transpose
286    def fn_transpose(fld):
287        exp = field.get_data(fld)
288        ityp = field.get_ty(fld)
289        if(field.is_Matrix(fld)):
290            [n, m] = fty.get_shape(ityp)
291            rtn = []
292            for i in range(n):
293                rtni = []
294                for j in range(m):
295                    rtni.append(exp[j][i])
296                rtn.append(rtni)
297            return rtn
298        else:
299            raise "unsupported transpose"
300
301    #evaluate det
302    def fn_det(fld):
303        exp = field.get_data(fld)
304        ityp = field.get_ty(fld)
305        rtn=[]
306        if(field.is_Matrix(fld)):
307            [n, m] = fty.get_shape(ityp)
308            if (n!=m):
309                raise Exception("matrix is not identitical")
310            a = exp[0][0]
311            d = exp[1][1]
312            c = exp[1][0]
313            b = exp[0][1]
314            if(n==2):
315                return a*d-b*c
316            elif(n==3):
317                a = exp[0][0]
318                b = exp[0][1]
319                c = exp[0][2]
320                d = exp[1][0]
321                e = exp[1][1]
322                f = exp[1][2]
323                g = exp[2][0]
324                h = exp[2][1]
325                i = exp[2][2]
326                return a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)
327            else:
328                raise "unsupported matrix field"
329        else:
330            raise "unsupported trace"
331
332
333    #evaluate outer product
334    def fn_outer(fld1, fld2):
335        exp1 = field.get_data(fld1)
336        ityp1 = field.get_ty(fld1)
337        exp2 = field.get_data(fld2)
338        ityp2 = field.get_ty(fld2)
339        rtn=[]
340        #print "exp1",exp1,"ityp1",ityp1.name,"-length",len(exp1)
341        #print "exp2",exp2,"ityp2",ityp2.name,"-length",len(exp2)
342
343        if(fty.is_Vector(ityp1)):
344            n1= fty.get_vecLength(ityp1)
345            if(fty.is_Vector(ityp2)):
346                #both vectors
347                n2= fty.get_vecLength(ityp2)
348                for i in  range(n1):
349                    tmpI = []
350                    for j in range(n2):
351                        tmpI.append(exp1[i]*exp2[j])
352                    rtn.append(tmpI)
353                return rtn
354            elif(fty.is_Matrix(ityp2)):
355                [n2,n3] = fty.get_shape(ityp2)
356                for i in  range(n1):
357                    tmpI = []
358                    for j in range(n2):
359                        tmpJ = []
360                        for k in range(n3):
361                            tmpJ.append(exp1[i]*exp2[j][k])
362                        tmpI.append(tmpJ)
363                    rtn.append(tmpI)
364                return rtn
365        elif(fty.is_Matrix(ityp1)):
366            [n1,n2] = fty.get_shape(ityp1)
367            if(fty.is_Vector(ityp2)):
368                n3= fty.get_vecLength(ityp2)
369                for i in  range(n1):
370                    tmpI = []
371                    for j in range(n2):
372                        tmpJ = []
373                        for k in range(n3):
374                            tmpJ.append(exp1[i][j]*exp2[k])
375                        tmpI.append(tmpJ)
376                    rtn.append(tmpI)
377                return rtn
378            elif(fty.is_Matrix(ityp2)):
379                [n3,n4] = fty.get_shape(ityp2)
380                for i in  range(n1):
381                    tmpI = []
382                    for j in range(n2):
383                        tmpJ = []
384                        for k in range(n3):
385                            tmpK = []
386                            for l in range(n4):
387                                tmpK.append(exp1[i][j]*exp2[k][l])
388                            tmpJ.append(tmpK)
389                        tmpI.append(tmpJ)
390                    rtn.append(tmpI)
391                return rtn
392            else:
393                raise "outer product is not supported"
394        else:
395            raise "outer product is not supported"
396    #evaluate inner product
397    def fn_inner(fld1, fld2):
398        exp1 = field.get_data(fld1)
399        ityp1 = field.get_ty(fld1)
400        exp2 = field.get_data(fld2)
401        ityp2 = field.get_ty(fld2)
402        #print " exp1: ",exp1," exp2: ",exp2
403        # vectors
404        if(fty.is_Vector(ityp1)):
405            n1 = fty.get_vecLength(ityp1) #vector
406            if(fty.is_Vector(ityp2)):
407                #length of vetors
408                rtn=0
409                n2 = fty.get_vecLength(ityp2)
410                for s in  range(n1):
411                    curr = exp1[s]*exp2[s]
412                    #print (" exp1[s]: ",exp1[s]," exp2[s]: ",exp2[s],"cur",curr)
413                    rtn += curr
414                return rtn
415            elif(fty.is_Matrix(ityp2)):
416                [n2] = fty.drop_last(ityp2)  #matrix
417                rtn=[]
418                for i in  range(n2):
419                    sumrtn=0
420                    for s in  range(n1):
421                        curr = exp1[s]*exp2[s][i]
422                        sumrtn += curr
423                    rtn.append(sumrtn)
424                return rtn
425            elif(fty.is_Ten3(ityp2)):
426                [n2,n3] = fty.drop_last(ityp2)
427                rtn = []
428                for i in  range(n2):
429                    tmpJ = []
430                    for j in  range(n3):
431                        sumrtn=0
432                        for s in  range(n1):
433                            curr = exp1[s]*exp2[s][i][j]
434                            sumrtn += curr
435                        tmpJ.append(sumrtn)
436                    rtn.append(tmpJ)
437                return rtn
438            else:
439                raise "inner product is not supported"
440        elif(fty.is_Matrix(ityp1)):
441            n2 = fty.get_first_ix(ityp1)  #matrix
442            if(fty.is_Vector(ityp2)):
443                ns = fty.get_vecLength(ityp2) #vector
444                rtn=[]
445                for i in  range(n2):
446                    sumrtn=0
447                    for s in  range(ns):
448                        curr = exp1[i][s]*exp2[s]
449                        sumrtn += curr
450                    rtn.append(sumrtn)
451                return rtn
452            else:
453                raise "inner product is not supported"
454        elif(fty.is_Ten3(ityp1)):
455            [n1,n2] = fty.drop_first(ityp1)
456            if(fty.is_Vector(ityp2)):
457                ns = fty.get_vecLength(ityp2)
458                rtn=[]
459                for i in  range(n1):
460                    tmpI=[]
461                    for j in  range(n2):
462                        sumrtn=0
463                        for s in  range(ns):
464                            curr = exp1[i][j][s]*exp2[s]
465                            sumrtn += curr
466                        tmpI.append(sumrtn)
467                    rtn.append(tmpI)
468                return rtn
469            else:
470                raise "inner product is not supported"
471        else:
472            raise "inner product is not supported"
473
474
475    # ***************************  generic apply operators ***************************
476    #unary operator on a vector
477    def applyToVector(vec, unary):
478        rtn = []
479        for v in vec:
480            rtn.append(unary(v))
481        return rtn
482    #binary operator on a vector
483    def applyToVectors(vecA, vecB,  binary):
484        rtn = []
485        for (a,b) in zip(vecA,vecB):
486            x= binary(a,b)
487            rtn.append(x)
488        return rtn
489
490    def applyToM(vec, unary):
491        rtn = []
492        for i in vec:
493            tmpI = []
494            for v in i:
495                tmpI.append(unary(v))
496            rtn.append(tmpI)
497        return rtn
498
499    def applyToT3(vec, unary):
500        rtn = []
501        for i in vec:
502            tmpI = []
503            for j in i:
504                tmpJ = []
505                for v in j:
506                    tmpJ.append(unary(v))
507                tmpI.append(tmpJ)
508            rtn.append(tmpI)
509        return rtn
510
511    # ***************************  apply to scalars or vectors ***************************
512    #apply operators to expression
513    # return output types and expression
514    # unary operator
515    # exp: scalar types
516
517    def applyToExp_U_S(fn_name, fld):
518        exp = field.get_data(fld)
519        dim = field.get_dim(fld)
520        #print fn_name
521        if(op_probe==fn_name): #probing
522            return  exp
523        elif(op_negation==fn_name): #negation
524            return fn_negation(exp)
527        else:
528            raise Exception("unsupported unary operator on scalar field:"+ fn_name.name)
529
530    # unary operator
531    # exp: vector  types
532    def applyToExp_U_V(fn_name, fld):
533        exp = field.get_data(fld)
534        if(op_probe==fn_name): #probing
535            return exp
536        elif(op_negation==fn_name): #negation
537            return applyToVector(exp, fn_negation)
538        elif(op_divergence==fn_name):
539            return fn_divergence(fld)
540        elif(op_curl==fn_name):
541            return fn_curl(fld)
542        elif(op_jacob==fn_name): #jacob
543            return fn_jacob(fld)
544        else:
545            raise Exception("unsupported unary operator:"+ fn_name.name)
546
547    def applyToExp_U_M(fn_name, fld):
548        exp = field.get_data(fld)
549        if(op_probe==fn_name): #probing
550            return exp
551        elif(op_negation==fn_name): #negation
552            return applyToM(exp, fn_negation)
553        elif(op_jacob==fn_name): #jacob
554            return fn_jacob(fld)
555        elif(op_slice==fn_name):
556            return fn_slice(fld)
557        elif(op_trace == fn_name):
558            return fn_trace(fld)
559        elif(op_transpose==fn_name):
560            return fn_transpose(fld)
561        elif(op_det==fn_name):
562            return fn_det(fld)
563        else:
564            raise Exception("unsupported unary operator:"+ fn_name.name)
565
566    def applyToExp_U_T3(fn_name, fld):
567        exp = field.get_data(fld)
568        if(op_probe==fn_name): #probing
569            return exp
570        elif(op_negation==fn_name): #negation
571
572            return applyToT3(exp, fn_negation)
573        elif(op_jacob==fn_name): #jacob
574            return fn_jacob(fld)
575        else:
576            raise Exception("unsupported unary operator:"+ fn_name.name)
577
578    # binary operator
579    # exp: scalar types
580    def applyToExp_B_S(e):
581        fn_name=e.opr
582        (fld1,fld2) =  apply.get_binary(e)
583        exp1 = field.get_data(fld1)
584        exp2 = field.get_data(fld2)
585        #print fn_name
588        elif(op_subtract==fn_name):#subtract
589            return fn_subtract(exp1,exp2)
590        elif(op_scale==fn_name): #scaling
591            return fn_scaling(fld1,fld2)
592        elif(op_division==fn_name): #division
593            return fn_division(fld1,fld2)
594        else:
595            raise Exception("unsupported binary operator on scalar fields:"+ fn_name.name)
596
597    # binary, args do not need to have the same shape
598    def applyToExp_B_uneven(e):
599        fn_name=e.opr
600        (fld1,fld2) =  apply.get_binary(e)
601        exp1 = field.get_data(fld1)
602        exp2 = field.get_data(fld2)
603        if(op_outer==fn_name):
604            return fn_outer(fld1, fld2)
605        elif(op_inner==fn_name):
606            return fn_inner(fld1, fld2)
607        elif(op_scale==fn_name): #scaling
608            return fn_scaling(fld1,fld2)
609        else:
610            raise Exception("unsupported unary operator:",op_name)
611
612
613    # binary operator
614    # args have the same shape
615    def applyToExp_B_V(e):
616        fn_name=e.opr
617        (fld1,fld2) =  apply.get_binary(e)
618        exp1 = field.get_data(fld1)
619        exp2 = field.get_data(fld2)
622        elif(op_subtract==fn_name):#subtract
623            return  applyToVectors(exp1, exp2, fn_subtract)
624        elif(op_cross==fn_name):
625            return fn_cross(fld1, fld2)
626        else:
627           return applyToExp_B_uneven(e)
628
629
630    # ***************************  unary/binary operators ***************************
631    def unary(e):
632        #apply.toStr(e)
633        fld =apply.get_unary(e)
634        fn_name=e.opr
635        exp = field.get_data(fld)
636        dim = field.get_dim(fld)
637        if(op_norm==fn_name):#norm
638            return fn_norm(fld, dim)
639        if(op_normalize==fn_name):#normalize
640            x= fn_normalize(fld, dim)
641            return x
642        elif (field.is_Scalar(fld)): # input is a scalar field
643            return applyToExp_U_S(fn_name, fld)
644        elif(field.is_Vector(fld)): # input is a vector field
645            return applyToExp_U_V(fn_name, fld)
646        elif(field.is_Matrix(fld)): # input is a vector field
647            return applyToExp_U_M(fn_name, fld)
648        else:
649            return applyToExp_U_T3(fn_name, fld)
650
651    def binary(e):
652        (f, g) =apply.get_binary(e)
653        fn_name = e.opr
654        # type is checked elsewhere or does not matter
655        if(op_division==fn_name): #division
656            return fn_division(f, g)
657        elif (field.is_Scalar(f) and field.is_Scalar(g)): # input is a scalar field
658            return applyToExp_B_S(e)
659        elif (field.is_Vector(f)):# input is a vector field
660            if(field.is_Vector(g)):
661                return applyToExp_B_V(e)
662            else: # input is a vector field, _
663                return applyToExp_B_uneven(e)
664        else:
665             return applyToExp_B_uneven(e)
666
667    def applyUnaryOnce(oexp_inner,app_inner,app_outer):
668        #print "applyUnaryOnce"
669        #apply.toStr(app_inner)
670        oty_inner = apply.get_oty(app_inner)
671        oty_outer = apply.get_oty(app_outer)
672        opr_outer = app_outer.opr
673        #print "oexp_inner",oexp_inner,"opr_outer",opr_outer.name
674        lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
675        app_tmp = apply("tmp", opr_outer, lhs_tmp, None, oty_outer, true, true)
676        oexp_tmp =unary(app_tmp)
677        #print " oexp_tmp", oexp_tmp
678        return (oty_outer, oexp_tmp)
679
680    def applyBinaryOnce(oexp_inner,app_inner,app_outer,rhs):
681        oty_inner = apply.get_oty(app_inner)
682        oty_outer = apply.get_oty(app_outer)
683        opr_outer = app_outer.opr
684
685        lhs_tmp = field(true, "tmp", oty_inner, "", oexp_inner, "")
686
687        app_tmp = apply("tmp", opr_outer, lhs_tmp, rhs, oty_outer, true,true)
688        oexp_tmp =binary(app_tmp)
689        return (oty_outer, oexp_tmp)
690
691
692    # operators with scalar field and vector field
693    def sort(e):
694        #apply.toStr(e)
695        arity = apply.get_arity(e)
696        if(e.isrootlhs): # is root
697            #print "sort is a root"
698            oty = apply.get_oty(e)
699            if (arity ==1):
700                return (oty, unary(e))
701            elif (arity ==2): # binary operator
702                return (oty, binary(e))
703            else:
704                raise Exception ("arity is not supported: "+str(arity_outer))
705        else:
706            app_outer = e
707            arity_outer = arity
708            #print "app_outer",app_outer.opr.name
709            if (arity_outer==1):  #assuming both arity
710                app_inner = apply.get_unary(app_outer)
711                #print "outer(1) app_inner:",app_inner.opr.name
712                arity_inner=  app_inner.opr.arity
713                if (arity_inner==1):
714                    oexp_inner = unary(app_inner)
715
716                    (oty_outer, oexp_tmp) =  applyUnaryOnce(oexp_inner ,app_inner, app_outer)
717
718                    return (oty_outer, oexp_tmp)
719                elif(arity_inner==2):
720                    oexp_inner = binary(app_inner)
721                    (oty_outer, oexp_tmp) =  applyUnaryOnce(oexp_inner ,app_inner, app_outer)
722
723                    return (oty_outer, oexp_tmp)
724                else:
725                    raise Exception ("arity is not supported: "+str(arity_outer))
726            elif (arity_outer==2):
727                (app_inner, G) = apply.get_binary(app_outer)
728                arity_inner=  app_inner.opr.arity
729                #print "outer(2) app_inner",app_inner.opr.name
730                if (arity_inner==1):
731                    oexp_inner = unary(app_inner)
732                    rhs = G
733                    (oty_outer, oexp_tmp) =  applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
734                    return (oty_outer, oexp_tmp)
735                elif(arity_inner==2):
736                    oexp_inner = binary(app_inner)
737                    #print "applying binary frst time", oexp_inner
738                    rhs = G
739                    (oty_outer, oexp_tmp) =  applyBinaryOnce(oexp_inner, app_inner, app_outer, rhs)
740                    #print "applying binary second time",  oexp_tmp
741                    return (oty_outer, oexp_tmp)
742                else:
743                    raise Exception ("arity is not supported: "+str(arity_outer))
744            else:
745                raise Exception ("arity is not supported: "+str(arity_outer))
746
747    # ***************************  evaluate at positions ***************************
748    #evaluate scalar field exp
749    def eval_d1(pos0, exp):
750        #print "exp:",exp
751  #evaluate field defined by coefficients at position  #evaluate field defined by coefficients at position
752  def eval(x,y,a,b,c,d,e,f,g,h,i):      exp = exp.subs(x,pos0)
753        return exp
754
755    def eval_d2(pos0, pos1, exp):
756        #print "exp:",exp
757        #evaluate field defined by coefficients at position
758        exp = exp.subs(x,pos0)
759        exp = exp.subs(y,pos1)
760        return exp
761
762    def eval_d3(pos0, pos1, pos2, exp):
763      #evaluate field defined by coefficients at position      #evaluate field defined by coefficients at position
764      t0 = b*x + c*y      exp = exp.subs(x,pos0)
765      t1 = d*x*x + g*y*x*x      exp = exp.subs(y,pos1)
766      t2 = f*y*y + h*x*y*y      exp = exp.subs(z,pos2)
767      t3 = a + e*x*y + i*(y*y)*(x*x)      return exp
768      t4= t0+t1+t2+t3
769      return t4  #evaluate vector field [exp]
770    def eval_vec_d1(pos0, vec):
771  def negation(x,y,a,b,c,d,e,f,g,h,i):      rtn = []
772      return (-1*eval(x,y,a,b,c,d,e,f,g,h,i))      for v in vec:
773            rtn.append(eval_d1(pos0, v))
774  #gradient of field defined by coefficients at position      return rtn
776      x0 = b  #evaluate vector field [exp,exp]
777      x1 = d*2*x + g*y*2*x  def eval_vec_d2(pos0, pos1, vec):
778      x2 = h*y*y      #print "eval_vec_d2 vec:",vec
779      x3 = e*y + i*(y*y)*(2*x)      rtn = []
780      x4= x0+x1+x2+x3 #xaxis      for v in vec:
781      y0 = c          rtn.append(eval_d2(pos0, pos1, v))
782      y1 = g*x*x      return rtn
783      y2 = f*2*y + h*x*2*y
784      y3 = e*x + i*(y*2)*(x*x)  def eval_ten3_d1(pos0,ten3):
785      y4= y0+y1+y2+y3 #yaxis      rtn = []
786      return [x4,y4]      for i in ten3:
787            for j in i:
788                for v in j:
789                    rtn.append(eval_d1(pos0, v))
790        return rtn
791
792
793
794    #evaluate vector field [exp,exp]
795    def eval_mat_d1(pos0, mat):
796        #print "eval_vec_d2 vec:",vec
797        rtn = []
798        for i in mat:
799            for v in i:
800                rtn.append(eval_d1(pos0, v))
801        return rtn
802
803    #evaluate vector field [exp,exp]
804    def eval_mat_d2(pos0, pos1, mat):
805        #print "eval_vec_d2 vec:",vec
806        rtn = []
807        print "eval_mat_d2 mat",mat
808        for i in mat:
809            for v in i:
810                rtn.append(eval_d2(pos0, pos1, v))
811        return rtn
812
813    def eval_ten3_d2(pos0, pos1, ten3):
814        #print "eval_vec_d2 vec:",vec
815        rtn = []
816        for i in ten3:
817            for j in i:
818                for v in j:
819                    rtn.append(eval_d2(pos0, pos1, v))
820        return rtn
821
822
823
824    #evaluate vector field [exp,exp]
825    def eval_vec_d3(pos0, pos1, pos2, vec):
826        rtn = []
827        for v in vec:
828            rtn.append(eval_d3(pos0, pos1, pos2, v))
829        return rtn
830
831
832    #evaluate vector field [exp,exp]
833    def eval_mat_d3(pos0, pos1, pos2, mat):
834        #print "eval_vec_d2 vec:",vec
835        rtn = []
836        for i in mat:
837            for v in i:
838                rtn.append(eval_d3(pos0, pos1, pos2, v))
839        return rtn
840
841    def eval_ten3_d3(pos0, pos1, pos2,ten3):
842        rtn = []
843        for i in ten3:
844            for j in i:
845                for v in j:
846                    rtn.append(eval_d3(pos0, pos1, pos2, v))
847        return rtn
848
849
850
851    def iter_d1(k, pos, exp):
852        corr = []
853        for x in pos:
854            val = k(x, exp)
855            corr.append(val)
856        return corr
857
858  def iter(k, pos, base, xsq, ysq, diag):  def iter_d2(k, pos, exp):
859      corr = []      corr = []
860      b = base[0]      #print "iter expr:", exp
861      c = base[1]      #print "pos", pos
d = xsq[0]
g = xsq[1]
f = ysq[0]
h = ysq[1]
a = diag[0]
e = diag[1]
i = diag[2]
862      for p in pos:      for p in pos:
863            #print "p", p
864          x=p[0]          x=p[0]
865          y=p[1]          y=p[1]
866          val = k(x,y,a,b,c,d,e,f,g,h,i)          val = k(x,y,exp)
867          corr.append(val)          corr.append(val)
868      return corr      return corr
869
870    def iter_d3(k, pos, exp):
871        corr = []
872        #print "iter exp:", exp
873        for p in pos:
874            x=p[0]
875            y=p[1]
876            z=p[2]
877            val = k(x,y,z, exp)
878            #print "pos: ",x,y,z, " val:", val
879            corr.append(val)
880        return corr
881
882  def unary(ex, pos,base,xsq,ysq,diag):  def probeField(otyp1,pos, ortn):
883      if(ex==op_probe): #probing      dim = fty.get_dim(otyp1)
884          return iter(eval,pos,base,xsq,ysq,diag)      #print "output type"+otyp1.name
885      elif(ex==op_negation): #negation      if (dim==1):
886          return iter(negation,pos,base,xsq,ysq,diag)          def get_k():
889                elif (fty.is_VectorField(otyp1)):
890                    return eval_vec_d1
891                elif (fty.is_Matrix(otyp1)):
892                    return eval_mat_d1
893                elif(fty.is_Ten3(otyp1)):
894                    return eval_ten3_d1
895      else:      else:
896          raise ("no operator value for "+str(ex))                  raise "error"+otyp1.name
897            return iter_d1(get_k(), pos, ortn)
898        elif (dim==2):
899            def get_k():
900                if (fty.is_ScalarField(otyp1)): # output is a scalar field
901                    return eval_d2
902                elif (fty.is_VectorField(otyp1)):
903                    return eval_vec_d2
904                elif (fty.is_Matrix(otyp1)):
905                    return eval_mat_d2
906                elif(fty.is_Ten3(otyp1)):
907                    return eval_ten3_d2
908                else:
909                    raise "error"+otyp1.name
910            return iter_d2(get_k(), pos, ortn)
911        elif (dim==3):
912            def get_k():
913                if (fty.is_ScalarField(otyp1)): # output is a scalar field
914                    return eval_d3
915                elif (fty.is_VectorField(otyp1)):
916                    return eval_vec_d3
917                elif (fty.is_Matrix(otyp1)):
918                    return eval_mat_d3
919                elif(fty.is_Ten3(otyp1)):
920                    return eval_ten3_d3
921                else:
922                    raise "error"+otyp1.name
923            return iter_d3(get_k(), pos, ortn)
924        else:
925            raise "unsupported dimension"
926
927    # ***************************  main  ***************************
928
929    # operators with scalar field and vector field
930    def eval(app, pos):
931        #print "evalname",app.name
932        # print "inside eval exp:", field.toStr(app.exp[0])
933        #apply.toStr(app)
934        (otyp1, ortn) = sort(app) #apply operations to expressions
935        print "ortn",ortn
936        rtn = probeField(otyp1, pos, ortn) #evaluate expression at positions
937        print "rtn", rtn
938        return rtn

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