Python sympy.simplify() Examples

def get_equations(self):
        """
        :return: Functions to calculate A, B and f given state x and input u
        """
        f = sp.zeros(3, 1)

        x = sp.Matrix(sp.symbols('x y theta', real=True))
        u = sp.Matrix(sp.symbols('v w', real=True))

        f[0, 0] = u[0, 0] * sp.cos(x[2, 0])
        f[1, 0] = u[0, 0] * sp.sin(x[2, 0])
        f[2, 0] = u[1, 0]

        f = sp.simplify(f)
        A = sp.simplify(f.jacobian(x))
        B = sp.simplify(f.jacobian(u))

        f_func = sp.lambdify((x, u), f, 'numpy')
        A_func = sp.lambdify((x, u), A, 'numpy')
        B_func = sp.lambdify((x, u), B, 'numpy')

        return f_func, A_func, B_func 

def testAlgebraInverse(self):
    dataset_objects = algorithmic_math.math_dataset_init(26)
    counter = 0
    for d in algorithmic_math.algebra_inverse(26, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      solve_var, expression = decoded_input.split(":")
      lhs, rhs = expression.split("=")

      # Solve for the solve-var.
      result = sympy.solve("%s-(%s)" % (lhs, rhs), solve_var)
      target_expression = dataset_objects.int_decoder(d["targets"])

      # Check that the target and sympy's solutions are equivalent.
      self.assertEqual(
          0, sympy.simplify(str(result[0]) + "-(%s)" % target_expression))
    self.assertEqual(counter, 10) 

def simplify_surd(value, sample_args, context=None):
  """E.g., "Simplify (2 + 5*sqrt(3))**2."."""
  del value  # unused
  if context is None:
    context = composition.Context()

  entropy, sample_args = sample_args.peel()

  while True:
    base = random.randint(2, 20)
    if sympy.Integer(base).is_prime:
      break
  num_primes_less_than_20 = 8
  entropy -= math.log10(num_primes_less_than_20)
  exp = _sample_surd(base, entropy, max_power=2, multiples_only=False)
  simplified = sympy.expand(sympy.simplify(exp))

  template = random.choice([
      'Simplify {exp}.',
  ])
  return example.Problem(
      question=example.question(context, template, exp=exp),
      answer=simplified) 

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 generate_algebra_simplify_sample(vlist, ops, min_depth, max_depth):
  """Randomly generate an algebra simplify dataset sample.

  Given an input expression, produce the simplified expression.

  Args:
    vlist: Variable list. List of chars that can be used in the expression.
    ops: List of ExprOp instances. The allowed operators for the expression.
    min_depth: Expression trees will not have a smaller depth than this. 0 means
        there is just a variable. 1 means there is one operation.
    max_depth: Expression trees will not have a larger depth than this. To make
        all trees have the same depth, set this equal to `min_depth`.

  Returns:
    sample: String representation of the input.
    target: String representation of the solution.
  """
  depth = random.randrange(min_depth, max_depth + 1)
  expr = random_expr(depth, vlist, ops)

  sample = str(expr)
  target = format_sympy_expr(sympy.simplify(sample))
  return sample, target 

def testAlgebraInverse(self):
    dataset_objects = algorithmic_math.math_dataset_init(26)
    counter = 0
    for d in algorithmic_math.algebra_inverse(26, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      solve_var, expression = decoded_input.split(":")
      lhs, rhs = expression.split("=")

      # Solve for the solve-var.
      result = sympy.solve("%s-(%s)" % (lhs, rhs), solve_var)
      target_expression = dataset_objects.int_decoder(d["targets"])

      # Check that the target and sympy's solutions are equivalent.
      self.assertEqual(
          0, sympy.simplify(str(result[0]) + "-(%s)" % target_expression))
    self.assertEqual(counter, 10) 

def generate_algebra_simplify_sample(vlist, ops, min_depth, max_depth):
  """Randomly generate an algebra simplify dataset sample.

  Given an input expression, produce the simplified expression.

  Args:
    vlist: Variable list. List of chars that can be used in the expression.
    ops: List of ExprOp instances. The allowed operators for the expression.
    min_depth: Expression trees will not have a smaller depth than this. 0 means
        there is just a variable. 1 means there is one operation.
    max_depth: Expression trees will not have a larger depth than this. To make
        all trees have the same depth, set this equal to `min_depth`.

  Returns:
    sample: String representation of the input.
    target: String representation of the solution.
  """
  depth = random.randrange(min_depth, max_depth + 1)
  expr = random_expr(depth, vlist, ops)

  sample = str(expr)
  target = format_sympy_expr(sympy.simplify(sample))
  return sample, target 

def testAlgebraInverse(self):
    dataset_objects = algorithmic_math.math_dataset_init(26)
    counter = 0
    for d in algorithmic_math.algebra_inverse(26, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      solve_var, expression = decoded_input.split(":")
      lhs, rhs = expression.split("=")

      # Solve for the solve-var.
      result = sympy.solve("%s-(%s)" % (lhs, rhs), solve_var)
      target_expression = dataset_objects.int_decoder(d["targets"])

      # Check that the target and sympy's solutions are equivalent.
      self.assertEqual(
          0, sympy.simplify(str(result[0]) + "-(%s)" % target_expression))
    self.assertEqual(counter, 10) 

def testCalculusIntegrate(self):
    dataset_objects = algorithmic_math.math_dataset_init(
        8, digits=5, functions={"log": "L"})
    counter = 0
    for d in algorithmic_math.calculus_integrate(8, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      var, expression = decoded_input.split(":")
      target = dataset_objects.int_decoder(d["targets"])

      for fn_name, fn_char in six.iteritems(dataset_objects.functions):
        target = target.replace(fn_char, fn_name)

      # Take the derivative of the target.
      derivative = str(sympy.diff(target, var))

      # Check that the derivative of the integral equals the input.
      self.assertEqual(0, sympy.simplify("%s-(%s)" % (expression, derivative)))
    self.assertEqual(counter, 10) 

def get_equations(self):
        """
        :return: Functions to calculate A, B and f given state x and input u
        """
        f = sp.zeros(6, 1)

        x = sp.Matrix(sp.symbols('rx ry vx vy t w', real=True))
        u = sp.Matrix(sp.symbols('gimbal T', real=True))

        f[0, 0] = x[2, 0]
        f[1, 0] = x[3, 0]
        f[2, 0] = 1 / self.m * sp.sin(x[4, 0] + u[0, 0]) * u[1, 0]
        f[3, 0] = 1 / self.m * (sp.cos(x[4, 0] + u[0, 0]) * u[1, 0] - self.m * self.g)
        f[4, 0] = x[5, 0]
        f[5, 0] = 1 / self.I * (-sp.sin(u[0, 0]) * u[1, 0] * self.r_T)

        f = sp.simplify(f)
        A = sp.simplify(f.jacobian(x))
        B = sp.simplify(f.jacobian(u))

        f_func = sp.lambdify((x, u), f, 'numpy')
        A_func = sp.lambdify((x, u), A, 'numpy')
        B_func = sp.lambdify((x, u), B, 'numpy')

        return f_func, A_func, B_func 

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 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)).simplify())
    print('grad^B =',grad^B) 

def test_check_generalized_BAC_CAB_formulas(self):

        a, b, c, d, e = Ga('a b c d e').mv()

        assert str(a|(b*c)) == '-(a.c)*b + (a.b)*c'
        assert str(a|(b^c)) == '-(a.c)*b + (a.b)*c'
        assert str(a|(b^c^d)) == '(a.d)*b^c - (a.c)*b^d + (a.b)*c^d'

        expr = (a|(b^c))+(c|(a^b))+(b|(c^a)) # = (a.b)*c - (b.c)*a - ((a.b)*c - (b.c)*a)
        assert str(expr.simplify()) == '0'

        assert str(a*(b^c)-b*(a^c)+c*(a^b)) == '3*a^b^c'
        assert str(a*(b^c^d)-b*(a^c^d)+c*(a^b^d)-d*(a^b^c)) == '4*a^b^c^d'
        assert str((a^b)|(c^d)) == '-(a.c)*(b.d) + (a.d)*(b.c)'
        assert str(((a^b)|c)|d) == '-(a.c)*(b.d) + (a.d)*(b.c)'
        assert str(Ga.com(a^b, c^d)) == '-(b.d)*a^c + (b.c)*a^d + (a.d)*b^c - (a.c)*b^d'
        assert str((a|(b^c))|(d^e)) == '(-(a.b)*(c.e) + (a.c)*(b.e))*d + ((a.b)*(c.d) - (a.c)*(b.d))*e' 

def test_derivatives_in_spherical_coordinates(self):

        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)

        assert str(f) == 'f'
        assert str(A) == 'A__r*e_r + A__theta*e_theta + A__phi*e_phi'
        assert str(B) == 'B__rtheta*e_r^e_theta + B__rphi*e_r^e_phi + B__thetaphi*e_theta^e_phi'

        assert str(grad*f) == 'D{r}f*e_r + D{theta}f*e_theta/r + D{phi}f*e_phi/(r*sin(theta))'
        assert str((grad|A).simplify()) == '(r*D{r}A__r + 2*A__r + A__theta/tan(theta) + D{theta}A__theta + D{phi}A__phi/sin(theta))/r'
        assert str(-s3d.I()*(grad^A)) == '(A__phi/tan(theta) + D{theta}A__phi - D{phi}A__theta/sin(theta))*e_r/r + (-r*D{r}A__phi - A__phi + D{phi}A__r/sin(theta))*e_theta/r + (r*D{r}A__theta + A__theta - D{theta}A__r)*e_phi/r'

        assert latex(grad) == r'\boldsymbol{e}_{r} \frac{\partial}{\partial r} + \boldsymbol{e}_{\theta } \frac{1}{r} \frac{\partial}{\partial \theta } + \boldsymbol{e}_{\phi } \frac{1}{r \sin{\left (\theta  \right )}} \frac{\partial}{\partial \phi }'
        assert latex(B|(eth^ephi)) == r'- B^{\theta \phi } {\left (r,\theta ,\phi  \right )}'

        assert str(grad^B) == '(r*D{r}B__thetaphi - B__rphi/tan(theta) + 2*B__thetaphi - D{theta}B__rphi + D{phi}B__rtheta/sin(theta))*e_r^e_theta^e_phi/r' 

def generate_algebra_simplify_sample(vlist, ops, min_depth, max_depth):
  """Randomly generate an algebra simplify dataset sample.

  Given an input expression, produce the simplified expression.

  Args:
    vlist: Variable list. List of chars that can be used in the expression.
    ops: List of ExprOp instances. The allowed operators for the expression.
    min_depth: Expression trees will not have a smaller depth than this. 0 means
        there is just a variable. 1 means there is one operation.
    max_depth: Expression trees will not have a larger depth than this. To make
        all trees have the same depth, set this equal to `min_depth`.

  Returns:
    sample: String representation of the input.
    target: String representation of the solution.
  """
  depth = random.randrange(min_depth, max_depth + 1)
  expr = random_expr(depth, vlist, ops)

  sample = str(expr)
  target = format_sympy_expr(sympy.simplify(sample))
  return sample, target 

def testCalculusIntegrate(self):
    dataset_objects = algorithmic_math.math_dataset_init(
        8, digits=5, functions={"log": "L"})
    counter = 0
    for d in algorithmic_math.calculus_integrate(8, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      var, expression = decoded_input.split(":")
      target = dataset_objects.int_decoder(d["targets"])

      for fn_name, fn_char in six.iteritems(dataset_objects.functions):
        target = target.replace(fn_char, fn_name)

      # Take the derivative of the target.
      derivative = str(sympy.diff(target, var))

      # Check that the derivative of the integral equals the input.
      self.assertEqual(0, sympy.simplify("%s-(%s)" % (expression, derivative)))
    self.assertEqual(counter, 10) 

def generate_algebra_simplify_sample(vlist, ops, min_depth, max_depth):
  """Randomly generate an algebra simplify dataset sample.

  Given an input expression, produce the simplified expression.

  See go/symbolic-math-dataset.

  Args:
    vlist: Variable list. List of chars that can be used in the expression.
    ops: List of ExprOp instances. The allowed operators for the expression.
    min_depth: Expression trees will not have a smaller depth than this. 0 means
        there is just a variable. 1 means there is one operation.
    max_depth: Expression trees will not have a larger depth than this. To make
        all trees have the same depth, set this equal to `min_depth`.

  Returns:
    sample: String representation of the input.
    target: String representation of the solution.
  """
  depth = random.randrange(min_depth, max_depth + 1)
  expr = random_expr(depth, vlist, ops)

  sample = str(expr)
  target = format_sympy_expr(sympy.simplify(sample))
  return sample, target 

def grad_log_norm_symbolic(self):
        import sympy
        D = self.dimension
        X = self.eigenvalues_symbol
        B = [1] * D
        for d in range(D):
            for dd in range(D):
                if d != dd:
                    B[d] = B[d] * (X[d] - X[dd])
        B = [1 / b for b in B]

        p_D = sympy.pi ** D

        tmp = [b * sympy.exp(x_) for x_, b in zip(X, B)]
        tmp = sum(tmp)
        symbolic_norm_for_bingham = 2 * p_D * tmp

        return [
            sympy.simplify(sympy.diff(
                sympy.log(symbolic_norm_for_bingham),
                x_
            ))
            for x_ in X
        ] 

def simplifyEquation(input_str):
    """
    Only simplify the equation if there is no obvious problem
    Equation is not simplified if it includes the following terms:
        exp
        log
    """
    s_list = list()
    # these are considered the "dangerous" operation that will
    # overflow/underflow in np
    s_list.append(len(input_str.atoms(exp)))
    s_list.append(len(input_str.atoms(log)))

    if np.sum(s_list) != 0:
        # it is dangerous to simplify!
        return input_str, True
    else:
        #TODO: Removed actual simplyify (do we need it?)
        return input_str, False 

def pureTransitionToOde(A):
    """
    Get the ode from a pure transition matrix

    Parameters
    ----------
    A: `sympy.Matrix`
        a transition matrix of size [n \times n]

    Returns
    -------
    b: `sympy.Matrix`
        a matrix of size [n \times 1] which is the ode
    """
    nrow, ncol = A.shape
    assert nrow == ncol, "Need a square matrix"
    B = [sum(A[:, i]) - sum(A[i, :]) for i in range(nrow)]
    return sympy.simplify(sympy.Matrix(B)) 

def test_hard(self):
        # the SLIARD model is considered to be hard because a state can
        # go to multiple state.  This is not as hard as the SEIHFR model
        # below.
        state_list = ['S', 'L', 'I', 'A', 'R', 'D']
        param_list = ['beta', 'p', 'kappa', 'alpha', 'f', 'delta', 'epsilon', 'N']
        ode_list = [
            Transition('S', '- beta * S/N * ( I + delta * A)', 'ODE'),
            Transition('L', 'beta * S/N * (I + delta * A) - kappa * L', 'ODE'),
            Transition('I', 'p * kappa * L - alpha * I', 'ODE'),
            Transition('A', '(1-p) * kappa * L - epsilon * A', 'ODE'),
            Transition('R', 'f * alpha * I + epsilon * A', 'ODE'),
            Transition('D', '(1-f) * alpha * I', 'ODE')
            ]

        ode = SimulateOde(state_list, param_list, ode=ode_list)

        ode2 = ode.get_unrolled_obj()
        diffEqZero = map(lambda x: x==0, sympy.simplify(ode.get_ode_eqn() - ode2.get_ode_eqn()))

        self.assertTrue(numpy.all(numpy.array(list(diffEqZero)))) 

def test_bd(self):
        state_list = ['S', 'I', 'R']
        param_list = ['beta', 'gamma', 'B', 'mu']
        ode_list = [
            Transition(origin='S',
                       equation='-beta * S * I + B - mu * S',
                       transition_type=TransitionType.ODE),
            Transition(origin='I',
                       equation='beta * S * I - gamma * I - mu * I',
                       transition_type=TransitionType.ODE),
            Transition(origin='R',
                       destination='R',
                       equation='gamma * I',
                       transition_type=TransitionType.ODE)
            ]

        ode = SimulateOde(state_list, param_list, ode=ode_list)

        ode2 = ode.get_unrolled_obj()
        diffEqZero = map(lambda x: x==0, sympy.simplify(ode.get_ode_eqn() - ode2.get_ode_eqn()))

        self.assertTrue(numpy.all(numpy.array(list(diffEqZero)))) 

def test_derived_param(self):
        # the derived parameters are treated separately when compared to the
        # normal parameters and the odes
        ode = common_models.Legrand_Ebola_SEIHFR()

        ode_list = [
            Transition('S', '-(beta_I*S*I + beta_H_Time*S*H + beta_F_Time*S*F)'),
            Transition('E', '(beta_I*S*I + beta_H_Time*S*H + beta_F_Time*S*F) - alpha*E'),
            Transition('I', '-gamma_I*(1 - theta_1)*(1 - delta_1)*I - gamma_D*(1 - theta_1)*delta_1*I - gamma_H*theta_1*I + alpha*E'),
            Transition('H', 'gamma_H*theta_1*I - gamma_DH*delta_2*H - gamma_IH*(1 - delta_2)*H'),
            Transition('F', '- gamma_F*F + gamma_DH*delta_2*H + gamma_D*(1 - theta_1)*delta_1*I'),
            Transition('R', 'gamma_I*(1 - theta_1)*(1 - delta_1)*I + gamma_F*F + gamma_IH*(1 - delta_2)*H'),
            Transition('tau', '1')
        ]

        ode1 = SimulateOde(ode.state_list, ode.param_list, ode._derivedParamEqn, ode=ode_list)

        ode2 = ode1.get_unrolled_obj()
        diffEqZero = map(lambda x: x==0, sympy.simplify(ode.get_ode_eqn() - ode2.get_ode_eqn()))

        self.assertTrue(numpy.all(numpy.array(list(diffEqZero)))) 

def test_simple(self):
        """
        This actually only test the internal consistency of the functions
        rather than the actual correctness of the result.  Nevertheless,
        it is a valid test and we will use this hack for now because
        testing the result of a sympy object against strings is now the
        simplest procedure.
        """
        ode = common_models.SIR_Birth_Death()
        disease_state = ['I']
        R0 = epi_analysis.R0(ode, ['I'])

        F, V = epi_analysis.disease_progression_matrices(ode, disease_state)
        e = epi_analysis.R0_from_matrix(F, V)
        dfe = epi_analysis.DFE(ode, ['I'])
        self.assertTrue(sympy.simplify(R0 - e[0].subs(dfe)) == 0) 

def _scheme_from_rc_sympy(alpha, beta):
    # Construct the triadiagonal matrix [sqrt(beta), alpha, sqrt(beta)]
    A = _sympy_tridiag(alpha, [sympy.sqrt(bta) for bta in beta])

    # Extract points and weights from eigenproblem
    x = []
    w = []
    for item in A.eigenvects():
        val, multiplicity, vec = item
        assert multiplicity == 1
        assert len(vec) == 1
        vec = vec[0]
        x.append(val)
        norm2 = sum([v ** 2 for v in vec])
        # simplifiction takes really long
        # w.append(sympy.simplify(beta[0] * vec[0]**2 / norm2))
        w.append(beta[0] * vec[0] ** 2 / norm2)
    # sort by x
    order = sorted(range(len(x)), key=lambda i: x[i])
    x = [x[i] for i in order]
    w = [w[i] for i in order]
    return x, w 

def test_integrate():
    moments = quadpy.tools.integrate(lambda x: [x ** k for k in range(5)], -1, +1)
    assert (moments == [2, 0, sympy.S(2) / 3, 0, sympy.S(2) / 5]).all()

    moments = quadpy.tools.integrate(
        lambda x: orthopy.line_segment.tree_legendre(x, 4, "monic", symbolic=True),
        -1,
        +1,
    )
    assert (moments == [2, 0, 0, 0, 0]).all()

    # Example from Gautschi's "How to and how not to" article
    moments = quadpy.tools.integrate(
        lambda x: [x ** k * sympy.exp(-(x ** 3) / 3) for k in range(5)], 0, sympy.oo
    )
    S = numpy.vectorize(sympy.S)
    gamma = numpy.vectorize(sympy.gamma)
    n = numpy.arange(5)
    reference = 3 ** (S(n - 2) / 3) * gamma(S(n + 1) / 3)
    assert numpy.all([sympy.simplify(m - r) == 0 for m, r in zip(moments, reference)]) 

def symbolize(flt: float) -> sympy.Symbol:
    """Attempt to convert a real number into a simpler symbolic
    representation.

    Returns:
        A sympy Symbol. (Convert to string with str(sym) or to latex with
            sympy.latex(sym)
    Raises:
        ValueError:     If cannot simplify float
    """
    try:
        ratio = rationalize(flt)
        res = sympy.simplify(ratio)
    except ValueError:
        ratio = rationalize(flt/np.pi)
        res = sympy.simplify(ratio) * sympy.pi
    return res

# fin 

def replace_pow(self, body, symexpr, powexprs, expr, exp):
        if exp == 1:
            return (symexpr, None)
        elif not (expr.base, exp) in powexprs:
            if exp == 2:
                operand = sp.simplify(expr.base)
                value = BinOp("Mult", self.sympyToAst.visit(operand), self.sympyToAst.visit(operand))
            elif exp % 2 == 1:
                _, operand = self.replace_pow(body, symexpr, powexprs, expr, exp - 1)
                value = BinOp("Mult", self.symbols[operand], self.sympyToAst.visit(sp.simplify(expr.base)))
            else:
                _, operand = self.replace_pow(body, symexpr, powexprs, expr, exp / 2)
                value = BinOp("Mult", self.symbols[operand], self.symbols[operand])
            name, self.next = "__sp{0}".format(self.next), self.next + 1
            self.symbols[name] = Variable(name, copy.deepcopy(value.shape), value.dtype)
            body.append(Assign(self.symbols[name], value))
            powexprs[(expr.base, exp)] = name
        if np.abs(expr.exp.p) == exp:
            symbol = sp.Symbol(powexprs[(expr.base, exp)])
            symexpr = symexpr.subs(expr, self.symbols[powexprs[(expr.base, exp)]].dtype(1) / symbol if expr.exp.is_negative else symbol)
        return (symexpr, powexprs[(expr.base, exp)]) 

def _build_logical_expression(self, grammar, terminal_component_names):
        terminal_component_symbols = eval("symbols('%s')"%(' '.join(terminal_component_names)))
        if isinstance(terminal_component_symbols, Symbol):
            terminal_component_symbols = [terminal_component_symbols]
        name_to_symbol = {terminal_component_names[i]:symbol for i, symbol in enumerate(terminal_component_symbols)}
        terminal_component_names = set(terminal_component_names)

        op_to_symbolic_operation = {'not':operator.invert, 'concat':operator.and_, 'gap':operator.and_, 'union':operator.or_, 'intersect':operator.and_}

        def logical_expression_builder(component):
            if component['id'] in terminal_component_names:
                return name_to_symbol[component['id']]
            else:
                children = component['components']
                return reduce(op_to_symbolic_operation[component['operation']],[logical_expression_builder(child) for child in children])

        return simplify(logical_expression_builder(grammar)) 

def test_cu3():
    E = eye(2)
    UPPER = Matrix([[1, 0], [0, 0]])
    LOWER = Matrix([[0, 0], [0, 1]])
    theta, phi, lambd = symbols("theta phi lambd")
    U = Circuit().rz(lambd)[0].ry(theta)[0].rz(phi)[0].run_with_sympy_unitary()

    actual_1 = Circuit().cu3(theta, phi, lambd)[0, 1].run(backend="sympy_unitary")
    expected_1 = reduce(TensorProduct, [E, UPPER]) + reduce(TensorProduct, [U, LOWER])
    print("actual")
    print(simplify(actual_1))
    print("expected")
    print(simplify(expected_1))
    print("diff")
    print(simplify(actual_1 - expected_1))
    assert simplify(actual_1 - expected_1) == zeros(4)