r"""
The symmetric positive definite cone `$S^{n}_{++}$` is the cone
consisting of all symmetric positive-definite matrices (as a subset of
$\mathbb{R}^{n \times n}$`. It is the interior of the symmetric positive
SEMI-definite cone.
"""

from sage.all import *
from mjo.cone.symmetric_psd import random_symmetric_psd

def is_symmetric_pd(A):
    """
    Determine whether or not the matrix ``A`` is symmetric
    positive-definite.

    INPUT:

    - ``A`` - The matrix in question.

    OUTPUT:

    Either ``True`` if ``A`` is symmetric positive-definite, or
    ``False`` otherwise.

    SETUP::

        sage: from mjo.cone.symmetric_pd import (is_symmetric_pd,
        ....:                                    random_symmetric_pd)

    EXAMPLES:

    The identity matrix is obviously symmetric and positive-definite::

        sage: A = identity_matrix(ZZ, ZZ.random_element(10))
        sage: is_symmetric_pd(A)
        True

    The following matrix is symmetric but not positive definite::

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

    This matrix isn't even symmetric::

        sage: A = matrix(ZZ, [[1, 2], [3, 4]])
        sage: is_symmetric_pd(A)
        False

    The trivial matrix in a trivial space is trivially symmetric and
    positive-definite::

        sage: A = matrix(QQ, 0,0)
        sage: is_symmetric_pd(A)
        True

    The implementation of "is_positive_definite" in Sage leaves a bit to
    be desired. This matrix is technically positive definite over the
    real numbers, but isn't symmetric::

        sage: A = matrix(RR,[[1,1],[-1,1]])
        sage: A.is_positive_definite()
        Traceback (most recent call last):
        ...
        ValueError: Could not see Real Field with 53 bits of precision as a
        subring of the real or complex numbers
        sage: is_symmetric_pd(A)
        False

    TESTS:

    Every gram matrix is positive-definite, and we can sum two
    positive-definite matrices (``A`` and its transpose) to get a new
    positive-definite matrix that happens to be symmetric::

        sage: set_random_seed()
        sage: V = VectorSpace(QQ, ZZ.random_element(5))
        sage: A = random_symmetric_pd(V)
        sage: is_symmetric_pd(A)
        True

    """

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

    # First make sure that ``A`` is symmetric.
    if not A.is_symmetric():
        return False

    # If ``A`` is symmetric, we only need to check that it is positive
    # definite. For that we can consult its minimum eigenvalue, which
    # should be greater than zero. Since ``A`` is symmetric, its
    # eigenvalues are guaranteed to be real.
    if A.is_zero():
        # A is trivial... so trivially positive-definite.
        return True
    else:
        return min(A.eigenvalues()) > 0


def random_symmetric_pd(V):
    r"""
    Generate a random symmetric positive-definite matrix over the
    vector space ``V``. That is, the returned matrix will be a linear
    transformation on ``V``, with the same base ring as ``V``.

    (We take a very loose interpretation of "random," here.)

    INPUT:

    - ``V`` -- The vector space on which the returned matrix will act.

    OUTPUT:

    A random symmetric positive-definite matrix, i.e. a linear
    transformation from ``V`` to itself.

    ALGORITHM:

    We request a full-rank positive-semidefinite matrix from the
    :func:`mjo.cone.symmetric_psd.random_symmetric_psd` function.

    SETUP::

        sage: from mjo.cone.symmetric_pd import random_symmetric_pd

    EXAMPLES:

    Well, it doesn't crash at least::

        sage: set_random_seed()
        sage: V = VectorSpace(QQ, 2)
        sage: A = random_symmetric_pd(V)
        sage: A.matrix_space()
        Full MatrixSpace of 2 by 2 dense matrices over Rational Field

    """
    # We accept_zero because the trivial matrix is (trivially)
    # symmetric and positive-definite, but it's also "zero." The rank
    # condition ensures that we don't get the zero matrix in a
    # nontrivial space.
    return random_symmetric_psd(V, rank=V.dimension())


def inverse_symmetric_pd(M):
    r"""
    Compute the inverse of a symmetric positive-definite matrix.

    The symmetric positive-definite matrix ``M`` has a Cholesky
    factorization, which can be used to invert ``M`` in a fast,
    numerically-stable way.

    INPUT:

    - ``M`` -- The matrix to invert.

    OUTPUT:

    The inverse of ``M``. The Cholesky factorization uses roots, so the
    returned matrix will have as its base ring the algebraic closure of
    the base ring of ``M``.

    SETUP::

        sage: from mjo.cone.symmetric_pd import inverse_symmetric_pd

    EXAMPLES:

    If we start with a rational matrix, its inverse will be over a field
    extension of that rationals that is necessarily contained in
    ``QQbar`` (which contails ALL roots)::

        sage: M = matrix(QQ, [[3,1],[1,5]])
        sage: M_inv = inverse_symmetric_pd(M)
        sage: M_inv.base_ring().is_subring(QQbar)
        True

    TESTS:

    This "inverse" should actually act like an inverse::

        sage: set_random_seed()
        sage: n = ZZ.random_element(5)
        sage: V = VectorSpace(QQ,n)
        sage: from mjo.cone.symmetric_pd import random_symmetric_pd
        sage: M = random_symmetric_pd(V)
        sage: M_inv = inverse_symmetric_pd(M)
        sage: I = identity_matrix(QQ,n)
        sage: M_inv * M == I
        True
        sage: M * M_inv == I
        True

    """
    # M = L * L.transpose()
    # The default "echelonize" inverse() method works just fine for
    # triangular matrices.
    L_inverse = M.cholesky().inverse()
    return L_inverse.transpose()*L_inverse