dyncmatx.h

#pragma once
#define TEXT_MODE

#ifdef TEXT_MODE
#include <iostream>
#include <iomanip>
#endif

#include <vector>
#include <cmath>
#include <random>
#include <string>
#include <iterator>
#include <algorithm>

using namespace std;

typedef unsigned short int uint;

// ESCALAR POR MATRIZ
template <typename S>
inline vector<vector<S>> multiplicarPorEscalar(const vector<vector<S>>& matriz,const S& escalar) {
    int filas = matriz.size();
    int columnas = matriz[0].size();

    // Inicializar la matriz resultado con el mismo tamaño que la matriz original
    vector<vector<S>> resultado(filas, vector<S>(columnas, 0.0));

    // Multiplicar cada elemento de la matriz por el escalar y almacenar el resultado
    for (uint i = 0; i < filas; ++i) {
        for (uint j = 0; j < columnas; ++j) {
            resultado[i][j] = matriz[i][j] * escalar;
        }
    }

    return resultado;
}
template <typename S>
inline vector<S> multiplicarPorEscalar(const vector<S>& vectorOriginal,const S& escalar) {
    int dimension = vectorOriginal.size();
    vector<S> resultado(dimension);

    for (uint i = 0; i < dimension; ++i) {
        resultado[i] = vectorOriginal[i] * escalar;
    }

    return resultado;
}

// Función para intercambiar dos filas de la matriz
template <typename S>
void swapRows(vector<vector<S>>& matrix,const uint& row1,const uint& row2) {
    vector<S> temp = matrix[row1];
    matrix[row1] = matrix[row2];
    matrix[row2] = temp;
}

// Función para multiplicar una fila por un escalar
template <typename S>
void scaleRow(vector<vector<S>>& matrix,const uint& row,const S& scalar) {
    for (uint i = 0; i != matrix[row].size(); ++i) {
        matrix[row][i] *= scalar;
    }
}

// Función para restar un múltiplo de una fila de otra fila
template <typename S>
void subtractRows(vector<vector<S>>& matrix,const uint& targetRow,const uint& sourceRow,const S& scalar) {
    for (uint i = 0; i != matrix[targetRow].size(); ++i) {
        matrix[targetRow][i] -= scalar * matrix[sourceRow][i];
    }
}

// Función para calcular la determinante de una matriz usando el método de Gauss-Jordan
template <typename S>
inline S detn(vector<vector<S>>& matrix) {
    uint n = matrix.size(); // Número de filas (y asumimos que es igual al número de columnas)

    S det = 1.0;

    for (uint i = 0; i < n; ++i) {
        // Encuentra la fila con el máximo valor absoluto en la columna actual
        uint maxRow = i;
        for (uint j = i + 1; j != n; ++j) {
            if (abs(matrix[j][i]) > abs(matrix[maxRow][i])) {
                maxRow = j;
            }
        }

        // Intercambia la fila actual con la fila con el máximo valor absoluto
        if (maxRow != i) {
            swapRows(matrix, i, maxRow);
            det *= -1.0; // Cambia el signo de la determinante por el intercambio de filas
        }

        // Hace que el elemento diagonal sea 1
        S diagonalElement = matrix[i][i];
        if (diagonalElement == 0.0) {
            return 0.0; // La matriz es singular, la determinante es cero
        }
        det *= diagonalElement; // Multiplica la determinante por el valor diagonal

        scaleRow(matrix, i, 1.0 / diagonalElement);

        // Hace ceros en los elementos por encima y por debajo del elemento diagonal actual
        for (uint j = 0; j != n; ++j) {
            if (j != i) {
                subtractRows(matrix, j, i, matrix[j][i]);
            }
        }
    }

    return det;
}

// Función para calcular la inversa de una matriz usando el método de Gauss-Jordan
template <typename S>
inline vector<vector<S>> invn(vector<vector<S>>& matrix) {
    uint n = matrix.size(); // Número de filas (y asumimos que es igual al número de columnas)

    // Crear una matriz identidad del mismo tamaño que la matriz original
    vector<vector<S>> identity(n, vector<S>(n, 0.0));
    for (uint i = 0; i < n; ++i) {
        identity[i][i] = 1.0;
    }

    // Aplicar el método de Gauss-Jordan simultáneamente a la matriz original y la matriz identidad
    for (uint i = 0; i < n; ++i) {
        // Hace que el elemento diagonal sea 1
        S diagonalElement = matrix[i][i];
        if (diagonalElement == 0.0) {
            cerr << "La matriz no tiene inversa (es singular)." << endl;
            exit(EXIT_FAILURE);
        }
        scaleRow(matrix, i, 1.0 / diagonalElement);
        scaleRow(identity, i, 1.0 / diagonalElement);

        // Hace ceros en los elementos por encima y por debajo del elemento diagonal actual
        for (uint j = 0; j < n; ++j) {
            if (j != i) {
                S scalar = matrix[j][i];
                subtractRows(matrix, j, i, scalar);
                subtractRows(identity, j, i, scalar);
            }
        }
    }

    return identity;
}

// TRANSPUESTA
template <typename S>
inline vector<vector<S>> transponerMatriz(const vector<vector<S>>& matriz) {
    uint filas = matriz.size();
    uint columnas = matriz[0].size();

    // Inicializar la matriz transpuesta con las dimensiones intercambiadas
    vector<vector<S>> transpuesta(columnas, vector<S>(filas, 0.0));

    // Llenar la matriz transpuesta intercambiando filas y columnas
    for (uint i = 0; i < filas; ++i) {
        for (uint j = 0; j < columnas; ++j) {
            transpuesta[j][i] = matriz[i][j];
        }
    }

    return transpuesta;
}

// SUMA CON CONSTANTE
template <typename S>
inline vector<vector<S>> sumarConstanteMatriz(const vector<vector<S>>& matriz,const S& constante) {
    uint filas = matriz.size();
    uint columnas = matriz[0].size();

    // Inicializar la matriz resultado con el mismo tamaño que la matriz original
    vector<vector<S>> resultado(filas, vector<S>(columnas, 0.0));

    // Sumar la constante a cada elemento de la matriz y almacenar el resultado
    for (int i = 0; i < filas; ++i) {
        for (int j = 0; j < columnas; ++j) {
            resultado[i][j] = matriz[i][j] + constante;
        }
    }

    return resultado;
}

// PRODUCTO DE MATRICES
template <typename S>
inline vector<vector<S>> multiplicarMatrices(const vector<vector<S>>& matriz1, const vector<vector<S>>& matriz2) {
    uint filas1 = matriz1.size();
    uint columnas1 = matriz1[0].size();
    uint filas2 = matriz2.size();
    uint columnas2 = matriz2[0].size();

    // Verificar si las matrices son multiplicables
    if (columnas1 != filas2) {
        cerr << "Error: Las matrices no son multiplicables." << endl;
        exit(EXIT_FAILURE);
    }

    // Inicializar la matriz resultado con el tamaño adecuado
    vector<vector<S>> resultado(filas1, vector<S>(columnas2, 0.0));

    // Realizar la multiplicación de matrices
    for (uint i = 0; i < filas1; ++i) {
        for (uint j = 0; j < columnas2; ++j) {
            for (uint k = 0; k < columnas1; ++k) {
                resultado[i][j] += matriz1[i][k] * matriz2[k][j];
            }
        }
    }

    return resultado;
}

// SUMA ENTRE MATRICES
template <typename S>
inline vector<vector<S>> sumarMatrices(const vector<vector<S>>& matriz1, const vector<vector<S>>& matriz2) {
    uint filas = matriz1.size();
    uint columnas = matriz1[0].size();

    // Verificar si las matrices tienen las mismas dimensiones
    if (matriz2.size() != filas || matriz2[0].size() != columnas) {
        cerr << "Error: Las matrices no tienen las mismas dimensiones para la suma." << endl;
        exit(EXIT_FAILURE);
    }

    // Inicializar la matriz resultado con el tamaño adecuado
    vector<vector<S>> resultado(filas, vector<S>(columnas, 0.0));

    // Realizar la suma de matrices
    for (uint i = 0; i < filas; ++i) {
        for (uint j = 0; j < columnas; ++j) {
            resultado[i][j] = matriz1[i][j] + matriz2[i][j];
        }
    }

    return resultado;
}

template <typename S>
inline vector<vector<S>> restarMatrices(const vector<vector<S>>& matriz1, const vector<vector<S>>& matriz2) {
    uint filas = matriz1.size();
    uint columnas = matriz1[0].size();

    // Verificar si las matrices tienen las mismas dimensiones
    if (matriz2.size() != filas || matriz2[0].size() != columnas) {
        cerr << "Error: Las matrices no tienen las mismas dimensiones para la suma." << endl;
        exit(EXIT_FAILURE);
    }

    // Inicializar la matriz resultado con el tamaño adecuado
    vector<vector<S>> resultado(filas, vector<S>(columnas, 0.0));

    // Realizar la suma de matrices
    for (uint i = 0; i < filas; ++i) {
        for (uint j = 0; j < columnas; ++j) {
            resultado[i][j] = matriz1[i][j] - matriz2[i][j];
        }
    }

    return resultado;
}

// MAIN CLASS

// ATENCIÓN Y CUIDADO

template <class T> class dyncmatx{
    private:
        uint d_r,d_c;
        vector<vector<T>> dyn;
        void initDyn(){
            dyn.reserve(d_r);
            for(uint i = 0;i < d_r;i++){
                vector<T> t1;
                t1.reserve(d_c);
                for(uint j = 0;j < d_c;j++){
                    t1.emplace_back(0);
                }
                dyn.emplace_back(t1);
                vector<T>().swap(t1);
            }
        }
    public:   
        dyncmatx(){} 
        dyncmatx(vector<vector<T>> dyn): dyn(dyn){ d_r = dyn.size(); d_c = dyn.at(0).size(); }
        dyncmatx(uint d_r,uint d_c): d_r(d_r),d_c(d_c){initDyn();}
        dyncmatx(vector<vector<T>> dyn,uint d_r,uint d_c): dyn(dyn),d_r(d_r),d_c(d_c){initDyn();}
        ~dyncmatx(){d_c = 0;d_r = 0; vector<vector<T>>().swap(dyn);}

        // CREACION Y COMPOSICIÓN

        inline vector<vector<T>> get_dyncmatx(){return dyn;}

        void resize(const uint& new_r,const uint& new_c){
            dyncmatx<T> change(new_r,new_c);
            vector<vector<T>> dyn_temp = dyn;
            dyn = change.get_dyncmatx();
            for(uint i = 0; i != d_r;i++){
                for(uint j = 0;j != d_c;j++){
                    dyn[i][j] = dyn_temp[i][j];
                }
            }
        }
        void set_dyncmatx(vector<vector<T>> in){dyn = in; d_r = in.size(); d_c = in.at(0).size();}

        void sd_element(const int& row, const int& col,const T& val){
            dyn.at( (row)%(d_r) + ( (row < 0)?1:0 ) * d_r ).at( (col)%(d_c) + ( (col < 0)?1:0 ) * d_c ) = val;
        }
        inline T gd_element(const int& row,const int& col){
            return dyn.at( (row)%(d_r) + ( (row < 0)?1:0 ) * d_r ).at( (col)%(d_c) + ( (col < 0)?1:0 ) * d_c );
        }
        
        uint r_len(){return d_r;}
        uint c_len(){return d_c;}

        void check_size(){
            d_r = dyn.size(); d_c = dyn.at(0).size();
        }

        void add_r_data(const vector<T>& data_in){
            dyn.emplace_back(data_in);
            check_size();
        }
        void add_c_data(const vector<T>& data_in){
            for(uint i = 0; i != data_in.size();i++){
                dyn[(i)%(d_r)].emplace_back(data_in[i]);
            }
            check_size();
        }

        void add_r_data(const vector<vector<T>>& data_in){
            for(uint i = 0; i != data_in.size();i++){
                dyn.emplace_back(data_in.at(i));
            }
            check_size();
        }
        void add_c_data(const vector<vector<T>>& data_in){
            for(uint i = 0; i != data_in.size();i++){
                add_c_data(data_in.at(i));
            }
            check_size();
        }
        void add_r_dyn(dyncmatx<T> data_m){
            add_r_data(data_m.get_dyncmatx());
            check_size();
        }

        void add_c_dyn(dyncmatx<T> data_m){
            add_c_data(data_m.get_dyncmatx());
            check_size();
        }
        void add_c_dyn_fixed(dyncmatx<T> data_m){
            ~data_m;
            add_c_data(data_m.get_dyncmatx());
            check_size();
        }

        void insert_r(const vector<T>& data_in,const uint& pos){
            dyn.emplace(dyn.begin()+pos,data_in);
            check_size();
        }
        void insert_c(const vector<T>& data_in,const uint& pos){
            for(uint i = 0; i != data_in.size();i++){
                dyn[i].emplace(dyn[i].begin()+pos,data_in.at(i));
            }
            check_size();
        }
        
        void rnd(const T& lim_min,const T& lim_max){
            random_device rd;
            mt19937 gen(rd());
            uniform_real_distribution<T> dist(lim_min,lim_max);

            for(uint i = 0;i != d_r;i++){
                for(uint j = 0;j != d_c;j++){
                    dyn[i][j] = dist(gen);
                }
            }
        }
        void uniform_val(const T& val){
            for(uint i = 0;i != d_r;i++){
                for(uint j = 0;j != d_c;j++){
                    dyn[i][j] = val;
                }
            }
        }

        #ifdef TEXT_MODE
        void show_dyn(){
        for (uint i = 0; i < dyn.size(); i++) {
        for (uint j = 0; j < dyn[i].size(); j++)
            cout << dyn[i][j] << " ";
        cout << endl;}
        } 

        void show_dyn(const uint& precision){
        for (uint i = 0; i < dyn.size(); i++) {
        for (uint j = 0; j < dyn[i].size(); j++)
            cout << fixed << setprecision(precision) << dyn[i][j] << " ";
        cout << endl;}
        }
        #endif

        void clear_dyn(){this->~dyncmatx();}
        

        // OPERACIONES LINEALES

        inline T det(){
            return detn<T>(dyn);
        }
        void _inv(){
            dyn = invn<T>(dyn);
        }
        void _escaleprod(const T& val){
            dyn = multiplicarPorEscalar<T>(dyn,val);
        }
        void trn(){
            dyn = transponerMatriz(dyn);
        }
        void operator~(){
            trn();
        }
        void _escalesum(const T& val){
            dyn = sumarConstanteMatriz(dyn,val);
        }
        void _selfprodmat(dyncmatx<T> mi){
            dyn = multiplicarMatrices(get_dyncmatx(),mi.get_dyncmatx());
        }
        void _selfsummat(dyncmatx<T> mi){
            dyn = sumarMatrices(get_dyncmatx(),mi.get_dyncmatx());
        }
        void _selfdifmat(dyncmatx<T> mi){
            dyn = restarMatrices(get_dyncmatx(),mi.get_dyncmatx());
        }
        void apply_func(T f(T)){
            for(uint i = 0; i != d_r;i++){
                for(uint j = 0;j != d_c;j++){
                    dyn[i][j] = f(dyn[i][j]);
                }
            }
        }
        dyncmatx(vector<T> in_vec){dyn = vector_to_dyn(in_vec).get_dyncmatx();}
};

// FUNCIONES PARTICULARES

template <typename S>
inline dyncmatx<S> prodmat(dyncmatx<S> m1,dyncmatx<S> m2){
    return multiplicarMatrices(m1.get_dyncmatx(),m2.get_dyncmatx());
}

template <typename S>
inline dyncmatx<S> prodmatesc(const S& val,dyncmatx<S> m2){
    return multiplicarPorEscalar(m2.get_dyncmatx(),val);
}

template <typename S>
inline dyncmatx<S> summat(dyncmatx<S> m1,dyncmatx<S> m2){
    return sumarMatrices(m1.get_dyncmatx(),m2.get_dyncmatx());
}

template <typename S>
inline dyncmatx<S> summatcons(const S& val,dyncmatx<S> m2){
    return sumarConstanteMatriz(m2.get_dyncmatx(),val);
}

template <typename S>
inline dyncmatx<S> difmat(dyncmatx<S> m1,dyncmatx<S> m2){
    return restarMatrices(m1.get_dyncmatx(),m2.get_dyncmatx());
}

template <typename S>
inline dyncmatx<S> vector_to_dyn(const vector<S>& in){
    dyncmatx<S> res(in.size(),1);
    for(uint k = 0; k != in.size();k++){
        res.sd_element(k,0,in.at(k));
    }
    return res;
}

template <typename S>
inline dyncmatx<S> mergemat_c(dyncmatx<S> m1,dyncmatx<S> m2){
    uint r_merged = min(m1.r_len(),m2.r_len());
    uint c_merged = m1.c_len() + m2.c_len();
    dyncmatx<S> merged(r_merged,c_merged);
    for(uint i = 0; i != r_merged;i++){
        for(uint j = 0; j != m1.c_len();j++){
            merged.sd_element(i,j,m1.gd_element(i,j));
        }
        for(uint j = m1.c_len(); j != c_merged;j++){
            merged.sd_element(i,j,m2.gd_element(i,j - m1.c_len()));
        }
    }
    return merged;
}

// FUNCIONES EXTRAS

template <typename T>
inline T CartesianArea(dyncmatx<T> arr){
    dyncmatx<T> temp(2,2);
    T res = 0;
    for(uint u = 0;u != arr.r_len();u++){
        if(u == arr.r_len() - 1){
            temp.sd_element(0,0,arr.gd_element(u,0));
            temp.sd_element(0,1,arr.gd_element(u,1));
            temp.sd_element(1,0,arr.gd_element(0,0));
            temp.sd_element(1,1,arr.gd_element(0,1));
        }
        else{
            temp.sd_element(0,0,arr.gd_element(u,0));
            temp.sd_element(0,1,arr.gd_element(u,1));
            temp.sd_element(1,0,arr.gd_element(u+1,0));
            temp.sd_element(1,1,arr.gd_element(u+1,1));
        }
        res = res + temp.det();
    }
    return abs(res/2);
}

template <typename T>
inline dyncmatx<T> hadamard_prod(dyncmatx<T> m1,dyncmatx<T> m2){
    if(m1.c_len() == m2.c_len() && m1.r_len() == m2.r_len()){
        dyncmatx<T> res(m1.r_len(),m2.c_len());

        for(uint i = 0;i != m1.r_len();i++){
            for(uint j = 0;j != m1.c_len();j++){
                res.sd_element(i,j,m1.gd_element(i,j)*m2.gd_element(i,j));
            }
        }
        return res;
    }
    else{
        cerr << "No se puede multiplicar" << endl;
        exit(EXIT_FAILURE);
    }
}

template <typename T>
inline T euclidean_norm(dyncmatx<T> v){
    T res = 0.0;
    if(v.c_len() == 1){
        for(uint i = 0;i != v.r_len();i++){
            res = res + v.gd_element(i,0)*v.gd_element(i,0);
        }
    }
    else if(v.r_len() == 1){
        for(uint i = 0;i != v.c_len();i++){
            res = res + v.gd_element(0,i)*v.gd_element(0,i);
        }
    }
    else{
        cerr << "No hay norma" << endl;
        exit(EXIT_SUCCESS);
    }
    return pow(res,0.5);
}

template<typename T>
inline dyncmatx<T> gen_rnd(const uint& rows,const uint& col,const T& lim_min,const T& lim_max){
    dyncmatx<T> res(rows,col);
    res.rnd(lim_min,lim_max);
    return res;
}

template<typename T>
inline dyncmatx<T> gen_unim(const uint& rows,const uint& col,const T& val){
    dyncmatx<T> res(rows,col);
    res.uniform_val(val);
    return res;
}

template <typename T>
inline T var_norm(dyncmatx<T> v){
    T res = 0.0;
    if(v.c_len() == 1){
        for(uint i = 0;i != v.r_len();i++){
            res = res + v.gd_element(i,0)*v.gd_element(i,0);
        }
    }
    else if(v.r_len() == 1){
        for(uint i = 0;i != v.c_len();i++){
            res = res + v.gd_element(0,i)*v.gd_element(0,i);
        }
    }
    else{
        cerr << "No hay norma" << endl;
        exit(EXIT_SUCCESS);
    }
    return res;
}

template <typename T>
inline dyncmatx<T> func_dyn(const dyncmatx<T>& in, T f(T)){
    dyncmatx<T> res = in;
    res.apply_func(f);
    return res;
}

// EXT LINEALES

template <typename T>
inline dyncmatx<T> trans(const dyncmatx<T>& in){
    dyncmatx<T> res = in;
    res.trn();
    return res;
}

template <typename T>
inline dyncmatx<T> inv(const dyncmatx<T>& in){
    dyncmatx<T> res = in;
    res._inv();
    return res;
}

template <typename T>
inline dyncmatx<T> scrap_r(dyncmatx<T> in,const uint& index){
    dyncmatx<T> res(in.c_len(),1);
    for(uint i = 0;i != in.c_len();i++){
        res.sd_element(i,0,in.gd_element(index,i));
    }
    return res;
}

template <typename T>
inline dyncmatx<T> scrap_c(dyncmatx<T> in,const uint& index){
    dyncmatx<T> res(in.r_len(),1);
    for(uint i = 0;i != in.r_len();i++){
        res.sd_element(i,0,in.gd_element(index,i));
    }
    return res;
}

template <typename T>
inline vector<T> scrap_r_vec(dyncmatx<T> in,const uint& index){
    vector<T> res;
    res.reserve(in.c_len());
    for(uint i = 0;i != in.c_len();i++){
        res.emplace_back(in.gd_element(index,i));
    }
    return res;
}
template <typename T>
inline vector<T> scrap_c_vec(dyncmatx<T> in,const uint& index){
    vector<T> res;
    res.reserve(in.r_len());
    for(uint i = 0;i != in.r_len();i++){
        res.emplace_back(in.gd_element(i,index));
    }
    return res;
}

// OPERADORES

template <typename T>
inline dyncmatx<T> operator+ (const dyncmatx<T>& lhs ,const dyncmatx<T>& rhs){
    return summat(lhs,rhs);
}
template <typename T>
inline dyncmatx<T> operator+ (const T& lhs ,const dyncmatx<T>& rhs){
    return summatcons(lhs,rhs);
}
template <typename T>
inline dyncmatx<T> operator+ (const dyncmatx<T>& lhs ,const T& rhs){
    return summatcons(rhs,lhs);
}

template <typename T>
inline dyncmatx<T> operator- (const dyncmatx<T>& lhs ,const dyncmatx<T>& rhs){
    return difmat(lhs,rhs);
}
template <typename T>
inline dyncmatx<T> operator- (const dyncmatx<T>& lhs ,const T& rhs){
    return summatcons(lhs, T(-1 * rhs));
}
template <typename T>
inline dyncmatx<T> operator- (const T& lhs ,const dyncmatx<T>& rhs){
    return summatcons(lhs,prodmatesc(T(-1),rhs));
}

template <typename T>
inline dyncmatx<T> operator* (const dyncmatx<T>& lhs ,const dyncmatx<T>& rhs){
    return prodmat(lhs,rhs);
}
template <typename T>
inline dyncmatx<T> operator* (const T& lhs ,const dyncmatx<T>& rhs){
    return prodmatesc(lhs,rhs);
}
template <typename T>
inline dyncmatx<T> operator* (const dyncmatx<T>& lhs ,const T& rhs){
    return prodmatesc(rhs,lhs);
}

template <typename T>
inline dyncmatx<T> operator% (const dyncmatx<T>& lhs ,const dyncmatx<T>& rhs){
    return hadamard_prod(lhs,rhs);
}

// POS-DEF

template <typename T>
inline dyncmatx<T> identity(const uint& ord){
    dyncmatx<T> res(ord,ord);
    for(uint i = 0; i != ord;i++){
        res.sd_element(i,i,T(1.0));
    }
    return res;
}

#ifdef TEXT_MODE
template <typename U>
inline ostream& operator << (ostream& out,dyncmatx<U> din){
        vector<vector<U>> outer = din.get_dyncmatx();
        for (uint i = 0; i < outer.size(); i++) {
            for (U & e: outer.at(i)){out << e << " ";}
            out << endl;
        }
        return out;
}
#endif

// definiciones

typedef dyncmatx<float> FMAT;
typedef dyncmatx<int> IMAT;
typedef dyncmatx<double> DMAT;
typedef dyncmatx<double long> DLMAT;
typedef dyncmatx<int long> ILMAT;
typedef dyncmatx<unsigned> UMAT;
typedef dyncmatx<unsigned int> UIMAT;

#define I(n) identity<double>(n)
#define N(n) dyncmatx<double>(n,n)
#define NULL_VECTOR(n) dyncmatx<double>(n,1)