Python sympy.symbols() Examples

def derivatives_in_paraboloidal_coordinates():
    #Print_Function()
    coords = (u,v,phi) = symbols('u v phi', real=True)
    (par3d,er,eth,ephi) = Ga.build('e_u e_v e_phi',X=[u*v*cos(phi),u*v*sin(phi),(u**2-v**2)/2],coords=coords,norm=True)
    grad = par3d.grad

    f = par3d.mv('f','scalar',f=True)
    A = par3d.mv('A','vector',f=True)
    B = par3d.mv('B','bivector',f=True)

    print('#Derivatives in Paraboloidal Coordinates')

    print('f =',f)
    print('A =',A)
    print('B =',B)

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    (-par3d.i*(grad^A)).Fmt(3,'grad\\times A = -I*(grad^A)')
    print('grad^B =',grad^B)

    return 

def test_construct_lazy(self):
        # adapted from https://gist.github.com/raddy/bd0e977dc8437a4f8276
        spot, strike, vol, dte, rate, cp = sy.symbols('spot strike vol dte rate cp')

        T = dte / 260.
        N = syNormal('N', 0.0, 1.0)

        d1 = (sy.ln(spot / strike) + (0.5 * vol ** 2) * T) / (vol * sy.sqrt(T))
        d2 = d1 - vol * sy.sqrt(T)

        TimeValueExpr = sy.exp(-rate * T) * (cp * spot * cdf(N)(cp * d1) - cp * strike * cdf(N)(cp * d2))

        PriceClass = ts.construct_lazy(TimeValueExpr)

        price = PriceClass(spot=210.59, strike=205, vol=14.04, dte=4, rate=.2175, cp=-1)

        x = price.evaluate()()

        assert price.evaluate()() == x

        price.strike = 210

        assert x != price.evaluate()() 

def associated_legendre_polynomials(k, zero_m_only=True):
    z = sym.symbols('z')
    P_l_m = [[0] * (j + 1) for j in range(k)]

    P_l_m[0][0] = 1
    if k > 0:
        P_l_m[1][0] = z

        for j in range(2, k):
            P_l_m[j][0] = sym.simplify(((2 * j - 1) * z * P_l_m[j - 1][0] -
                                        (j - 1) * P_l_m[j - 2][0]) / j)
        if not zero_m_only:
            for i in range(1, k):
                P_l_m[i][i] = sym.simplify((1 - 2 * i) * P_l_m[i - 1][i - 1])
                if i + 1 < k:
                    P_l_m[i + 1][i] = sym.simplify(
                        (2 * i + 1) * z * P_l_m[i][i])
                for j in range(i + 2, k):
                    P_l_m[j][i] = sym.simplify(
                        ((2 * j - 1) * z * P_l_m[j - 1][i] -
                         (i + j - 1) * P_l_m[j - 2][i]) / (j - i))

    return P_l_m 

def bessel_basis(n, k):
    zeros = Jn_zeros(n, k)
    normalizer = []
    for order in range(n):
        normalizer_tmp = []
        for i in range(k):
            normalizer_tmp += [0.5 * Jn(zeros[order, i], order + 1)**2]
        normalizer_tmp = 1 / np.array(normalizer_tmp)**0.5
        normalizer += [normalizer_tmp]

    f = spherical_bessel_formulas(n)
    x = sym.symbols('x')
    bess_basis = []
    for order in range(n):
        bess_basis_tmp = []
        for i in range(k):
            bess_basis_tmp += [
                sym.simplify(normalizer[order][i] *
                             f[order].subs(x, zeros[order, i] * x))
            ]
        bess_basis += [bess_basis_tmp]
    return bess_basis 

def diff_inference():
    
    r,s=symbols('r,s')
    
    la1,la2,lb1,lb2=symbols('l^a_1,l^a_2,l^b_1,l^b_2')
    
    la1=(1-s)/2
    la2=(1+s)/2
    lb1=(1-r)/2
    lb2=(1+r)/2
    N1=la1*lb1
    N2=la1*lb2
    N3=la2*lb1
    N4=la2*lb2
    
    N=Matrix([[N1,0,N2,0,N3,0,N4,0],
              [0,N1,0,N2,0,N3,0,N4]]) 

def test_jscode_Indexed():
    from sympy.tensor import IndexedBase, Idx
    from sympy import symbols
    i, j, k, n, m, o = symbols('i j k n m o', integer=True)

    p = JavascriptCodePrinter()
    p._not_c = set()

    x = IndexedBase('x')[Idx(j, n)]
    assert p._print_Indexed(x) == 'x[j]'
    A = IndexedBase('A')[Idx(i, m), Idx(j, n)]
    assert p._print_Indexed(A) == 'A[%s]' % str(j + n*i)
    B = IndexedBase('B')[Idx(i, m), Idx(j, n), Idx(k, o)]
    assert p._print_Indexed(B) == 'B[%s]' % str(k + i*n*o + j*o)

    assert p._not_c == set() 

def test_jscode_loops_matrix_vector():
    n, m = symbols('n m', integer=True)
    A = IndexedBase('A')
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx('i', m)
    j = Idx('j', n)

    s = (
        'for (var i=0; i<m; i++){\n'
        '   y[i] = 0;\n'
        '}\n'
        'for (var i=0; i<m; i++){\n'
        '   for (var j=0; j<n; j++){\n'
        '      y[i] = y[i] + A[i*n + j]*x[j];\n'
        '   }\n'
        '}'
    )
    c = jscode(A[i, j]*x[j], assign_to=y[i])
    assert c == s 

def test_dummy_loops():
    # the following line could also be
    # [Dummy(s, integer=True) for s in 'im']
    # or [Dummy(integer=True) for s in 'im']
    i, m = symbols('i m', integer=True, cls=Dummy)
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx(i, m)

    expected = (
        'for (var i_%(icount)i=0; i_%(icount)i<m_%(mcount)i; i_%(icount)i++){\n'
        '   y[i_%(icount)i] = x[i_%(icount)i];\n'
        '}'
    ) % {'icount': i.label.dummy_index, 'mcount': m.dummy_index}
    code = jscode(x[i], assign_to=y[i])
    assert code == expected 

def test_jscode_loops_add():
    from sympy.tensor import IndexedBase, Idx
    from sympy import symbols
    n, m = symbols('n m', integer=True)
    A = IndexedBase('A')
    x = IndexedBase('x')
    y = IndexedBase('y')
    z = IndexedBase('z')
    i = Idx('i', m)
    j = Idx('j', n)

    s = (
        'for (var i=0; i<m; i++){\n'
        '   y[i] = x[i] + z[i];\n'
        '}\n'
        'for (var i=0; i<m; i++){\n'
        '   for (var j=0; j<n; j++){\n'
        '      y[i] = y[i] + A[i*n + j]*x[j];\n'
        '   }\n'
        '}'
    )
    c = jscode(A[i, j]*x[j] + x[i] + z[i], assign_to=y[i])
    assert c == s 

def test_ccode_Indexed():
    from sympy.tensor import IndexedBase, Idx
    from sympy import symbols
    n, m, o = symbols('n m o', integer=True)
    i, j, k = Idx('i', n), Idx('j', m), Idx('k', o)
    p = CCodePrinter()
    p._not_c = set()

    x = IndexedBase('x')[j]
    assert p._print_Indexed(x) == 'x[j]'
    A = IndexedBase('A')[i, j]
    assert p._print_Indexed(A) == 'A[%s]' % (m*i+j)
    B = IndexedBase('B')[i, j, k]
    assert p._print_Indexed(B) == 'B[%s]' % (i*o*m+j*o+k)

    assert p._not_c == set() 

def test_dummy_loops():
    # the following line could also be
    # [Dummy(s, integer=True) for s in 'im']
    # or [Dummy(integer=True) for s in 'im']
    i, m = symbols('i m', integer=True, cls=Dummy)
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx(i, m)

    expected = (
        'for (int i_%(icount)i=0; i_%(icount)i<m_%(mcount)i; i_%(icount)i++){\n'
        '   y[i_%(icount)i] = x[i_%(icount)i];\n'
        '}'
    ) % {'icount': i.label.dummy_index, 'mcount': m.dummy_index}
    code = ccode(x[i], assign_to=y[i])
    assert code == expected 

def test_ccode_loops_add():
    from sympy.tensor import IndexedBase, Idx
    from sympy import symbols
    n, m = symbols('n m', integer=True)
    A = IndexedBase('A')
    x = IndexedBase('x')
    y = IndexedBase('y')
    z = IndexedBase('z')
    i = Idx('i', m)
    j = Idx('j', n)

    s = (
        'for (int i=0; i<m; i++){\n'
        '   y[i] = x[i] + z[i];\n'
        '}\n'
        'for (int i=0; i<m; i++){\n'
        '   for (int j=0; j<n; j++){\n'
        '      y[i] = y[i] + A[%s]*x[j];\n' % (i*n + j) +\
        '   }\n'
        '}'
    )
    c = ccode(A[i, j]*x[j] + x[i] + z[i], assign_to=y[i])
    assert c == s 

def __init__(self, in_vars, out_vars):
        """Initialize physics layer.

        Args:
          in_vars: str, a string of input variable names separated by space.
          E.g., 'x y t' for the three variables x, y and t.
          out_vars: str, a string of output variable names separated by space.
          E.g., 'u v p' for the three variables u, v and p.
        """
        self.in_vars = sympy.symbols(in_vars)
        self.out_vars = sympy.symbols(out_vars)
        if not isinstance(self.in_vars, tuple): self.in_vars = (self.in_vars,)
        if not isinstance(self.out_vars, tuple): self.out_vars = (self.out_vars,)
        self.n_in = len(self.in_vars)
        self.n_out = len(self.out_vars)
        self.all_vars = list(self.in_vars) + list(self.out_vars)
        self.eqns_raw = {}  # raw string equations
        self.eqns_fn = {}  # lambda function for the equations
        self.forward_method = None 

def _sympy_to_native(symbol):
        try:
            if symbol.is_Integer:
                return int(symbol)
            if symbol.is_Float:
                return float(symbol)
            if symbol.is_Complex:
                return complex(symbol)
            return symbol
        except Exception as ex:
            raise RuntimeError(
                "Could not convert sympy symbol to native type."
                "It may be due to misinterpretation of some symbols by sympy."
                "Try to use sympy expressions as gate parameters' values "
                "explicitly."
            ) from ex 

def derivatives_in_spherical_coordinates():
    Print_Function()
    X = (r,th,phi) = symbols('r theta phi')
    curv = [[r*cos(phi)*sin(th),r*sin(phi)*sin(th),r*cos(th)],[1,r,r*sin(th)]]
    (er,eth,ephi,grad) = MV.setup('e_r e_theta e_phi',metric='[1,1,1]',coords=X,curv=curv)

    f = MV('f','scalar',fct=True)
    A = MV('A','vector',fct=True)
    B = MV('B','grade2',fct=True)

    print('f =',f)
    print('A =',A)
    print('B =',B)

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',-MV.I*(grad^A))
    print('grad^B =',grad^B)
    return 

def derivatives_in_spherical_coordinates():

    X = (r,th,phi) = symbols('r theta phi')
    curv = [[r*cos(phi)*sin(th),r*sin(phi)*sin(th),r*cos(th)],[1,r,r*sin(th)]]
    (er,eth,ephi,grad) = MV.setup('e_r e_theta e_phi',metric='[1,1,1]',coords=X,curv=curv)

    f = MV('f','scalar',fct=True)
    A = MV('A','vector',fct=True)
    B = MV('B','grade2',fct=True)

    print('f =',f)
    print('A =',A)
    print('B =',B)

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',-MV.I*(grad^A))
    print('grad^B =',grad^B)
    return 

def MV_setup_options():
    (e1,e2,e3) = MV.setup('e_1 e_2 e_3','[1,1,1]')
    v = MV('v', 'vector')
    print(v)

    (e1,e2,e3) = MV.setup('e*1|2|3','[1,1,1]')
    v = MV('v', 'vector')
    print(v)

    (e1,e2,e3) = MV.setup('e*x|y|z','[1,1,1]')
    v = MV('v', 'vector')
    print(v)

    coords = symbols('x y z')
    (e1,e2,e3,grad) = MV.setup('e','[1,1,1]',coords=coords)
    v = MV('v', 'vector')
    print(v)

    return 

def derivatives_in_rectangular_coordinates():
    Print_Function()
    X = (x,y,z) = symbols('x y z')
    (ex,ey,ez,grad) = MV.setup('e_x e_y e_z',metric='[1,1,1]',coords=X)

    f = MV('f','scalar',fct=True)
    A = MV('A','vector',fct=True)
    B = MV('B','grade2',fct=True)
    C = MV('C','mv')
    print('f =',f)
    print('A =',A)
    print('B =',B)
    print('C =',C)

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('grad*A =',grad*A)

    print(-MV.I)

    print('-I*(grad^A) =',-MV.I*(grad^A))
    print('grad*B =',grad*B)
    print('grad^B =',grad^B)
    print('grad|B =',grad|B)
    return 

def derivatives_in_spherical_coordinates():
    Print_Function()
    X = (r,th,phi) = symbols('r theta phi')
    curv = [[r*cos(phi)*sin(th),r*sin(phi)*sin(th),r*cos(th)],[1,r,r*sin(th)]]
    (er,eth,ephi,grad) = MV.setup('e_r e_theta e_phi',metric='[1,1,1]',coords=X,curv=curv)

    f = MV('f','scalar',fct=True)
    A = MV('A','vector',fct=True)
    B = MV('B','grade2',fct=True)

    print('f =',f)
    print('A =',A)
    print('B =',B)

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',(-MV.I*(grad^A)).simplify())
    print('grad^B =',grad^B) 

def main():
    Format()
    a = Matrix ( 2, 2, ( 1, 2, 3, 4 ) )
    b = Matrix ( 2, 1, ( 5, 6 ) )
    c = a * b
    print(a,b,'=',c)

    x, y = symbols ( 'x, y' )

    d = Matrix ( 1, 2, ( x ** 3, y ** 3 ))
    e = Matrix ( 2, 2, ( x ** 2, 2 * x * y, 2 * x * y, y ** 2 ) )
    f = d * e

    print('%',d,e,'=',f)

    # xpdf()
    xpdf(pdfprog=None)
    return 

def Dirac_Equation_in_Geometric_Calculus():
    Print_Function()
    vars = symbols('t x y z')
    (g0,g1,g2,g3,grad) = MV.setup('gamma*t|x|y|z',metric='[1,-1,-1,-1]',coords=vars)
    I = MV.I

    (m,e) = symbols('m e')

    psi = MV('psi','spinor',fct=True)
    A = MV('A','vector',fct=True)
    sig_z = g3*g0

    print('\\text{4-Vector Potential\\;\\;}\\bm{A} =',A)
    print('\\text{8-component real spinor\\;\\;}\\bm{\\psi} =',psi)

    dirac_eq = (grad*psi)*I*sig_z-e*A*psi-m*psi*g0
    dirac_eq.simplify()

    dirac_eq.Fmt(3,r'%\text{Dirac Equation\;\;}\nabla \bm{\psi} I \sigma_{z}-e\bm{A}\bm{\psi}-m\bm{\psi}\gamma_{t} = 0')

    return 

def Simple_manifold_with_scalar_function_derivative():
    coords = (x,y,z) = symbols('x y z')
    basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3',metric='[1,1,1]',coords=coords)
    # Define surface
    mfvar = (u,v) = symbols('u v')
    X = u*e1+v*e2+(u**2+v**2)*e3
    print X
    MF = Manifold(X,mfvar)

    # Define field on the surface.
    g = (v+1)*log(u)

    # Method 1: Using old Manifold routines.
    VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g,u) + (MF.rbasis[1]/MF.E_sq)*diff(g,v)
    print 'Vector derivative =', VectorDerivative.subs({u:1,v:0})

    # Method 2: Using new Manifold routines.
    dg = MF.Grad(g)
    print 'Vector derivative =', dg.subs({u:1,v:0})
    return 

def Simple_manifold_with_scalar_function_derivative():
    Print_Function()
    coords = (x,y,z) = symbols('x y z')
    basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3',metric='[1,1,1]',coords=coords)
    # Define surface
    mfvar = (u,v) = symbols('u v')
    X = u*e1+v*e2+(u**2+v**2)*e3
    print '\\f{X}{u,v} =',X
    MF = Manifold(X,mfvar)
    (eu,ev) = MF.Basis()
    # Define field on the surface.
    g = (v+1)*log(u)

    print '\\f{g}{u,v} =',g

    # Method 1: Using old Manifold routines.
    VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g,u) + (MF.rbasis[1]/MF.E_sq)*diff(g,v)
    print '\\eval{\\nabla g}{u=1,v=0} =', VectorDerivative.subs({u:1,v:0})

    # Method 2: Using new Manifold routines.
    dg = MF.Grad(g)
    print '\\eval{\\f{Grad}{g}}{u=1,v=0} =', dg.subs({u:1,v:0})
    dg = MF.grad*g
    print '\\eval{\\nabla g}{u=1,v=0} =', dg.subs({u:1,v:0})
    return 

def derivatives_in_spherical_coordinates():
    Print_Function()

    X = (r, th, phi) = symbols('r theta phi')
    s3d = Ga('e_r e_theta e_phi', g=[1, r ** 2, r ** 2 * sin(th) ** 2], coords=X, norm=True)
    (er, eth, ephi) = s3d.mv()
    grad = s3d.grad

    f = s3d.mv('f', 'scalar', f=True)
    A = s3d.mv('A', 'vector', f=True)
    B = s3d.mv('B', 'bivector', f=True)

    print('f =', f)
    print('A =', A)
    print('B =', B)

    print('grad*f =', grad * f)
    print('grad|A =', grad | A)
    print('-I*(grad^A) =', -s3d.E() * (grad ^ A))
    print('grad^B =', grad ^ B)
    return 

def Mv_setup_options():
    Print_Function()

    (o3d,e1,e2,e3) = Ga.build('e_1 e_2 e_3',g=[1,1,1])
    v = o3d.mv('v', 'vector')
    print(v)

    (o3d,e1,e2,e3) = Ga.build('e*1|2|3',g=[1,1,1])
    v = o3d.mv('v', 'vector')
    print(v)

    (o3d,e1,e2,e3) = Ga.build('e*x|y|z',g=[1,1,1])
    v = o3d.mv('v', 'vector')
    print(v)

    coords = symbols('x y z',real=True)
    (o3d,e1,e2,e3) = Ga.build('e',g=[1,1,1],coords=coords)
    v = o3d.mv('v', 'vector')
    print(v)

    print(v.grade(2))
    print(v.i_grade)

    return 

def main():
    Eprint()

    X = (x,y,z) = symbols('x y z',real=True)
    (o3d,ex,ey,ez) = Ga.build('e_x e_y e_z',g=[1,1,1],coords=(x,y,z))

    A = x*(ey^ez) + y*(ez^ex) + z*(ex^ey)
    print('A =', A)
    print('grad^A =',(o3d.grad^A).simplify())
    print()

    f = o3d.mv(1/sqrt(x**2 + y**2 + z**2))
    print('f =', f)
    print('grad*f =',(o3d.grad*f).simplify())
    print()

    B = f*A
    print('B =', B)
    print()

    Curl_B = o3d.grad^B

    print('grad^B =', Curl_B.simplify())

    return 

def simplify_power(value, sample_args, context=None):
  """E.g., "Simplify ((x**2)**3/x**4)**2/x**3."."""
  del value  # unused
  if context is None:
    context = composition.Context()

  entropy, sample_args = sample_args.peel()

  variable = sympy.symbols(context.pop(), positive=True)
  unsimplified = polynomials.sample_messy_power(variable, entropy)
  answer = unsimplified.sympy()

  template = random.choice([
      'Simplify {unsimplified} assuming {variable} is positive.',
  ])
  return example.Problem(
      example.question(
          context, template, unsimplified=unsimplified, variable=variable),
      answer) 

def test_others(self):
        # adapted from https://gist.github.com/raddy/bd0e977dc8437a4f8276
        spot, strike, vol, dte, rate, cp = sy.symbols('spot strike vol dte rate cp')
        T = dte / 260.
        N = syNormal('N', 0.0, 1.0)
        d1 = (sy.ln(spot / strike) + (0.5 * vol ** 2) * T) / (vol * sy.sqrt(T))
        d2 = d1 - vol * sy.sqrt(T)
        TimeValueExpr = sy.exp(-rate * T) * (cp * spot * cdf(N)(cp * d1) - cp * strike * cdf(N)(cp * d2))
        PriceClass = ts.construct_lazy(TimeValueExpr)
        price = PriceClass(spot=210.59, strike=205, vol=14.04, dte=4, rate=.2175, cp=-1)
        price.evaluate()()
        ts.graphviz(TimeValueExpr)
        assert ts.traversal(TimeValueExpr)
        assert ts.symbols(TimeValueExpr) 

def test_construct_streaming(self):
        # adapted from https://gist.github.com/raddy/bd0e977dc8437a4f8276
        # spot, strike, vol, days till expiry, interest rate, call or put (1,-1)
        spot, strike, vol, dte, rate, cp = sy.symbols('spot strike vol dte rate cp')

        T = dte / 260.
        N = syNormal('N', 0.0, 1.0)

        d1 = (sy.ln(spot / strike) + (0.5 * vol ** 2) * T) / (vol * sy.sqrt(T))
        d2 = d1 - vol * sy.sqrt(T)

        TimeValueExpr = sy.exp(-rate * T) * (cp * spot * cdf(N)(cp * d1) - cp * strike * cdf(N)(cp * d2))

        PriceClass = ts.construct_streaming(TimeValueExpr)

        def strikes():
            strike = 205
            while strike <= 220:
                yield strike
                strike += 2.5

        price = PriceClass(spot=tss.Const(210.59),
                           #    strike=tss.Print(tss.Const(205), text='strike'),
                           strike=tss.Foo(strikes, interval=1),
                           vol=tss.Const(14.04),
                           dte=tss.Const(4),
                           rate=tss.Const(.2175),
                           cp=tss.Const(-1))

        ret = tss.run(tss.Print(price._node))
        time.sleep(2)
        print(ret)
        assert len(ret) == 7 

def Maxwells_Equations_in_Geometric_Calculus():
    Print_Function()
    X = symbols('t x y z')
    (g0,g1,g2,g3,grad) = MV.setup('gamma*t|x|y|z',metric='[1,-1,-1,-1]',coords=X)
    I = MV.I

    B = MV('B','vector',fct=True)
    E = MV('E','vector',fct=True)
    B.set_coef(1,0,0)
    E.set_coef(1,0,0)
    B *= g0
    E *= g0
    J = MV('J','vector',fct=True)
    F = E+I*B

    print(r'\text{Pseudo Scalar\;\;}I =',I)
    print('\\text{Magnetic Field Bi-Vector\\;\\;} B = \\bm{B\\gamma_{t}} =',B)
    print('\\text{Electric Field Bi-Vector\\;\\;} E = \\bm{E\\gamma_{t}} =',E)
    print('\\text{Electromagnetic Field Bi-Vector\\;\\;} F = E+IB =',F)
    print('%\\text{Four Current Density\\;\\;} J =',J)
    gradF = grad*F
    print('#Geometric Derivative of Electomagnetic Field Bi-Vector')
    gradF.Fmt(3,'grad*F')

    print('#Maxwell Equations')
    print('grad*F = J')
    print('#Div $E$ and Curl $H$ Equations')
    (gradF.grade(1)-J).Fmt(3,'%\\grade{\\nabla F}_{1} -J = 0')
    print('#Curl $E$ and Div $B$ equations')
    (gradF.grade(3)).Fmt(3,'%\\grade{\\nabla F}_{3} = 0')
    return