algebra.hpp

#ifndef ALGEBRA_HPP_
#define ALGEBRA_HPP_
// clang-format off
#include <iostream>
#include <vector>
#include <array>
#include <map>
#include <math.h>
#include <complex>
#include <string>
#include <fstream>
#include <random>

namespace algebra{

//@note Take the habit of using doxygenated comments.
//@note good idea enucleating the data for compressed matrix in a struct. It could have been made internal to the class, but is ok also as external object.
    //structure to store the data in compressed format
    template<typename T> struct CompressedData{
        std::vector<size_t> inner;
        std::vector<size_t> outer;
        std::vector<T> values;
    };

    //useful enums
    enum StorageOrder{
        row_wise,
        column_wise
    };


    enum NormType{
        one,
        infinity,
        frobenius
    };

    //declaration of compare operator (definition at bottom)
    template <StorageOrder S> struct CompareOperator;

    //declaration of operators for std::complex
    //@note The first two are provided already by the standard library, no need to redefine them.
    // As for the last two, the standard library provides the one where T1=T2, so in this case you already have the comparison with a scalar. 
    // Yours last two are an extension to the standard library, so they are fine, but maybe not really necessary, do you need to compare a complex<double> with an int?
    template<typename T> bool operator==(const std::complex<T> a, std::complex<T> b);
    template<typename T> bool operator!=(const std::complex<T> a, std::complex<T> b);
    template<typename T1, typename T2> bool operator==(const std::complex<T1> a, T2 b);
    template<typename T1, typename T2> bool operator!=(const std::complex<T1> a, T2 b);
    
    template<typename T, StorageOrder S> class Matrix{

        private:

        size_t rows=0; //number of rows
        size_t columns=0; //number of columns
        bool compressed=false; //current status
        std::map<std::array<std::size_t, 2>, T, CompareOperator<S>> data; //map with values in uncompressed format
        CompressedData<T> compr; //struct with values in compressed format

        //constructors
        public:

        Matrix(): rows(0), columns(0){} //default constructor
        Matrix(size_t r, size_t c): rows(r), columns(c){} //constructs an empty r*c matrix
        Matrix(std::vector<std::vector<T>> M){ //contructs the matrix from a non-sparse matrix defined as vector<vector<T>>
            size_t r=M.size();
            size_t c=M[0].size();
            for(size_t i=0; i<r; ++i){
                if(M[i].size() != c){
                    std::cerr << "Incompatible dimensions of rows, matrix has not been created row " << i << std::endl;
                    return;
                }
                for(size_t j=0; j<c; ++j){
                    if (M[i][j]!=0)
                        this -> data[{i,j}]=M[i][j];
                }
            }

            rows=r;
            columns=c;
            
        }
        //@note Ok, but it can be made nicer, and you have to ensure that the comparison operators are the same!!!
        // Otherwise it fails if you pass a map with a different comparison. 
        // If you want to accept giving in input a map with different comparison operator, things are getting trickier.
        // Here I show you a version with just the same comparison operator, but wich may exploit move semantic and is a bit more efficient
        // (you can make it parallel)
        /*
        template <typename MAP=std::map<std::array<std::size_t, 2>, T, CompareOperator<S>>>
        Matrix(MAP&& M): data(std::forward<MAP>(M)){ //contruct the matrix from a map
         rows=0;
         cols=0;
         std::for_each(data.begin(), data.end(), [&rows,&cols](const auto& el){
            rows=std::max(rows, el.first[0]);
            cols=std::max(cols, el.first[1]);
             });
            // If you prefer the ranges version
            // std::ranges::for_each(data, [&rows,&cols](const auto& el){
            // rows=std::max(rows, el.first[0]);
            // cols=std::max(cols, el.first[1]);
            //});
        };

        If you want a version that may take also a map with a different comparison operator, you have to do something like this:

        template <typename MAP>
        Matrix (const MAP&& M)
        {
            if constexpr (std::is_same_v<typename MAP::key_compare, CompareOperator<S>>)
            {
                data=std::forward<MAP>(M);

            }
            else
            {
                std::for_each(M.begin(), M.end(), [this](const auto& el){
                    data[el.first]=el.second;
                });
            }
            rows=0;
            cols=0;
            std::for_each(data.begin(), data.end(), [&rows,&cols](const auto& el){
            rows=std::max(rows, el.first[0]);
            cols=std::max(cols, el.first[1]);
            });
        }

        */
        Matrix(std::map<std::array<std::size_t, 2>, T> M): data(M){ //contruct the matrix from a map //@note fails if map M has different compop.
            rows=0;
            columns=0;
            for(auto it=M.begin(); it!=M.end(); ++it){
                if(it->first[0] > rows){
                    rows=it->first[0];
                }
                if(it->first[1] > columns){
                    columns=it->first[1];
                }
            }
        };

        //METHODS:
        //@note why not returning a pair or an array of 2 ints? 
        size_t get_dimension(size_t i){ //method to get dimension i of the matrix (dimension 0 = rows, dimension 1 = columns)
            if(i==0){
                return rows;
            }
            if(i==1){
                return columns;
            }
            else{
                std::cerr << "matrix has 2 dimensions, insert 0 or 1" << std::endl;
            }
        }

        bool is_compressed() const{
            return compressed;
        };

        void read(std::string filename){ //method to read a matrix in matrix market format
            std::ifstream file{filename};
            if (!file.is_open()) {
                std::cerr << "Error in opening the file" << std::endl;
                return;
            }
            data.clear();
            compr.inner.clear();
            compr.outer.clear();
            compr.values.clear();
            compressed=false;
            std::string line;
            size_t n;
            while(std::getline(file, line)){ //this portion of code skips all the text lines before the data
                if(line[0]!='%'){
                    std::istringstream iss(line);
                    iss >> rows >> columns >> n; //first line with nnumbers will contain the dimensions and the number of non-zeros
                    break;
                }
            }
            T value;
            size_t i, j;
            for (size_t k=0; k<n; ++k){
                std::getline(file, line);
                std::istringstream iss(line);
                iss >> i >> j >> value;
                data[{i-1,j-1}] = value; //indexes are 1-based, we need to subtract 1
            }
            file.close();
        }
 
        void compress(){ //fills in the compr structure, delete the uncompressed data and set the variable compressed to true
            if (is_compressed())
                return;
            //@note with a little effort you can reduce the duplication of code at a minimum!
            // for intance, the first two lines may be taken out of the if constexpr blocks.
            // you can just define two ints that take the value 0 or 1 according to the storage order and then the while
            // loop is the same for both cases.
            // Your code is fine but if you have a lot of duplication maintenance is made more difficult.
            if constexpr (S==row_wise){
                auto it= data.begin();
                size_t j=0;
                for (size_t i=0; i<rows; ++i){
                    compr.inner.push_back(j);
                    while (it!=data.end() and it -> first[0] == i){
                        compr.values.push_back(it -> second);
                        compr.outer.push_back(it -> first[1]);
                        ++j;
                        ++it;
                    }
                }
                compr.inner.push_back(j);
            }

           else{
                auto it= data.begin();
                size_t j=0;
                for (size_t i=0; i<columns; ++i){
                    compr.inner.push_back(j);
                    while (it!=data.end() and it -> first[1] == i){
                        compr.values.push_back(it -> second);
                        compr.outer.push_back(it -> first[0]);
                        ++j;
                        ++it;
                    }
                }
                compr.inner.push_back(j);
            }

            data.clear();
            compressed=true;

        };

        void uncompress(){ //fills in the uncompressed map, delete the compressed data and set the variable compressed to false
            if(!is_compressed())
                return;
            //@note also here duplication can be reduced
            if constexpr (S==row_wise){
                auto it=compr.inner.begin()+1;
                size_t i=0;
                size_t k=0;
                while(k < compr.values.size()){
                    if(k < *it){
                        data[{i, compr.outer[k]}]=compr.values[k];
                        ++k;
                    }
                    else{
                        ++i;
                        ++it;
                    }
                }
            }

            else{ //column-wise
                auto it=compr.inner.begin()+1;
                size_t i=0;
                size_t k=0;
                while(k < compr.values.size()){
                    if(k < *it){
                        data[{compr.outer[k], i}] = compr.values[k];
                        ++k;
                    }
                    else{
                        ++i;
                        ++it;
                    }
                }
            }

            compr.inner.clear();
            compr.outer.clear();
            compr.values.clear();
            //@note remember that in case of vectors clear sets size to 0 but leaves capacity untouched. If you
            // want to release memory you need to use shrink_to_fit after clearing.
            compressed=false;
        };

        void resize(const size_t r, const size_t c){ //if we resize to a smaller dimension it erases the values out of new bounds
            if (is_compressed()){
                //@note an alternative is decompress, resize and compress again. But of course it may be made also ouside of the function.
                std::cerr << "Cannot resize a compressed matrix" << std::endl;
                return;
            }
            if constexpr(S==row_wise){
                if(r<rows){
                    data.erase(data.lower_bound({r, 0}), data.end());
                }
                if(c<columns){
                    for(auto it=data.begin(); it!=data.end(); ++it){
                        if(it->first[1]>=c){
                            data.erase(it);
                        }
                    }
                }
            }
            else{
                if(c<columns){
                    data.erase(data.lower_bound({0, c}), data.end());
                }
                if(r<rows){
                    for(auto it=data.begin(); it!=data.end(); ++it){
                        if(it->first[0]>=r){
                            data.erase(it);
                        }
                    }
                }
            }
            rows=r;
            columns=c;

        };
        
        T operator()(const size_t i, const size_t j) const{ //Access operator
            if(i>=rows or j>=columns){
                std::cerr << "indexes out of bounds" << std::endl;
                return T(); //@note Better return T{}.
            }
            if(is_compressed()){
                std::cerr << "cannot access to elements of a compressed matrix" << std::endl;
                //@note You can still search for the element and  return it, if not found return 0. 
                return T();
            }
            auto it = data.find({i, j});
            if (it!= data.end()) {
                return it->second;
            } 
            else {
                return T();
            }
        }

        T& operator()(const size_t i, const size_t j){ //Access operator
        //@note I admit the the request of the challange was unclear, but this solution is not working.
        // The problem is that if you return by reference youo are allowing to chenge the object. Here you are allowing
        // the user to change a static variable. This is not what is expected and also probebly not what you wanted!
        // A safer cade would work as followe. If the element is not present in the case of compressed matrix you issue an error, with a compressed
        // matrix you are allowed only to change an existing value. Id the matrix is uncompresses you may decide that if
        // the element is not contained you add it and return the reference to the added element.
            static T default_value=0;
            if(i>=rows or j>=columns){
                std::cerr << "indexes out of bounds" << std::endl;
                default_value=0;
                return default_value;
            }
            if(is_compressed()){
                std::cerr << "cannot access to elements of a compressed matrix" << std::endl;
                default_value=0;
                return default_value;
            }
            return data[{i,j}];
        }

        std::map<size_t, T> extract(const size_t i) const{  //method to extract a row or a column (depending on the storage order)
            
            std::map<size_t, T> result;
            size_t n;
            size_t m;

            //check that the index exists
            if constexpr (S==row_wise){
                n=rows;
                m=columns;
            }
            else{
                n=columns;
                m=rows;
            }

            if (i >= n){
                std::cerr << "index out of bounds" << std::endl;
                return result;
            }

            if constexpr(S==row_wise){
                for(auto it=data.lower_bound({i,0}); it!=data.upper_bound({i,m}); ++it){
                    result[it->first[1]]=it->second;
                }
            }
            else{
                for(auto it=data.lower_bound({0,i}); it!=data.upper_bound({m,i}); ++it){
                    result[it->first[0]]=it->second;
                }
            }
            return result;
        }

        template<NormType N> double norm() const{
            //@note Nice jobs. Yet, normally for high level functions/methods one assumes that the matrix is in compressed form, whcih is 
            // the normal case for sparse matrices when you have methods that only require reading the values.
            double result=0;
            if constexpr (N==one){
                if constexpr (S==row_wise){
                    std::vector<double> sum_columns(columns, 0);
                    if(!is_compressed()){
                        for(auto it=data.begin(); it!=data.end(); ++it){
                            sum_columns[it->first[1]] += std::abs(it->second);
                        }
                    }
                    else{
                        //@note you can use std::reduce to make the code more compact and readable (and maybe parallel)
                        for(size_t i=0; i<compr.values.size(); ++i){
                            sum_columns[compr.outer[i]] += std::abs(compr.values[i]);
                        }
                    }
                    for(auto it=sum_columns.begin(); it!=sum_columns.end(); ++it){
                        if(*it > result){
                            result=*it;
                        }
                    }
                }
                else{ //column_wise
                    if(!is_compressed()){
                        for(size_t i=0; i<columns; ++i){
                            double sum_column=0;
                            for(auto it=data.lower_bound({0,i}); it!=data.upper_bound({rows,i}); ++it){
                                sum_column+=std::abs(it->second);
                            }
                            if(sum_column>result){
                                result=sum_column;
                            }
                        }
                    }
                    else{
                        for (size_t i=0; i<columns; ++i){
                            double sum_column=0;
                            for(size_t j=compr.inner[i]; j<compr.inner[i+1]; ++j){
                                sum_column += std::abs(compr.values[j]);
                            }
                            if(sum_column>result){
                                result=sum_column;
                            }
                        }   
                    }
                }
            }

            else if constexpr (N==infinity){
               if constexpr (S==row_wise){
                    if(!is_compressed()){
                        for(size_t i=0; i<rows; ++i){
                            double sum_row=0;
                            for(auto it=data.lower_bound({i,0}); it!=data.upper_bound({i,columns}); ++it){
                                sum_row+=std::abs(it->second);
                            }
                            if(sum_row>result){
                                result=sum_row;
                            }
                        }
                    }
                    else{
                        for (size_t i=0; i<rows; ++i){
                            double sum_row=0;
                            for(size_t j=compr.inner[i]; j<compr.inner[i+1]; ++j){
                                sum_row += std::abs(compr.values[j]);
                            }
                            if(sum_row>result){
                                result=sum_row;
                            }
                        }   
                    }
                }
                else{ //column_wise
                    std::vector<double> sum_rows(rows, 0);
                    if(!is_compressed()){
                        for(auto it=data.begin(); it!=data.end(); ++it){
                            sum_rows[it->first[0]] += std::abs(it->second);
                        }
                    }
                    else{
                        for(size_t i=0; i<compr.values.size(); ++i){
                            sum_rows[compr.outer[i]] += std::abs(compr.values[i]);
                        }
                    }
                    for(auto it=sum_rows.begin(); it!=sum_rows.end(); ++it){
                        if(*it > result){
                            result=*it;
                        }
                    }
                }
            }

            else{ //frobenius (both row-wise and column-wise cases)
                if(!is_compressed()){
                    for(auto it=data.begin(); it!=data.end(); ++it){
                        result += pow(abs(it->second), 2);
                    }
                    result = sqrt(result);
                }
                else{
                    for(auto it=compr.values.begin(); it!=compr.values.end(); ++it){
                        result += pow(abs(*it), 2);
                    }
                    result = sqrt(result);
                }
            }

            return result;

        };

        void print() const{ //method to print the matrix in compressed or uncompressed form depending on its current status
            if(!is_compressed()){
                std::cout << rows << " x " << columns << " matrix:" << std::endl;
                for (auto it = data.begin(); it != data.end(); ++it) {
                    std::cout << "(" << it->first[0] << ", " << it->first[1] << ") --> " << it->second << std::endl;
                }
            }
            else{
                std::cout << "inner: [";
                for (auto it = compr.inner.begin(); it != compr.inner.end(); ++it){
                    std::cout << *it << ", ";
                }
                std::cout << "]" << std::endl;

                std::cout << "outer: [";
                for (auto it = compr.outer.begin(); it != compr.outer.end(); ++it){
                    std::cout << *it << ", ";
                }
                std::cout << "]" << std::endl;

                std::cout << "values: [";
                for (auto it = compr.values.begin(); it != compr.values.end(); ++it){
                    std::cout << *it << ", ";
                }
                std::cout << "]" << std::endl;
            }
        }

        //friend functions:
        template<typename Ty, StorageOrder St> friend std::vector<Ty> operator*(const Matrix<Ty,St> M, const std::vector<Ty> v);
        template<typename Ty, StorageOrder St1, StorageOrder St2> friend std::vector<Ty> operator*(const Matrix<Ty,St1> M, const Matrix<Ty, St2> N);
    
    };

    //matrix*vector operators
    template<typename T, StorageOrder S> std::vector<T> operator*(const Matrix<T,S> M, const std::vector<T> v){ // operator * if vector is passed as a std::vector

        std::vector<T> result;

        if(M.columns != v.size()){ //checking the dimensions
            std::cerr << "Incorrect dimensions for matrix-vector multiplication" << std::endl;
            return result;
        }
//@note reserve capacity for the result vector, it is a good practice.
// std::vector<T> result(M.rows, 0);
        if constexpr(S==row_wise){
            if (!M.is_compressed()){
                for (size_t i=0; i<M.rows; ++i){
                    T sum=0;
                    for(auto it = M.data.lower_bound({i, 0}); it != M.data.upper_bound({i, M.columns}); ++it){//@note good idea to use lower_bound and upper_bound
                        sum += v[it->first[1]] * it->second;
                    }
                    result.push_back(sum);
                }
            }
            else{
                for (size_t i=0; i<M.rows; ++i){
                    T sum=0;
                    for(size_t j=M.compr.inner[i]; j<M.compr.inner[i+1]; ++j){
                        sum += M.compr.values[j] * v[M.compr.outer[j]];
                    }
                    result.push_back(sum);
                }
            }
        }
        else{
            result.resize(M.columns);
            if(!M.is_compressed()){
                for (auto it=M.data.begin(); it!=M.data.end(); ++it){
                    result[it->first[0]] += it->second * v[it->first[1]];
                }
            }
            else{
                for (size_t i=0; i<M.columns; ++i){
                    for(size_t j=M.compr.inner[i]; j<M.compr.inner[i+1]; ++j){
                        result[M.compr.outer[j]] += M.compr.values[j] * v[i];
                    }
                } 
            }
        }
        return result;
    }

    template<typename T, StorageOrder S1, StorageOrder S2> std::vector<T> operator*(const Matrix<T,S1> M, const Matrix<T, S2> N){//operator * if vector is passed as a n*1 matrix
        std::vector<T> result;
        if (N.columns!=1 or M.columns!=N.rows){
            std::cerr <<  "Incorrect dimensions for matrix-vector multiplication" << std::endl;
            return result;
        }
        std::vector<T> v(N.rows, 0);
        if (!N.is_compressed()){
            for(auto it=N.data.begin(); it!=N.data.end(); ++it){
                v[it->first[0]] = it->second;
            }
        }
        else{
            if constexpr (S2==row_wise){
                for(size_t i=0; i<N.rows; ++i){
                    if(N.compr.inner[i]!=N.compr.inner[i+1]){
                        v[i]=N.compr.values[N.compr.inner[i]];
                    }
                }
            }
            else{
                for(size_t i=0; i<N.compr.values.size(); ++i){
                    v[N.compr.outer[i]]=N.compr.values[i];
                }
            }     
        }
        result=M*v;
        return result;
    }
    
    //definition of == and != operators for std::complex
    template<typename T> bool operator==(const std::complex<T> a, std::complex<T> b){
        return (a.real()==b.real() and a.imag()==b.imag());
    }

    template<typename T> bool operator!=(const std::complex<T> a, std::complex<T> b){
        return !(a==b);
    }

    template<typename T1, typename T2> bool operator==(const std::complex<T1> a, T2 b){
        if constexpr(std::is_integral<T2>::value or std::is_floating_point<T2>::value){
            if(a.imag()==0){
                return a.real()==b;
            }
            else return false;
        }
        else return false;
    }

    template<typename T1, typename T2> bool operator!=(const std::complex<T1> a, T2 b){
        return !(a==b);
    }

    //definition of compare operator for the map
    template <StorageOrder S> struct CompareOperator{
//@note better use const references not const values. You are uselessy copying the values.
        bool operator()(const std::array<std::size_t, 2> a, const std::array<std::size_t, 2> b) const{

            if constexpr(S==row_wise){
                return a<b;
            }
            else{
                //@note if you want to exploit some nice features of C++17 you can use std::tie
                // return std::tie(a[1], a[0]) < std::tie(b[1], b[0]);
                if(a[1] != b[1]){
                    return a[1] < b[1];
                }
                else{
                    return a[0] < b[0];
                }
            }
        }
    };

    //OTHER FUNCTIONS TO BE USED IN THE MAIN
    //function to turn a vector into a n*1 matrix (default row-wise storage order)
    template <typename T, StorageOrder S> Matrix<T, S> vector_to_matrix(std::vector<T> v){
        Matrix<T,S> M(v.size(), 1);
        for(size_t i=0; i<v.size(); ++i){
            T a=v[i];
            if(a!=0){
                M(i, 0)=a;
            }
        }
        return M;
    }

    //function to define random vector of n elements with value between a and b (integers or real numbers)
    template <typename T> std::vector<T> generate_vector(const size_t n, const T a, const T b){


        std::vector<T> v;
        std::random_device rd;

        // @note This way you are everytime regenerating the random number generator, which is not necessary and costly.
        // Here a possible alternative:
        /*
        std::mt19937 gen(rd());
        if constexpr(std::is_integral<T>::value){

            std::uniform_int_distribution<T> distr(a, b);
            for(size_t i=0; i<n; ++i){
                v.push_back(distr(gen));
            }
        }
        else if constexpr(std::is_floating_point<T>::value){
            std::uniform_real_distribution<T> distr(a, b);
            for(size_t i=0; i<n; ++i){
                v.push_back(distr(gen));
            }
        }
*/

        for(size_t i=0; i<n; ++i){
            std::mt19937 gen(rd());
            if constexpr(std::is_integral<T>::value){
                std::uniform_int_distribution<T> distr(a, b);
                v.push_back(distr(gen));
            } 
            else if constexpr(std::is_floating_point<T>::value){
                std::uniform_real_distribution<T> distr(a, b);
                v.push_back(distr(gen));
            }
        }
        return v;
    }

    //variant to define the vector as a n*1 matrix
    template <typename T, StorageOrder S> Matrix<T,S> generate_column_vector(const size_t n, const T a, const T b){
        std::vector<T> v=generate_vector<T>(n,a,b);
        Matrix<T, S> M= vector_to_matrix<T,S>(v);
        return M;
    }

    //function to print a vector
    template <typename T> void print(const std::vector<T> v){
        size_t i=0;
        for(auto it=v.begin(); it!=v.end(); ++it){
            std::cout << "[" << i << "]: " << *it << std::endl;
            ++i;
        }
    }

}

#endif // ALGEBRA_HPP_