"""
Utility functions for working with CVXOPT matrices (instances of the
``cvxopt.base.matrix`` class).
"""

from math import sqrt
from cvxopt import matrix
from cvxopt.lapack import syev

def append_col(left, right):
    """
    Append the matrix ``right`` to the right side of the matrix ``left``.

    EXAMPLES:

        >>> A = matrix([1,2,3,4], (2,2))
        >>> B = matrix([5,6,7,8,9,10], (2,3))
        >>> print(append_col(A,B))
        [  1   3   5   7   9]
        [  2   4   6   8  10]
        <BLANKLINE>

    """
    return matrix([left.trans(), right.trans()]).trans()

def append_row(top, bottom):
    """
    Append the matrix ``bottom`` to the bottom of the matrix ``top``.

    EXAMPLES:

        >>> A = matrix([1,2,3,4], (2,2))
        >>> B = matrix([5,6,7,8,9,10], (3,2))
        >>> print(append_row(A,B))
        [  1   3]
        [  2   4]
        [  5   8]
        [  6   9]
        [  7  10]
        <BLANKLINE>

    """
    return matrix([top, bottom])


def eigenvalues(real_matrix):
    """
    Return the eigenvalues of the given ``real_matrix``.

    EXAMPLES:

        >>> A = matrix([[2,1],[1,2]], tc='d')
        >>> eigenvalues(A)
        [1.0, 3.0]

    """
    domain_dim = real_matrix.size[0] # Assume ``real_matrix`` is square.
    eigs = matrix(0, (domain_dim, 1), tc='d')
    syev(real_matrix, eigs)
    return list(eigs)


def identity(domain_dim):
    """
    Return a ``domain_dim``-by-``domain_dim`` dense integer identity
    matrix.

    EXAMPLES:

       >>> print(identity(3))
       [ 1  0  0]
       [ 0  1  0]
       [ 0  0  1]
       <BLANKLINE>

    """
    if domain_dim <= 0:
        raise ValueError('domain dimension must be positive')

    entries = [int(i == j)
               for i in range(domain_dim)
               for j in range(domain_dim)]
    return matrix(entries, (domain_dim, domain_dim))


def inner_product(vec1, vec2):
    """
    Compute the (Euclidean) inner product of the two vectors ``vec1``
    and ``vec2``.

    EXAMPLES:

        >>> x = [1,2,3]
        >>> y = [3,4,1]
        >>> inner_product(x,y)
        14

        >>> x = matrix([1,1,1])
        >>> y = matrix([2,3,4], (1,3))
        >>> inner_product(x,y)
        9

        >>> x = [1,2,3]
        >>> y = [1,1]
        >>> inner_product(x,y)
        Traceback (most recent call last):
        ...
        TypeError: the lengths of vec1 and vec2 must match

    """
    if not len(vec1) == len(vec2):
        raise TypeError('the lengths of vec1 and vec2 must match')

    return sum([x*y for (x,y) in zip(vec1,vec2)])


def norm(matrix_or_vector):
    """
    Return the Frobenius norm of ``matrix_or_vector``, which is the same
    thing as its Euclidean norm when it's a vector (when one of its
    dimensions is unity).

    EXAMPLES:

        >>> v = matrix([1,1])
        >>> print('{:.5f}'.format(norm(v)))
        1.41421

        >>> A = matrix([1,1,1,1], (2,2))
        >>> norm(A)
        2.0

    """
    return sqrt(inner_product(matrix_or_vector,matrix_or_vector))


def vec(real_matrix):
    """
    Create a long vector in column-major order from ``real_matrix``.

    EXAMPLES:

        >>> A = matrix([[1,2],[3,4]])
        >>> print(A)
        [ 1  3]
        [ 2  4]
        <BLANKLINE>

        >>> print(vec(A))
        [ 1]
        [ 2]
        [ 3]
        [ 4]
        <BLANKLINE>

    Note that if ``real_matrix`` is a vector, this function is a no-op:

        >>> v = matrix([1,2,3,4], (4,1))
        >>> print(v)
        [ 1]
        [ 2]
        [ 3]
        [ 4]
        <BLANKLINE>
        >>> print(vec(v))
        [ 1]
        [ 2]
        [ 3]
        [ 4]
        <BLANKLINE>

    """
    return matrix(real_matrix, (len(real_matrix), 1))