r"""
The completely positive cone `$\mathcal{K}$` over `\mathbb{R}^{n}$` is
the set of all matrices `$A$`of the form `$\sum uu^{T}$` for `$u \in
\mathbb{R}^{n}_{+}$`. Equivalently, `$A = XX{T}$` where all entries of
`$X$` are nonnegative.
"""

from sage.all import *
from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd
from mjo.cone.doubly_nonnegative import (is_doubly_nonnegative,
                                         is_extreme_doubly_nonnegative)

def is_completely_positive(A):
    """
    Determine whether or not the matrix ``A`` is completely
    positive. This is a known-hard problem, and we'll just give up and
    throw a ``ValueError`` if we can't come to a decision one way or
    another.

    INPUT:

    - ``A`` - The matrix in question

    OUTPUT:

    Either ``True`` if ``A`` is completely positive, or ``False``
    otherwise.

    SETUP::

        sage: from mjo.cone.completely_positive import is_completely_positive

    EXAMPLES:

    Generate an extreme completely positive matrix and check that we
    identify it correctly::

        sage: v = vector([1,2,3,4])
        sage: A = v.column() * v.row()
        sage: A = A.change_ring(QQ)
        sage: is_completely_positive(A)
        True

    Generate an intractible matrix (which does happen to be completely
    positive) and explode::

        sage: v = vector(map(abs, random_vector(ZZ, 10)))
        sage: A = v.column() * v.row()
        sage: A = A.change_ring(QQ)
        sage: is_completely_positive(A)
        Traceback (most recent call last):
        ...
        ValueError: Unable to determine extremity of ``A``.

    Generate a *non*-extreme completely positive matrix and check we
    we identify it correctly. This doesn't work so well if we generate
    random vectors because our algorithm can give up::

        sage: v1 = vector([1,2,3])
        sage: v2 = vector([4,5,6])
        sage: A = v1.column()*v1.row() + v2.column()*v2.row()
        sage: A = A.change_ring(QQ)
        sage: is_completely_positive(A)
        True

    The following matrix isn't positive semidefinite, so it can't be
    completely positive::

        sage: A = matrix(QQ, [[1, 2], [2, 1]])
        sage: is_completely_positive(A)
        False

    This matrix isn't even symmetric, so it can't be completely
    positive::

        sage: A = matrix(QQ, [[1, 2], [3, 4]])
        sage: is_completely_positive(A)
        False

    This one is symmetric but has full rank::

        sage: A = matrix(QQ, [[1, 0], [0, 4]])
        sage: is_completely_positive(A)
        True

    """

    if A.base_ring() == SR:
        msg = 'The matrix ``A`` cannot be symbolic.'
        raise ValueError.new(msg)

    if not is_symmetric_psd(A):
        return False

    n = A.nrows() # Makes sense since ``A`` is symmetric.

    if n <= 4:
        # The two sets are the same for n <= 4.
        return is_doubly_nonnegative(A)

    # No luck.
    raise ValueError('Unable to determine extremity of ``A``.')



def is_extreme_completely_positive(A):
    """
    Determine whether or not the matrix ``A`` is an extreme vector of
    the completely positive cone. This is in general a known-hard
    problem, so our algorithm will give up and throw a ``ValueError`` if
    a decision cannot be made one way or another.

    INPUT:

    - ``A`` - The matrix whose extremity we want to discover

    OUTPUT:

    Either ``True`` if ``A`` is an extreme vector of the completely
    positive cone, or ``False`` otherwise.

    REFERENCES:

    1. Berman, Abraham and Shaked-Monderer, Naomi. Completely Positive
       Matrices. World Scientific, 2003.

    SETUP::

        sage: from mjo.cone.completely_positive import is_extreme_completely_positive

    EXAMPLES:

    Generate an extreme completely positive matrix and check that we
    identify it correctly::

        sage: v = vector(map(abs, random_vector(ZZ, 10)))
        sage: A = v.column() * v.row()
        sage: A = A.change_ring(QQ)
        sage: is_extreme_completely_positive(A)
        True

    Generate a *non*-extreme completely positive matrix and check we
    we identify it correctly. This doesn't work so well if we generate
    random vectors because our algorithm can give up::

        sage: v1 = vector([1,2,3])
        sage: v2 = vector([4,5,6])
        sage: A = v1.column()*v1.row() + v2.column()*v2.row()
        sage: A = A.change_ring(QQ)
        sage: is_extreme_completely_positive(A)
        False

    The following matrix isn't positive semidefinite, so it can't be
    completely positive::

        sage: A = matrix(QQ, [[1, 2], [2, 1]])
        sage: is_extreme_completely_positive(A)
        False

    This matrix isn't even symmetric, so it can't be completely
    positive::

        sage: A = matrix(QQ, [[1, 2], [3, 4]])
        sage: is_extreme_completely_positive(A)
        False

    This one is symmetric but has full rank::

        sage: A = matrix(QQ, [[1, 0], [0, 4]])
        sage: is_extreme_completely_positive(A)
        False

    """

    if A.base_ring() == SR:
        msg = 'The matrix ``A`` cannot be symbolic.'
        raise ValueError(msg)

    if not is_symmetric_psd(A):
        return False

    n = A.nrows() # Makes sense since ``A`` is symmetric.

    # Easy case, see the reference (Remark 2.4).
    if n <= 4:
        return is_extreme_doubly_nonnegative(A)

    # If n > 5, we don't really know. But we might get lucky? See the
    # reference again for this claim.
    if A.rank() == 1 and is_doubly_nonnegative(A):
        return True

    # We didn't get lucky. We have a characterization of the extreme
    # vectors, but it isn't useful if we start with ``A`` because the
    # factorization into `$XX^{T}$` may not be unique!
    raise ValueError('Unable to determine extremity of ``A``.')