"""
There is an explicit isomorphism between all finite-dimensional vector
spaces. In particular, there is an isomorphism between the m-by-n
matrices and `$R^(m \times n)$`. Since most vector operations are not
available on Sage matrices, we have to go back and forth between these
two vector spaces often.
"""

from sage.all import *

def isomorphism(matrix_space):
    """
    Create isomorphism (i.e. the function) that converts elements
    of a matrix space into those of the corresponding finite-dimensional
    vector space.

    INPUT:

    - matrix_space: A finite-dimensional ``MatrixSpace`` object.

    OUTPUT:

    - (phi, phi_inverse): If ``matrix_space`` has dimension m*n, then
                          ``phi`` will map m-by-n matrices to R^(m*n).
                          The inverse mapping ``phi_inverse`` will go
                          the other way.

    SETUP::

        sage: from mjo.matrix_vector import isomorphism

    EXAMPLES::

        sage: M = MatrixSpace(QQ,4,4)
        sage: (p, p_inv) = isomorphism(M)
        sage: m = M(xrange(16))
        sage: p_inv(p(m)) == m
        True

    """
    from sage.matrix.matrix_space import is_MatrixSpace
    if not is_MatrixSpace(matrix_space):
        raise TypeError('argument must be a matrix space')

    base_ring = matrix_space.base_ring()
    vector_space = VectorSpace(base_ring, matrix_space.dimension())

    def phi(m):
        return vector_space(m.list())

    def phi_inverse(v):
        return matrix_space(v.list())

    return (phi, phi_inverse)



def matrix_of_transformation(T, V):
    """
    Compute the matrix of a linear transformation ``T``, `$T : V
    \rightarrow V$` with domain/range ``V``. This essentially uses the
    Riesz representation theorem to represent the entries of the matrix
    of ``T`` in terms of inner products.

    INPUT:

    - ``T`` -- The linear transformation whose matrix we should
               compute. This should be a callable function that
               takes as its single argument an element of ``V``.

    - ``V`` -- The vector or matrix space on which ``T`` is defined.

    OUTPUT:

    If the dimension of ``V`` is `$n$`, we return an `$n \times n$`
    matrix that represents ``T`` with respect to the standard basis of
    ``V``.

    SETUP::

        sage: from mjo.matrix_vector import isomorphism, matrix_of_transformation

    EXAMPLES:

    The matrix of a transformation on a simple vector space should be
    the expected matrix::

        sage: V = VectorSpace(QQ, 3)
        sage: def f(x):
        ....:     return 3*x
        ....:
        sage: matrix_of_transformation(f, V)
        [3 0 0]
        [0 3 0]
        [0 0 3]

    A more complicated example confirms that we get a matrix consistent
    with our ``matrix_to_vector`` function::

        sage: M = MatrixSpace(QQ,3,3)
        sage: Q = M([[0,1,0],[1,0,0],[0,0,1]])
        sage: def f(x):
        ....:     return Q*x*Q.inverse()
        ....:
        sage: F = matrix_of_transformation(f, M)
        sage: F
        [0 0 0 0 1 0 0 0 0]
        [0 0 0 1 0 0 0 0 0]
        [0 0 0 0 0 1 0 0 0]
        [0 1 0 0 0 0 0 0 0]
        [1 0 0 0 0 0 0 0 0]
        [0 0 1 0 0 0 0 0 0]
        [0 0 0 0 0 0 0 1 0]
        [0 0 0 0 0 0 1 0 0]
        [0 0 0 0 0 0 0 0 1]
        sage: phi, phi_inv = isomorphism(M)
        sage: X = M([[1,2,3],[4,5,6],[7,8,9]])
        sage: F*phi(X)
        (5, 4, 6, 2, 1, 3, 8, 7, 9)
        sage: phi(f(X))
        (5, 4, 6, 2, 1, 3, 8, 7, 9)
        sage: F*phi(X) == phi(f(X))
        True

    """
    n = V.dimension()
    B = list(V.basis())

    def inner_product(v, w):
        # An inner product function that works for both matrices and
        # vectors.
        if callable(getattr(v, 'inner_product', None)):
            return v.inner_product(w)
        elif callable(getattr(v, 'matrix_space', None)):
            # V must be a matrix space?
            return (v*w.transpose()).trace()
        else:
            raise ValueError('inner_product only works on vectors and matrices')

    def apply(L, x):
        # A "call" that works for both matrices and functions.
        if callable(getattr(L, 'matrix_space', None)):
            # L is a matrix, and we need to use "multiply" to call it.
            return L*x
        else:
            # If L isn't a matrix, try this. It works for python
            # functions at least.
            return L(x)

    entries = []
    for j in xrange(n):
        for i in xrange(n):
            entry = inner_product(apply(T,B[i]), B[j])
            entries.append(entry)

    # Construct the matrix space in which our return value will lie.
    W = MatrixSpace(V.base_ring(), n, n)

    # And make a matrix out of our list of entries.
    A = W(entries)

    return A