// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
// License: Boost Software License   See LICENSE.txt for the full license.
#undef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_
#ifdef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_ 

#include "matrix_abstract.h"
#include <complex>

namespace dlib
{

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
//                             Global linear algebra functions 
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    const matrix_exp::matrix_type inv (
        const matrix_exp& m
    );
    /*!
        requires
            - m is a square matrix
        ensures
            - returns the inverse of m 
              (Note that if m is singular or so close to being singular that there
              is a lot of numerical error then the returned matrix will be bogus.  
              You can check by seeing if m*inv(m) is an identity matrix)
    !*/

// ----------------------------------------------------------------------------------------

    const matrix pinv (
        const matrix_exp& m
    );
    /*!
        ensures
            - returns the Moore-Penrose pseudoinverse of m.
            - The returned matrix has m.nc() rows and m.nr() columns.
    !*/

// ----------------------------------------------------------------------------------------

    void svd (
        const matrix_exp& m,
        matrix<matrix_exp::type>& u,
        matrix<matrix_exp::type>& w,
        matrix<matrix_exp::type>& v
    );
    /*!
        ensures
            - computes the singular value decomposition of m
            - m == #u*#w*trans(#v)
            - trans(#u)*#u == identity matrix
            - trans(#v)*#v == identity matrix
            - diag(#w) == the singular values of the matrix m in no 
              particular order.  All non-diagonal elements of #w are
              set to 0.
            - #u.nr() == m.nr()
            - #u.nc() == m.nc()
            - #w.nr() == m.nc()
            - #w.nc() == m.nc()
            - #v.nr() == m.nc()
            - #v.nc() == m.nc()
    !*/

// ----------------------------------------------------------------------------------------

    long svd2 (
        bool withu, 
        bool withv, 
        const matrix_exp& m,
        matrix<matrix_exp::type>& u,
        matrix<matrix_exp::type>& w,
        matrix<matrix_exp::type>& v
    );
    /*!
        requires
            - m.nr() >= m.nc()
        ensures
            - computes the singular value decomposition of matrix m
            - m == subm(#u,get_rect(m))*diagm(#w)*trans(#v)
            - trans(#u)*#u == identity matrix
            - trans(#v)*#v == identity matrix
            - #w == the singular values of the matrix m in no 
              particular order.  
            - #u.nr() == m.nr()
            - #u.nc() == m.nr()
            - #w.nr() == m.nc()
            - #w.nc() == 1 
            - #v.nr() == m.nc()
            - #v.nc() == m.nc()
            - if (widthu == false) then
                - ignore the above regarding #u, it isn't computed and its
                  output state is undefined.
            - if (widthv == false) then
                - ignore the above regarding #v, it isn't computed and its
                  output state is undefined.
            - returns an error code of 0, if no errors and 'k' if we fail to
              converge at the 'kth' singular value.
    !*/

// ----------------------------------------------------------------------------------------

    void svd3 (
        const matrix_exp& m,
        matrix<matrix_exp::type>& u,
        matrix<matrix_exp::type>& w,
        matrix<matrix_exp::type>& v
    );
    /*!
        ensures
            - computes the singular value decomposition of m
            - m == #u*diagm(#w)*trans(#v)
            - trans(#u)*#u == identity matrix
            - trans(#v)*#v == identity matrix
            - #w == the singular values of the matrix m in no 
              particular order.  
            - #u.nr() == m.nr()
            - #u.nc() == m.nc()
            - #w.nr() == m.nc()
            - #w.nc() == 1 
            - #v.nr() == m.nc()
            - #v.nc() == m.nc()
    !*/

// ----------------------------------------------------------------------------------------

    const matrix real_eigenvalues (
        const matrix_exp& m
    );
    /*!
        requires
            - m.nr() == m.nc()
            - matrix_exp::type == float or double
        ensures
            - returns a matrix E such that:
                - E.nr() == m.nr()
                - E.nc() == 1
                - E contains the real part of all eigenvalues of the matrix m.
                  (note that the eigenvalues are not sorted)
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::type det (
        const matrix_exp& m
    );
    /*!
        requires
            - m is a square matrix
        ensures
            - returns the determinant of m
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::type trace (
        const matrix_exp& m
    );
    /*!
        requires
            - m is a square matrix
        ensures
            - returns the trace of m
              (i.e. returns sum(diag(m)))
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::matrix_type chol (
        const matrix_exp& A
    );
    /*!
        requires
            - A is a square matrix
        ensures
            - if (A has a Cholesky Decomposition) then
                - returns the decomposition of A.  That is, returns a matrix L
                  such that L*trans(L) == A.  L will also be lower triangular.
            - else
                - returns a matrix with the same dimensions as A but it 
                  will have a bogus value.  I.e. it won't be a decomposition.
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::matrix_type inv_lower_triangular (
        const matrix_exp& A
    );
    /*!
        requires
            - A is a square matrix
        ensures
            - if (A is lower triangular) then
                - returns the inverse of A. 
            - else
                - returns a matrix with the same dimensions as A but it 
                  will have a bogus value.  I.e. it won't be an inverse.
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::matrix_type inv_upper_triangular (
        const matrix_exp& A
    );
    /*!
        requires
            - A is a square matrix
        ensures
            - if (A is upper triangular) then
                - returns the inverse of A. 
            - else
                - returns a matrix with the same dimensions as A but it 
                  will have a bogus value.  I.e. it won't be an inverse.
    !*/

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
//                             Matrix decomposition classes 
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    template <
        typename matrix_exp_type
        >
    class lu_decomposition
    {
        /*!
            REQUIREMENTS ON matrix_exp_type
                must be some kind of matrix expression as defined in the 
                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
                The matrix type must also contain float or double values.

            WHAT THIS OBJECT REPRESENTS
                This object represents something that can compute an LU 
                decomposition of a real valued matrix.  That is, for any 
                matrix A it computes matrices L, U, and a pivot vector P such 
                that rowm(A,P) == L*U.

                The LU decomposition with pivoting always exists, even if the matrix is
                singular, so the constructor will never fail.  The primary use of the
                LU decomposition is in the solution of square systems of simultaneous
                linear equations.  This will fail if is_singular() returns true (or
                if A is very nearly singular).
        !*/

    public:

        const static long NR = matrix_exp_type::NR;
        const static long NC = matrix_exp_type::NC;
        typedef typename matrix_exp_type::type type;
        typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
        typedef typename matrix_exp_type::layout_type layout_type;

        typedef matrix<type,0,0,mem_manager_type,layout_type>  matrix_type;
        typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;
        typedef matrix<long,NR,1,mem_manager_type,layout_type> pivot_column_vector_type;

        template <typename EXP>
        lu_decomposition (
            const matrix_exp<EXP> &A
        );
        /*!
            requires
                - EXP::type == lu_decomposition::type 
                - A.size() > 0
            ensures
                - #nr() == A.nr()
                - #nc() == A.nc()
                - #is_square() == (A.nr() == A.nc())
                - computes the LU factorization of the given A matrix.
        !*/

        bool is_square (
        ) const;
        /*!
            ensures
                - if (the input A matrix was a square matrix) then
                    - returns true
                - else
                    - returns false
        !*/

        bool is_singular (
        ) const;
        /*!
            requires
                - is_square() == true
            ensures
                - if (the input A matrix is singular) then
                    - returns true
                - else
                    - returns false
        !*/

        long nr(
        ) const;
        /*!
            ensures
                - returns the number of rows in the input matrix
        !*/

        long nc(
        ) const;
        /*!
            ensures
                - returns the number of columns in the input matrix
        !*/

        const matrix_type get_l (
        ) const; 
        /*!
            ensures
                - returns the lower triangular L factor of the LU factorization.  
                - L.nr() == nr()
                - L.nc() == min(nr(),nc())
        !*/

        const matrix_type get_u (
        ) const;
        /*!
            ensures
                - returns the upper triangular U factor of the LU factorization.  
                - U.nr() == min(nr(),nc())
                - U.nc() == nc()
        !*/

        const pivot_column_vector_type& get_pivot (
        ) const;
        /*!
            ensures
                - returns the pivot permutation vector.  That is,
                  if A is the input matrix then this function 
                  returns a vector P such that:
                    - rowm(A,P) == get_l()*get_u() 
                    - P.nr() == A.nr()
        !*/

        type det (
        ) const;
        /*!
            requires
                - is_square() == true
            ensures
                - computes and returns the determinant of the input 
                  matrix using LU factors.
        !*/

        template <typename EXP>
        const matrix_type solve (
            const matrix_exp<EXP> &B
        ) const;
        /*!
            requires
                - EXP::type == lu_decomposition::type
                - is_square() == true
                - B.nr() == nr()
            ensures
                - Let A denote the input matrix to this class's constructor.  
                  Then this function solves A*X == B for X and returns X.  
                - Note that if A is singular (or very close to singular) then
                  the X returned by this function won't fit A*X == B very well (if at all).
        !*/

    }; 

// ----------------------------------------------------------------------------------------

    template <
        typename matrix_exp_type
        >
    class cholesky_decomposition
    {
        /*! 
            REQUIREMENTS ON matrix_exp_type
                must be some kind of matrix expression as defined in the 
                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
                The matrix type must also contain float or double values.

            WHAT THIS OBJECT REPRESENTS
                This object represents something that can compute a cholesky 
                decomposition of a real valued matrix.  That is, for any 
                symmetric, positive definite matrix A, it computes a lower 
                triangular matrix L such that A == L*trans(L).
                
                If the matrix is not symmetric or positive definite, the function
                computes only a partial decomposition.  This can be tested with
                the is_spd() flag.
        !*/

    public:

        const static long NR = matrix_exp_type::NR;
        const static long NC = matrix_exp_type::NC;
        typedef typename matrix_exp_type::type type;
        typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
        typedef typename matrix_exp_type::layout_type layout_type;

        typedef typename matrix_exp_type::matrix_type matrix_type;
        typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;

        template <typename EXP>
        cholesky_decomposition(
            const matrix_exp<EXP>& A
        );
        /*!
            requires
                - EXP::type == cholesky_decomposition::type 
                - A.size() > 0
                - A.nr() == A.nc() 
                  (i.e. A must be a square matrix)
            ensures
                - if (A is symmetric positive-definite) then
                    - #is_spd() == true 
                    - Constructs a lower triangular matrix L, such that L*trans(L) == A.
                      and #get_l() == L
                - else
                    - #is_spd() == false
        !*/

        bool is_spd(
        ) const;
        /*!
            ensures
                - if (the input matrix was symmetric positive-definite) then
                    - returns true
                - else
                    - returns false
        !*/

        const matrix_type& get_l(
        ) const;
        /*!
            ensures
                - returns the lower triangular factor, L, such that L*trans(L) == A
                  (where A is the input matrix to this class's constructor)
                - Note that if A is not symmetric positive definite or positive semi-definite
                  then the equation L*trans(L) == A won't hold.  
        !*/

        template <typename EXP>
        const matrix solve (
            const matrix_exp<EXP>& B
        ) const;
        /*!
            requires
                - EXP::type == cholesky_decomposition::type
                - B.nr() == get_l().nr()
                  (i.e. the number of rows in B must match the number of rows in the
                  input matrix A)
            ensures
                - Let A denote the input matrix to this class's constructor.  Then 
                  this function solves A*X = B for X and returns X.  
                - Note that if is_spd() == false or A was really close to being
                  non-SPD then the solver will fail to find an accurate solution.
        !*/

    };

// ----------------------------------------------------------------------------------------

    template <
        typename matrix_exp_type
        >
    class qr_decomposition 
    {
        /*! 
            REQUIREMENTS ON matrix_exp_type
                must be some kind of matrix expression as defined in the 
                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
                The matrix type must also contain float or double values.

            WHAT THIS OBJECT REPRESENTS
                This object represents something that can compute a classical
                QR decomposition of an m-by-n real valued matrix A with m >= n.  

                The QR decomposition is an m-by-n orthogonal matrix Q and an 
                n-by-n upper triangular matrix R so that A == Q*R. The QR decomposition 
                always exists, even if the matrix does not have full rank, so the 
                constructor will never fail.  The primary use of the QR decomposition 
                is in the least squares solution of non-square systems of simultaneous 
                linear equations.  This will fail if is_full_rank() returns false or
                A is very nearly not full rank.

                The Q and R factors can be retrieved via the get_q() and get_r()
                methods. Furthermore, a solve() method is provided to find the
                least squares solution of Ax=b using the QR factors.  
        !*/

    public:

        const static long NR = matrix_exp_type::NR;
        const static long NC = matrix_exp_type::NC;
        typedef typename matrix_exp_type::type type;
        typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
        typedef typename matrix_exp_type::layout_type layout_type;

        typedef matrix<type,0,0,mem_manager_type,layout_type> matrix_type;
        typedef matrix<type,0,1,mem_manager_type,layout_type> column_vector_type;

        template <typename EXP>
        qr_decomposition(
            const matrix_exp<EXP>& A
        );
        /*!
            requires
                - EXP::type == qr_decomposition::type
                - A.nr() >= A.nc()
                - A.size() > 0
            ensures
                - #nr() == A.nr()
                - #nc() == A.nc()
                - computes the QR decomposition of the given A matrix.
        !*/

        bool is_full_rank(
        ) const;
        /*!
            ensures
                - if (the input A matrix had full rank) then
                    - returns true
                - else
                    - returns false
        !*/

        long nr(
        ) const;
        /*!
            ensures
                - returns the number of rows in the input matrix
        !*/

        long nc(
        ) const;
        /*!
            ensures
                - returns the number of columns in the input matrix
        !*/

        const matrix_type get_householder (
        )  const;
        /*!
            ensures
                - returns a matrix H such that:
                    - H is the lower trapezoidal matrix whose columns define the 
                      Householder reflection vectors from QR factorization 
                    - H.nr() == nr()
                    - H.nc() == nc()
        !*/

        const matrix_type get_r (
        ) const;
        /*!
            ensures
                - returns a matrix R such that: 
                    - R is the upper triangular factor, R, of the QR factorization
                    - get_q()*R == input matrix A
                    - R.nr() == nc()
                    - R.nc() == nc()
        !*/

        const matrix_type get_q (
        ) const;
        /*!
            ensures
                - returns a matrix Q such that:
                    - Q is the economy-sized orthogonal factor Q from the QR 
                      factorization.  
                    - trans(Q)*Q == identity matrix
                    - Q*get_r() == input matrix A 
                    - Q.nr() == nr()
                    - Q.nc() == nc()
        !*/

        template <typename EXP>
        const matrix_type solve (
            const matrix_exp<EXP>& B
        ) const;
        /*!
            requires
                - EXP::type == qr_decomposition::type
                - B.nr() == nr()
            ensures
                - Let A denote the input matrix to this class's constructor.  
                  Then this function finds the least squares solution to the equation A*X = B 
                  and returns X.  X has the following properties: 
                    - X is the matrix that minimizes the two norm of A*X-B.  That is, it
                      minimizes sum(squared(A*X - B)).
                    - X.nr() == nc()
                    - X.nc() == B.nc()
                - Note that this function will fail to output a good solution if is_full_rank() == false
                  or the A matrix is close to not being full rank.
        !*/

    };

// ----------------------------------------------------------------------------------------

    template <
        typename matrix_exp_type
        >
    class eigenvalue_decomposition
    {
        /*!
            REQUIREMENTS ON matrix_exp_type
                must be some kind of matrix expression as defined in the 
                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
                The matrix type must also contain float or double values.

            WHAT THIS OBJECT REPRESENTS
                This object represents something that can compute an eigenvalue 
                decomposition of a real valued matrix.   So it gives 
                you the set of eigenvalues and eigenvectors for a matrix.   

                Let A denote the input matrix to this object's constructor.  Then 
                what this object does is it finds two matrices, D and V, such that
                    - A*V == V*D
                Where V is a square matrix that contains all the eigenvectors
                of the A matrix (each column of V is an eigenvector) and
                D is a diagonal matrix containing the eigenvalues of A.


                It is important to note that if A is symmetric or non-symmetric you
                get somewhat different results.  If A is a symmetric matrix (i.e. A == trans(A))
                then:
                    - All the eigenvalues and eigenvectors of A are real numbers. 
                        - Because of this there isn't really any point in using the
                          part of this class's interface that returns complex matrices.
                          All you need are the get_real_eigenvalues() and
                          get_pseudo_v() functions.  
                    - V*trans(V) should be equal to the identity matrix.  That is, all the
                      eigenvectors in V should be orthonormal. 
                        - So A == V*D*trans(V)

                On the other hand, if A is not symmetric then:
                    - Some of the eigenvalues and eigenvectors might be complex numbers.  
                        - An eigenvalue is complex if and only if its corresponding eigenvector 
                          is complex.  So you can check for this case by just checking 
                          get_imag_eigenvalues() to see if any values are non-zero.  You don't 
                          have to check the V matrix as well.
                    - V*trans(V) won't be equal to the identity matrix but it is usually
                      invertible.  So A == V*D*inv(V) is usually a valid statement but
                      A == V*D*trans(V) won't be.
        !*/

    public:

        const static long NR = matrix_exp_type::NR;
        const static long NC = matrix_exp_type::NC;
        typedef typename matrix_exp_type::type type;
        typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
        typedef typename matrix_exp_type::layout_type layout_type;

        typedef typename matrix_exp_type::matrix_type matrix_type;
        typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;

        typedef matrix<std::complex<type>,0,0,mem_manager_type,layout_type> complex_matrix_type;
        typedef matrix<std::complex<type>,NR,1,mem_manager_type,layout_type> complex_column_vector_type;


        template <typename EXP>
        eigenvalue_decomposition(
            const matrix_exp<EXP>& A
        ); 
        /*!
            requires
                - A.nr() == A.nc() 
                - A.size() > 0
                - EXP::type == eigenvalue_decomposition::type 
            ensures
                - #dim() == A.nr()
                - computes the eigenvalue decomposition of A.  
                - #get_eigenvalues() == the eigenvalues of A
                - #get_v() == all the eigenvectors of A
        !*/

        long dim (
        ) const;
        /*!
            ensures
                - dim() == the number of rows/columns in the input matrix A 
        !*/

        const complex_column_vector_type get_eigenvalues (
        ) const;
        /*!
            ensures
                - returns diag(get_d()).  That is, returns a 
                  vector that contains the eigenvalues of the input 
                  matrix.
                - the returned vector has dim() rows
                - the eigenvalues are not sorted in any particular way
        !*/

        const column_vector_type& get_real_eigenvalues (
        ) const;
        /*! 
            ensures
                - returns the real parts of the eigenvalues.  That is,
                  returns real(get_eigenvalues()) 
                - the returned vector has dim() rows
                - the eigenvalues are not sorted in any particular way
        !*/

        const column_vector_type& get_imag_eigenvalues (
        ) const;
        /*! 
            ensures
                - returns the imaginary parts of the eigenvalues.  That is,
                  returns imag(get_eigenvalues()) 
                - the returned vector has dim() rows
                - the eigenvalues are not sorted in any particular way
        !*/

        const complex_matrix_type get_v (
        ) const;
        /*!
            ensures
                - returns the eigenvector matrix V that is 
                  dim() rows by dim() columns
                - Each column in V is one of the eigenvectors of the input 
                  matrix
        !*/

        const complex_matrix_type get_d (
        ) const; 
        /*!
            ensures
                - returns a matrix D such that:
                    - D.nr() == dim()
                    - D.nc() == dim()
                    - diag(D) == get_eigenvalues()
                      (i.e. the diagonal of D contains all the eigenvalues in the input matrix)
                    - all off diagonal elements of D are set to 0
        !*/

        const matrix_type& get_pseudo_v (
        ) const;
        /*!
            ensures
                - returns a matrix that is dim() rows by dim() columns
                - Let A denote the input matrix given to this object's constructor.
                - if (A has any imaginary eigenvalues) then
                    - returns the pseudo-eigenvector matrix V  
                    - The matrix V returned by this function is structured such that:
                        - A*V == V*get_pseudo_d()
                - else
                    - returns the eigenvector matrix V with A's eigenvectors as
                      the columns of V
                    - A*V == V*diagm(get_real_eigenvalues())
        !*/

        const matrix_type get_pseudo_d (
        ) const; 
        /*!
            ensures
                - The returned matrix is dim() rows by dim() columns
                - Computes and returns the block diagonal eigenvalue matrix.
                  If the original matrix A is not symmetric, then the eigenvalue 
                  matrix D is block diagonal with the real eigenvalues in 1-by-1 
                  blocks and any complex eigenvalues,
                  a + i*b, in 2-by-2 blocks, (a, b; -b, a).  That is, if the complex
                  eigenvalues look like

                      u + iv     .        .          .      .    .
                        .      u - iv     .          .      .    .
                        .        .      a + ib       .      .    .
                        .        .        .        a - ib   .    .
                        .        .        .          .      x    .
                        .        .        .          .      .    y

                  Then D looks like

                        u        v        .          .      .    .
                       -v        u        .          .      .    . 
                        .        .        a          b      .    .
                        .        .       -b          a      .    .
                        .        .        .          .      x    .
                        .        .        .          .      .    y

                  This keeps V (The V you get from get_pseudo_v()) a real matrix in both 
                  symmetric and non-symmetric cases, and A*V = V*D.
                - the eigenvalues are not sorted in any particular way
        !*/

    };

// ----------------------------------------------------------------------------------------

}

#endif // DLIB_MATRIx_LA_FUNCTS_ABSTRACT_

