from sage.all import *
from sage.geometry.cone import is_Cone

def is_positive_on(L,K):
    r"""
    Determine whether or not ``L`` is positive on ``K``.

    We say that ``L`` is positive on a closed convex cone ``K`` if
    `L\left\lparen x \right\rparen` belongs to ``K`` for all `x` in
    ``K``. This property need only be checked for generators of ``K``.

    To reliably check whether or not ``L`` is positive, its base ring
    must be either exact (for example, the rationals) or ``SR``. An
    exact ring is more reliable, but in some cases a matrix whose
    entries contain symbolic constants like ``e`` and ``pi`` will work.

    INPUT:

    - ``L`` -- A matrix over either an exact ring or ``SR``.

    - ``K`` -- A polyhedral closed convex cone.

    OUTPUT:

    If the base ring of ``L`` is exact, then ``True`` will be returned if
    and only if ``L`` is positive on ``K``.

    If the base ring of ``L`` is ``SR``, then the situation is more
    complicated:

    - ``True`` will be returned if it can be proven that ``L``
      is positive on ``K``.
    - ``False`` will be returned if it can be proven that ``L``
      is not positive on ``K``.
    - ``False`` will also be returned if we can't decide; specifically
      if we arrive at a symbolic inequality that cannot be resolved.

    .. SEEALSO::

          :func:`is_cross_positive_on`,
          :func:`is_Z_operator_on`,
          :func:`is_lyapunov_like_on`

    EXAMPLES:

    Nonnegative matrices are positive operators on the nonnegative
    orthant::

        sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
        sage: L = random_matrix(QQ,3).apply_map(abs)
        sage: is_positive_on(L,K)
        True

    Your matrix can be over any exact ring, but you may get unexpected
    answers with weirder rings. For example, any rational matrix is
    positive on the plane, but if your matrix contains polynomial
    variables, the answer will be negative::

        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
        sage: K.is_full_space()
        True
        sage: L = matrix(QQ[x], [[x,0],[0,1]])
        sage: is_positive_on(L,K)
        False

    The previous example is "unexpected" because it depends on how we
    check whether or not ``L`` is positive. For exact base rings, we
    check whether or not ``L*z`` belongs to ``K`` for each ``z in K``.
    If ``K`` is closed, then an equally-valid test would be to check
    whether the inner product of ``L*z`` and ``s`` is nonnegative for
    every ``z`` in ``K`` and ``s`` in ``K.dual()``. In fact, that is
    what we do over inexact rings. In the previous example, that test
    would return an affirmative answer::

        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
        sage: L = matrix(QQ[x], [[x,0],[0,1]])
        sage: all([ (L*z).inner_product(s) for z in K for s in K.dual() ])
        True
        sage: is_positive_on(L.change_ring(SR), K)
        True

    TESTS:

    The identity operator is always positive::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=8)
        sage: L = identity_matrix(K.lattice_dim())
        sage: is_positive_on(L,K)
        True

    The "zero" operator is always positive::

        sage: K = random_cone(max_ambient_dim=8)
        sage: R = K.lattice().vector_space().base_ring()
        sage: L = zero_matrix(R, K.lattice_dim())
        sage: is_positive_on(L,K)
        True

    Everything in ``K.positive_operators_gens()`` should be
    positive on ``K``::

        sage: K = random_cone(max_ambient_dim=5)
        sage: all([ is_positive_on(L,K)                     # long time
        ....:       for L in K.positive_operators_gens() ]) # long time
        True
        sage: all([ is_positive_on(L.change_ring(SR),K)     # long time
        ....:       for L in K.positive_operators_gens() ]) # long time
        True

    Technically we could test this, but for now only closed convex cones
    are supported as our ``K`` argument::

        sage: K = [ vector([1,2,3]), vector([5,-1,7]) ]
        sage: L = identity_matrix(3)
        sage: is_positive_on(L,K)
        Traceback (most recent call last):
        ...
        TypeError: K must be a Cone.

    We can't give reliable answers over inexact rings::

        sage: K = Cone([(1,2,3), (4,5,6)])
        sage: L = identity_matrix(RR,3)
        sage: is_positive_on(L,K)
        Traceback (most recent call last):
        ...
        ValueError: The base ring of L is neither SR nor exact.

    """

    if not is_Cone(K):
        raise TypeError('K must be a Cone.')
    if not L.base_ring().is_exact() and not L.base_ring() is SR:
        raise ValueError('The base ring of L is neither SR nor exact.')

    if L.base_ring().is_exact():
        # This should be way faster than computing the dual and
        # checking a bunch of inequalities, but it doesn't work if
        # ``L*x`` is symbolic. For example, ``e in Cone([(1,)])``
        # is true, but returns ``False``.
        return all([ L*x in K for x in K ])
    else:
        # Fall back to inequality-checking when the entries of ``L``
        # might be symbolic.
        return all([ s*(L*x) >= 0 for x in K for s in K.dual() ])


def is_cross_positive_on(L,K):
    r"""
    Determine whether or not ``L`` is cross-positive on ``K``.

    We say that ``L`` is cross-positive on a closed convex cone``K`` if
    `\left\langle L\left\lparenx\right\rparen,s\right\rangle \ge 0` for
    all pairs `\left\langle x,s \right\rangle` in the complementarity
    set of ``K``. This property need only be checked for generators of
    ``K`` and its dual.

    To reliably check whether or not ``L`` is cross-positive, its base
    ring must be either exact (for example, the rationals) or ``SR``. An
    exact ring is more reliable, but in some cases a matrix whose
    entries contain symbolic constants like ``e`` and ``pi`` will work.

    INPUT:

    - ``L`` -- A matrix over either an exact ring or ``SR``.

    - ``K`` -- A polyhedral closed convex cone.

    OUTPUT:

    If the base ring of ``L`` is exact, then ``True`` will be returned if
    and only if ``L`` is cross-positive on ``K``.

    If the base ring of ``L`` is ``SR``, then the situation is more
    complicated:

    - ``True`` will be returned if it can be proven that ``L``
      is cross-positive on ``K``.
    - ``False`` will be returned if it can be proven that ``L``
      is not cross-positive on ``K``.
    - ``False`` will also be returned if we can't decide; specifically
      if we arrive at a symbolic inequality that cannot be resolved.

    .. SEEALSO::

          :func:`is_positive_on`,
          :func:`is_Z_operator_on`,
          :func:`is_lyapunov_like_on`

    EXAMPLES:

    The identity operator is always cross-positive::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=8)
        sage: L = identity_matrix(K.lattice_dim())
        sage: is_cross_positive_on(L,K)
        True

    The "zero" operator is always cross-positive::

        sage: K = random_cone(max_ambient_dim=8)
        sage: R = K.lattice().vector_space().base_ring()
        sage: L = zero_matrix(R, K.lattice_dim())
        sage: is_cross_positive_on(L,K)
        True

    TESTS:

    Everything in ``K.cross_positive_operators_gens()`` should be
    cross-positive on ``K``::

        sage: K = random_cone(max_ambient_dim=5)
        sage: all([ is_cross_positive_on(L,K)                     # long time
        ....:       for L in K.cross_positive_operators_gens() ]) # long time
        True
        sage: all([ is_cross_positive_on(L.change_ring(SR),K)     # long time
        ....:       for L in K.cross_positive_operators_gens() ]) # long time
        True

    Technically we could test this, but for now only closed convex cones
    are supported as our ``K`` argument::

        sage: L = identity_matrix(3)
        sage: K = [ vector([8,2,-8]), vector([5,-5,7]) ]
        sage: is_cross_positive_on(L,K)
        Traceback (most recent call last):
        ...
        TypeError: K must be a Cone.

    We can't give reliable answers over inexact rings::

        sage: K = Cone([(1,2,3), (4,5,6)])
        sage: L = identity_matrix(RR,3)
        sage: is_cross_positive_on(L,K)
        Traceback (most recent call last):
        ...
        ValueError: The base ring of L is neither SR nor exact.

    """
    if not is_Cone(K):
        raise TypeError('K must be a Cone.')
    if not L.base_ring().is_exact() and not L.base_ring() is SR:
        raise ValueError('The base ring of L is neither SR nor exact.')

    return all([ s*(L*x) >= 0
                 for (x,s) in K.discrete_complementarity_set() ])

def is_Z_operator_on(L,K):
    r"""
    Determine whether or not ``L`` is a Z-operator on ``K``.

    We say that ``L`` is a Z-operator on a closed convex cone``K`` if
    `\left\langle L\left\lparenx\right\rparen,s\right\rangle \le 0` for
    all pairs `\left\langle x,s \right\rangle` in the complementarity
    set of ``K``. It is known that this property need only be checked
    for generators of ``K`` and its dual.

    A matrix is a Z-operator on ``K`` if and only if its negation is a
    cross-positive operator on ``K``.

    To reliably check whether or not ``L`` is a Z operator, its base
    ring must be either exact (for example, the rationals) or ``SR``. An
    exact ring is more reliable, but in some cases a matrix whose
    entries contain symbolic constants like ``e`` and ``pi`` will work.

    INPUT:

    - ``L`` -- A matrix over either an exact ring or ``SR``.

    - ``K`` -- A polyhedral closed convex cone.

    OUTPUT:

    If the base ring of ``L`` is exact, then ``True`` will be returned if
    and only if ``L`` is a Z-operator on ``K``.

    If the base ring of ``L`` is ``SR``, then the situation is more
    complicated:

    - ``True`` will be returned if it can be proven that ``L``
      is a Z-operator on ``K``.
    - ``False`` will be returned if it can be proven that ``L``
      is not a Z-operator on ``K``.
    - ``False`` will also be returned if we can't decide; specifically
      if we arrive at a symbolic inequality that cannot be resolved.

    .. SEEALSO::

          :func:`is_positive_on`,
          :func:`is_cross_positive_on`,
          :func:`is_lyapunov_like_on`

    EXAMPLES:

    The identity operator is always a Z-operator::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=8)
        sage: L = identity_matrix(K.lattice_dim())
        sage: is_Z_operator_on(L,K)
        True

    The "zero" operator is always a Z-operator::

        sage: K = random_cone(max_ambient_dim=8)
        sage: R = K.lattice().vector_space().base_ring()
        sage: L = zero_matrix(R, K.lattice_dim())
        sage: is_Z_operator_on(L,K)
        True

    TESTS:

    Everything in ``K.Z_operators_gens()`` should be a Z-operator
    on ``K``::

        sage: K = random_cone(max_ambient_dim=5)
        sage: all([ is_Z_operator_on(L,K)            # long time
        ....:       for L in K.Z_operators_gens() ]) # long time
        True
        sage: all([ is_Z_operator_on(L.change_ring(SR),K) # long time
        ....:       for L in K.Z_operators_gens() ])      # long time
        True

    Technically we could test this, but for now only closed convex cones
    are supported as our ``K`` argument::

        sage: L = identity_matrix(3)
        sage: K = [ vector([-4,20,3]), vector([1,-5,2]) ]
        sage: is_Z_operator_on(L,K)
        Traceback (most recent call last):
        ...
        TypeError: K must be a Cone.


    We can't give reliable answers over inexact rings::

        sage: K = Cone([(1,2,3), (4,5,6)])
        sage: L = identity_matrix(RR,3)
        sage: is_Z_operator_on(L,K)
        Traceback (most recent call last):
        ...
        ValueError: The base ring of L is neither SR nor exact.

    """
    return is_cross_positive_on(-L,K)


def is_lyapunov_like_on(L,K):
    r"""
    Determine whether or not ``L`` is Lyapunov-like on ``K``.

    We say that ``L`` is Lyapunov-like on a closed convex cone ``K`` if
    `\left\langle L\left\lparenx\right\rparen,s\right\rangle = 0` for
    all pairs `\left\langle x,s \right\rangle` in the complementarity
    set of ``K``. This property need only be checked for generators of
    ``K`` and its dual.

    An operator is Lyapunov-like on ``K`` if and only if both the
    operator itself and its negation are cross-positive on ``K``.

    To reliably check whether or not ``L`` is Lyapunov-like, its base
    ring must be either exact (for example, the rationals) or ``SR``. An
    exact ring is more reliable, but in some cases a matrix whose
    entries contain symbolic constants like ``e`` and ``pi`` will work.

    INPUT:

    - ``L`` -- A matrix over either an exact ring or ``SR``.

    - ``K`` -- A polyhedral closed convex cone.

    OUTPUT:

    If the base ring of ``L`` is exact, then ``True`` will be returned if
    and only if ``L`` is Lyapunov-like on ``K``.

    If the base ring of ``L`` is ``SR``, then the situation is more
    complicated:

    - ``True`` will be returned if it can be proven that ``L``
      is Lyapunov-like on ``K``.
    - ``False`` will be returned if it can be proven that ``L``
      is not Lyapunov-like on ``K``.
    - ``False`` will also be returned if we can't decide; specifically
      if we arrive at a symbolic inequality that cannot be resolved.

    .. SEEALSO::

          :func:`is_positive_on`,
          :func:`is_cross_positive_on`,
          :func:`is_Z_operator_on`

    EXAMPLES:

    Diagonal matrices are Lyapunov-like operators on the nonnegative
    orthant::

        sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
        sage: L = diagonal_matrix(random_vector(QQ,3))
        sage: is_lyapunov_like_on(L,K)
        True

    TESTS:

    The identity operator is always Lyapunov-like::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=8)
        sage: L = identity_matrix(K.lattice_dim())
        sage: is_lyapunov_like_on(L,K)
        True

    The "zero" operator is always Lyapunov-like::

        sage: K = random_cone(max_ambient_dim=8)
        sage: R = K.lattice().vector_space().base_ring()
        sage: L = zero_matrix(R, K.lattice_dim())
        sage: is_lyapunov_like_on(L,K)
        True

    Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
    on ``K``::

        sage: K = random_cone(max_ambient_dim=5)
        sage: all([ is_lyapunov_like_on(L,K)            # long time
        ....:       for L in K.lyapunov_like_basis() ]) # long time
        True
        sage: all([ is_lyapunov_like_on(L.change_ring(SR),K) # long time
        ....:       for L in K.lyapunov_like_basis() ])      # long time
        True

    Technically we could test this, but for now only closed convex cones
    are supported as our ``K`` argument::

        sage: L = identity_matrix(3)
        sage: K = [ vector([2,2,-1]), vector([5,4,-3]) ]
        sage: is_lyapunov_like_on(L,K)
        Traceback (most recent call last):
        ...
        TypeError: K must be a Cone.

    We can't give reliable answers over inexact rings::

        sage: K = Cone([(1,2,3), (4,5,6)])
        sage: L = identity_matrix(RR,3)
        sage: is_lyapunov_like_on(L,K)
        Traceback (most recent call last):
        ...
        ValueError: The base ring of L is neither SR nor exact.

    An operator is Lyapunov-like on a cone if and only if both the
    operator and its negation are cross-positive on the cone::

        sage: K = random_cone(max_ambient_dim=5)
        sage: R = K.lattice().vector_space().base_ring()
        sage: L = random_matrix(R, K.lattice_dim())
        sage: actual = is_lyapunov_like_on(L,K)          # long time
        sage: expected = (is_cross_positive_on(L,K) and  # long time
        ....:             is_cross_positive_on(-L,K))    # long time
        sage: actual == expected                         # long time
        True

    """
    if not is_Cone(K):
        raise TypeError('K must be a Cone.')
    if not L.base_ring().is_exact() and not L.base_ring() is SR:
        raise ValueError('The base ring of L is neither SR nor exact.')

    # Even though ``discrete_complementarity_set`` is a cached method
    # of cones, this is faster than calling ``is_cross_positive_on``
    # twice: doing so checks twice as many inequalities as the number
    # of equalities that we're about to check.
    return all([ s*(L*x) == 0
                 for (x,s) in K.discrete_complementarity_set() ])


def LL_cone(K):
    gens = K.lyapunov_like_basis()
    L = ToricLattice(K.lattice_dim()**2)
    return Cone([ g.list() for g in gens ], lattice=L, check=False)

def Sigma_cone(K):
    gens = K.cross_positive_operators_gens()
    L = ToricLattice(K.lattice_dim()**2)
    return Cone([ g.list() for g in gens ], lattice=L, check=False)

def Z_cone(K):
    gens = K.Z_operators_gens()
    L = ToricLattice(K.lattice_dim()**2)
    return Cone([ g.list() for g in gens ], lattice=L, check=False)

def pi_cone(K1, K2=None):
    if K2 is None:
        K2 = K1
    gens = K1.positive_operators_gens(K2)
    L = ToricLattice(K1.lattice_dim()*K2.lattice_dim())
    return Cone([ g.list() for g in gens ], lattice=L, check=False)