# ifndef TMATRIXLIB_HPP
# define TMATRIXLIB_HPP

#include "TLib.hpp"

// Transpose a matrix formed by any STL container of the form container<container<T>>

/**
Usage:
deque<deque<bool>> tdata = transpose<bool, deque<bool>>(datanull);
also:
printcol(transpose<bool,deque<bool>>(sequences)[snp]);


**/
template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
inline void transpose(const Tcolofcol& data, Tcolofcol& tdata) {
    // Clean if necessary
    if(!tdata.empty()){
        for(auto& element:tdata){
            element.clear();
            Tcol().swap(element);
        }
        tdata.clear();
        Tcolofcol().swap(tdata);
    }

    // This assumes that all inner vectors have the same size and
    // allocates space for the complete result in advance
    tdata.resize(data[0].size(), Tcol(data.size()));

    for (typename Tcol::size_type i = 0; i < data[0].size(); ++i)
        for (typename Tcol::size_type j = 0; j < data.size(); ++j) {
            tdata[i][j] = data[j][i];
        }
}


// Transpose in situ a matrix formed by any STL container of the form container<container<T>>
// the matrix passed as argument is transposed
/**
Usage:

vector<vector<double>>data(100,vector<double>(200000,0));
int id=1;
for(auto& d:data)
    for(auto& dd:d){
        dd+=id;
        ++id;
    }

cout<<data.size()<<" x " << data[0].size()<<endl;

//for(auto d:data)
//    printcol(d);

cout<<endl;
transposeinsitu(data);

cout<<data.size()<<" x " << data[0].size()<<endl;


**/
template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
inline void transposeinsitu(Tcolofcol& data) {

    Tcolofcol tdata;

    // This assumes that all inner vectors have the same size and
    // allocates space for the complete result in advance
    tdata.resize(data[0].size(), Tcol(data.size()));

    for (typename Tcol::size_type i = 0; i < data[0].size(); ++i)
        for (typename Tcol::size_type j = 0; j < data.size(); ++j) {
            tdata[i][j] = data[j][i];
        }
    data = tdata;

}

// Another way that do not copy the data. This way is useful just to access or write the trasposed elements as in the example but not for storing and managing the elements i.e. wrapper do not has the .begin() etc
/**
Usage:
typedef std::vector<std::vector<int> > TwoDVector;
TwoDVector  data(3,std::vector<int>(2,0));
TwoDPivotWrapper<TwoDVector> wrapper(data);

cout<<data[2][1]<<"\t"<<wrapper[1][2]<<endl;

**/
template<typename T>
class TwoDPivotWrapper
{
    public:
        // These two typedef's were done with std::vector
        // in mind. But with a small amount of effort I am
        // sure they can be generalized. They are solely to define
        // value_type (the data stored in the 2-D array).
   // In the STL containers value_type is the type of the elements in the container
        typedef typename T::value_type OneDType; // OneDType is the type (might be a container) of the elements in container T
        typedef typename OneDType::value_type value_type; // value_type is the type of the elements in container OneDType

        // A constructor that wraps a 2-D structure.
        TwoDPivotWrapper(T& o)
            : object(o)
        {}

        // A helper class used to store the row after the first array accesses.
        class Row
        {
            friend class TwoDPivotWrapper;
            Row(TwoDPivotWrapper& w, size_t r)
                : wrapper(w)
                , row(r)
            {}
            TwoDPivotWrapper&    wrapper;
            size_t               row;

            public:
                value_type operator[](size_t col)
                {
                    return wrapper.get(row,col);
                }
        };

        // The operator [] returns a Row object that overloads
        // the operator[] for the next dimension.
        Row operator[](size_t row)              {return Row(*this, row);}

        // Generic get function used to access elements.
        // Notice we swap the row/col around when accessing
        // the underlying object.
        value_type get(size_t row, size_t col)  {return object[col][row];}

    private:
        T&  object;
};


template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
inline void mproduct(const Tcolofcol& mat1, const Tcolofcol& mat2, Tcolofcol& matP)
{

    Tcolofcol Tmat2;
    transpose(mat2,Tmat2);

    matP.resize(mat1.size(), Tcol(Tmat2.size()));

    size_t r=0, c=0;
    for(auto rowm1: mat1){
        for(auto rowm2: Tmat2){
            matP[r][c] = inner_product(rowm1.begin(), rowm1.end(),rowm2.begin(),T(0));
            ++c;
        }
        c=0;
        ++r;
    }
}

template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
inline void escalarbyM(const T& scalar, const Tcolofcol& mat, Tcolofcol& matP){

    matP =mat;
    for(auto& row: matP)
        for(auto& element: row)
            element*=scalar;

}

//// Kronecker product: mat1 X mat2. Implementation by blocks
//template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
//inline void kronecker(const Tcolofcol& mat1, const Tcolofcol& mat2, Tcolofcol& matK)
//{
//    size_t m = mat1.size(), n= mat1[0].size(), p = mat2.size(), q=mat2[0].size();
//    size_t nr = m*p;
//    size_t nc = n*q;
//    matK.resize(nr,Tcol(nc));
//
//    size_t rini=0,cini=0,r=0, c=0;
//    for(const auto& row1: mat1){
//        for(const auto& element1: row1){
//
//            for(const auto& row2: mat2){
//                for(const auto& element2: row2){
//                    matK[r][c] = element1*element2;
//                    ++c;
//                }
//                ++r;
//                c=cini;
//            }
//
//
//            r=rini;
//            cini+=q;
//            c=cini;
//
//        }
//        cini=0;
//        c=cini;
//        rini+=p;
//        r=rini;
//
//    }
//
//}

//// Kronecker product: mat1 X mat2. Implementation by rows
template <typename T= double, typename Tcol= deque<T>, typename Tcolofcol = deque<Tcol> >
inline void kronecker(const Tcolofcol& mat1, const Tcolofcol& mat2, Tcolofcol& matK)
{
    size_t m = mat1.size(), n= mat1[0].size(), p = mat2.size(), q=mat2[0].size();
    size_t nr = m*p;
    size_t nc = n*q;

    Tcol temp;
    if (nc<Tcol().max_size())
        temp.resize(nc);
    else {cout << "kronecker: That size exceeds the limit.\n";exit(3);}

    size_t r=0, c=0;

    for(const auto& row1: mat1){
        for(const auto& row2: mat2){
            for(const auto& element1: row1){
                for(const auto& element2: row2){
                    temp[c] = element1*element2;
                    ++c;
                }

            }
            matK.push_back(temp);
            ++r;
            c=0;
        }
    }

}

//template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
//inline void slowmproduct(const Tcolofcol& mat1, const Tcolofcol& mat2, Tcolofcol& matP)
//{
//
//    Tcolofcol Tmat2;
//    size_t row_a=mat1.size();
//    size_t row_b=mat2.size();
//    size_t col_b= mat2[0].size();
//
//    matP.resize(row_a, Tcol(col_b,T(0)));
//
//    size_t r=0, c=0;
//
//    for(int i = 0; i < row_a; i++)
//	  for(int j = 0; j < col_b; j++)
//		for(int k = 0; k < row_b; k++)
//			matP[i][j] += mat1[i][k]*mat2[k][j];
//}

//
//

/// LU descomposition and inverse form NRC++ p. 49

template <typename Tcol1= vector<int>, typename Tcol2= vector<double>, typename Tcolofcol = vector<Tcol2> >
inline void ludcmp( Tcolofcol& a, Tcol1& index, double& d){
    const double TINY=1.0e-20;
    int imax;
    double big,dum,sum,temp;

    int n=a.size();
    Tcol2 vv(n);
    d=1.0;
    for(int i=0; i<n; ++i){
        big=0.0;
        for(int j=0;j<n;++j){
            if((temp=abs(a[i][j]))> big) big=temp;
        }
        if(big==0.0) {printval<int>(1,"Singular matrix in ludcmp","Errfact.log");exit(1);}
        vv[i]=1.0/big;
    } // close for i
    for(int j=0;j<n;++j){
        for(int i=0; i<j; ++i){
            sum = a[i][j];
            for(int k=0; k<i; ++k) sum-= a[i][k]*a[k][j];
            a[i][j] = sum;
        }
        big=0.0;
        for(int i=j; i<n; ++i){
            sum=a[i][j];
            for(int k=0; k<j; ++k) sum-= a[i][k]*a[k][j];
            a[i][j] = sum;
            if((dum=vv[i]*abs(sum))>=big){
                big=dum;
                imax=i;
            }
        } // close for i
        if(j!=imax){
            for(int k=0; k<n; ++k){
                swap(a[imax][k],a[j][k]);
            }
            d=-d;
            vv[imax]=vv[j];
        } // close if
        index[j]=imax;
        if(a[j][j]==0) a[j][j]=TINY;
        if(j != n-1){
            dum = 1.0/a[j][j];
            for(int i = j+1; i<n; ++i) a[i][j]*=dum;
        } // close if
    } // close for j
}

template <typename Tcol1= vector<int>, typename Tcol2= vector<double>, typename Tcolofcol = vector<Tcol2> >
inline void lubksb(const Tcolofcol& a, const Tcol1& index, Tcol2& b){
    int ii=0,ip;
    double sum;
    int n=a.size();
    for(int i=0; i<n; ++i){
        ip=index[i];
        sum=b[ip];
        b[ip]=b[i];
        if(ii!=0)
            for(int j=ii-1; j<i; ++j) sum -= a[i][j]*b[j];
        else if(sum != 0.0)
            ii=i+1;
        b[i]=sum;
    } // close for i
    for(int i = n-1; i>=0; --i){
        sum=b[i];
        for(int j = i+1; j<n; ++j) sum -= a[i][j]*b[j];
        b[i] = sum/a[i][i];
    } // close for i
}
///Inverse of a matrix: NRC++ p. 51
template <typename Tcol1= vector<int>, typename Tcol2= vector<double>, typename Tcolofcol = vector<Tcol2> >
inline void inverseM(const Tcolofcol& M, Tcolofcol& InvM){
    int N = M.size();
    Tcol2 col(N);
    Tcol1 index(N);
    double d;
    Tcolofcol a = M; // because a will be destroyed be do not want to use M
    InvM.resize(N,Tcol2(N));
// Decompose the matrix
    ludcmp(a, index,d);
// Find inverse by columns
    for(int j=0; j < N; ++j){
        for(int i=0; i < N; ++i) col[i]=0;
        col[j]=1;
        lubksb(a, index, col);
        for(int i=0; i < N; ++i) InvM[i][j] = col[i];
    }

}


# endif
