from sage.all import *

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

    We say that ``L`` is positive on ``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``.

    INPUT:

    - ``L`` -- A linear transformation or matrix.

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

    OUTPUT:

    ``True`` if it can be proven that ``L`` is positive on ``K``,
    and ``False`` otherwise.

    .. WARNING::

        If this function returns ``True``, then ``L`` is positive
        on ``K``. However, if ``False`` is returned, that could mean one
        of two things. The first is that ``L`` is definitely not
        positive on ``K``. The second is more of an "I don't know"
        answer, returned (for example) if we cannot prove that an inner
        product is nonnegative.

    EXAMPLES:

    The identity is always positive in a nontrivial space::

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

    As is the "zero" transformation::

        sage: K = random_cone(min_ambient_dim=1, 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

    TESTS:

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

        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
        sage: all([ is_positive_on(L,K)
        ....:       for L in K.positive_operators_gens() ])
        True
        sage: all([ is_positive_on(L.change_ring(SR),K)
        ....:       for L in K.positive_operators_gens() ])
        True

    """
    if L.base_ring().is_exact():
        # This could potentially be extended to other types of ``K``...
        return all([ L*x in K for x in K ])
    elif L.base_ring() is SR:
        # 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 ])
    else:
        # The only inexact ring that we're willing to work with is SR,
        # since it can still be exact when working with symbolic
        # constants like pi and e.
        raise ValueError('base ring of operator L is neither SR nor exact')


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 ``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.

    INPUT:

    - ``L`` -- A linear transformation or matrix.

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

    OUTPUT:

    ``True`` if it can be proven that ``L`` is cross-positive on ``K``,
    and ``False`` otherwise.

    .. WARNING::

        If this function returns ``True``, then ``L`` is cross-positive
        on ``K``. However, if ``False`` is returned, that could mean one
        of two things. The first is that ``L`` is definitely not
        cross-positive on ``K``. The second is more of an "I don't know"
        answer, returned (for example) if we cannot prove that an inner
        product is nonnegative.

    EXAMPLES:

    The identity is always cross-positive in a nontrivial space::

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

    As is the "zero" transformation::

        sage: K = random_cone(min_ambient_dim=1, 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(min_ambient_dim=1, max_ambient_dim=6)
        sage: all([ is_cross_positive_on(L,K)
        ....:       for L in K.cross_positive_operators_gens() ])
        True
        sage: all([ is_cross_positive_on(L.change_ring(SR),K)
        ....:       for L in K.cross_positive_operators_gens() ])
        True

    """
    if L.base_ring().is_exact() or L.base_ring() is SR:
        return all([ s*(L*x) >= 0
                     for (x,s) in K.discrete_complementarity_set() ])
    else:
        # The only inexact ring that we're willing to work with is SR,
        # since it can still be exact when working with symbolic
        # constants like pi and e.
        raise ValueError('base ring of operator L is neither SR nor exact')


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

    We say that ``L`` is a Z-operator on ``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``.

    INPUT:

    - ``L`` -- A linear transformation or matrix.

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

    OUTPUT:

    ``True`` if it can be proven that ``L`` is a Z-operator on ``K``,
    and ``False`` otherwise.

    .. WARNING::

        If this function returns ``True``, then ``L`` is a Z-operator
        on ``K``. However, if ``False`` is returned, that could mean one
        of two things. The first is that ``L`` is definitely not
        a Z-operator on ``K``. The second is more of an "I don't know"
        answer, returned (for example) if we cannot prove that an inner
        product is nonnegative.

    EXAMPLES:

    The identity is always a Z-operator in a nontrivial space::

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

    As is the "zero" transformation::

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

    TESTS:

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

        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
        sage: all([ is_Z_on(L,K)
        ....:       for L in K.Z_operators_gens() ])
        True
        sage: all([ is_Z_on(L.change_ring(SR),K)
        ....:       for L in K.Z_operators_gens() ])
        True

    """
    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 ``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.

    INPUT:

    - ``L`` -- A linear transformation or matrix.

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

    OUTPUT:

    ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
    and ``False`` otherwise.

    .. WARNING::

        If this function returns ``True``, then ``L`` is Lyapunov-like
        on ``K``. However, if ``False`` is returned, that could mean one
        of two things. The first is that ``L`` is definitely not
        Lyapunov-like on ``K``. The second is more of an "I don't know"
        answer, returned (for example) if we cannot prove that an inner
        product is zero.

    EXAMPLES:

    The identity is always Lyapunov-like in a nontrivial space::

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

    As is the "zero" transformation::

        sage: K = random_cone(min_ambient_dim=1, 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

    TESTS:

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

        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
        sage: all([ is_lyapunov_like_on(L,K)
        ....:       for L in K.lyapunov_like_basis() ])
        True
        sage: all([ is_lyapunov_like_on(L.change_ring(SR),K)
        ....:       for L in K.lyapunov_like_basis() ])
        True

    """
    if L.base_ring().is_exact() or L.base_ring() is SR:
        # The "fast method" of creating a vector space based on a
        # ``lyapunov_like_basis`` is actually slower than this.
        return all([ s*(L*x) == 0
                     for (x,s) in K.discrete_complementarity_set() ])
    else:
        # The only inexact ring that we're willing to work with is SR,
        # since it can still be exact when working with symbolic
        # constants like pi and e.
        raise ValueError('base ring of operator L is neither SR nor exact')

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)