Created
April 24, 2023 12:54
-
-
Save ATMI/d70a0ea752acd43e1445287d673105f4 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// Artemii Miasoedov | |
// a.miasoedov@innopolis.university | |
#include <iostream> | |
#include <memory> | |
#include <functional> | |
#include <type_traits> | |
#include <sstream> | |
#include <iomanip> | |
#include <limits> | |
#include <cmath> | |
#include <fstream> | |
namespace la { | |
/** | |
* Exception, thrown, when invalid matrix size is encountered | |
*/ | |
class dimensions_error : public std::runtime_error { | |
public: | |
explicit dimensions_error() | |
: runtime_error("Error: the dimensional problem occurred") { | |
} | |
explicit dimensions_error(const std::string &message) | |
: runtime_error(message) { | |
} | |
}; | |
/** | |
* Steps of elimination process | |
*/ | |
enum EliminationStep { | |
ELIMINATION, | |
PERMUTATION, | |
BACK_ELIMINATION, | |
DIAG_NORM | |
}; | |
// epsilon values for float and double types, used to compare with 0 and other values | |
static thread_local float f_epsilon = std::numeric_limits<float>::epsilon(); | |
static thread_local double d_epsilon = std::numeric_limits<double>::epsilon(); | |
template<typename T> | |
requires std::is_arithmetic_v<T> | |
inline T epsilon(); | |
// returns epsilon for float type | |
template<> | |
inline float epsilon<float>() { | |
return f_epsilon; | |
} | |
// returns epsilon for double type | |
template<> | |
inline double epsilon<double>() { | |
return d_epsilon; | |
} | |
template<typename T> | |
requires std::is_arithmetic_v<T> | |
inline bool is_zero(T t); | |
// compares int value with 0 | |
template<> | |
inline bool is_zero<int>(int t) { | |
return t == 0; | |
} | |
// compares float value with 0 | |
template<> | |
inline bool is_zero<float>(float t) { | |
return std::fabs(t) < epsilon<float>(); | |
} | |
// compares double value with 0 | |
template<> | |
inline bool is_zero<double>(double t) { | |
return std::fabs(t) < epsilon<double>(); | |
} | |
template<typename T> | |
requires std::is_floating_point_v<T> | |
bool is_less_than_10_nn(T value, int n) { | |
auto t = std::abs(value) * std::pow(10, n); | |
if (t > 1) return false; | |
return t * 10 < 5; | |
// return std::abs(value) < std::pow(10, (T) -n); | |
} | |
template<typename T> requires std::is_arithmetic_v<T> | |
class Matrix; | |
template<typename T> requires std::is_arithmetic_v<T> | |
class ColumnVector; | |
/** | |
* Generic Matrix class, supporting different operations such as: multiplication, summation, and others | |
* @tparam T type of Matrix values | |
*/ | |
template<typename T> requires std::is_arithmetic_v<T> | |
class Matrix { | |
protected: | |
using dim_t = int; | |
dim_t rows_ = 0, cols_ = 0; | |
std::unique_ptr<std::unique_ptr<T[]>[]> row_vectors_; | |
public: | |
Matrix() | |
: row_vectors_(nullptr) { | |
} | |
/** | |
* Creates new matrix with specified width and height filled with 0s | |
* @param height number of rows in matrix | |
* @param width number of cols in matrix | |
*/ | |
Matrix(dim_t height, dim_t width) | |
: cols_(width), rows_(height), row_vectors_(std::make_unique<std::unique_ptr<T[]>[]>(rows_)) { | |
for (dim_t i = 0; i < rows_; ++i) { | |
std::unique_ptr<T[]> &row = (row_vectors_[i] = std::make_unique<T[]>(cols_)); | |
std::fill(&row[0], &row[cols_], 0); | |
} | |
} | |
/** | |
* Creates new square matrix filled with 0s | |
* @param size number of rows and cols in the square matrix | |
*/ | |
explicit Matrix(dim_t size) | |
: Matrix(size, size) { | |
} | |
/** | |
* Creates matrix from initializer list | |
* @param m initializer list | |
*/ | |
Matrix(const std::initializer_list<std::initializer_list<T>> &m) | |
: Matrix(m.size(), m.begin()->size()) { | |
dim_t y = 0, x = 0; | |
for (const auto &row: m) { | |
if (row.size() != cols_) { | |
throw dimensions_error("All rows must have the same number of columns"); | |
} | |
for (const auto &item: row) { | |
at(y, x) = item; | |
x++; | |
} | |
x = 0; | |
y++; | |
} | |
} | |
/** | |
* Copy constructor | |
* @param m matrix to copy | |
*/ | |
Matrix(const Matrix<T> &m) | |
: Matrix(m.rows_, m.cols_) { | |
for (dim_t i = 0; i < rows_; ++i) { | |
std::copy(&m.row_vectors_[i][0], &m.row_vectors_[i][cols_], &row_vectors_[i][0]); | |
} | |
} | |
/** | |
* Move constructor | |
* @param m matrix to move | |
*/ | |
Matrix(Matrix<T> &&m) noexcept | |
: rows_(m.rows_), cols_(m.cols_) { | |
row_vectors_ = std::move(m.row_vectors_); | |
} | |
/** | |
* Creates identity matrix with specified size | |
* @param size number of rows and cols in the matrix | |
* @return identity matrix with specified size | |
*/ | |
static Matrix<T> identity(dim_t size) { | |
Matrix<T> result(size); | |
for (int i = 0; i < size; ++i) { | |
result[i][i] = 1; | |
} | |
return result; | |
} | |
/** | |
* Creates permutation matrix for specified matrix | |
* @param m matrix, which rows/cols should be swapped | |
* @param a first row/column to swap | |
* @param b second row/column to swap | |
* @return new permutation matrix | |
*/ | |
static Matrix<T> permutation(const Matrix<T> &m, dim_t a, dim_t b) { | |
Matrix<T> result = Matrix<T>::identity(m.rows_); | |
result[a][a] = result[b][b] = 0; | |
result[a][b] = result[b][a] = 1; | |
return result; | |
} | |
/** | |
* Finds pivot in the specified column | |
* @param col where to find a pivot | |
* @return row where pivot is located and pivot itself | |
*/ | |
std::tuple<dim_t, T> pivot(dim_t col) { | |
dim_t pivot_row = col; | |
T pivot = std::abs(at(col, col)); | |
for (dim_t row = col + 1; row < rows_; ++row) { | |
// todo: | |
T e = std::abs(at(row, col)); | |
if (pivot < e) { | |
pivot = e; | |
pivot_row = row; | |
} | |
} | |
return std::make_tuple(pivot_row, at(pivot_row, col)); | |
} | |
using elimination_callback = std::function<void(EliminationStep, dim_t, dim_t, T)>; | |
/** | |
* Subtracts one row from another | |
* @param row_a subtrahend row | |
* @param row_b minuend row | |
* @param scalar factor of a subtrahend row | |
*/ | |
void eliminate(dim_t row_a, dim_t row_b, T scalar) { | |
for (int i = 0; i < cols_; ++i) { | |
at(row_a, i) -= scalar * at(row_b, i); | |
} | |
} | |
/** | |
* Swaps two rows in the matrix | |
* @param row_a first row index to swap | |
* @param row_b second row index to swap | |
*/ | |
void swap(dim_t row_a, dim_t row_b) { | |
std::swap(row_vectors_[row_a], row_vectors_[row_b]); | |
} | |
/** | |
* Divides row by the specified factor | |
* @param row row to divide | |
* @param scalar quotient | |
*/ | |
void div_row(dim_t row, T scalar) { | |
auto &r = row_vectors_[row]; | |
for (int i = 0; i < cols_; ++i) { | |
r[i] /= scalar; | |
} | |
} | |
// todo: check cols >= rows | |
/** | |
* Performs row echelon form elimination | |
* @param callback callback function for each step of REF elimination | |
*/ | |
void ref(const elimination_callback &callback = elimination_callback()) { | |
for (dim_t col = 0; col < rows_ - 1; ++col) { | |
auto [p_row, p] = pivot(col); | |
if (la::is_zero(p)) continue; | |
if (p_row != col) { | |
std::swap(row_vectors_[col], row_vectors_[p_row]); | |
callback(PERMUTATION, col, p_row, 0); | |
p_row = col; | |
} | |
for (dim_t e_row = col + 1; e_row < rows_; ++e_row) { | |
T scalar = at(e_row, col) / p; | |
at(e_row, col) = 0; | |
for (dim_t e_col = col + 1; e_col < cols_; ++e_col) { | |
at(e_row, e_col) -= (scalar * at(p_row, e_col)); | |
} | |
callback(ELIMINATION, e_row, p_row, scalar); | |
} | |
} | |
} | |
/** | |
* Performs reduced row echelon form elimination | |
* @param callback callback function for each step of RREF elimination | |
*/ | |
void rref(const elimination_callback &callback = elimination_callback()) { | |
auto max_pivots = std::min(rows_, cols_); | |
for (dim_t col = 0; col < max_pivots; ++col) { | |
auto [p_row, p] = pivot(col); | |
if (is_zero(p)) continue; | |
if (p_row != col) { | |
std::swap(row_vectors_[col], row_vectors_[p_row]); | |
callback(PERMUTATION, col, p_row, 0); | |
p_row = col; | |
} | |
for (dim_t e_row = col + 1; e_row < rows_; ++e_row) { | |
T scalar = at(e_row, col) / p; | |
if (is_zero(scalar)) continue; | |
at(e_row, col) = 0; | |
for (dim_t e_col = col; e_col < cols_; ++e_col) { | |
at(e_row, e_col) -= (scalar * at(p_row, e_col)); | |
} | |
callback(ELIMINATION, e_row, p_row, scalar); | |
} | |
} | |
callback(BACK_ELIMINATION, 0, 0, 0); | |
for (dim_t col = max_pivots - 1; col >= 0; --col) { | |
T p = at(col, col); | |
if (is_zero(p)) continue; | |
for (dim_t e_row = col - 1; e_row >= 0; --e_row) { | |
T scalar = at(e_row, col) / p; | |
if (is_zero(scalar)) continue; | |
at(e_row, col) = 0; | |
for (dim_t e_col = col; e_col < cols_; ++e_col) { | |
at(e_row, e_col) -= (scalar * at(col, e_col)); | |
} | |
callback(ELIMINATION, e_row, col, scalar); | |
} | |
} | |
for (int i = 0; i < diag(); ++i) { | |
T scalar = at(i, i); | |
if (is_zero(scalar) || is_zero(1 - scalar)) continue; | |
div_row(i, scalar); | |
callback(DIAG_NORM, i, 0, scalar); | |
} | |
} | |
/** | |
* Assignment operator | |
* @param m value to be copied | |
* @return copied value | |
*/ | |
Matrix<T> &operator =(const Matrix<T> &m) { | |
cols_ = m.cols_; | |
rows_ = m.rows_; | |
row_vectors_ = std::make_unique<std::unique_ptr<T[]>[]>(rows_); | |
for (dim_t i = 0; i < rows_; ++i) { | |
std::unique_ptr<T[]> &row = (row_vectors_[i] = std::make_unique<T[]>(cols_)); | |
std::copy(&row[0], &row[cols_], &row_vectors_[0]); | |
} | |
return *this; | |
} | |
/** | |
* Assignment operator | |
* @param m value to moved | |
* @return moved value | |
*/ | |
Matrix<T> &operator =(Matrix<T> &&m) noexcept { | |
cols_ = m.cols_; | |
rows_ = m.rows_; | |
row_vectors_ = std::move(m.row_vectors_); | |
return *this; | |
} | |
/** | |
* Returns element at specified position | |
* @param row index of value row | |
* @param col index of value column | |
* @return value at specified position | |
*/ | |
[[nodiscard]] virtual const T &at(dim_t row, dim_t col) const { | |
return row_vectors_[row][col]; | |
} | |
/** | |
* Returns element at specified position | |
* @param row index of value row | |
* @param col index of value column | |
* @return value at specified position | |
*/ | |
virtual T &at(dim_t row, dim_t col) { | |
return row_vectors_[row][col]; | |
} | |
/** | |
* Returns row with specified index | |
* @param row index of row to return | |
* @return row with specified index | |
*/ | |
const T *operator [](dim_t row) const { | |
return &row_vectors_[row][0]; | |
} | |
/** | |
* Returns row with specified index | |
* @param row index of row to return | |
* @return row with specified index | |
*/ | |
virtual T *operator [](dim_t row) { | |
return &row_vectors_[row][0]; | |
} | |
/** | |
* Sums 2 matrices | |
* @param b summand matrix | |
* @return sum of 2 matrices | |
*/ | |
virtual Matrix<T> operator +(const Matrix<T> &b) const { | |
if (b.rows_ != rows_ || b.cols_ != cols_) { | |
throw dimensions_error(); | |
} | |
Matrix<T> result(rows_, cols_); | |
for (dim_t i = 0; i < rows_; ++i) { | |
for (dim_t j = 0; j < cols_; ++j) { | |
result[i][j] = at(i, j) + b[i][j]; | |
} | |
} | |
return result; | |
} | |
/** | |
* Subtracts 2 matrices | |
* @param b subtrahend matrix | |
* @return difference of 2 matrices | |
*/ | |
virtual Matrix<T> operator -(const Matrix<T> &b) const { | |
if (b.rows_ != rows_ || b.cols_ != cols_) { | |
throw dimensions_error(); | |
} | |
Matrix<T> result(rows_, cols_); | |
for (dim_t i = 0; i < rows_; ++i) { | |
for (dim_t j = 0; j < cols_; ++j) { | |
result[i][j] = at(i, j) - b[i][j]; | |
} | |
} | |
return result; | |
} | |
/** | |
* Multiplies 2 matrices | |
* @param b multiplier matrix | |
* @return the result of matrix multiplication | |
*/ | |
virtual Matrix<T> operator *(const Matrix<T> &b) const { | |
if (b.rows_ != cols_) { | |
throw dimensions_error(); | |
} | |
Matrix<T> result(rows_, b.cols_); | |
for (dim_t i = 0; i < rows_; ++i) { | |
for (dim_t j = 0; j < b.cols_; ++j) { | |
T e_ij = 0; | |
for (dim_t k = 0; k < cols_; ++k) { | |
e_ij += at(i, k) * b[k][j]; | |
} | |
result[i][j] = e_ij; | |
} | |
} | |
return result; | |
} | |
/** | |
* Multiplies matrix by scalar | |
* @param scalar scalar to be multiplied by | |
* @return the result of multiplication of matrix by scalar | |
*/ | |
virtual Matrix<T> operator *(const T scalar) const { | |
Matrix<T> result(rows_, cols_); | |
for (dim_t i = 0; i < rows_; ++i) { | |
for (dim_t j = 0; j < cols_; ++j) { | |
result[i][j] = scalar * at(i, j); | |
} | |
} | |
return result; | |
} | |
virtual ColumnVector<T> operator *(const ColumnVector<T> &vector) const; | |
/** | |
* Transposes the matrix | |
* @return transposition result | |
*/ | |
[[nodiscard]] virtual Matrix<T> transpose() const { | |
T transposed_width = rows_; | |
T transposed_height = cols_; | |
Matrix<T> result(transposed_height, transposed_width); | |
for (dim_t i = 0; i < rows_; ++i) { | |
for (dim_t j = 0; j < cols_; ++j) { | |
result[j][i] = at(i, j); | |
} | |
} | |
return result; | |
} | |
[[nodiscard]] virtual Matrix<T> inverse() const { | |
if (cols_ != rows_) { | |
throw dimensions_error(); | |
} | |
Matrix<T> tmp = *this; | |
Matrix<T> result = Matrix<T>::identity(rows_); | |
tmp.rref( | |
[&result](la::EliminationStep step, int row_a, int row_b, double value) { | |
switch (step) { | |
case la::ELIMINATION: | |
result.eliminate(row_a, row_b, value); | |
break; | |
case la::PERMUTATION: | |
result.swap(row_a, row_b); | |
break; | |
case la::DIAG_NORM: | |
result.div_row(row_a, value); | |
break; | |
default: | |
break; | |
} | |
} | |
); | |
return result; | |
} | |
[[nodiscard]] dim_t rows() const { | |
return rows_; | |
} | |
[[nodiscard]] dim_t cols() const { | |
return cols_; | |
} | |
[[nodiscard]] dim_t diag() const { | |
return std::min(cols_, rows_); | |
} | |
[[nodiscard]] dim_t order() const { | |
return rows_ * cols_; | |
} | |
/** | |
* Reads matrix from stream | |
* @param is stream to read matrix from | |
* @param m matrix, where values will be read | |
* @return passed in stream | |
*/ | |
friend std::istream &operator >>(std::istream &is, Matrix<T> &m) { | |
for (dim_t i = 0; i < m.rows_; ++i) { | |
for (dim_t j = 0; j < m.cols_; ++j) { | |
is >> m[i][j]; | |
} | |
} | |
return is; | |
} | |
/** | |
* Prints the row from matrix | |
* @param os stream to print to | |
* @param row index of row to be printed | |
* @return passed in stream | |
*/ | |
std::ostream &print_row(std::ostream &os, dim_t row) const { | |
auto precision = os.precision(); | |
for (dim_t j = 0; j < cols_; ++j) { | |
const auto value = at(row, j); | |
os << (la::is_less_than_10_nn(value, precision) ? 0 : value); | |
if (j < cols_ - 1) { | |
os << ' '; | |
} | |
} | |
return os; | |
} | |
/** | |
* Prints the matrix to the stream | |
* @param os stream to print matrix to | |
* @param m matrix to be printed | |
* @return passed in stream | |
*/ | |
friend std::ostream &operator <<(std::ostream &os, const Matrix<T> &m) { | |
for (dim_t i = 0; i < m.rows_; ++i) { | |
m.print_row(os, i); | |
os << std::endl; | |
} | |
return os; | |
} | |
}; | |
template<> | |
std::ostream &Matrix<int>::print_row(std::ostream &os, la::Matrix<int>::dim_t row) const { | |
for (dim_t j = 0; j < cols_; ++j) { | |
os << at(row, j); | |
if (j < cols_ - 1) { | |
os << ' '; | |
} | |
} | |
return os; | |
} | |
/** | |
* Column vector | |
* @tparam T type of elements in the ColumnVector | |
*/ | |
template<typename T> requires std::is_arithmetic_v<T> | |
class ColumnVector { | |
using dim_t = unsigned int; | |
dim_t rows_ = 0; | |
std::unique_ptr<T[]> column_; | |
public: | |
/** | |
* Creates new ColumnVector with specified number of rows | |
* @param rows number of rows in ColumnVector | |
*/ | |
explicit ColumnVector(dim_t rows) | |
: rows_(rows), column_(std::make_unique<T[]>(rows)) { | |
} | |
ColumnVector(const ColumnVector<T> &m) | |
: ColumnVector(m.rows_) { | |
for (dim_t i = 0; i < rows_; ++i) { | |
this->column_[i] = m.column_[i]; | |
} | |
} | |
ColumnVector<T> &operator =(const ColumnVector<T> &m) { | |
rows_ = m.rows_; | |
column_ = std::make_unique<T[]>(rows_); | |
for (dim_t i = 0; i < rows_; ++i) { | |
this->column_[i] = m.column_[i]; | |
} | |
return *this; | |
} | |
[[nodiscard]] dim_t rows() const { | |
return rows_; | |
} | |
T &operator [](dim_t row) { | |
return column_[row]; | |
} | |
const T &operator [](dim_t row) const { | |
return column_[row]; | |
} | |
[[nodiscard]] T norm() const { | |
T r = 0; | |
for (int i = 0; i < rows_; ++i) { | |
r += (column_[i] * column_[i]); | |
} | |
return std::sqrt(r); | |
} | |
/** | |
* Sums 2 ColumnVectors | |
* @param a first summand | |
* @param b second summand | |
* @return sum of 2 vectors | |
*/ | |
friend ColumnVector<T> operator +(const ColumnVector<T> &a, const ColumnVector<T> &b) { | |
if (a.rows_ != b.rows_) { | |
throw dimensions_error("Vectors should have the same number of rows"); | |
} | |
ColumnVector<T> r(a.rows_); | |
for (int i = 0; i < a.rows_; ++i) { | |
r[i] = a[i] + b[i]; | |
} | |
return r; | |
} | |
friend ColumnVector<T> operator -(const ColumnVector<T> &a, const ColumnVector<T> &b) { | |
if (a.rows_ != b.rows_) { | |
throw dimensions_error("Vectors should have the same number of rows"); | |
} | |
ColumnVector<T> r(a.rows_); | |
for (int i = 0; i < a.rows_; ++i) { | |
r[i] = a[i] - b[i]; | |
} | |
return r; | |
} | |
/** | |
* Calculates the dot-product of 2 vectors | |
* @param a first vector | |
* @param b second vector | |
* @return dot product of 2 vectors | |
*/ | |
friend T operator *(const ColumnVector<T> &a, const ColumnVector<T> &b) { | |
if (a.rows_ != b.rows_) { | |
throw dimensions_error("Vectors should have the same number of rows"); | |
} | |
T r = 0; | |
for (int i = 0; i < a.rows_; ++i) { | |
r += a[i] * b[i]; | |
} | |
return r; | |
} | |
/** | |
* Reads ColumnVector from stream | |
* @param is stream to read from | |
* @param v ColumnVector to read to | |
* @return passed in stream | |
*/ | |
friend std::istream &operator >>(std::istream &is, ColumnVector<T> &v) { | |
for (dim_t i = 0; i < v.rows_; ++i) { | |
is >> v[i]; | |
} | |
return is; | |
} | |
/** | |
* Prints ColumnVector to stream | |
* @param os stream to print to | |
* @param v ColumnVector to print | |
* @return passed in stream | |
*/ | |
friend std::ostream &operator <<(std::ostream &os, const ColumnVector<T> &v) { | |
auto precision = os.precision(); | |
for (dim_t i = 0; i < v.rows_; ++i) { | |
const auto value = v[i]; | |
os << (la::is_less_than_10_nn(value, precision) ? 0 : value) << std::endl; | |
} | |
return os; | |
} | |
}; | |
using mati = Matrix<int>; | |
using matd = Matrix<double>; | |
using cvecd = ColumnVector<double>; | |
} | |
int main() { | |
la::d_epsilon = 10e-9; | |
std::istream &in = std::cin; | |
std::ofstream out("points.dat"); | |
out << std::fixed << std::setprecision(4); | |
int dataset_size, degree; | |
in >> dataset_size; | |
double d[dataset_size]; | |
la::cvecd b(dataset_size); | |
for (int i = 0; i < dataset_size; ++i) { | |
double x, y; | |
in >> x >> y, out << x << ' ' << y << std::endl; | |
d[i] = x, b[i] = y; | |
} | |
in >> degree; | |
la::matd a(dataset_size, degree + 1); | |
for (int i = 0; i < a.rows(); ++i) { | |
for (int j = 0; j < a.cols(); ++j) { | |
a[i][j] = std::pow(d[i], j); | |
} | |
} | |
la::matd a_t = a.transpose(); | |
b = (a_t * a).inverse() * (a_t * b); | |
FILE *gnuplot = popen("gnuplot -persistent", "w"); | |
fprintf(gnuplot, "set title 'Least squares approximation'\n"); | |
fprintf(gnuplot, "plot 'points.dat' with points title 'Data', \\\n"); | |
for (int i = 0; i < b.rows(); ++i) { | |
fprintf(gnuplot, "%f*x**%d", b[i], i); | |
if (i < b.rows() - 1) fprintf(gnuplot, "+"); | |
} | |
fprintf(gnuplot, " title 'Approximation'\n"); | |
pclose(gnuplot); | |
return 0; | |
} | |
template<typename T> | |
requires std::is_arithmetic_v<T> | |
la::ColumnVector<T> la::Matrix<T>::operator *(const ColumnVector<T> &vector) const { | |
if (cols_ != vector.rows()) { | |
throw la::dimensions_error(); | |
} | |
ColumnVector<T> result(rows_); | |
for (dim_t i = 0; i < rows_; ++i) { | |
for (dim_t j = 0; j < cols_; ++j) { | |
result[i] += vector[j] * at(i, j); | |
} | |
} | |
return result; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment