In [17]:

    

def generateBinv():
    from sympy import symbols, Rational, Matrix
    from sympy.vector import Vector

    x,y,z = symbols('x y z')
    wholeFunc = Rational(0)
    alphaVec = []
    i = 0
    while i <= 3:
        j = 0
        while j <= 3:
            k = 0
            while k <= 3:
                a = symbols('a_{0}'.format(k+4*(j + 4*i)))
                alphaVec.append(a)
                wholeFunc += a*x**Rational(i)*y**Rational(j)*z**Rational(k)
                k += 1
            j += 1
        i += 1

    #print(wholeFunc)

    #order per corner f,fx,fy,fz,fxy,fxz,fyz,fxyz
    cornerVec = Matrix([wholeFunc,
                 wholeFunc.diff(x),
                 wholeFunc.diff(y),
                 wholeFunc.diff(z),
                 wholeFunc.diff(x,y),
                 wholeFunc.diff(x,z),
                 wholeFunc.diff(y,z),
                 wholeFunc.diff(x,y,z)])
    #print(cornerVec)
    #build B suchthat b = B.a

    def getBRows(vec,alpha):
        '''get the rows from cornervec evaluated at the specific points'''
        B = []
        for el in vec:
            row = []
            for a in alpha:
                k = Rational(el.diff(a))
                row.append(k)
            B.append(row)
        return Matrix(B)

    #x=0,y=0,z=0
    row0 = getBRows(cornerVec.subs({x:0,y:0,z:0}),alphaVec)
    row1 = getBRows(cornerVec.subs({x:0,y:0,z:1}),alphaVec)
    row2 = getBRows(cornerVec.subs({x:0,y:1,z:0}),alphaVec)
    row3 = getBRows(cornerVec.subs({x:0,y:1,z:1}),alphaVec)
    row4 = getBRows(cornerVec.subs({x:1,y:0,z:0}),alphaVec)
    row5 = getBRows(cornerVec.subs({x:1,y:0,z:1}),alphaVec)
    row6 = getBRows(cornerVec.subs({x:1,y:1,z:0}),alphaVec)
    row7 = getBRows(cornerVec.subs({x:1,y:1,z:1}),alphaVec)
    Binv = Matrix([row0,row1,row2,row3,row4,row5,row6,row7]).inv()
    string = "["
    for i in range(64):
        string += '['+('{:},'*63).format(*Binv[i,:64])+'{:}'.format(Binv[i,-1])+'],\n'
    string += ']'
    print(string)
    
def optimizeBvecFormation(order=1):
    from sympy import symbols, Matrix
    from sympy import cse, Function, Rational
    from scipy.special import binom
    
    def centralDiff(order, inFunc,var):
        #args = inFunc.args
        #func = inFunc.func
        #terms = []
        x,y,z,h = symbols('x y z h')
        outFunc = 0
        k = 0
        while k <= order:
            if var == 'x':
                outFunc += Rational(-1)**Rational(k) * Rational(int(binom(order,k))) * inFunc.subs({x:x + (Rational(order,2) - Rational(k))*h})
            if var == 'y':
                outFunc += Rational(-1)**Rational(k) * Rational(int(binom(order,k))) * inFunc.subs({y:y + (Rational(order,2) - Rational(k))*h})
            if var == 'z':
                outFunc += Rational(-1)**Rational(k) * Rational(int(binom(order,k))) * inFunc.subs({z:z + (Rational(order,2) - Rational(k))*h})
            k += 1
        return outFunc
    
    def centralDiff(order, inFunc,var):
        #args = inFunc.args
        #func = inFunc.func
        #terms = []
        x,y,z = symbols('xi yi zi', intergers=True)
        h = symbols('h')
        #print(var)
        if var == 'x':
            var = x
        if var == 'y':
            var = y
        if var == 'z':
            var = z
        outFunc = Rational(8)*(inFunc.subs({var:(var + h)}) - inFunc.subs({var:(var - h)}))/Rational(12)/h + (inFunc.subs({var:(var - Rational(2)*h)}) - inFunc.subs({var:(var + Rational(2)*h)}))/Rational(12)/h
        #print(inFunc)
        #print(Rational(8)*inFunc.subs({var:(var + h)})/Rational(12)/h)
        #print(-Rational(8)*inFunc.subs({var:(var - h)})/Rational(12)/h)
        #print(inFunc.subs({var:(var - Rational(2)*h)})/Rational(12)/h)
        #print(-inFunc.subs({var:(var + Rational(2)*h)})/Rational(12)/h)
        #outFunc = (inFunc.subs({var:var + h}) - inFunc.subs({var:var}))/h
        return outFunc
    
    vec = []
    
    xi,yi,zi,nz,ny = symbols('xi yi zi nz ny', intergers=True)
    func = Function('f')
    xf,yf,zf = Function('xf'),Function('yf'),Function('zf')
    vec = []
    #order
    #000,001,010,011,100,101,110,111
    i = 0
    while i <= 1:
        j = 0
        while j <= 1:
            k = 0
            while k <= 1:
                
                dx = (xf(xi+Rational(i+1)) - xf(xi+Rational(i-1)))/Rational(2)
                dy = (yf(yi+Rational(j+1)) - yf(yi+Rational(j-1)))/Rational(2)
                dz = (zf(zi+Rational(k+1)) - zf(zi+Rational(k-1)))/Rational(2)
                #f 
                f = func((zi+Rational(k)) + nz*((yi+Rational(j)) + ny*(xi+Rational(i))))#x+Rational(i),y+Rational(j),z+Rational(k))
                vec.append(f)
                
                #fx
                fx = centralDiff(order,f,'x')/dx
                vec.append(fx)
                
                #fy
                fy = centralDiff(order,f,'y')/dy
                vec.append(fy)
                
                #fz
                fz = centralDiff(order,f,'z')/dz
                vec.append(fz)
                
                #fxy
                fxy = centralDiff(order,fx,'y')/dy
                vec.append(fxy)
                
                #fxz
                fxz = centralDiff(order,fx,'z')/dz
                vec.append(fxz)
                
                #fyz
                fyz = centralDiff(order,fy,'z')/dz
                vec.append(fyz)
                
                #fxyz
                fxyz = centralDiff(order,fxy,'z')/dz
                vec.append(fxyz)
                
                k += 1
            j += 1
        i += 1
    vec = Matrix(vec).subs({'h':1})
    #print(vec)
    cseFunc = cse(vec,optimizations='basic')
    lines = []
    for optLine in cseFunc[0]:
        line = "{0} = {1}".format(optLine[0],optLine[1].evalf())
        line = line.replace("xf","self.get_xvec").replace("yf","self.get_yvec").replace("zf","self.get_zvec").replace("f(","self.get_m(")
        line = line.replace("xi","i").replace("yi","j").replace("zi","k")
        #line = line.replace("(","[").replace(")","]")
        lines.append(line)
        print(line)
    
    #print(lines)
    out = "{0}".format(cseFunc[1][0].transpose()[0,:].evalf())
    out = out.replace("xf","self.get_xvec").replace("yf","self.get_yvec").replace("zf","self.get_zvec").replace("f(","self.get_m(")
    out = out.replace("xi","i").replace("yi","j").replace("zi","k")
    print(out)
    
    def index1(i):
        str = ''
        if i == -1:
            str += 'm'
        if i == 0:
            str += 'z'
        if i == 1:
            str += 'p'
        if i == 2:
            str += 'P'
        return str
    
    def index2(i):
        str = ''
        if i == -1:
            str += 'm'
        if i == 0:
            str += 'z'
        if i == 1:
            str += 'p'
        if i == 2:
            str += 'P'
        return str
    
    def index(i,j,k):
        str = ''
        str += index1(i) + index1(j) + index1(k)
        return str

    vec = []

    i = 0
    while i <= 1:
        j = 0
        while j <= 1:
            k = 0
            while k <= 1:
                #f
                vec.append(symbols('f_{0}'.format(index(i,j,k))))
                
                #x10 = symbols('x_{0}'.format(index1(i))) - symbols('x_{0}'.format(index1(i-1)))
                #x02 = symbols('x_{0}'.format(index1(i+1))) - symbols('x_{0}'.format(index1(i)))
                #y10 = symbols('y_{0}'.format(index1(j))) - symbols('y_{0}'.format(index1(j-1)))
                #y02 = symbols('y_{0}'.format(index1(j+1))) - symbols('y_{0}'.format(index1(j)))
                #z10 = symbols('z_{0}'.format(index1(k))) - symbols('z_{0}'.format(index1(k-1)))
                #z02 = symbols('z_{0}'.format(index1(k+1))) - symbols('z_{0}'.format(index1(k)))
                
                x12 = symbols('x_{0}'.format(index1(i+1))) - symbols('x_{0}'.format(index1(i-1)))
                y12 = symbols('y_{0}'.format(index1(j+1))) - symbols('y_{0}'.format(index1(j-1)))
                z12 = symbols('z_{0}'.format(index1(k+1))) - symbols('z_{0}'.format(index1(k-1)))
                
                #fx,fy,fz
                #f0 = symbols('f_{0}'.format(index(i,j,k)))
                #fmzz = symbols('f_{0}'.format(index(i-1,j,k)))
                #fpzz = symbols('f_{0}'.format(index(i+1,j,k))) 
                #vec.append(((fmzz - f0)*x02 - (fpzz - f0)*x10)/(2*x10*x02))
                
                #fzmz = symbols('f_{0}'.format(index(i,j-1,k)))
                #fzpz = symbols('f_{0}'.format(index(i,j+1,k))) 
                #vec.append(((fzmz - f0)*y02 - (fzpz - f0)*y10)/(2*y10*y02))
                
                #fzzm = symbols('f_{0}'.format(index(i,j,k-1)))
                #fzzp = symbols('f_{0}'.format(index(i,j,k+1))) 
                #vec.append(((fzzm - f0)*z02 - (fzzp - f0)*z10)/(2*z10*z02))
                
                vec.append((symbols('f_{0}'.format(index(i+1,j,k))) - symbols('f_{0}'.format(index(i-1,j,k))) )/x12)
                vec.append((symbols('f_{0}'.format(index(i,j+1,k))) - symbols('f_{0}'.format(index(i,j-1,k))) )/y12)
                vec.append((symbols('f_{0}'.format(index(i,j,k+1))) - symbols('f_{0}'.format(index(i,j,k-1))) )/z12)

                #fxy,fxz,fyz
                vec.append((((symbols('f_{0}'.format(index(i+1,j+1,k))) - symbols('f_{0}'.format(index(i-1,j+1,k))) )/x12)-((symbols('f_{0}'.format(index(i+1,j-1,k))) - symbols('f_{0}'.format(index(i-1,j-1,k))) )/x12))/y12)
                vec.append((((symbols('f_{0}'.format(index(i+1,j,k+1))) - symbols('f_{0}'.format(index(i-1,j,k+1))) )/x12)-((symbols('f_{0}'.format(index(i+1,j,k-1))) - symbols('f_{0}'.format(index(i-1,j,k-1))) )/x12))/z12)
                vec.append((((symbols('f_{0}'.format(index(i,j+1,k+1))) - symbols('f_{0}'.format(index(i,j-1,k+1))) )/y12)-((symbols('f_{0}'.format(index(i,j+1,k-1))) - symbols('f_{0}'.format(index(i,j-1,k-1))) )/y12))/z12)
                
                #fxyz
                vec.append((((((symbols('f_{0}'.format(index(i+1,j+1,k+1))) - symbols('f_{0}'.format(index(i-1,j+1,k+1))) )/x12)-((symbols('f_{0}'.format(index(i+1,j-1,k+1))) - symbols('f_{0}'.format(index(i-1,j-1,k+1))) )/x12))/y12)-((((symbols('f_{0}'.format(index(i+1,j+1,k-1))) - symbols('f_{0}'.format(index(i-1,j+1,k-1))) )/x12)-((symbols('f_{0}'.format(index(i+1,j-1,k-1))) - symbols('f_{0}'.format(index(i-1,j-1,k-1))) )/x12))/y12))/z12)         

                k += 1
            j += 1
        i += 1
    vec = Matrix(vec)

    cseFunc = cse(vec,optimizations='basic')
    #generate the indices
    lines = ['im = i - 1','iz = i','ip = i + 1','iP = i + 2',
            'jm = j - 1','jz = j','jp = j + 1','jP = j + 2',
            'km = k - 1','kz = k','kp = k + 1','kP = k + 2']
    i = -1
    while i <= 2:
        j = -1
        while j <= 2:
            k = -1
            while k <= 2:
                var = index(i,j,k)
                line = "{0} = self.index(i{1},j{2},k{3})".format(index(i,j,k),index2(i),index2(j),index2(k))
                lines.append(line)
                k += 1
            j += 1
        i += 1
    def replaceIndices(f):
        i = -1
        while i <= 2:
            j = -1
            while j <= 2:
                k = -1
                while k <= 2:
                    var = index(i,j,k)
                    f = f.replace(var,'[{0}]'.format(var))
                    k += 1
                j += 1
            i += 1
        return f
    def replaceIndices2(f):
        f = f.replace('x_m','self.xvec[im]')
        f = f.replace('x_z','self.xvec[iz]')
        f = f.replace('x_p','self.xvec[ip]')
        f = f.replace('x_P','self.xvec[iP]')
        f = f.replace('y_m','self.yvec[jm]')
        f = f.replace('y_z','self.yvec[jz]')
        f = f.replace('y_p','self.yvec[jp]')
        f = f.replace('y_P','self.yvec[jP]')
        f = f.replace('z_m','self.zvec[km]')
        f = f.replace('z_z','self.zvec[kz]')
        f = f.replace('z_p','self.zvec[kp]')
        f = f.replace('z_P','self.zvec[kP]')
        return f
        
    for expr in cseFunc[0]:
        f = str(expr[1])
        f = replaceIndices(f)
        f = replaceIndices2(f)
        f = f.replace('f_','self.m')
        line = '{0} = {1}'.format(expr[0], f)
        lines.append(line)
        
    bvec = str(cseFunc[1][0].transpose())
    bvec = bvec.replace('Matrix([','bvec = np.array(')
    bvec = bvec.replace('])',')')
    bvec = bvec.replace(',',',\n')
    bvec = replaceIndices(bvec)
    bvec = replaceIndices2(bvec)
    bvec = bvec.replace('f_','self.m')
    lines.append(bvec)
    code = ''
    for line in lines:
        code += line+'\n'
    print(code)

if __name__=='__main__':
    genFuncCalls()
    #generateBinv()
    #optimizeBvecFormation()
    #testResult()


    

    




---------------------------------------------------------------------------
ShapeError                                Traceback (most recent call last)
<ipython-input-17-847da4192470> in <module>()
    364 
    365 if __name__=='__main__':
--> 366     genFuncCalls()
    367     #generateBinv()
    368     #optimizeBvecFormation()

<ipython-input-17-847da4192470> in genFuncCalls()
     33 
     34     cseFunc = cse(cornerVec,optimizations='basic')
---> 35     print(cornerVec.dot(alphaVec.transpose()))
     36 
     37 

C:\Users\josh\Anaconda3\envs\mayavi_env\lib\site-packages\sympy\matrices\matrices.pyc in dot(self, b)
   1746                 mat = mat.T
   1747                 b = b.T
-> 1748             prod = flatten((mat*b).tolist())
   1749             if len(prod) == 1:
   1750                 return prod[0]

C:\Users\josh\Anaconda3\envs\mayavi_env\lib\site-packages\sympy\core\decorators.pyc in binary_op_wrapper(self, other)
    116                     else:
    117                         return f(self)
--> 118             return func(self, other)
    119         return binary_op_wrapper
    120     return priority_decorator

C:\Users\josh\Anaconda3\envs\mayavi_env\lib\site-packages\sympy\matrices\dense.pyc in __mul__(self, other)
    553     @call_highest_priority('__rmul__')
    554     def __mul__(self, other):
--> 555         return super(DenseMatrix, self).__mul__(_force_mutable(other))
    556 
    557     @call_highest_priority('__rmul__')

C:\Users\josh\Anaconda3\envs\mayavi_env\lib\site-packages\sympy\matrices\matrices.pyc in __mul__(self, other)
    502             B = other
    503             if A.cols != B.rows:
--> 504                 raise ShapeError("Matrices size mismatch.")
    505             if A.cols == 0:
    506                 return classof(A, B)._new(A.rows, B.cols, lambda i, j: 0)

ShapeError: Matrices size mismatch.

In [52]:

    

from sympy import *

x,y,z = symbols('x y z')
f = Function('f')
f1 = f(x+1,y,z)
f2 = f(f1.args[0] +1,y,z) + f1
e = f1.func
h = f1 - f(x-3,y,z)
g = Function('g')
g = f1 + f(x,y,z)


    

In [54]:

    

g.subs({'x':x+6})


    

    
Out[54]:





f(x + 6, y, z) + f(x + 7, y, z)

In [116]:

    

a,o = symbols('a o')
from sympy.solvers import solve
ne = a*exp(-x**2/(o*log(2)*Rational(7,2))**2)
c = cse(ne.diff(x,x,x,x,x).subs({o:10,a:1e13}))
print(c)
for k in c[0]:
    print("{0} = {1}".format(k[0],k[1]))


    

    




([(x0, log(2)), (x1, x**2/x0**2)], [18.1254818541641*x*(-4*x**4/(625*x0**4) + 196*x1/5 - 36015)*exp(-x1/1225)/x0**6])
x0 = log(2)
x1 = x**2/x0**2

In [102]:

    

import numpy as np
f=lambdify(x,ne.diff(x,x,x,x,x).subs({o:10,a:1e13}),"numpy")
import pylab as plt
plt.plot(np.linspace(0,100,1000),f(np.linspace(0,100,1000)))
plt.show()
ne.diff(x,x,x,x,x).subs({x:o})


    

    
Out[102]:





-0.387942463307711*a/o**5

In [46]:

    

def genFuncCalls():
    from sympy import symbols, Rational, Matrix,Function
    from sympy.vector import Vector
    from sympy import cse

    x,y,z = symbols('x y z')
    a = Function('a')
    wholeFunc = []
    alphaVec = []
    i = 0
    while i <= 3:
        j = 0
        while j <= 3:
            k = 0
            while k <= 3:
                afunc = a(Rational(k)+Rational(4)*(Rational(j) + Rational(4)*Rational(i)))
                alphaVec.append(afunc)
                wholeFunc.append(x**Rational(i)*y**Rational(j)*z**Rational(k))
                k += 1
            j += 1
        i += 1
    
    alphaVec = Matrix(alphaVec)
    wholeFunc = Matrix(wholeFunc)
    cornerVec = Matrix([wholeFunc.transpose()[:],
                 wholeFunc.diff(x).transpose()[:],
                 wholeFunc.diff(y).transpose()[:],
                 wholeFunc.diff(z).transpose()[:],
                 wholeFunc.diff(x,y).transpose()[:],
                 wholeFunc.diff(x,z).transpose()[:],
                 wholeFunc.diff(y,z).transpose()[:],
                 wholeFunc.diff(x,y,z).transpose()[:]])
    
    cseFunc = cse(cornerVec,optimizations='basic')
    #print(cornerVec.shape,alphaVec.shape)
    #print(Matrix(cornerVec.dot(alphaVec)))
    
    #print (cseFunc)
    lines = []
    for optLine in cseFunc[0]:
        line = "{0} = {1}".format(optLine[0],optLine[1])
        #line = line.replace("xf","self.get_xvec").replace("yf","self.get_yvec").replace("zf","self.get_zvec").replace("f(","self.get_m(")
        #line = line.replace("xi","i").replace("yi","j").replace("zi","k")
        #line = line.replace("(","[").replace(")","]")
        lines.append(line)
        print(line)
    
    #print(lines)
    out = "np.array({0},dtype=np.double)".format(cseFunc[1][0])
    #out = out.replace("xf","self.get_xvec").replace("yf","self.get_yvec").replace("zf","self.get_zvec").replace("f(","self.get_m(")
    #out = out.replace("xi","i").replace("yi","j").replace("zi","k")
    print(out)
    
    

genFuncCalls()


    

    




x0 = z**2
x1 = z**3
x2 = y*z
x3 = x0*y
x4 = x1*y
x5 = y**2
x6 = x5*z
x7 = x0*x5
x8 = x1*x5
x9 = y**3
x10 = x9*z
x11 = x0*x9
x12 = x1*x9
x13 = x*z
x14 = x*x0
x15 = x*x1
x16 = x*y
x17 = x16*z
x18 = x0*x16
x19 = x*x5
x20 = x13*x5
x21 = x0*x19
x22 = x*x9
x23 = x0*x22
x24 = x**2
x25 = x24*z
x26 = x0*x24
x27 = x1*x24
x28 = x24*y
x29 = x2*x24
x30 = x0*x28
x31 = x24*x5
x32 = x25*x5
x33 = x0*x31
x34 = x24*x9
x35 = x0*x34
x36 = x**3
x37 = x36*z
x38 = x0*x36
x39 = x1*x36
x40 = x36*y
x41 = x0*x40
x42 = x36*x5
x43 = x37*x5
x44 = x0*x42
x45 = x36*x9
x46 = 2*x
x47 = x46*z
x48 = x0*x46
x49 = x1*x46
x50 = x46*y
x51 = x2*x46
x52 = x0*x50
x53 = x1*x50
x54 = x46*x5
x55 = x47*x5
x56 = x46*x9
x57 = x47*x9
x58 = 3*x24
x59 = x58*z
x60 = x0*x58
x61 = x1*x58
x62 = x58*y
x63 = x3*x58
x64 = x5*x58
x65 = x58*x6
x66 = x58*x7
x67 = x58*x8
x68 = x58*x9
x69 = x11*x58
x70 = 2*y
x71 = x70*z
x72 = x0*x70
x73 = x1*x70
x74 = 3*x5
x75 = x74*z
x76 = x0*x74
x77 = x1*x74
x78 = x*x74
x79 = x14*x74
x80 = x24*x70
x81 = x24*x71
x82 = x36*x70
x83 = x36*x71
x84 = x36*x74
x85 = x38*x74
x86 = 2*z
x87 = 3*x0
x88 = x87*y
x89 = x5*x86
x90 = x86*x9
x91 = x87*x9
x92 = x*x87
x93 = x24*x86
x94 = x36*x86
x95 = x36*x87
x96 = 4*x16
x97 = x96*z
x98 = 6*x19
x99 = 6*x20
x100 = x0*x98
x101 = 6*x28
x102 = 6*x29
x103 = x0*x101
x104 = 9*x31
x105 = x0*x104
x106 = 4*x13
x107 = 6*x14
x108 = 6*x18
x109 = 6*x25
x110 = 9*x26
x111 = x109*x5
x112 = 4*x2
x113 = 6*x3
x114 = 6*x6
x115 = 9*x7
np.array(Matrix([[1, z, x0, x1, y, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x, x13, x14, x15, x16, x17, x18, x1*x16, x19, x20, x21, x1*x19, x22, x13*x9, x23, x1*x22, x24, x25, x26, x27, x28, x29, x30, x1*x28, x31, x32, x33, x1*x31, x34, x25*x9, x35, x1*x34, x36, x37, x38, x39, x40, x2*x36, x41, x1*x40, x42, x43, x44, x1*x42, x45, x37*x9, x0*x45, x1*x45], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, z, x0, x1, y, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x46, x47, x48, x49, x50, x51, x52, x53, x54, x55, x0*x54, x1*x54, x56, x57, x0*x56, x1*x56, x58, x59, x60, x61, x62, x2*x58, x63, x4*x58, x64, x65, x66, x67, x68, x10*x58, x69, x12*x58], [0, 0, 0, 0, 1, z, x0, x1, x70, x71, x72, x73, x74, x75, x76, x77, 0, 0, 0, 0, x, x13, x14, x15, x50, x51, x52, x53, x78, x13*x74, x79, x15*x74, 0, 0, 0, 0, x24, x25, x26, x27, x80, x81, x0*x80, x1*x80, x64, x65, x66, x67, 0, 0, 0, 0, x36, x37, x38, x39, x82, x83, x0*x82, x1*x82, x84, x37*x74, x85, x39*x74], [0, 1, x86, x87, 0, y, x71, x88, 0, x5, x89, x76, 0, x9, x90, x91, 0, x, x47, x92, 0, x16, x51, x16*x87, 0, x19, x55, x79, 0, x22, x57, x22*x87, 0, x24, x93, x60, 0, x28, x81, x63, 0, x31, x5*x93, x66, 0, x34, x9*x93, x69, 0, x36, x94, x95, 0, x40, x83, x40*x87, 0, x42, x5*x94, x85, 0, x45, x9*x94, x45*x87], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, z, x0, x1, x70, x71, x72, x73, x74, x75, x76, x77, 0, 0, 0, 0, x46, x47, x48, x49, x96, x97, x0*x96, x1*x96, x98, x99, x100, x1*x98, 0, 0, 0, 0, x58, x59, x60, x61, x101, x102, x103, x1*x101, x104, 9*x32, x105, x1*x104], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, x86, x87, 0, y, x71, x88, 0, x5, x89, x76, 0, x9, x90, x91, 0, x46, x106, x107, 0, x50, x97, x108, 0, x54, x106*x5, x100, 0, x56, x106*x9, 6*x23, 0, x58, x109, x110, 0, x62, x102, 9*x30, 0, x64, x111, x105, 0, x68, x109*x9, 9*x35], [0, 0, 0, 0, 0, 1, x86, x87, 0, x70, x112, x113, 0, x74, x114, x115, 0, 0, 0, 0, 0, x, x47, x92, 0, x50, x97, x108, 0, x78, x99, 9*x21, 0, 0, 0, 0, 0, x24, x93, x60, 0, x80, x112*x24, x103, 0, x64, x111, x105, 0, 0, 0, 0, 0, x36, x94, x95, 0, x82, x112*x36, 6*x41, 0, x84, 6*x43, 9*x44], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, x86, x87, 0, x70, x112, x113, 0, x74, x114, x115, 0, 0, 0, 0, 0, x46, x106, x107, 0, x96, 8*x17, 12*x18, 0, x98, 12*x20, 18*x21, 0, 0, 0, 0, 0, x58, x109, x110, 0, x101, 12*x29, 18*x30, 0, x104, 18*x32, 27*x33]]),dtype=np.double)

In [ ]: