def test_assumptions():
    n = Symbol('n')
    A = MatrixSymbol('A', 1, n)
    P = PermutationMatrix(A)
    assert ask(Q.integer_elements(P))
    assert ask(Q.real_elements(P))
    assert ask(Q.complex_elements(P))

def test_POSV():
    X = MatrixSymbol('X', n, n)
    Y = MatrixSymbol('Y', n, m)
    posv = POSV(X, Y)
    assert posv.outputs[0] == X.I*Y
    assert not POSV.valid(posv.inputs, True)
    assert POSV.valid(posv.inputs, Q.symmetric(X) & Q.positive_definite(X))

def test_pow2():
    # powers of (-1)
    assert refine((-1)**(x + y), Q.even(x)) == (-1)**y
    assert refine((-1)**(x + y + z), Q.odd(x) & Q.odd(z)) == (-1)**y
    assert refine((-1)**(x + y + 1), Q.odd(x)) == (-1)**y
    assert refine((-1)**(x + y + 2), Q.odd(x)) == (-1)**(y + 1)
    assert refine((-1)**(x + 3)) == (-1)**(x + 1)

def test_reduce_poly_inequalities_real_relational():
    global_assumptions.add(Q.real(x))
    global_assumptions.add(Q.real(y))

    assert reduce_poly_inequalities([[Eq(x ** 2, 0)]], x, relational=True) == Eq(x, 0)
    assert reduce_poly_inequalities([[Le(x ** 2, 0)]], x, relational=True) == Eq(x, 0)
    assert reduce_poly_inequalities([[Lt(x ** 2, 0)]], x, relational=True) == False
    assert reduce_poly_inequalities([[Ge(x ** 2, 0)]], x, relational=True) == True
    assert reduce_poly_inequalities([[Gt(x ** 2, 0)]], x, relational=True) == Or(Lt(x, 0), Lt(0, x))
    assert reduce_poly_inequalities([[Ne(x ** 2, 0)]], x, relational=True) == Or(Lt(x, 0), Lt(0, x))

    assert reduce_poly_inequalities([[Eq(x ** 2, 1)]], x, relational=True) == Or(Eq(x, -1), Eq(x, 1))
    assert reduce_poly_inequalities([[Le(x ** 2, 1)]], x, relational=True) == And(Le(-1, x), Le(x, 1))
    assert reduce_poly_inequalities([[Lt(x ** 2, 1)]], x, relational=True) == And(Lt(-1, x), Lt(x, 1))
    assert reduce_poly_inequalities([[Ge(x ** 2, 1)]], x, relational=True) == Or(Le(x, -1), Le(1, x))
    assert reduce_poly_inequalities([[Gt(x ** 2, 1)]], x, relational=True) == Or(Lt(x, -1), Lt(1, x))
    assert reduce_poly_inequalities([[Ne(x ** 2, 1)]], x, relational=True) == Or(
        Lt(x, -1), And(Lt(-1, x), Lt(x, 1)), Lt(1, x)
    )

    assert (
        reduce_poly_inequalities([[Eq(x ** 2, 1.0)]], x, relational=True).evalf() == Or(Eq(x, -1.0), Eq(x, 1.0)).evalf()
    )
    assert reduce_poly_inequalities([[Le(x ** 2, 1.0)]], x, relational=True) == And(Le(-1.0, x), Le(x, 1.0))
    assert reduce_poly_inequalities([[Lt(x ** 2, 1.0)]], x, relational=True) == And(Lt(-1.0, x), Lt(x, 1.0))
    assert reduce_poly_inequalities([[Ge(x ** 2, 1.0)]], x, relational=True) == Or(Le(x, -1.0), Le(1.0, x))
    assert reduce_poly_inequalities([[Gt(x ** 2, 1.0)]], x, relational=True) == Or(Lt(x, -1.0), Lt(1.0, x))
    assert reduce_poly_inequalities([[Ne(x ** 2, 1.0)]], x, relational=True) == Or(
        Lt(x, -1.0), And(Lt(-1.0, x), Lt(x, 1.0)), Lt(1.0, x)
    )

    global_assumptions.remove(Q.real(x))
    global_assumptions.remove(Q.real(y))

def test_reduce_rational_inequalities_real_relational():
    def OpenInterval(a, b):
        return Interval(a, b, True, True)
    def LeftOpenInterval(a, b):
        return Interval(a, b, True, False)
    def RightOpenInterval(a, b):
        return Interval(a, b, False, True)

    with assuming(Q.real(x)):
        assert reduce_rational_inequalities([[(x**2 + 3*x + 2)/(x**2 - 16) >= 0]], x, relational=False) == \
            Union(OpenInterval(-oo, -4), Interval(-2, -1), OpenInterval(4, oo))

        assert reduce_rational_inequalities([[((-2*x - 10)*(3 - x))/((x**2 + 5)*(x - 2)**2) < 0]], x, relational=False) == \
            Union(OpenInterval(-5, 2), OpenInterval(2, 3))

        assert reduce_rational_inequalities([[(x + 1)/(x - 5) <= 0]], x, assume=Q.real(x), relational=False) == \
            RightOpenInterval(-1, 5)

        assert reduce_rational_inequalities([[(x**2 + 4*x + 3)/(x - 1) > 0]], x, assume=Q.real(x), relational=False) == \
            Union(OpenInterval(-3, -1), OpenInterval(1, oo))

        assert reduce_rational_inequalities([[(x**2 - 16)/(x - 1)**2 < 0]], x, assume=Q.real(x), relational=False) == \
            Union(OpenInterval(-4, 1), OpenInterval(1, 4))

        assert reduce_rational_inequalities([[(3*x + 1)/(x + 4) >= 1]], x, assume=Q.real(x), relational=False) == \
            Union(OpenInterval(-oo, -4), RightOpenInterval(S(3)/2, oo))

        assert reduce_rational_inequalities([[(x - 8)/x <= 3 - x]], x, assume=Q.real(x), relational=False) == \
            Union(LeftOpenInterval(-oo, -2), LeftOpenInterval(0, 4))

def test_zero_pow():
    assert satask(Q.zero(x**y), Q.zero(x) & Q.positive(y)) is True
    assert satask(Q.zero(x**y), Q.nonzero(x) & Q.zero(y)) is False

    assert satask(Q.zero(x), Q.zero(x**y)) is True

    assert satask(Q.zero(x**y), Q.zero(x)) is None

def test_refine_issue_12724():
    expr1 = refine(Abs(x * y), Q.positive(x))
    expr2 = refine(Abs(x * y * z), Q.positive(x))
    assert expr1 == x * Abs(y)
    assert expr2 == x * Abs(y * z)
    y1 = Symbol('y1', real = True)
    expr3 = refine(Abs(x * y1**2 * z), Q.positive(x))
    assert expr3 == x * y1**2 * Abs(z)

def test_dtype_of():
    X = MatrixSymbol('X', n, n)
    assert 'integer' in dtype_of(X, Q.integer_elements(X))
    assert 'real' in dtype_of(X, Q.real_elements(X))
    assert 'complex' in dtype_of(X, Q.complex_elements(X))

    alpha = Symbol('alpha', integer=True)
    assert 'integer' in dtype_of(alpha)

def test_Abs():
    assert refine(Abs(x), Q.positive(x)) == x
    assert refine(1 + Abs(x), Q.positive(x)) == 1 + x
    assert refine(Abs(x), Q.negative(x)) == -x
    assert refine(1 + Abs(x), Q.negative(x)) == 1 - x

    assert refine(Abs(x ** 2)) != x ** 2
    assert refine(Abs(x ** 2), Q.real(x)) == x ** 2

def test_pow4():
    assert refine((-1)**((-1)**x/2 - 7*S.Half), Q.integer(x)) == (-1)**(x + 1)
    assert refine((-1)**((-1)**x/2 - 9*S.Half), Q.integer(x)) == (-1)**x

    # powers of Abs
    assert refine(Abs(x)**2, Q.real(x)) == x**2
    assert refine(Abs(x)**3, Q.real(x)) == Abs(x)**3
    assert refine(Abs(x)**2) == Abs(x)**2

def test_numerics():
    import numpy as np
    with assuming(Q.real_elements(X), Q.real_elements(y)):
        f = build(ic, inputs, outputs, filename='numerics.f90',
                modname='numerics')
    nX, ny = np.ones((5, 5)), np.ones(5)
    result = np.matrix(nX) * np.matrix(ny).T
    assert np.allclose(f(nX, ny), result)

def test_valid():
    A = MatrixSymbol('A', n, n)
    B = MatrixSymbol('B', n, n)
    C = MatrixSymbol('C', n, n)
    assert GEMM.valid((1, A, B, 2, C), True)
    assert not SYMM.valid((1, A, B, 2, C), True)
    assert SYMM.valid((1, A, B, 2, C), Q.symmetric(A))
    assert SYMM.valid((1, A, B, 2, C), Q.symmetric(B))

def test_triangular():
    assert ask(Q.upper_triangular(X + Z.T + Identity(2)), Q.upper_triangular(X) &
            Q.lower_triangular(Z)) is True
    assert ask(Q.upper_triangular(X*Z.T), Q.upper_triangular(X) &
            Q.lower_triangular(Z)) is True
    assert ask(Q.lower_triangular(Identity(3))) is True
    assert ask(Q.lower_triangular(ZeroMatrix(3, 3))) is True
    assert ask(Q.triangular(X), Q.unit_triangular(X))

def test_POSV_codemap():
    A = MatrixSymbol('A', n, n)
    B = MatrixSymbol('B', n, m)
    assumptions = Q.positive_definite(A) & Q.symmetric(A)
    codemap = POSV(A, B).codemap('A B INFO'.split(), assumptions)
    call = POSV.fortran_template % codemap
    assert "('U', n, m, A, n, B, n, INFO)" in call
    assert 'dposv' in call.lower()

def test_ElemProd_code():
    from computations.matrices.fortran.core import generate
    ic = inplace_compile(c)
    with assuming(Q.real_elements(x), Q.real_elements(y)):
        s = generate(ic, [x,y], [HadamardProduct(x,y)])
    with open('tmp/elem_test.f90','w') as f:
        f.write(s)
    assert "= X * Y" in s

 def _eval_determinant(self):
     if self.blockshape == (2, 2):
         [[A, B],
          [C, D]] = self.blocks.tolist()
         if ask(Q.invertible(A)):
             return det(A)*det(D - C*A.I*B)
         elif ask(Q.invertible(D)):
             return det(D)*det(A - B*D.I*C)
     return Determinant(self)

def test_execution():
    with assuming(Q.real_elements(X), Q.real_elements(Y)):
        f = build(pgemm, [X, Y], [pgemm.duration], filename='tmp/profile.f90',
                modname='profile')
    assert callable(f)
    nX, nY = np.random.rand(500, 500), np.random.rand(500, 500)
    result = f(nX, nY)
    assert isinstance(result, float)
    assert result > 0

def test_is_literal():
    assert is_literal(True) is True
    assert is_literal(False) is True
    assert is_literal(A) is True
    assert is_literal(~A) is True
    assert is_literal(Or(A, B)) is False
    assert is_literal(Q.zero(A)) is True
    assert is_literal(Not(Q.zero(A))) is True
    assert is_literal(Or(A, B)) is False
    assert is_literal(And(Q.zero(A), Q.zero(B))) is False

def test_GEMM():
    with assuming(Q.real_elements(A), Q.real_elements(X)):
        f = build(GEMM(1.0, A, X, 0.0, ZeroMatrix(A.rows, X.cols)),
                [A, X], [A*X], modname='gemmtest', filename='tmp/gemmtest.f90')

    nA = np.asarray([[1, 2], [3, 4]], dtype=np.float64, order='F')
    nX = np.asarray([[1, 1], [1, 1]], dtype=np.float64, order='F')
    expected = np.asarray([[3., 3.], [7., 7.]])
    result = f(nA, nX)
    assert np.allclose(expected, result)

def test_Join_build():
    c = Join(B)
    with assuming(Q.real_elements(X), Q.real_elements(Y)):
        f = build(c, [X, Y], [B], filename='tmp/join.f90', modname='join')
    nX = np.ones((2, 3), dtype=np.float64)
    nY = np.ones((2, 3), dtype=np.float64)
    result   = f(nX, nY)
    expected = np.ones((2, 6), dtype=np.float64)

    assert np.allclose(result, expected)