Program Listing for File Matrix.h¶
↰ Return to documentation for file (necsim/Matrix.h)
//This file is part of necsim project which is released under MIT license.
//See file **LICENSE.txt** or visit https://opensource.org/licenses/MIT) for full license details.
#ifndef MATRIX
#define MATRIX
#define null 0
#include <cstdio>
#include <iostream>
#include <sstream>
#include <fstream>
#include <cstdlib>
#include <cstring>
#include <stdexcept>
#include <vector>
#ifdef use_csv
#include<cmath>
#include <stdexcept>
#include "fast-cpp-csv-parser/csv.h"
#endif
#include <cstdint>
#include "Logging.h"
using namespace std;
namespace necsim
{
// Array of data sizes for importing tif files.
const int gdal_data_sizes[] = {0, 8, 16, 16, 32, 32, 32, 64};
template<class T>
class Matrix
{
protected:
// number of rows and columns
unsigned long num_cols{};
unsigned long num_rows{};
// a matrix is an array of rows
vector<T> matrix;
public:
explicit Matrix(unsigned long rows = 0, unsigned long cols = 0) : matrix(rows * cols, T())
{
}
Matrix(const Matrix &m) : matrix()
{
this = m;
}
virtual ~Matrix()
{
}
void setSize(unsigned long rows, unsigned long cols)
{
if(!matrix.empty())
{
matrix.clear();
}
matrix.resize(rows * cols);
num_cols = cols;
num_rows = rows;
}
unsigned long getCols() const
{
return num_cols;
}
unsigned long getRows() const
{
return num_rows;
}
void fill(T val)
{
std::fill(matrix.begin(), matrix.end(), val);
}
unsigned long index(const unsigned long &row, const unsigned long &col) const
{
#ifdef DEBUG
if(num_cols == 0 || num_rows == 0)
{
throw out_of_range("Matrix has 0 rows and columns for indexing from.");
}
if(col + num_cols * row > matrix.size())
{
stringstream ss;
ss << "Index of " << col + num_cols * row << ", (" << row << ", " << col << ")";
ss << " is out of range of matrix vector with size " << matrix.size() << endl;
throw out_of_range(ss.str());
}
#endif // DEBUG
return col + num_cols * row;
}
T &get(const unsigned long &row, const unsigned long &col)
{
#ifdef DEBUG
if(row < 0 || row >= num_rows || col < 0 || col >= num_cols)
{
stringstream ss;
ss << "Index of " << row << ", " << col << " is out of range of matrix with size " << num_rows;
ss << ", " << num_cols << endl;
throw out_of_range(ss.str());
}
#endif
return matrix[index(row, col)];
}
const T &get(const unsigned long &row, const unsigned long &col) const
{
#ifdef DEBUG
if(row < 0 || row >= num_rows || col < 0 || col >= num_cols)
{
stringstream ss;
ss << "Index of " << row << ", " << col << " is out of range of matrix with size " << num_rows;
ss << ", " << num_cols << endl;
throw out_of_range(ss.str());
}
#endif
return matrix[index(row, col)];
}
T getCopy(const unsigned long &row, const unsigned long &col) const
{
#ifdef DEBUG
if(row < 0 || row >= num_rows || col < 0 || col >= num_cols)
{
stringstream ss;
ss << "Index of " << row << ", " << col << " is out of range of matrix with size " << num_rows;
ss << ", " << num_cols << endl;
throw out_of_range(ss.str());
}
#endif
return matrix[index(row, col)];
}
typename vector<T>::iterator begin()
{
return matrix.begin();
}
typename vector<T>::iterator end()
{
return matrix.end();
}
typename vector<T>::const_iterator begin() const
{
return matrix.begin();
}
typename vector<T>::const_iterator end() const
{
return matrix.end();
}
double getMean() const
{
if(matrix.empty())
{
return 0.0;
}
T total = sum();
return double(total)/matrix.size();
}
T sum() const
{
if(matrix.empty())
{
return T();
}
else
{
T total = 0;
for(const auto &item : matrix)
{
total += item;
}
return total;
}
}
Matrix &operator=(const Matrix &m)
{
this->matrix = m.matrix;
this->num_cols = m.num_cols;
this->num_rows = m.num_rows;
return *this;
}
Matrix operator+(const Matrix &m) const
{
//Since addition creates a new matrix, we don't want to return a reference, but an actual matrix object.
unsigned long new_num_cols = findMinCols(this, m);
unsigned long new_num_rows = findMinRows(this, m);
Matrix result(new_num_rows, new_num_cols);
for(unsigned long r = 0; r < new_num_rows; r++)
{
for(unsigned long c = 0; c < new_num_cols; c++)
{
result.get(r, c) = get(r, c) + m.get(r, c);
}
}
return result;
}
Matrix operator-(const Matrix &m) const
{
unsigned long new_num_cols = findMinCols(this, m);
unsigned long new_num_rows = findMinRows(this, m);
Matrix result(new_num_rows, new_num_cols);
for(unsigned long r = 0; r < new_num_rows; r++)
{
for(unsigned long c = 0; c < new_num_cols; c++)
{
result.get(r, c) = get(r, c) - m.get(r, c);
}
}
return result;
}
Matrix &operator+=(const Matrix &m)
{
unsigned long new_num_cols = findMinCols(this, m);
unsigned long new_num_rows = findMinRows(this, m);
for(unsigned long r = 0; r < new_num_rows; r++)
{
for(unsigned long c = 0; c < new_num_cols; c++)
{
get(r, c) += m.get(r, c);
}
}
return *this;
}
Matrix &operator-=(const Matrix &m)
{
unsigned long new_num_cols = findMinCols(this, m);
unsigned long new_num_rows = findMinRows(this, m);
for(unsigned long r = 0; r < new_num_rows; r++)
{
for(unsigned long c = 0; c < new_num_cols; c++)
{
matrix.get(r, c) -= m.get(r, c);
}
}
return *this;
}
Matrix operator*(const double s) const
{
Matrix result(num_rows, num_cols);
for(unsigned long r = 0; r < num_rows; r++)
{
for(unsigned long c = 0; c < num_cols; c++)
{
result.get(r, c) = get(r, c) * s;
}
}
return result;
}
Matrix operator*(Matrix &m) const
{
unsigned long new_num_cols = findMinCols(this, m);
unsigned long new_num_rows = findMinRows(this, m);
Matrix result(num_rows, m.num_cols);
for(unsigned long r = 0; r < new_num_rows; r++)
{
for(unsigned long c = 0; c < new_num_cols; c++)
{
result.get(r, c) = get(r, c) * m.get(r, c);
}
}
return result;
}
Matrix &operator*=(const double s)
{
for(unsigned long r = 0; r < num_rows; r++)
{
for(unsigned long c = 0; c < num_cols; c++)
{
get(r, c) *= s;
}
}
return *this;
}
Matrix &operator*=(const Matrix &m)
{
unsigned long new_num_cols = findMinCols(this, m);
unsigned long new_num_rows = findMinRows(this, m);
for(unsigned long r = 0; r < new_num_rows; r++)
{
for(unsigned long c = 0; c < new_num_cols; c++)
{
get(r, c) *= m.get(r, c);
}
}
return *this;
}
Matrix operator/(const double s) const
{
Matrix result(num_rows, num_cols);
for(unsigned long r = 0; r < num_rows; r++)
{
for(unsigned long c = 0; c < num_cols; c++)
{
result.get(r, c) = get(r, c) / s;
}
}
return result;
}
Matrix &operator/=(const double s)
{
for(unsigned long r = 0; r < num_rows; r++)
{
for(unsigned long c = 0; c < num_cols; c++)
{
get(r, c) /= s;
}
}
return *this;
}
Matrix &operator/=(const Matrix &m)
{
unsigned long new_num_cols = findMinCols(this, m);
unsigned long new_num_rows = findMinRows(this, m);
for(unsigned long r = 0; r < new_num_rows; r++)
{
for(unsigned long c = 0; c < new_num_cols; c++)
{
get(r, c) /= m.get(r, c);
}
}
return *this;
}
friend ostream &writeOut(ostream &os, const Matrix &m)
{
for(unsigned long r = 0; r < m.num_rows; r++)
{
for(unsigned long c = 0; c < m.num_cols; c++)
{
os << m.getCopy(r, c) << ",";
}
os << "\n";
}
return os;
}
friend istream &readIn(istream &is, Matrix &m)
{
char delim;
for(unsigned long r = 0; r < m.num_rows; r++)
{
for(unsigned long c = 0; c < m.num_cols; c++)
{
is >> m.get(r, c);
is >> delim;
}
}
return is;
}
friend ostream &operator<<(ostream &os, const Matrix &m)
{
return writeOut(os, m);
}
friend istream &operator>>(istream &is, Matrix &m)
{
return readIn(is, m);
}
void setValue(const unsigned long &row, const unsigned long &col, const char* value)
{
matrix[index(row, col)] = static_cast<T>(*value);
}
void setValue(const unsigned long &row, const unsigned long &col, const T &value)
{
matrix[index(row, col)] = value;
}
virtual void import(const string &filename)
{
if(!importCsv(filename))
{
string s = "Type detection failed for " + filename + ". Check file_name is correct.";
throw runtime_error(s);
}
}
#ifdef use_csv
bool importCsv(const string &file_name)
{
if(file_name.find(".csv") != string::npos)
{
stringstream os;
os << "Importing " << file_name << " " << flush;
writeInfo(os.str());
// LineReader option
io::LineReader in(file_name);
// Keep track of whether we've printed to terminal or not.
bool bPrint = false;
// Initialies empty variable so that the setValue operator overloading works properly.
unsigned int number_printed = 0;
for(unsigned long i =0; i<num_rows; i++)
{
char* line = in.next_line();
if(line == nullptr)
{
if(!bPrint)
{
writeError("Input dimensions incorrect - read past end of file.");
bPrint = true;
}
break;
}
else
{
char *dToken;
dToken = strtok(line,",");
for(unsigned long j = 0; j<num_cols; j++)
{
if(dToken == nullptr)
{
if(!bPrint)
{
writeError("Input dimensions incorrect - read past end of file.");
bPrint = true;
}
break;
}
else
{
// This function is overloaded to correctly determine the type of the template
setValue(i,j,dToken);
dToken = strtok(NULL,",");
}
}
// output the percentage complete
double dComplete = ((double)i/(double)num_rows)*20;
if( number_printed < dComplete)
{
stringstream os;
os << "\rImporting " << file_name << " ";
number_printed = 0;
while(number_printed < dComplete)
{
os << ".";
number_printed ++;
}
os << flush;
writeInfo(os.str());
}
}
}
writeInfo("done.\n");
return true;
}
return false;
}
#endif
#ifndef use_csv
bool importCsv(const string &filename)
{
if(filename.find(".csv") != string::npos)
{
stringstream os;
os << "Importing" << filename << " " << flush;
ifstream inputstream;
inputstream.open(filename.c_str());
unsigned long number_printed = 0;
for(uint32_t j = 0; j < num_rows; j++)
{
string line;
getline(inputstream, line);
istringstream iss(line);
for(uint32_t i = 0; i < num_cols; i++)
{
char delim;
T val;
iss >> val >> delim;
this->setValue(j, i, val);
}
double dComplete = ((double) j / (double) num_rows) * 5;
if(number_printed < dComplete)
{
os << "\rImporting " << filename << " " << flush;
while(number_printed < dComplete)
{
os << ".";
number_printed++;
}
os << flush;
writeInfo(os.str());
}
}
stringstream os2;
os2 << "\rImporting" << filename << "..." << "done." << " " << endl;
inputstream.close();
writeInfo(os2.str());
return true;
}
return false;
}
#endif // use_csv
};
template<typename T> const unsigned long findMinCols(const Matrix<T> &matrix1, const Matrix<T> &matrix2)
{
if(matrix1.getCols() < matrix2.getCols())
{
return matrix1.getCols();
}
return matrix2.getCols();
}
template<typename T> const unsigned long findMinRows(const Matrix<T> &matrix1, const Matrix<T> &matrix2)
{
if(matrix1.getRows() < matrix2.getRows())
{
return matrix1.getRows();
}
return matrix2.getRows();
}
}
#endif // MATRIX