"""
The doubly-nonnegative cone in `S^{n}` is the set of all such matrices
that both,

  a) are positive semidefinite

  b) have only nonnegative entries

It is represented typically by either `\mathcal{D}^{n}` or
`\mathcal{DNN}`.

"""

from sage.all import *

# Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
# have to explicitly mangle our sitedir here so that "mjo.cone"
# resolves.
from os.path import abspath
from site import addsitedir
addsitedir(abspath('../../'))
from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd



def is_doubly_nonnegative(A):
    """
    Determine whether or not the matrix ``A`` is doubly-nonnegative.

    INPUT:

    - ``A`` - The matrix in question

    OUTPUT:

    Either ``True`` if ``A`` is doubly-nonnegative, or ``False``
    otherwise.

    EXAMPLES:

    Every completely positive matrix is doubly-nonnegative::

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

    The following matrix is nonnegative but non positive semidefinite::

        sage: A = matrix(ZZ, [[1, 2], [2, 1]])
        sage: is_doubly_nonnegative(A)
        False

    """

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

    # Check that all of the entries of ``A`` are nonnegative.
    if not all([ a >= 0 for a in A.list() ]):
        return False

    # It's nonnegative, so all we need to do is check that it's
    # symmetric positive-semidefinite.
    return is_symmetric_psd(A)



def has_admissible_extreme_rank(A):
    """
    The extreme matrices of the doubly-nonnegative cone have some
    restrictions on their ranks. This function checks to see whether or
    not ``A`` could be extreme based on its rank.

    INPUT:

    - ``A`` - The matrix in question

    OUTPUT:

    ``False`` if the rank of ``A`` precludes it from being an extreme
    matrix of the doubly-nonnegative cone, ``True`` otherwise.

    REFERENCE:

    Hamilton-Jester, Crista Lee; Li, Chi-Kwong. Extreme Vectors of
    Doubly Nonnegative Matrices. Rocky Mountain Journal of Mathematics
    26 (1996), no. 4, 1371--1383. doi:10.1216/rmjm/1181071993.
    http://projecteuclid.org/euclid.rmjm/1181071993.

    EXAMPLES:

    The zero matrix has rank zero, which is admissible::

    sage: A = zero_matrix(QQ, 5, 5)
    sage: has_admissible_extreme_rank(A)
    True

    """
    if not A.is_symmetric():
        # This function is more or less internal, so blow up if passed
        # something unexpected.
        raise ValueError('The matrix ``A`` must be symmetric.')

    r = rank(A)
    n = ZZ(A.nrows()) # Columns would work, too, since ``A`` is symmetric.

    if r == 0:
        # Zero is in the doubly-nonnegative cone.
        return True

    # See Theorem 3.1 in the cited reference.
    if r == 2:
        return False

    if n.mod(2) == 0:
        # n is even
        return r <= max(1, n-3)
    else:
        # n is odd
        return r <= max(1, n-2)


def E(matrix_space, i,j):
    """
    Return the ``i``,``j``th element of the standard basis in
    ``matrix_space``.

    INPUT:

    - ``matrix_space`` - The underlying matrix space of whose basis
                         the returned matrix is an element

    - ``i`` - The row index of the single nonzero entry

    - ``j`` - The column index of the single nonzero entry

    OUTPUT:

    A basis element of ``matrix_space``. It has a single \"1\" in the
    ``i``,``j`` row,column and zeros elsewhere.

    EXAMPLES::

        sage: M = MatrixSpace(ZZ, 2, 2)
        sage: E(M,0,0)
        [1 0]
        [0 0]
        sage: E(M,0,1)
        [0 1]
        [0 0]
        sage: E(M,1,0)
        [0 0]
        [1 0]
        sage: E(M,1,1)
        [0 0]
        [0 1]
        sage: E(M,2,1)
        Traceback (most recent call last):
        ...
        IndexError: Index `i` is out of bounds.
        sage: E(M,1,2)
        Traceback (most recent call last):
        ...
        IndexError: Index `j` is out of bounds.

    """
    # We need to check these ourselves, see below.
    if i >= matrix_space.nrows():
        raise IndexError('Index `i` is out of bounds.')
    if j >= matrix_space.ncols():
        raise IndexError('Index `j` is out of bounds.')

    # The basis here is returned as a one-dimensional list, so we need
    # to compute the offset into it based on ``i`` and ``j``. Since we
    # compute the index ourselves, we need to do bounds-checking
    # manually. Otherwise for e.g. a 2x2 matrix space, the index (0,2)
    # would be computed as offset 3 into a four-element list and we
    # would succeed incorrectly.
    idx = matrix_space.ncols()*i + j
    return matrix_space.basis()[idx]



def is_extreme_doubly_nonnegative(A):
    """
    Returns ``True`` if the given matrix is an extreme matrix of the
    doubly-nonnegative cone, and ``False`` otherwise.

    EXAMPLES:

    The zero matrix is an extreme matrix::

        sage: A = zero_matrix(QQ, 5, 5)
        sage: is_extreme_doubly_nonnegative(A)
        True

    """

    r = A.rank()

    if r == 0:
        # Short circuit, we know the zero matrix is extreme.
        return True

    if not is_admissible_extreme_rank(r):
        return False

    raise NotImplementedError()