#******************************************************************************
#   Copyright (C) 2013 Kenneth L. Ho
#
#   Redistribution and use in source and binary forms, with or without
#   modification, are permitted provided that the following conditions are met:
#
#   Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer. Redistributions in binary
#   form must reproduce the above copyright notice, this list of conditions and
#   the following disclaimer in the documentation and/or other materials
#   provided with the distribution.
#
#   None of the names of the copyright holders may be used to endorse or
#   promote products derived from this software without specific prior written
#   permission.
#
#   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
#   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
#   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
#   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
#   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
#   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
#   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
#   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
#   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
#   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
#   POSSIBILITY OF SUCH DAMAGE.
#******************************************************************************

"""
Direct wrappers for Fortran `id_dist` backend.
"""

import scipy.linalg._interpolative as _id
import numpy as np

_RETCODE_ERROR = RuntimeError("nonzero return code")


#------------------------------------------------------------------------------
# id_rand.f
#------------------------------------------------------------------------------

def id_srand(n):
    """
    Generate standard uniform pseudorandom numbers via a very efficient lagged
    Fibonacci method.

    :param n:
        Number of pseudorandom numbers to generate.
    :type n: int

    :return:
        Pseudorandom numbers.
    :rtype: :class:`numpy.ndarray`
    """
    return _id.id_srand(n)


def id_srandi(t):
    """
    Initialize seed values for :func:`id_srand` (any appropriately random
    numbers will do).

    :param t:
        Array of 55 seed values.
    :type t: :class:`numpy.ndarray`
    """
    t = np.asfortranarray(t)
    _id.id_srandi(t)


def id_srando():
    """
    Reset seed values to their original values.
    """
    _id.id_srando()


#------------------------------------------------------------------------------
# idd_frm.f
#------------------------------------------------------------------------------

def idd_frm(n, w, x):
    """
    Transform real vector via a composition of Rokhlin's random transform,
    random subselection, and an FFT.

    In contrast to :func:`idd_sfrm`, this routine works best when the length of
    the transformed vector is the power-of-two integer output by
    :func:`idd_frmi`, or when the length is not specified but instead
    determined a posteriori from the output. The returned transformed vector is
    randomly permuted.

    :param n:
        Greatest power-of-two integer satisfying `n <= x.size` as obtained from
        :func:`idd_frmi`; `n` is also the length of the output vector.
    :type n: int
    :param w:
        Initialization array constructed by :func:`idd_frmi`.
    :type w: :class:`numpy.ndarray`
    :param x:
        Vector to be transformed.
    :type x: :class:`numpy.ndarray`

    :return:
        Transformed vector.
    :rtype: :class:`numpy.ndarray`
    """
    return _id.idd_frm(n, w, x)


def idd_sfrm(l, n, w, x):
    """
    Transform real vector via a composition of Rokhlin's random transform,
    random subselection, and an FFT.

    In contrast to :func:`idd_frm`, this routine works best when the length of
    the transformed vector is known a priori.

    :param l:
        Length of transformed vector, satisfying `l <= n`.
    :type l: int
    :param n:
        Greatest power-of-two integer satisfying `n <= x.size` as obtained from
        :func:`idd_sfrmi`.
    :type n: int
    :param w:
        Initialization array constructed by :func:`idd_sfrmi`.
    :type w: :class:`numpy.ndarray`
    :param x:
        Vector to be transformed.
    :type x: :class:`numpy.ndarray`

    :return:
        Transformed vector.
    :rtype: :class:`numpy.ndarray`
    """
    return _id.idd_sfrm(l, n, w, x)


def idd_frmi(m):
    """
    Initialize data for :func:`idd_frm`.

    :param m:
        Length of vector to be transformed.
    :type m: int

    :return:
        Greatest power-of-two integer `n` satisfying `n <= m`.
    :rtype: int
    :return:
        Initialization array to be used by :func:`idd_frm`.
    :rtype: :class:`numpy.ndarray`
    """
    return _id.idd_frmi(m)


def idd_sfrmi(l, m):
    """
    Initialize data for :func:`idd_sfrm`.

    :param l:
        Length of output transformed vector.
    :type l: int
    :param m:
        Length of the vector to be transformed.
    :type m: int

    :return:
        Greatest power-of-two integer `n` satisfying `n <= m`.
    :rtype: int
    :return:
        Initialization array to be used by :func:`idd_sfrm`.
    :rtype: :class:`numpy.ndarray`
    """
    return _id.idd_sfrmi(l, m)


#------------------------------------------------------------------------------
# idd_id.f
#------------------------------------------------------------------------------

def iddp_id(eps, A):
    """
    Compute ID of a real matrix to a specified relative precision.

    :param eps:
        Relative precision.
    :type eps: float
    :param A:
        Matrix.
    :type A: :class:`numpy.ndarray`

    :return:
        Rank of ID.
    :rtype: int
    :return:
        Column index array.
    :rtype: :class:`numpy.ndarray`
    :return:
        Interpolation coefficients.
    :rtype: :class:`numpy.ndarray`
    """
    A = np.asfortranarray(A)
    k, idx, rnorms = _id.iddp_id(eps, A)
    n = A.shape[1]
    proj = A.T.ravel()[:k*(n-k)].reshape((k, n-k), order='F')
    return k, idx, proj


def iddr_id(A, k):
    """
    Compute ID of a real matrix to a specified rank.

    :param A:
        Matrix.
    :type A: :class:`numpy.ndarray`
    :param k:
        Rank of ID.
    :type k: int

    :return:
        Column index array.
    :rtype: :class:`numpy.ndarray`
    :return:
        Interpolation coefficients.
    :rtype: :class:`numpy.ndarray`
    """
    A = np.asfortranarray(A)
    idx, rnorms = _id.iddr_id(A, k)
    n = A.shape[1]
    proj = A.T.ravel()[:k*(n-k)].reshape((k, n-k), order='F')
    return idx, proj


def idd_reconid(B, idx, proj):
    """
    Reconstruct matrix from real ID.

    :param B:
        Skeleton matrix.
    :type B: :class:`numpy.ndarray`
    :param idx:
        Column index array.
    :type idx: :class:`numpy.ndarray`
    :param proj:
        Interpolation coefficients.
    :type proj: :class:`numpy.ndarray`

    :return:
        Reconstructed matrix.
    :rtype: :class:`numpy.ndarray`
    """
    B = np.asfortranarray(B)
    if proj.size > 0:
        return _id.idd_reconid(B, idx, proj)
    else:
        return B[:, np.argsort(idx)]


def idd_reconint(idx, proj):
    """
    Reconstruct interpolation matrix from real ID.

    :param idx:
        Column index array.
    :type idx: :class:`numpy.ndarray`
    :param proj:
        Interpolation coefficients.
    :type proj: :class:`numpy.ndarray`

    :return:
        Interpolation matrix.
    :rtype: :class:`numpy.ndarray`
    """
    return _id.idd_reconint(idx, proj)


def idd_copycols(A, k, idx):
    """
    Reconstruct skeleton matrix from real ID.

    :param A:
        Original matrix.
    :type A: :class:`numpy.ndarray`
    :param k:
        Rank of ID.
    :type k: int
    :param idx:
        Column index array.
    :type idx: :class:`numpy.ndarray`

    :return:
        Skeleton matrix.
    :rtype: :class:`numpy.ndarray`
    """
    A = np.asfortranarray(A)
    return _id.idd_copycols(A, k, idx)


#------------------------------------------------------------------------------
# idd_id2svd.f
#------------------------------------------------------------------------------

def idd_id2svd(B, idx, proj):
    """
    Convert real ID to SVD.

    :param B:
        Skeleton matrix.
    :type B: :class:`numpy.ndarray`
    :param idx:
        Column index array.
    :type idx: :class:`numpy.ndarray`
    :param proj:
        Interpolation coefficients.
    :type proj: :class:`numpy.ndarray`

    :return:
        Left singular vectors.
    :rtype: :class:`numpy.ndarray`
    :return:
        Right singular vectors.
    :rtype: :class:`numpy.ndarray`
    :return:
        Singular values.
    :rtype: :class:`numpy.ndarray`
    """
    B = np.asfortranarray(B)
    U, V, S, ier = _id.idd_id2svd(B, idx, proj)
    if ier:
        raise _RETCODE_ERROR
    return U, V, S


#------------------------------------------------------------------------------
# idd_snorm.f
#------------------------------------------------------------------------------

def idd_snorm(m, n, matvect, matvec, its=20):
    """
    Estimate spectral norm of a real matrix by the randomized power method.