254 lines
7.0 KiB
C++
254 lines
7.0 KiB
C++
#ifndef _matrix_n_
|
||
#define _matrix_n_
|
||
|
||
#include "iostream"
|
||
#include <vector>
|
||
#include <random>
|
||
|
||
#include "./utils/vector.h"
|
||
|
||
|
||
namespace utils{
|
||
//#######################################
|
||
//# MATRIX TYPE #
|
||
//#######################################
|
||
template <typename T>
|
||
struct Matrix{
|
||
utils::Vector<T> m;
|
||
|
||
|
||
using vector_type = typename decltype(std::declval<T>().v)::value_type;
|
||
|
||
// Constructor to initialize matrix with rows × cols and a fill value
|
||
Matrix(uint64_t rows, uint64_t cols, typename T::value_type value = {}) {
|
||
fill(rows, cols, value);
|
||
}
|
||
|
||
Matrix() = default; // Default constructor
|
||
|
||
T& operator[](uint64_t idx) { return m[idx]; } // Makes it able to do matr[1][1]
|
||
const T& operator[](uint64_t idx) const { return m[idx]; } // Makes it able to do matr[1][1]
|
||
|
||
|
||
void fill(uint64_t rows, uint64_t cols, const vector_type num=0){
|
||
m.clear();
|
||
for (uint64_t i = 0; i < rows; i++){
|
||
T temp_vec;
|
||
|
||
for (uint64_t j = 0; j < cols; j++){
|
||
temp_vec.v.push_back(num);
|
||
}
|
||
m.push_back(temp_vec);
|
||
}
|
||
}
|
||
|
||
void fill_RNG(const uint64_t rows, const uint64_t cols, const vector_type min = 0, const vector_type max = 1){
|
||
m.clear();
|
||
|
||
std::mt19937_64 rng{};
|
||
rng.seed( std::random_device{}());
|
||
|
||
for (uint64_t i = 0; i < rows; i++){
|
||
T temp_vec;
|
||
for (uint64_t j = 0; j < cols; j++){
|
||
temp_vec.v.push_back(std::uniform_real_distribution<>{min, max}(rng));
|
||
}
|
||
m.push_back(temp_vec);
|
||
}
|
||
}
|
||
|
||
|
||
inline friend std::ostream& operator << (std::ostream& out, const Matrix& mat){
|
||
out << "[";
|
||
for (uint64_t i = 0; i < mat.m.size(); i++){
|
||
out << "[";
|
||
for (uint64_t j = 0; j < mat.m[i].v.size(); j++){
|
||
if (j % mat.m[i].v.size() == mat.m[i].v.size() -1 && i == mat.m.size()-1){
|
||
out << mat.m[i].v[j] << "]";
|
||
}
|
||
else if ((j % mat.m[i].v.size() == mat.m[i].v.size() -1)){
|
||
out << mat.m[i].v[j] << "]," << std::endl;
|
||
}
|
||
else{
|
||
out << mat.m[i].v[j] << ", ";
|
||
}
|
||
}
|
||
}
|
||
out << "]";
|
||
return out;
|
||
}
|
||
void print() const{
|
||
std::cout << *this << std::endl;
|
||
}
|
||
|
||
void inplace_transpose(){
|
||
utils::Vector<T> temp_m = m;
|
||
m.clear();
|
||
|
||
uint64_t rows = temp_m.size();
|
||
uint64_t cols = temp_m[0].v.size();
|
||
|
||
for (uint64_t i = 0; i < cols; i++){
|
||
T temp_vec;
|
||
for (uint64_t j = 0; j < rows; j++){
|
||
temp_vec.v.push_back(temp_m[j].v[i]);
|
||
}
|
||
m.push_back(temp_vec);
|
||
}
|
||
}
|
||
Matrix<T> transpose()const{
|
||
Matrix<T> copy = *this;
|
||
copy.inplace_transpose();
|
||
return copy;
|
||
}
|
||
|
||
|
||
void inplace_inverse(std::string method = "Gauss-Jordan"){
|
||
//Matrix<T> temp_m = *this; // Copies the m into temp_m correctly (Before: utils::Vector<T> temp_m = m;)
|
||
if (method == "Gauss-Jordan"){
|
||
Matrix<T> temp_m(m.v.size(),m[0].v.size(),0);
|
||
|
||
//std::cout << temp_m.m.v[0].size() << std::endl;
|
||
//std::cout << m.v.size() << std::endl;
|
||
|
||
uint64_t icol,irow,N=m.v.size(),M=temp_m.m.v[0].size();
|
||
double big,dum,pivinv;
|
||
Vi indxc(N,0),indxr(N,0),ipiv(N,0);
|
||
|
||
//for (uint64_t j = 0; j < N; ++j){ ipiv[j] = 0;}
|
||
for (uint64_t i = 0; i < N; i++){
|
||
big=0.0;
|
||
for (uint64_t j = 0; j < N; j++){
|
||
if (ipiv[j] != 1){
|
||
for (uint64_t k = 0; k < N; k++){
|
||
if (ipiv[k] == 0){
|
||
if (abs(m[j].v[k]) >= big){
|
||
big = abs(m[j].v[k]);
|
||
irow = j;
|
||
icol = k;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
ipiv[icol]++;
|
||
if (irow != icol){
|
||
for (uint64_t l = 0; l < N; l++){ // SWAP
|
||
double temp = m[irow].v[l];
|
||
m[irow].v[l] = m[icol].v[l];
|
||
m[icol].v[l] = temp;
|
||
}
|
||
for (uint64_t l = 0; l < M; l++){ // SWAP temp matrix
|
||
double temp = temp_m.m[irow].v[l];
|
||
temp_m.m[irow].v[l] = temp_m.m[icol].v[l];
|
||
temp_m.m[icol].v[l] = temp;
|
||
}
|
||
}
|
||
|
||
indxr[i] = irow;
|
||
indxc[i] = icol;
|
||
if (m[icol].v[icol] == 0.0){
|
||
throw std::runtime_error("utill:Matrix.Gauss-Jordan - Singular Matrix");
|
||
}
|
||
pivinv= 1.0/m[icol].v[icol];
|
||
m[icol].v[icol]=1.0;
|
||
for (uint64_t l = 0; l < N; l++){
|
||
m[icol].v[l] *= pivinv;
|
||
}
|
||
for (uint64_t l = 0; l < M; l++){
|
||
temp_m.m[icol].v[l] *= pivinv;
|
||
}
|
||
for (uint64_t ll = 0; ll < N; ll++){
|
||
if (ll != icol){
|
||
dum = m[ll].v[icol];
|
||
m[ll].v[icol] = 0;
|
||
for (uint64_t l = 0; l < N; l++){
|
||
m[ll].v[l] -= m[icol].v[l]*dum;
|
||
}
|
||
for (uint64_t l = 0; l < N; l++){
|
||
temp_m.m[ll].v[l] -= temp_m.m[icol].v[l]*dum;
|
||
}
|
||
}
|
||
}
|
||
|
||
}
|
||
//m = temp_m;
|
||
for (int64_t l = N-1; l >= 0; l--){
|
||
if (indxr[l] != indxc[l]){
|
||
for (uint64_t k = 0; k < N; k++){
|
||
double temp = m[k].v[indxr[l]];
|
||
m[k].v[indxr[l]] = m[k].v[indxc[l]];
|
||
m[k].v[indxc[l]] = temp;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
else{
|
||
throw std::runtime_error("utill:Matrix." + method + " - Not implemented yet \r \nImplemented: 'Gauss-Jordan',");
|
||
}
|
||
}
|
||
Matrix<T> inverse(std::string method = "Gauss-Jordan")const{
|
||
Matrix<T> copy = *this;
|
||
copy.inplace_inverse(method);
|
||
return copy;
|
||
}
|
||
|
||
utils::Vector<vector_type> vecmult(const utils::Vector<vector_type>& Vec)const{
|
||
|
||
if (m[0].size() != Vec.size()){
|
||
throw std::runtime_error("utill:Matrix.vecmult - Dimentions does not fit");
|
||
}
|
||
|
||
// Create a temporary result vector
|
||
utils::Vector<vector_type> copy(Vec.size(), 0);
|
||
|
||
|
||
for (uint64_t i = 0; i < m.size(); ++i) {
|
||
for (uint64_t j = 0; j < m[0].size(); ++j) {
|
||
copy[i] += m[i][j] * Vec[j];
|
||
}
|
||
}
|
||
return copy;
|
||
}
|
||
|
||
|
||
void inplace_matmult(const Matrix<T>& Mat){
|
||
|
||
if (m.v[0].size() != Mat.m.v.size()){
|
||
throw std::runtime_error("utill:Matrix.matmult - Dimentions does not fit");
|
||
}
|
||
|
||
// Dimensions of the result
|
||
uint64_t rows = m.v.size(); // rows in *this
|
||
uint64_t cols = Mat.m[0].v.size(); // columns in Mat
|
||
uint64_t inner = m.v[0].size(); // shared dimension
|
||
|
||
// Create a temporary result matrix
|
||
Matrix<T> temp_m(rows, cols, 0);
|
||
|
||
// Perform matrix multiplication
|
||
for (uint64_t i = 0; i < rows; i++){
|
||
for (uint64_t j = 0; j < cols; j++){
|
||
for (uint64_t k = 0; k < inner; k++){
|
||
temp_m.m[i].v[j] += m[i].v[k] * Mat.m[k].v[j];
|
||
}
|
||
}
|
||
}
|
||
*this = temp_m;
|
||
}
|
||
Matrix<T> matmult(const Matrix<T>& Mat)const{
|
||
Matrix<T> copy = *this;
|
||
copy.inplace_matmult(Mat);
|
||
return copy;
|
||
}
|
||
|
||
|
||
};
|
||
typedef Matrix<Vi> Mi;
|
||
typedef Matrix<Vf> Mf;
|
||
typedef Matrix<Vd> Md;
|
||
|
||
}
|
||
|
||
|
||
#endif // _numerics_n_
|