Python sympy 模块，simplify() 实例源码

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 variance(X, condition=None, **kwargs):
    """
    Variance of a random expression

    Expectation of (X-E(X))**2

    Examples
    ========

    >>> from sympy.stats import Die, E, Bernoulli, variance
    >>> from sympy import simplify, Symbol

    >>> X = Die('X', 6)
    >>> p = Symbol('p')
    >>> B = Bernoulli('B', p, 1, 0)

    >>> variance(2*X)
    35/3

    >>> simplify(variance(B))
    p*(-p + 1)
    """
    return cmoment(X, 2, condition, **kwargs)

def standard_deviation(X, condition=None, **kwargs):
    """
    Standard Deviation of a random expression

    Square root of the Expectation of (X-E(X))**2

    Examples
    ========

    >>> from sympy.stats import Bernoulli, std
    >>> from sympy import Symbol, simplify

    >>> p = Symbol('p')
    >>> B = Bernoulli('B', p, 1, 0)

    >>> simplify(std(B))
    sqrt(p*(-p + 1))
    """
    return sqrt(variance(X, condition, **kwargs))

def test_karr_proposition_2a():
    # Test Karr, page 309, proposition 2, part a
    i = Symbol("i", integer=True)
    u = Symbol("u", integer=True)
    v = Symbol("v", integer=True)

    def test_the_product(m, n):
        # g
        g = i**3 + 2*i**2 - 3*i
        # f = Delta g
        f = simplify(g.subs(i, i+1) / g)
        # The product
        a = m
        b = n - 1
        P = Product(f, (i, a, b)).doit()
        # Test if Product_{m <= i < n} f(i) = g(n) / g(m)
        assert simplify(P / (g.subs(i, n) / g.subs(i, m))) == 1

    # m < n
    test_the_product(u,   u+v)
    # m = n
    test_the_product(u,   u  )
    # m > n
    test_the_product(u+v, u  )

def test_simplify():
    y, t, b, c = symbols('y, t, b, c', integer = True)

    assert simplify(Product(x*y, (x, n, m), (y, a, k)) * \
        Product(y, (x, n, m), (y, a, k))) == \
            Product(x*y**2, (x, n, m), (y, a, k))
    assert simplify(3 * y* Product(x, (x, n, m)) * Product(x, (x, m + 1, a))) \
        == 3 * y * Product(x, (x, n, a))
    assert simplify(Product(x, (x, k + 1, a)) * Product(x, (x, n, k))) == \
        Product(x, (x, n, a))
    assert simplify(Product(x, (x, k + 1, a)) * Product(x + 1, (x, n, k))) == \
        Product(x, (x, k + 1, a)) * Product(x + 1, (x, n, k))
    assert simplify(Product(x, (t, a, b)) * Product(y, (t, a, b)) * \
        Product(x, (t, b+1, c))) == Product(x*y, (t, a, b)) * \
            Product(x, (t, b+1, c))
    assert simplify(Product(x, (t, a, b)) * Product(x, (t, b+1, c)) * \
        Product(y, (t, a, b))) == Product(x*y, (t, a, b)) * \
            Product(x, (t, b+1, c))

def test_karr_proposition_2a():
    # Test Karr, page 309, proposition 2, part a
    i = Symbol("i", integer=True)
    u = Symbol("u", integer=True)
    v = Symbol("v", integer=True)

    def test_the_sum(m, n):
        # g
        g = i**3 + 2*i**2 - 3*i
        # f = Delta g
        f = simplify(g.subs(i, i+1) - g)
        # The sum
        a = m
        b = n - 1
        S = Sum(f, (i, a, b)).doit()
        # Test if Sum_{m <= i < n} f(i) = g(n) - g(m)
        assert simplify(S - (g.subs(i, n) - g.subs(i, m))) == 0

    # m < n
    test_the_sum(u,   u+v)
    # m = n
    test_the_sum(u,   u  )
    # m > n
    test_the_sum(u+v, u  )

def test_wavefunction():
    a = 1/Z
    R = {
        (1, 0): 2*sqrt(1/a**3) * exp(-r/a),
        (2, 0): sqrt(1/(2*a**3)) * exp(-r/(2*a)) * (1 - r/(2*a)),
        (2, 1): S(1)/2 * sqrt(1/(6*a**3)) * exp(-r/(2*a)) * r/a,
        (3, 0): S(2)/3 * sqrt(1/(3*a**3)) * exp(-r/(3*a)) *
        (1 - 2*r/(3*a) + S(2)/27 * (r/a)**2),
        (3, 1): S(4)/27 * sqrt(2/(3*a**3)) * exp(-r/(3*a)) *
        (1 - r/(6*a)) * r/a,
        (3, 2): S(2)/81 * sqrt(2/(15*a**3)) * exp(-r/(3*a)) * (r/a)**2,
        (4, 0): S(1)/4 * sqrt(1/a**3) * exp(-r/(4*a)) *
        (1 - 3*r/(4*a) + S(1)/8 * (r/a)**2 - S(1)/192 * (r/a)**3),
        (4, 1): S(1)/16 * sqrt(5/(3*a**3)) * exp(-r/(4*a)) *
        (1 - r/(4*a) + S(1)/80 * (r/a)**2) * (r/a),
        (4, 2): S(1)/64 * sqrt(1/(5*a**3)) * exp(-r/(4*a)) *
        (1 - r/(12*a)) * (r/a)**2,
        (4, 3): S(1)/768 * sqrt(1/(35*a**3)) * exp(-r/(4*a)) * (r/a)**3,
    }
    for n, l in R:
        assert simplify(R_nl(n, l, r, Z) - R[(n, l)]) == 0

def test_simplify():
    f, n = symbols('f, n')

    m = Matrix([[1, x], [x + 1/x, x - 1]])
    m = m.row_join(eye(m.cols))
    raw = m.rref(simplify=lambda x: x)[0]
    assert raw != m.rref(simplify=True)[0]

    M = Matrix([[            1/x + 1/y,                 (x + x*y) / x  ],
                [ (f(x) + y*f(x))/f(x), 2 * (1/n - cos(n * pi)/n) / pi ]])
    M.simplify()
    assert M == Matrix([[ (x + y)/(x * y),                        1 + y ],
                        [           1 + y, 2*((1 - 1*cos(pi*n))/(pi*n)) ]])
    eq = (1 + x)**2
    M = Matrix([[eq]])
    M.simplify()
    assert M == Matrix([[eq]])
    M.simplify(ratio=oo) == M
    assert M == Matrix([[eq.simplify(ratio=oo)]])

def test_invertible_check():
    # sometimes a singular matrix will have a pivot vector shorter than
    # the number of rows in a matrix...
    assert Matrix([[1, 2], [1, 2]]).rref() == (Matrix([[1, 2], [0, 0]]), [0])
    raises(ValueError, lambda: Matrix([[1, 2], [1, 2]]).inv())
    # ... but sometimes it won't, so that is an insufficient test of
    # whether something is invertible.
    m = Matrix([
        [-1, -1,  0],
        [ x,  1,  1],
        [ 1,  x, -1],
    ])
    assert len(m.rref()[1]) == m.rows
    # in addition, unless simplify=True in the call to rref, the identity
    # matrix will be returned even though m is not invertible
    assert m.rref()[0] == eye(3)
    assert m.rref(simplify=signsimp)[0] != eye(3)
    raises(ValueError, lambda: m.inv(method="ADJ"))
    raises(ValueError, lambda: m.inv(method="GE"))
    raises(ValueError, lambda: m.inv(method="LU"))

def test_anti_symmetric():
    assert Matrix([1, 2]).is_anti_symmetric() is False
    m = Matrix(3, 3, [0, x**2 + 2*x + 1, y, -(x + 1)**2, 0, x*y, -y, -x*y, 0])
    assert m.is_anti_symmetric() is True
    assert m.is_anti_symmetric(simplify=False) is False
    assert m.is_anti_symmetric(simplify=lambda x: x) is False

    # tweak to fail
    m[2, 1] = -m[2, 1]
    assert m.is_anti_symmetric() is False
    # untweak
    m[2, 1] = -m[2, 1]

    m = m.expand()
    assert m.is_anti_symmetric(simplify=False) is True
    m[0, 0] = 1
    assert m.is_anti_symmetric() is False

def test_pinv():
    from sympy.abc import a, b, c, d, e, f
    # Pseudoinverse of an invertible matrix is the inverse.
    A1 = Matrix([[a, b], [c, d]])
    assert simplify(A1.pinv()) == simplify(A1.inv())
    # Test the four properties of the pseudoinverse for various matrices.
    As = [Matrix([[13, 104], [2212, 3], [-3, 5]]),
          Matrix([[1, 7, 9], [11, 17, 19]]),
          Matrix([a, b])]
    for A in As:
        A_pinv = A.pinv()
        AAp = A * A_pinv
        ApA = A_pinv * A
        assert simplify(AAp * A) == A
        assert simplify(ApA * A_pinv) == A_pinv
        assert AAp.H == AAp
        assert ApA.H == ApA

def power_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if symbol not in exp.free_symbols and isinstance(base, sympy.Symbol):
        if sympy.simplify(exp + 1) == 0:
            return ReciprocalRule(base, integrand, symbol)
        return PowerRule(base, exp, integrand, symbol)
    elif symbol not in base.free_symbols and isinstance(exp, sympy.Symbol):
        rule = ExpRule(base, exp, integrand, symbol)

        if sympy.ask(~sympy.Q.zero(sympy.log(base))):
            return rule
        elif sympy.ask(sympy.Q.zero(sympy.log(base))):
            return ConstantRule(1, 1, symbol)

        return PiecewiseRule([
            (ConstantRule(1, 1, symbol), sympy.Eq(sympy.log(base), 0)),
            (rule, True)
        ], integrand, symbol)

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():
    enhance_print()

    (ex, ey, ez) = MV.setup('e*x|y|z', metric='[1,1,1]')

    u = MV('u', 'vector')
    v = MV('v', 'vector')
    w = MV('w', 'vector')
    print(u)
    print(v)
    print(w)

    uv = u ^ v
    print(uv)
    print(uv.is_blade())
    uvw = u ^ v ^ w
    print(uvw)
    print(uvw.is_blade())

    print(simplify((uv*uv).scalar()))
    return

def compose_u3(theta1, phi1, lambda1, theta2, phi2, lambda2):
    """Return a triple theta, phi, lambda for the product.

    u3(theta, phi, lambda)
       = u3(theta1, phi1, lambda1).u3(theta2, phi2, lambda2)
       = Rz(phi1).Ry(theta1).Rz(lambda1+phi2).Ry(theta2).Rz(lambda2)
       = Rz(phi1).Rz(phi').Ry(theta').Rz(lambda').Rz(lambda2)
       = u3(theta', phi1 + phi', lambda2 + lambda')

    Return theta, phi, lambda.
    """
    # Careful with the factor of two in yzy_to_zyz
    thetap, phip, lambdap = yzy_to_zyz((lambda1 + phi2) / 2,
                                       theta1 / 2, theta2 / 2)
    (theta, phi, lamb) = (2 * thetap, phi1 + 2 * phip, lambda2 + 2 * lambdap)
    return (theta.simplify(), phi.simplify(), lamb.simplify())

def linearity_index_inverse_depth():
    """Linearity index of Inverse Depth Parameterization"""
    D, rho, rho_0, d_1, sigma_rho = sympy.symbols("D,rho,rho_0,d_1,sigma_rho")
    alpha = sympy.symbols("alpha")
    u = (rho * sympy.sin(alpha)) / (rho_0 * d_1 * (rho_0 - rho) + rho * sympy.cos(alpha))  # NOQA

    # first order derivative of u
    u_p = sympy.diff(u, rho)
    u_p = sympy.simplify(u_p)

    # second order derivative of u
    u_pp = sympy.diff(u_p, rho)
    u_pp = sympy.simplify(u_pp)

    # Linearity index
    L = (u_pp * 2 * sigma_rho) / (u_p)
    L = sympy.simplify(L)

    print()
    print("u: ", u)
    print("u': ", u_p)
    print("u'': ", u_pp)
    # print("L = ", L)
    print("L = ", L.subs(rho, 0))
    print()

def compute_stencil(order, n, x_value, h_value, x0=0.):
    """
    computes a stencil of Order order
    """
    h,x = symbols('h x')
    xs = [x+h*cos(i*pi/n) for i in range(n,-1,-1)]
    cs = finite_diff_weights(order, xs, x0)[order][n]
    cs = [simplify(c) for c in cs]

    cs = [simplify(expr.subs(h, h_value)) for expr in cs]
    xs = [simplify(expr.subs(h, h_value)) for expr in xs]

    cs = [simplify(expr.subs(x, x_value)) for expr in cs]
    xs = [simplify(expr.subs(x, x_value)) for expr in xs]

    return xs, cs

##############################################

def test_meijerg_formulae():
    from sympy.simplify.hyperexpand import MeijerFormulaCollection
    formulae = MeijerFormulaCollection().formulae
    for sig in formulae:
        for formula in formulae[sig]:
            g = meijerg(formula.func.an, formula.func.ap,
                        formula.func.bm, formula.func.bq,
                        formula.z)
            rep = {}
            for sym in formula.symbols:
                rep[sym] = randcplx()

            # first test if the closed-form is actually correct
            g = g.subs(rep)
            closed_form = formula.closed_form.subs(rep)
            z = formula.z
            assert tn(g, closed_form, z)

            # now test the computed matrix
            cl = (formula.C * formula.B)[0].subs(rep)
            assert tn(closed_form, cl, z)
            deriv1 = z*formula.B.diff(z)
            deriv2 = formula.M * formula.B
            for d1, d2 in zip(deriv1, deriv2):
                assert tn(d1.subs(rep), d2.subs(rep), z)

def test_simplify():
    y, t, b, c = symbols('y, t, b, c', integer = True)

    assert simplify(Product(x*y, (x, n, m), (y, a, k)) * \
        Product(y, (x, n, m), (y, a, k))) == \
            Product(x*y**2, (x, n, m), (y, a, k))
    assert simplify(3 * y* Product(x, (x, n, m)) * Product(x, (x, m + 1, a))) \
        == 3 * y * Product(x, (x, n, a))
    assert simplify(Product(x, (x, k + 1, a)) * Product(x, (x, n, k))) == \
        Product(x, (x, n, a))
    assert simplify(Product(x, (x, k + 1, a)) * Product(x + 1, (x, n, k))) == \
        Product(x, (x, k + 1, a)) * Product(x + 1, (x, n, k))
    assert simplify(Product(x, (t, a, b)) * Product(y, (t, a, b)) * \
        Product(x, (t, b+1, c))) == Product(x*y, (t, a, b)) * \
            Product(x, (t, b+1, c))
    assert simplify(Product(x, (t, a, b)) * Product(x, (t, b+1, c)) * \
        Product(y, (t, a, b))) == Product(x*y, (t, a, b)) * \
            Product(x, (t, b+1, c))

def test_gosper_sum_AeqB_part2():
    f2a = n**2*a**n
    f2b = (n - r/2)*binomial(r, n)
    f2c = factorial(n - 1)**2/(factorial(n - x)*factorial(n + x))

    g2a = -a*(a + 1)/(a - 1)**3 + a**(
        m + 1)*(a**2*m**2 - 2*a*m**2 + m**2 - 2*a*m + 2*m + a + 1)/(a - 1)**3
    g2b = (m - r)*binomial(r, m)/2
    ff = factorial(1 - x)*factorial(1 + x)
    g2c = 1/ff*(
        1 - 1/x**2) + factorial(m)**2/(x**2*factorial(m - x)*factorial(m + x))

    g = gosper_sum(f2a, (n, 0, m))
    assert g is not None and simplify(g - g2a) == 0
    g = gosper_sum(f2b, (n, 0, m))
    assert g is not None and simplify(g - g2b) == 0
    g = gosper_sum(f2c, (n, 1, m))
    assert g is not None and simplify(g - g2c) == 0

def test_karr_proposition_2a():
    # Test Karr, page 309, proposition 2, part a
    i = Symbol("i", integer=True)
    u = Symbol("u", integer=True)
    v = Symbol("v", integer=True)

    def test_the_sum(m, n):
        # g
        g = i**3 + 2*i**2 - 3*i
        # f = Delta g
        f = simplify(g.subs(i, i+1) - g)
        # The sum
        a = m
        b = n - 1
        S = Sum(f, (i, a, b)).doit()
        # Test if Sum_{m <= i < n} f(i) = g(n) - g(m)
        assert simplify(S - (g.subs(i, n) - g.subs(i, m))) == 0

    # m < n
    test_the_sum(u,   u+v)
    # m = n
    test_the_sum(u,   u  )
    # m > n
    test_the_sum(u+v, u  )

def test_hypersum():
    from sympy import sin
    assert simplify(summation(x**n/fac(n), (n, 1, oo))) == -1 + exp(x)
    assert summation((-1)**n * x**(2*n) / fac(2*n), (n, 0, oo)) == cos(x)
    assert simplify(summation((-1)**n*x**(2*n + 1) /
        factorial(2*n + 1), (n, 3, oo))) == -x + sin(x) + x**3/6 - x**5/120

    assert summation(1/(n + 2)**3, (n, 1, oo)) == -S(9)/8 + zeta(3)
    assert summation(1/n**4, (n, 1, oo)) == pi**4/90

    s = summation(x**n*n, (n, -oo, 0))
    assert s.is_Piecewise
    assert s.args[0].args[0] == -1/(x*(1 - 1/x)**2)
    assert s.args[0].args[1] == (abs(1/x) < 1)

    m = Symbol('n', integer=True, positive=True)
    assert summation(binomial(m, k), (k, 0, m)) == 2**m

def test_wavefunction():
    a = 1/Z
    R = {
        (1, 0): 2*sqrt(1/a**3) * exp(-r/a),
        (2, 0): sqrt(1/(2*a**3)) * exp(-r/(2*a)) * (1 - r/(2*a)),
        (2, 1): S(1)/2 * sqrt(1/(6*a**3)) * exp(-r/(2*a)) * r/a,
        (3, 0): S(2)/3 * sqrt(1/(3*a**3)) * exp(-r/(3*a)) *
        (1 - 2*r/(3*a) + S(2)/27 * (r/a)**2),
        (3, 1): S(4)/27 * sqrt(2/(3*a**3)) * exp(-r/(3*a)) *
        (1 - r/(6*a)) * r/a,
        (3, 2): S(2)/81 * sqrt(2/(15*a**3)) * exp(-r/(3*a)) * (r/a)**2,
        (4, 0): S(1)/4 * sqrt(1/a**3) * exp(-r/(4*a)) *
        (1 - 3*r/(4*a) + S(1)/8 * (r/a)**2 - S(1)/192 * (r/a)**3),
        (4, 1): S(1)/16 * sqrt(5/(3*a**3)) * exp(-r/(4*a)) *
        (1 - r/(4*a) + S(1)/80 * (r/a)**2) * (r/a),
        (4, 2): S(1)/64 * sqrt(1/(5*a**3)) * exp(-r/(4*a)) *
        (1 - r/(12*a)) * (r/a)**2,
        (4, 3): S(1)/768 * sqrt(1/(35*a**3)) * exp(-r/(4*a)) * (r/a)**3,
    }
    for n, l in R:
        assert simplify(R_nl(n, l, r, Z) - R[(n, l)]) == 0

def test_n_link_pendulum_on_cart_higher_order():
    l0, m0 = symbols("l0 m0")
    l1, m1 = symbols("l1 m1")
    m2 = symbols("m2")
    g = symbols("g")
    q0, q1, q2 = dynamicsymbols("q0 q1 q2")
    u0, u1, u2 = dynamicsymbols("u0 u1 u2")
    F, T1 = dynamicsymbols("F T1")

    kane1 = models.n_link_pendulum_on_cart(2)
    massmatrix1 = Matrix([[m0 + m1 + m2, -l0*m1*cos(q1) - l0*m2*cos(q1),
                           -l1*m2*cos(q2)],
                          [-l0*m1*cos(q1) - l0*m2*cos(q1), l0**2*m1 + l0**2*m2,
                           l0*l1*m2*(sin(q1)*sin(q2) + cos(q1)*cos(q2))],
                          [-l1*m2*cos(q2),
                           l0*l1*m2*(sin(q1)*sin(q2) + cos(q1)*cos(q2)),
                           l1**2*m2]])
    forcing1 = Matrix([[-l0*m1*u1**2*sin(q1) - l0*m2*u1**2*sin(q1) -
                        l1*m2*u2**2*sin(q2) + F],
                       [g*l0*m1*sin(q1) + g*l0*m2*sin(q1) -
                        l0*l1*m2*(sin(q1)*cos(q2) - sin(q2)*cos(q1))*u2**2],
                       [g*l1*m2*sin(q2) - l0*l1*m2*(-sin(q1)*cos(q2) +
                                                    sin(q2)*cos(q1))*u1**2]])
    assert simplify(massmatrix1 - kane1.mass_matrix) == zeros(3)
    assert simplify(forcing1 - kane1.forcing) == Matrix([0, 0, 0])

def test_pend():
    q, u = dynamicsymbols('q u')
    qd, ud = dynamicsymbols('q u', 1)
    m, l, g = symbols('m l g')
    N = ReferenceFrame('N')
    P = Point('P')
    P.set_vel(N, -l * u * sin(q) * N.x + l * u * cos(q) * N.y)
    kd = [qd - u]

    FL = [(P, m * g * N.x)]
    pa = Particle('pa', P, m)
    BL = [pa]

    KM = KanesMethod(N, [q], [u], kd)
    with warnings.catch_warnings():
        warnings.filterwarnings("ignore", category=SymPyDeprecationWarning)
        KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    rhs.simplify()
    assert expand(rhs[0]) == expand(-g / l * sin(q))
    assert simplify(KM.rhs() -
                    KM.mass_matrix_full.LUsolve(KM.forcing_full)) == zeros(2, 1)

def test_form_1():
    symsystem1 = SymbolicSystem(states, comb_explicit_rhs,
                                alg_con=alg_con_full, output_eqns=out_eqns,
                                coord_idxs=coord_idxs, speed_idxs=speed_idxs,
                                bodies=bodies, loads=loads)

    assert symsystem1.coordinates == Matrix([x, y])
    assert symsystem1.speeds == Matrix([u, v])
    assert symsystem1.states == Matrix([x, y, u, v, lam])

    assert symsystem1.alg_con == [4]

    inter = comb_explicit_rhs
    assert simplify(symsystem1.comb_explicit_rhs - inter) == zeros(5, 1)

    assert set(symsystem1.dynamic_symbols()) == set([y, v, lam, u, x])
    assert type(symsystem1.dynamic_symbols()) == tuple
    assert set(symsystem1.constant_symbols()) == set([l, g, m])
    assert type(symsystem1.constant_symbols()) == tuple

    assert symsystem1.output_eqns == out_eqns

    assert symsystem1.bodies == (Pa,)
    assert symsystem1.loads == ((P, g * m * N.x),)

def test_simplify():
    f, n = symbols('f, n')

    m = Matrix([[1, x], [x + 1/x, x - 1]])
    m = m.row_join(eye(m.cols))
    raw = m.rref(simplify=lambda x: x)[0]
    assert raw != m.rref(simplify=True)[0]

    M = Matrix([[            1/x + 1/y,                 (x + x*y) / x  ],
                [ (f(x) + y*f(x))/f(x), 2 * (1/n - cos(n * pi)/n) / pi ]])
    M.simplify()
    assert M == Matrix([[ (x + y)/(x * y),                        1 + y ],
                        [           1 + y, 2*((1 - 1*cos(pi*n))/(pi*n)) ]])
    eq = (1 + x)**2
    M = Matrix([[eq]])
    M.simplify()
    assert M == Matrix([[eq]])
    M.simplify(ratio=oo) == M
    assert M == Matrix([[eq.simplify(ratio=oo)]])

def test_jordan_form_complex_issue_9274():
    A = Matrix([[ 2,  4,  1,  0],
                [-4,  2,  0,  1],
                [ 0,  0,  2,  4],
                [ 0,  0, -4,  2]])
    p = 2 - 4*I;
    q = 2 + 4*I;
    Jmust1 = Matrix([[p, 1, 0, 0],
                     [0, p, 0, 0],
                     [0, 0, q, 1],
                     [0, 0, 0, q]])
    Jmust2 = Matrix([[q, 1, 0, 0],
                     [0, q, 0, 0],
                     [0, 0, p, 1],
                     [0, 0, 0, p]])
    P, J = A.jordan_form()
    assert J == Jmust1 or J == Jmust2
    assert simplify(P*J*P.inv()) == A

def test_anti_symmetric():
    assert Matrix([1, 2]).is_anti_symmetric() is False
    m = Matrix(3, 3, [0, x**2 + 2*x + 1, y, -(x + 1)**2, 0, x*y, -y, -x*y, 0])
    assert m.is_anti_symmetric() is True
    assert m.is_anti_symmetric(simplify=False) is False
    assert m.is_anti_symmetric(simplify=lambda x: x) is False

    # tweak to fail
    m[2, 1] = -m[2, 1]
    assert m.is_anti_symmetric() is False
    # untweak
    m[2, 1] = -m[2, 1]

    m = m.expand()
    assert m.is_anti_symmetric(simplify=False) is True
    m[0, 0] = 1
    assert m.is_anti_symmetric() is False

def test_rank_regression_from_so():
    # see:
    # http://stackoverflow.com/questions/19072700/why-does-sympy-give-me-the-wrong-answer-when-i-row-reduce-a-symbolic-matrix

    nu, lamb = symbols('nu, lambda')
    A = Matrix([[-3*nu,         1,                  0,  0],
                [ 3*nu, -2*nu - 1,                  2,  0],
                [    0,      2*nu, (-1*nu) - lamb - 2,  3],
                [    0,         0,          nu + lamb, -3]])
    expected_reduced = Matrix([[1, 0, 0, 1/(nu**2*(-lamb - nu))],
                               [0, 1, 0,    3/(nu*(-lamb - nu))],
                               [0, 0, 1,         3/(-lamb - nu)],
                               [0, 0, 0,                      0]])
    expected_pivots = [0, 1, 2]

    reduced, pivots = A.rref()

    assert simplify(expected_reduced - reduced) == zeros(*A.shape)
    assert pivots == expected_pivots

def test_pinv():
    # Pseudoinverse of an invertible matrix is the inverse.
    A1 = Matrix([[a, b], [c, d]])
    assert simplify(A1.pinv()) == simplify(A1.inv())
    # Test the four properties of the pseudoinverse for various matrices.
    As = [Matrix([[13, 104], [2212, 3], [-3, 5]]),
          Matrix([[1, 7, 9], [11, 17, 19]]),
          Matrix([a, b])]
    for A in As:
        A_pinv = A.pinv()
        AAp = A * A_pinv
        ApA = A_pinv * A
        assert simplify(AAp * A) == A
        assert simplify(ApA * A_pinv) == A_pinv
        assert AAp.H == AAp
        assert ApA.H == ApA

def power_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if symbol not in exp.free_symbols and isinstance(base, sympy.Symbol):
        if sympy.simplify(exp + 1) == 0:
            return ReciprocalRule(base, integrand, symbol)
        return PowerRule(base, exp, integrand, symbol)
    elif symbol not in base.free_symbols and isinstance(exp, sympy.Symbol):
        rule = ExpRule(base, exp, integrand, symbol)

        if fuzzy_not(sympy.log(base).is_zero):
            return rule
        elif sympy.log(base).is_zero:
            return ConstantRule(1, 1, symbol)

        return PiecewiseRule([
            (ConstantRule(1, 1, symbol), sympy.Eq(sympy.log(base), 0)),
            (rule, True)
        ], integrand, symbol)

def equals(self, o):
        """
        Returns True if self and o are the same mathematical entities.

        Examples
        ========

        >>> from sympy import Plane, Point3D
        >>> a = Plane(Point3D(1, 2, 3), normal_vector=(1, 1, 1))
        >>> b = Plane(Point3D(1, 2, 3), normal_vector=(2, 2, 2))
        >>> c = Plane(Point3D(1, 2, 3), normal_vector=(-1, 4, 6))
        >>> a.equals(a)
        True
        >>> a.equals(b)
        True
        >>> a.equals(c)
        False
        """
        if isinstance(o, Plane):
            a = self.equation()
            b = o.equation()
            return simplify(a / b).is_constant()
        else:
            return False

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 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 _gauss_from_coefficients_sympy(alpha, beta):
    assert isinstance(alpha[0], sympy.Rational)
    # 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_normality(n=3):
    '''Make sure that the polynomials are normal.
    '''
    polar = sympy.Symbol('theta', real=True)
    azimuthal = sympy.Symbol('phi', real=True)
    tree = numpy.concatenate(
            orthopy.sphere.sph_tree(
                n, polar, azimuthal, normalization='quantum mechanic',
                symbolic=True
                ))

    for val in tree:
        integrand = sympy.simplify(
            val * sympy.conjugate(val) * sympy.sin(polar)
            )
        assert sympy.integrate(
            integrand,
            (azimuthal, 0, 2*pi), (polar, 0, pi)
            ) == 1
    return

def test_orthogonality(normalization, n=4):
    polar = sympy.Symbol('theta', real=True)
    azimuthal = sympy.Symbol('phi', real=True)
    tree = numpy.concatenate(
            orthopy.sphere.sph_tree(
                n, polar, azimuthal, normalization=normalization,
                symbolic=True
                ))
    vals = tree * sympy.conjugate(numpy.roll(tree, 1, axis=0))

    for val in vals:
        integrand = sympy.simplify(val * sympy.sin(polar))
        assert sympy.integrate(
            integrand,
            (azimuthal, 0, 2*pi), (polar, 0, pi)
            ) == 0
    return

def eqs_match(origeqs, origspecies, removed, neweqs, newspecies, debug = False):
    """Check that the original equations match the new equations after
    replacing the species with the expression in removed."""
    origeqsdict = dict((origspecies[s], origeqs[s]) for s in range(len(origspecies)))
    replaceeqs = {}
    for s in origeqsdict:
        e = origeqsdict[s]
        for i in range(len(removed)): e = e.subs(removed[i][0], removed[i][1])
        e.simplify()
        replaceeqs[s] = e
    neweqsdict = dict((newspecies[s], neweqs[s]) for s in range(len(newspecies)))

    if debug:
        for s in newspecies: print("{}: {}".format(s, replaceeqs[s] - neweqsdict[s]).simplify())

    return sum((replaceeqs[s] - neweqsdict[s]).simplify() != 0 for s in newspecies)