// This file is part of QuadProg++: a C++ library implementing // the algorithm of Goldfarb and Idnani for the solution of a (convex) // Quadratic Programming problem by means of an active-set dual method. // Copyright (C) 2007-2009 Luca Di Gaspero. // Copyright (C) 2009 Eric Moyer. // // QuadProg++ is free software: you can redistribute it and/or modify // it under the terms of the GNU Lesser General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // QuadProg++ is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Lesser General Public License for more details. // // You should have received a copy of the GNU Lesser General Public License // along with QuadProg++. If not, see . #ifndef _ARRAY_HH #define _ARRAY_HH #include #include #include #include #include #include namespace QuadProgPP{ enum MType { DIAG }; template class Vector { public: Vector(); Vector(const unsigned int n); Vector(const T& a, const unsigned int n); //initialize to constant value Vector(const T* a, const unsigned int n); // Initialize to array Vector(const Vector &rhs); // copy constructor ~Vector(); // destructor inline void set(const T* a, const unsigned int n); Vector extract(const std::set& indexes) const; inline T& operator[](const unsigned int& i); //i-th element inline const T& operator[](const unsigned int& i) const; inline unsigned int size() const; inline void resize(const unsigned int n); inline void resize(const T& a, const unsigned int n); Vector& operator=(const Vector& rhs); //assignment Vector& operator=(const T& a); //assign a to every element inline Vector& operator+=(const Vector& rhs); inline Vector& operator-=(const Vector& rhs); inline Vector& operator*=(const Vector& rhs); inline Vector& operator/=(const Vector& rhs); inline Vector& operator^=(const Vector& rhs); inline Vector& operator+=(const T& a); inline Vector& operator-=(const T& a); inline Vector& operator*=(const T& a); inline Vector& operator/=(const T& a); inline Vector& operator^=(const T& a); private: unsigned int n; // size of array. upper index is n-1 T* v; // storage for data }; template Vector::Vector() : n(0), v(0) {} template Vector::Vector(const unsigned int n) : v(new T[n]) { this->n = n; } template Vector::Vector(const T& a, const unsigned int n) : v(new T[n]) { this->n = n; for (unsigned int i = 0; i < n; i++) v[i] = a; } template Vector::Vector(const T* a, const unsigned int n) : v(new T[n]) { this->n = n; for (unsigned int i = 0; i < n; i++) v[i] = *a++; } template Vector::Vector(const Vector& rhs) : v(new T[rhs.n]) { this->n = rhs.n; for (unsigned int i = 0; i < n; i++) v[i] = rhs[i]; } template Vector::~Vector() { if (v != 0) delete[] (v); } template void Vector::resize(const unsigned int n) { if (n == this->n) return; if (v != 0) delete[] (v); v = new T[n]; this->n = n; } template void Vector::resize(const T& a, const unsigned int n) { resize(n); for (unsigned int i = 0; i < n; i++) v[i] = a; } template inline Vector& Vector::operator=(const Vector& rhs) // postcondition: normal assignment via copying has been performed; // if vector and rhs were different sizes, vector // has been resized to match the size of rhs { if (this != &rhs) { resize(rhs.n); for (unsigned int i = 0; i < n; i++) v[i] = rhs[i]; } return *this; } template inline Vector & Vector::operator=(const T& a) //assign a to every element { for (unsigned int i = 0; i < n; i++) v[i] = a; return *this; } template inline T & Vector::operator[](const unsigned int& i) //subscripting { return v[i]; } template inline const T& Vector::operator[](const unsigned int& i) const //subscripting { return v[i]; } template inline unsigned int Vector::size() const { return n; } template inline void Vector::set(const T* a, unsigned int n) { resize(n); for (unsigned int i = 0; i < n; i++) v[i] = a[i]; } template inline Vector Vector::extract(const std::set& indexes) const { Vector tmp(indexes.size()); unsigned int i = 0; for (std::set::const_iterator el = indexes.begin(); el != indexes.end(); el++) { if (*el >= n) throw std::logic_error("Error extracting subvector: the indexes are out of vector bounds"); tmp[i++] = v[*el]; } return tmp; } template inline Vector& Vector::operator+=(const Vector& rhs) { if (this->size() != rhs.size()) throw std::logic_error("Operator+=: vectors have different sizes"); for (unsigned int i = 0; i < n; i++) v[i] += rhs[i]; return *this; } template inline Vector& Vector::operator+=(const T& a) { for (unsigned int i = 0; i < n; i++) v[i] += a; return *this; } template inline Vector operator+(const Vector& rhs) { return rhs; } template inline Vector operator+(const Vector& lhs, const Vector& rhs) { if (lhs.size() != rhs.size()) throw std::logic_error("Operator+: vectors have different sizes"); Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = lhs[i] + rhs[i]; return tmp; } template inline Vector operator+(const Vector& lhs, const T& a) { Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = lhs[i] + a; return tmp; } template inline Vector operator+(const T& a, const Vector& rhs) { Vector tmp(rhs.size()); for (unsigned int i = 0; i < rhs.size(); i++) tmp[i] = a + rhs[i]; return tmp; } template inline Vector& Vector::operator-=(const Vector& rhs) { if (this->size() != rhs.size()) throw std::logic_error("Operator-=: vectors have different sizes"); for (unsigned int i = 0; i < n; i++) v[i] -= rhs[i]; return *this; } template inline Vector& Vector::operator-=(const T& a) { for (unsigned int i = 0; i < n; i++) v[i] -= a; return *this; } template inline Vector operator-(const Vector& rhs) { return (T)(-1) * rhs; } template inline Vector operator-(const Vector& lhs, const Vector& rhs) { if (lhs.size() != rhs.size()) throw std::logic_error("Operator-: vectors have different sizes"); Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = lhs[i] - rhs[i]; return tmp; } template inline Vector operator-(const Vector& lhs, const T& a) { Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = lhs[i] - a; return tmp; } template inline Vector operator-(const T& a, const Vector& rhs) { Vector tmp(rhs.size()); for (unsigned int i = 0; i < rhs.size(); i++) tmp[i] = a - rhs[i]; return tmp; } template inline Vector& Vector::operator*=(const Vector& rhs) { if (this->size() != rhs.size()) throw std::logic_error("Operator*=: vectors have different sizes"); for (unsigned int i = 0; i < n; i++) v[i] *= rhs[i]; return *this; } template inline Vector& Vector::operator*=(const T& a) { for (unsigned int i = 0; i < n; i++) v[i] *= a; return *this; } template inline Vector operator*(const Vector& lhs, const Vector& rhs) { if (lhs.size() != rhs.size()) throw std::logic_error("Operator*: vectors have different sizes"); Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = lhs[i] * rhs[i]; return tmp; } template inline Vector operator*(const Vector& lhs, const T& a) { Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = lhs[i] * a; return tmp; } template inline Vector operator*(const T& a, const Vector& rhs) { Vector tmp(rhs.size()); for (unsigned int i = 0; i < rhs.size(); i++) tmp[i] = a * rhs[i]; return tmp; } template inline Vector& Vector::operator/=(const Vector& rhs) { if (this->size() != rhs.size()) throw std::logic_error("Operator/=: vectors have different sizes"); for (unsigned int i = 0; i < n; i++) v[i] /= rhs[i]; return *this; } template inline Vector& Vector::operator/=(const T& a) { for (unsigned int i = 0; i < n; i++) v[i] /= a; return *this; } template inline Vector operator/(const Vector& lhs, const Vector& rhs) { if (lhs.size() != rhs.size()) throw std::logic_error("Operator/: vectors have different sizes"); Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = lhs[i] / rhs[i]; return tmp; } template inline Vector operator/(const Vector& lhs, const T& a) { Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = lhs[i] / a; return tmp; } template inline Vector operator/(const T& a, const Vector& rhs) { Vector tmp(rhs.size()); for (unsigned int i = 0; i < rhs.size(); i++) tmp[i] = a / rhs[i]; return tmp; } template inline Vector operator^(const Vector& lhs, const Vector& rhs) { if (lhs.size() != rhs.size()) throw std::logic_error("Operator^: vectors have different sizes"); Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = pow(lhs[i], rhs[i]); return tmp; } template inline Vector operator^(const Vector& lhs, const T& a) { Vector tmp(lhs.size()); for (unsigned int i = 0; i < lhs.size(); i++) tmp[i] = pow(lhs[i], a); return tmp; } template inline Vector operator^(const T& a, const Vector& rhs) { Vector tmp(rhs.size()); for (unsigned int i = 0; i < rhs.size(); i++) tmp[i] = pow(a, rhs[i]); return tmp; } template inline Vector& Vector::operator^=(const Vector& rhs) { if (this->size() != rhs.size()) throw std::logic_error("Operator^=: vectors have different sizes"); for (unsigned int i = 0; i < n; i++) v[i] = pow(v[i], rhs[i]); return *this; } template inline Vector& Vector::operator^=(const T& a) { for (unsigned int i = 0; i < n; i++) v[i] = pow(v[i], a); return *this; } template inline bool operator==(const Vector& v, const Vector& w) { if (v.size() != w.size()) throw std::logic_error("Vectors of different size are not confrontable"); for (unsigned i = 0; i < v.size(); i++) if (v[i] != w[i]) return false; return true; } template inline bool operator!=(const Vector& v, const Vector& w) { if (v.size() != w.size()) throw std::logic_error("Vectors of different size are not confrontable"); for (unsigned i = 0; i < v.size(); i++) if (v[i] != w[i]) return true; return false; } template inline bool operator<(const Vector& v, const Vector& w) { if (v.size() != w.size()) throw std::logic_error("Vectors of different size are not confrontable"); for (unsigned i = 0; i < v.size(); i++) if (v[i] >= w[i]) return false; return true; } template inline bool operator<=(const Vector& v, const Vector& w) { if (v.size() != w.size()) throw std::logic_error("Vectors of different size are not confrontable"); for (unsigned i = 0; i < v.size(); i++) if (v[i] > w[i]) return false; return true; } template inline bool operator>(const Vector& v, const Vector& w) { if (v.size() != w.size()) throw std::logic_error("Vectors of different size are not confrontable"); for (unsigned i = 0; i < v.size(); i++) if (v[i] <= w[i]) return false; return true; } template inline bool operator>=(const Vector& v, const Vector& w) { if (v.size() != w.size()) throw std::logic_error("Vectors of different size are not confrontable"); for (unsigned i = 0; i < v.size(); i++) if (v[i] < w[i]) return false; return true; } /** Input/Output */ template inline std::ostream& operator<<(std::ostream& os, const Vector& v) { os << std::endl << v.size() << std::endl; for (unsigned int i = 0; i < v.size() - 1; i++) os << std::setw(20) << std::setprecision(16) << v[i] << ", "; os << std::setw(20) << std::setprecision(16) << v[v.size() - 1] << std::endl; return os; } template std::istream& operator>>(std::istream& is, Vector& v) { int elements; char comma; is >> elements; v.resize(elements); for (unsigned int i = 0; i < elements; i++) is >> v[i] >> comma; return is; } /** Index utilities */ std::set seq(unsigned int s, unsigned int e); std::set singleton(unsigned int i); template class CanonicalBaseVector : public Vector { public: CanonicalBaseVector(unsigned int i, unsigned int n); inline void reset(unsigned int i); private: unsigned int e; }; template CanonicalBaseVector::CanonicalBaseVector(unsigned int i, unsigned int n) : Vector((T)0, n), e(i) { (*this)[e] = (T)1; } template inline void CanonicalBaseVector::reset(unsigned int i) { (*this)[e] = (T)0; e = i; (*this)[e] = (T)1; } #include template inline T sum(const Vector& v) { T tmp = (T)0; for (unsigned int i = 0; i < v.size(); i++) tmp += v[i]; return tmp; } template inline T prod(const Vector& v) { T tmp = (T)1; for (unsigned int i = 0; i < v.size(); i++) tmp *= v[i]; return tmp; } template inline T mean(const Vector& v) { T sum = (T)0; for (unsigned int i = 0; i < v.size(); i++) sum += v[i]; return sum / v.size(); } template inline T median(const Vector& v) { Vector tmp = sort(v); if (v.size() % 2 == 1) // it is an odd-sized vector return tmp[v.size() / 2]; else return 0.5 * (tmp[v.size() / 2 - 1] + tmp[v.size() / 2]); } template inline T stdev(const Vector& v, bool sample_correction = false) { return sqrt(var(v, sample_correction)); } template inline T var(const Vector& v, bool sample_correction = false) { T sum = (T)0, ssum = (T)0; unsigned int n = v.size(); for (unsigned int i = 0; i < n; i++) { sum += v[i]; ssum += (v[i] * v[i]); } if (!sample_correction) return (ssum / n) - (sum / n) * (sum / n); else return n * ((ssum / n) - (sum / n) * (sum / n)) / (n - 1); } template inline T max(const Vector& v) { T value = v[0]; for (unsigned int i = 1; i < v.size(); i++) value = std::max(v[i], value); return value; } template inline T min(const Vector& v) { T value = v[0]; for (unsigned int i = 1; i < v.size(); i++) value = std::min(v[i], value); return value; } template inline unsigned int index_max(const Vector& v) { unsigned int max = 0; for (unsigned int i = 1; i < v.size(); i++) if (v[i] > v[max]) max = i; return max; } template inline unsigned int index_min(const Vector& v) { unsigned int min = 0; for (unsigned int i = 1; i < v.size(); i++) if (v[i] < v[min]) min = i; return min; } template inline T dot_prod(const Vector& a, const Vector& b) { T sum = (T)0; if (a.size() != b.size()) throw std::logic_error("Dotprod error: the vectors are not the same size"); for (unsigned int i = 0; i < a.size(); i++) sum += a[i] * b[i]; return sum; } /** Single element mathematical functions */ template inline Vector exp(const Vector& v) { Vector tmp(v.size()); for (unsigned int i = 0; i < v.size(); i++) tmp[i] = exp(v[i]); return tmp; } template inline Vector log(const Vector& v) { Vector tmp(v.size()); for (unsigned int i = 0; i < v.size(); i++) tmp[i] = log(v[i]); return tmp; } template inline Vector sqrt(const Vector& v) { Vector tmp(v.size()); for (unsigned int i = 0; i < v.size(); i++) tmp[i] = sqrt(v[i]); return tmp; } template inline Vector pow(const Vector& v, double a) { Vector tmp(v.size()); for (unsigned int i = 0; i < v.size(); i++) tmp[i] = pow(v[i], a); return tmp; } template inline Vector abs(const Vector& v) { Vector tmp(v.size()); for (unsigned int i = 0; i < v.size(); i++) tmp[i] = (T)fabs(v[i]); return tmp; } template inline Vector sign(const Vector& v) { Vector tmp(v.size()); for (unsigned int i = 0; i < v.size(); i++) tmp[i] = v[i] > 0 ? +1 : v[i] == 0 ? 0 : -1; return tmp; } template inline unsigned int partition(Vector& v, unsigned int begin, unsigned int end) { unsigned int i = begin + 1, j = begin + 1; T pivot = v[begin]; while (j <= end) { if (v[j] < pivot) { std::swap(v[i], v[j]); i++; } j++; } v[begin] = v[i - 1]; v[i - 1] = pivot; return i - 2; } template inline void quicksort(Vector& v, unsigned int begin, unsigned int end) { if (end > begin) { unsigned int index = partition(v, begin, end); quicksort(v, begin, index); quicksort(v, index + 2, end); } } template inline Vector sort(const Vector& v) { Vector tmp(v); quicksort(tmp, 0, tmp.size() - 1); return tmp; } template inline Vector rank(const Vector& v) { Vector tmp(v); Vector tmp_rank(0.0, v.size()); for (unsigned int i = 0; i < tmp.size(); i++) { unsigned int smaller = 0, equal = 0; for (unsigned int j = 0; j < tmp.size(); j++) if (i == j) continue; else if (tmp[j] < tmp[i]) smaller++; else if (tmp[j] == tmp[i]) equal++; tmp_rank[i] = smaller + 1; if (equal > 0) { for (unsigned int j = 1; j <= equal; j++) tmp_rank[i] += smaller + 1 + j; tmp_rank[i] /= (double)(equal + 1); } } return tmp_rank; } //enum MType { DIAG }; template class Matrix { public: Matrix(); // Default constructor Matrix(const unsigned int n, const unsigned int m); // Construct a n x m matrix Matrix(const T& a, const unsigned int n, const unsigned int m); // Initialize the content to constant a Matrix(MType t, const T& a, const T& o, const unsigned int n, const unsigned int m); Matrix(MType t, const Vector& v, const T& o, const unsigned int n, const unsigned int m); Matrix(const T* a, const unsigned int n, const unsigned int m); // Initialize to array Matrix(const Matrix& rhs); // Copy constructor ~Matrix(); // destructor inline T* operator[](const unsigned int& i) { return v[i]; } // Subscripting: row i inline const T* operator[](const unsigned int& i) const { return v[i]; }; // const subsctipting inline void resize(const unsigned int n, const unsigned int m); inline void resize(const T& a, const unsigned int n, const unsigned int m); inline Vector extractRow(const unsigned int i) const; inline Vector extractColumn(const unsigned int j) const; inline Vector extractDiag() const; inline Matrix extractRows(const std::set& indexes) const; inline Matrix extractColumns(const std::set& indexes) const; inline Matrix extract(const std::set& r_indexes, const std::set& c_indexes) const; inline void set(const T* a, unsigned int n, unsigned int m); inline void set(const std::set& r_indexes, const std::set& c_indexes, const Matrix& m); inline void setRow(const unsigned int index, const Vector& v); inline void setRow(const unsigned int index, const Matrix& v); inline void setRows(const std::set& indexes, const Matrix& m); inline void setColumn(const unsigned int index, const Vector& v); inline void setColumn(const unsigned int index, const Matrix& v); inline void setColumns(const std::set& indexes, const Matrix& m); inline unsigned int nrows() const { return n; } // number of rows inline unsigned int ncols() const { return m; } // number of columns inline Matrix& operator=(const Matrix& rhs); // Assignment operator inline Matrix& operator=(const T& a); // Assign to every element value a inline Matrix& operator+=(const Matrix& rhs); inline Matrix& operator-=(const Matrix& rhs); inline Matrix& operator*=(const Matrix& rhs); inline Matrix& operator/=(const Matrix& rhs); inline Matrix& operator^=(const Matrix& rhs); inline Matrix& operator+=(const T& a); inline Matrix& operator-=(const T& a); inline Matrix& operator*=(const T& a); inline Matrix& operator/=(const T& a); inline Matrix& operator^=(const T& a); inline operator Vector(); private: unsigned int n; // number of rows unsigned int m; // number of columns T **v; // storage for data }; template Matrix::Matrix() : n(0), m(0), v(0) {} template Matrix::Matrix(unsigned int n, unsigned int m) : v(new T*[n]) { register unsigned int i; this->n = n; this->m = m; v[0] = new T[m * n]; for (i = 1; i < n; i++) v[i] = v[i - 1] + m; } template Matrix::Matrix(const T& a, unsigned int n, unsigned int m) : v(new T*[n]) { register unsigned int i, j; this->n = n; this->m = m; v[0] = new T[m * n]; for (i = 1; i < n; i++) v[i] = v[i - 1] + m; for (i = 0; i < n; i++) for (j = 0; j < m; j++) v[i][j] = a; } template Matrix::Matrix(const T* a, unsigned int n, unsigned int m) : v(new T*[n]) { register unsigned int i, j; this->n = n; this->m = m; v[0] = new T[m * n]; for (i = 1; i < n; i++) v[i] = v[i - 1] + m; for (i = 0; i < n; i++) for (j = 0; j < m; j++) v[i][j] = *a++; } template Matrix::Matrix(MType t, const T& a, const T& o, unsigned int n, unsigned int m) : v(new T*[n]) { register unsigned int i, j; this->n = n; this->m = m; v[0] = new T[m * n]; for (i = 1; i < n; i++) v[i] = v[i - 1] + m; switch (t) { case DIAG: for (i = 0; i < n; i++) for (j = 0; j < m; j++) if (i != j) v[i][j] = o; else v[i][j] = a; break; default: throw std::logic_error("Matrix type not supported"); } } template Matrix::Matrix(MType t, const Vector& a, const T& o, unsigned int n, unsigned int m) : v(new T*[n]) { register unsigned int i, j; this->n = n; this->m = m; v[0] = new T[m * n]; for (i = 1; i < n; i++) v[i] = v[i - 1] + m; switch (t) { case DIAG: for (i = 0; i < n; i++) for (j = 0; j < m; j++) if (i != j) v[i][j] = o; else v[i][j] = a[i]; break; default: throw std::logic_error("Matrix type not supported"); } } template Matrix::Matrix(const Matrix& rhs) : v(new T*[rhs.n]) { register unsigned int i, j; n = rhs.n; m = rhs.m; v[0] = new T[m * n]; for (i = 1; i < n; i++) v[i] = v[i - 1] + m; for (i = 0; i < n; i++) for (j = 0; j < m; j++) v[i][j] = rhs[i][j]; } template Matrix::~Matrix() { if (v != 0) { delete[] (v[0]); delete[] (v); } } template inline Matrix& Matrix::operator=(const Matrix &rhs) // postcondition: normal assignment via copying has been performed; // if matrix and rhs were different sizes, matrix // has been resized to match the size of rhs { if (this != &rhs) { register unsigned int i, j; resize(rhs.n, rhs.m); for (i = 0; i < n; i++) for (j = 0; j < m; j++) v[i][j] = rhs[i][j]; } return *this; } template inline Matrix& Matrix::operator=(const T& a) // assign a to every element { register unsigned int i, j; for (i = 0; i < n; i++) for (j = 0; j < m; j++) v[i][j] = a; return *this; } template inline void Matrix::resize(const unsigned int n, const unsigned int m) { register unsigned int i; if (n == this->n && m == this->m) return; if (v != 0) { delete[] (v[0]); delete[] (v); } this->n = n; this->m = m; v = new T*[n]; v[0] = new T[m * n]; for (i = 1; i < n; i++) v[i] = v[i - 1] + m; } template inline void Matrix::resize(const T& a, const unsigned int n, const unsigned int m) { register unsigned int i, j; resize(n, m); for (i = 0; i < n; i++) for (j = 0; j < m; j++) v[i][j] = a; } template inline Vector Matrix::extractRow(const unsigned int i) const { if (i >= n) throw std::logic_error("Error in extractRow: trying to extract a row out of matrix bounds"); Vector tmp(v[i], m); return tmp; } template inline Vector Matrix::extractColumn(const unsigned int j) const { register unsigned int i; if (j >= m) throw std::logic_error("Error in extractRow: trying to extract a row out of matrix bounds"); Vector tmp(n); for (i = 0; i < n; i++) tmp[i] = v[i][j]; return tmp; } template inline Vector Matrix::extractDiag() const { register unsigned int d = std::min(n, m), i; Vector tmp(d); for (i = 0; i < d; i++) tmp[i] = v[i][i]; return tmp; } template inline Matrix Matrix::extractRows(const std::set& indexes) const { Matrix tmp(indexes.size(), m); register unsigned int i = 0, j; for (std::set::const_iterator el = indexes.begin(); el != indexes.end(); el++) { for (j = 0; j < m; j++) { if (*el >= n) throw std::logic_error("Error extracting rows: the indexes are out of matrix bounds"); tmp[i][j] = v[*el][j]; } i++; } return tmp; } template inline Matrix Matrix::extractColumns(const std::set& indexes) const { Matrix tmp(n, indexes.size()); register unsigned int i, j = 0; for (std::set::const_iterator el = indexes.begin(); el != indexes.end(); el++) { for (i = 0; i < n; i++) { if (*el >= m) throw std::logic_error("Error extracting columns: the indexes are out of matrix bounds"); tmp[i][j] = v[i][*el]; } j++; } return tmp; } template inline Matrix Matrix::extract(const std::set& r_indexes, const std::set& c_indexes) const { Matrix tmp(r_indexes.size(), c_indexes.size()); register unsigned int i = 0, j; for (std::set::const_iterator r_el = r_indexes.begin(); r_el != r_indexes.end(); r_el++) { if (*r_el >= n) throw std::logic_error("Error extracting submatrix: the indexes are out of matrix bounds"); j = 0; for (std::set::const_iterator c_el = c_indexes.begin(); c_el != c_indexes.end(); c_el++) { if (*c_el >= m) throw std::logic_error("Error extracting rows: the indexes are out of matrix bounds"); tmp[i][j] = v[*r_el][*c_el]; j++; } i++; } return tmp; } template inline void Matrix::setRow(unsigned int i, const Vector& a) { if (i >= n) throw std::logic_error("Error in setRow: trying to set a row out of matrix bounds"); if (this->m != a.size()) throw std::logic_error("Error setting matrix row: ranges are not compatible"); for (unsigned int j = 0; j < ncols(); j++) v[i][j] = a[j]; } template inline void Matrix::setRow(unsigned int i, const Matrix& a) { if (i >= n) throw std::logic_error("Error in setRow: trying to set a row out of matrix bounds"); if (this->m != a.ncols()) throw std::logic_error("Error setting matrix column: ranges are not compatible"); if (a.nrows() != 1) throw std::logic_error("Error setting matrix column with a non-row matrix"); for (unsigned int j = 0; j < ncols(); j++) v[i][j] = a[0][j]; } template inline void Matrix::setRows(const std::set& indexes, const Matrix& m) { unsigned int i = 0; if (indexes.size() != m.nrows() || this->m != m.ncols()) throw std::logic_error("Error setting matrix rows: ranges are not compatible"); for (std::set::const_iterator el = indexes.begin(); el != indexes.end(); el++) { for (unsigned int j = 0; j < ncols(); j++) { if (*el >= n) throw std::logic_error("Error in setRows: trying to set a row out of matrix bounds"); v[*el][j] = m[i][j]; } i++; } } template inline void Matrix::setColumn(unsigned int j, const Vector& a) { if (j >= m) throw std::logic_error("Error in setColumn: trying to set a column out of matrix bounds"); if (this->n != a.size()) throw std::logic_error("Error setting matrix column: ranges are not compatible"); for (unsigned int i = 0; i < nrows(); i++) v[i][j] = a[i]; } template inline void Matrix::setColumn(unsigned int j, const Matrix& a) { if (j >= m) throw std::logic_error("Error in setColumn: trying to set a column out of matrix bounds"); if (this->n != a.nrows()) throw std::logic_error("Error setting matrix column: ranges are not compatible"); if (a.ncols() != 1) throw std::logic_error("Error setting matrix column with a non-column matrix"); for (unsigned int i = 0; i < nrows(); i++) v[i][j] = a[i][0]; } template inline void Matrix::setColumns(const std::set& indexes, const Matrix& a) { unsigned int j = 0; if (indexes.size() != a.ncols() || this->n != a.nrows()) throw std::logic_error("Error setting matrix columns: ranges are not compatible"); for (std::set::const_iterator el = indexes.begin(); el != indexes.end(); el++) { for (unsigned int i = 0; i < nrows(); i++) { if (*el >= m) throw std::logic_error("Error in setColumns: trying to set a column out of matrix bounds"); v[i][*el] = a[i][j]; } j++; } } template inline void Matrix::set(const std::set& r_indexes, const std::set& c_indexes, const Matrix& a) { unsigned int i = 0, j; if (c_indexes.size() != a.ncols() || r_indexes.size() != a.nrows()) throw std::logic_error("Error setting matrix elements: ranges are not compatible"); for (std::set::const_iterator r_el = r_indexes.begin(); r_el != r_indexes.end(); r_el++) { if (*r_el >= n) throw std::logic_error("Error in set: trying to set a row out of matrix bounds"); j = 0; for (std::set::const_iterator c_el = c_indexes.begin(); c_el != c_indexes.end(); c_el++) { if (*c_el >= m) throw std::logic_error("Error in set: trying to set a column out of matrix bounds"); v[*r_el][*c_el] = a[i][j]; j++; } i++; } } template inline void Matrix::set(const T* a, unsigned int n, unsigned int m) { if (this->n != n || this->m != m) resize(n, m); unsigned int k = 0; for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] = a[k++]; } template Matrix operator+(const Matrix& rhs) { return rhs; } template Matrix operator+(const Matrix& lhs, const Matrix& rhs) { if (lhs.ncols() != rhs.ncols() || lhs.nrows() != rhs.nrows()) throw std::logic_error("Operator+: matrices have different sizes"); Matrix tmp(lhs.nrows(), lhs.ncols()); for (unsigned int i = 0; i < lhs.nrows(); i++) for (unsigned int j = 0; j < lhs.ncols(); j++) tmp[i][j] = lhs[i][j] + rhs[i][j]; return tmp; } template Matrix operator+(const Matrix& lhs, const T& a) { Matrix tmp(lhs.nrows(), lhs.ncols()); for (unsigned int i = 0; i < lhs.nrows(); i++) for (unsigned int j = 0; j < lhs.ncols(); j++) tmp[i][j] = lhs[i][j] + a; return tmp; } template Matrix operator+(const T& a, const Matrix& rhs) { Matrix tmp(rhs.nrows(), rhs.ncols()); for (unsigned int i = 0; i < rhs.nrows(); i++) for (unsigned int j = 0; j < rhs.ncols(); j++) tmp[i][j] = a + rhs[i][j]; return tmp; } template inline Matrix& Matrix::operator+=(const Matrix& rhs) { if (m != rhs.ncols() || n != rhs.nrows()) throw std::logic_error("Operator+=: matrices have different sizes"); for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] += rhs[i][j]; return *this; } template inline Matrix& Matrix::operator+=(const T& a) { for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] += a; return *this; } template Matrix operator-(const Matrix& rhs) { return (T)(-1) * rhs; } template Matrix operator-(const Matrix& lhs, const Matrix& rhs) { if (lhs.ncols() != rhs.ncols() || lhs.nrows() != rhs.nrows()) throw std::logic_error("Operator-: matrices have different sizes"); Matrix tmp(lhs.nrows(), lhs.ncols()); for (unsigned int i = 0; i < lhs.nrows(); i++) for (unsigned int j = 0; j < lhs.ncols(); j++) tmp[i][j] = lhs[i][j] - rhs[i][j]; return tmp; } template Matrix operator-(const Matrix& lhs, const T& a) { Matrix tmp(lhs.nrows(), lhs.ncols()); for (unsigned int i = 0; i < lhs.nrows(); i++) for (unsigned int j = 0; j < lhs.ncols(); j++) tmp[i][j] = lhs[i][j] - a; return tmp; } template Matrix operator-(const T& a, const Matrix& rhs) { Matrix tmp(rhs.nrows(), rhs.ncols()); for (unsigned int i = 0; i < rhs.nrows(); i++) for (unsigned int j = 0; j < rhs.ncols(); j++) tmp[i][j] = a - rhs[i][j]; return tmp; } template inline Matrix& Matrix::operator-=(const Matrix& rhs) { if (m != rhs.ncols() || n != rhs.nrows()) throw std::logic_error("Operator-=: matrices have different sizes"); for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] -= rhs[i][j]; return *this; } template inline Matrix& Matrix::operator-=(const T& a) { for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] -= a; return *this; } template Matrix operator*(const Matrix& lhs, const Matrix& rhs) { if (lhs.ncols() != rhs.ncols() || lhs.nrows() != rhs.nrows()) throw std::logic_error("Operator*: matrices have different sizes"); Matrix tmp(lhs.nrows(), lhs.ncols()); for (unsigned int i = 0; i < lhs.nrows(); i++) for (unsigned int j = 0; j < lhs.ncols(); j++) tmp[i][j] = lhs[i][j] * rhs[i][j]; return tmp; } template Matrix operator*(const Matrix& lhs, const T& a) { Matrix tmp(lhs.nrows(), lhs.ncols()); for (unsigned int i = 0; i < lhs.nrows(); i++) for (unsigned int j = 0; j < lhs.ncols(); j++) tmp[i][j] = lhs[i][j] * a; return tmp; } template Matrix operator*(const T& a, const Matrix& rhs) { Matrix tmp(rhs.nrows(), rhs.ncols()); for (unsigned int i = 0; i < rhs.nrows(); i++) for (unsigned int j = 0; j < rhs.ncols(); j++) tmp[i][j] = a * rhs[i][j]; return tmp; } template inline Matrix& Matrix::operator*=(const Matrix& rhs) { if (m != rhs.ncols() || n != rhs.nrows()) throw std::logic_error("Operator*=: matrices have different sizes"); for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] *= rhs[i][j]; return *this; } template inline Matrix& Matrix::operator*=(const T& a) { for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] *= a; return *this; } template Matrix operator/(const Matrix& lhs, const Matrix& rhs) { if (lhs.ncols() != rhs.ncols() || lhs.nrows() != rhs.nrows()) throw std::logic_error("Operator+: matrices have different sizes"); Matrix tmp(lhs.nrows(), lhs.ncols()); for (unsigned int i = 0; i < lhs.nrows(); i++) for (unsigned int j = 0; j < lhs.ncols(); j++) tmp[i][j] = lhs[i][j] / rhs[i][j]; return tmp; } template Matrix operator/(const Matrix& lhs, const T& a) { Matrix tmp(lhs.nrows(), lhs.ncols()); for (unsigned int i = 0; i < lhs.nrows(); i++) for (unsigned int j = 0; j < lhs.ncols(); j++) tmp[i][j] = lhs[i][j] / a; return tmp; } template Matrix operator/(const T& a, const Matrix& rhs) { Matrix tmp(rhs.nrows(), rhs.ncols()); for (unsigned int i = 0; i < rhs.nrows(); i++) for (unsigned int j = 0; j < rhs.ncols(); j++) tmp[i][j] = a / rhs[i][j]; return tmp; } template inline Matrix& Matrix::operator/=(const Matrix& rhs) { if (m != rhs.ncols() || n != rhs.nrows()) throw std::logic_error("Operator+=: matrices have different sizes"); for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] /= rhs[i][j]; return *this; } template inline Matrix& Matrix::operator/=(const T& a) { for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] /= a; return *this; } template Matrix operator^(const Matrix& lhs, const T& a) { Matrix tmp(lhs.nrows(), lhs.ncols()); for (unsigned int i = 0; i < lhs.nrows(); i++) for (unsigned int j = 0; j < lhs.ncols(); j++) tmp[i][j] = pow(lhs[i][j], a); return tmp; } template inline Matrix& Matrix::operator^=(const Matrix& rhs) { if (m != rhs.ncols() || n != rhs.nrows()) throw std::logic_error("Operator^=: matrices have different sizes"); for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] = pow(v[i][j], rhs[i][j]); return *this; } template inline Matrix& Matrix::operator^=(const T& a) { for (unsigned int i = 0; i < n; i++) for (unsigned int j = 0; j < m; j++) v[i][j] = pow(v[i][j], a); return *this; } template inline Matrix::operator Vector() { if (n > 1 && m > 1) throw std::logic_error("Error matrix cast to vector: trying to cast a multi-dimensional matrix"); if (n == 1) return extractRow(0); else return extractColumn(0); } template inline bool operator==(const Matrix& a, const Matrix& b) { if (a.nrows() != b.nrows() || a.ncols() != b.ncols()) throw std::logic_error("Matrices of different size are not confrontable"); for (unsigned i = 0; i < a.nrows(); i++) for (unsigned j = 0; j < a.ncols(); j++) if (a[i][j] != b[i][j]) return false; return true; } template inline bool operator!=(const Matrix& a, const Matrix& b) { if (a.nrows() != b.nrows() || a.ncols() != b.ncols()) throw std::logic_error("Matrices of different size are not confrontable"); for (unsigned i = 0; i < a.nrows(); i++) for (unsigned j = 0; j < a.ncols(); j++) if (a[i][j] != b[i][j]) return true; return false; } /** Input/Output */ template std::ostream& operator<<(std::ostream& os, const Matrix& m) { os << std::endl << m.nrows() << " " << m.ncols() << std::endl; for (unsigned int i = 0; i < m.nrows(); i++) { for (unsigned int j = 0; j < m.ncols() - 1; j++) os << std::setw(20) << std::setprecision(16) << m[i][j] << ", "; os << std::setw(20) << std::setprecision(16) << m[i][m.ncols() - 1] << std::endl; } return os; } template std::istream& operator>>(std::istream& is, Matrix& m) { int rows, cols; char comma; is >> rows >> cols; m.resize(rows, cols); for (unsigned int i = 0; i < rows; i++) for (unsigned int j = 0; j < cols; j++) is >> m[i][j] >> comma; return is; } template T sign(const T& v) { if (v >= (T)0.0) return (T)1.0; else return (T)-1.0; } template T dist(const T& a, const T& b) { T abs_a = (T)fabs(a), abs_b = (T)fabs(b); if (abs_a > abs_b) return abs_a * sqrt((T)1.0 + (abs_b / abs_a) * (abs_b / abs_a)); else return (abs_b == (T)0.0 ? (T)0.0 : abs_b * sqrt((T)1.0 + (abs_a / abs_b) * (abs_a / abs_b))); } template void svd(const Matrix& A, Matrix& U, Vector& W, Matrix& V) { int m = A.nrows(), n = A.ncols(), i, j, k, l, jj, nm; const unsigned int max_its = 30; bool flag; Vector rv1(n); U = A; W.resize(n); V.resize(n, n); T anorm, c, f, g, h, s, scale, x, y, z; g = scale = anorm = (T)0.0; // Householder reduction to bidiagonal form for (i = 0; i < n; i++) { l = i + 1; rv1[i] = scale * g; g = s = scale = (T)0.0; if (i < m) { for (k = i; k < m; k++) scale += fabs(U[k][i]); if (scale != (T)0.0) { for (k = i; k < m; k++) { U[k][i] /= scale; s += U[k][i] * U[k][i]; } f = U[i][i]; g = -sign(f) * sqrt(s); h = f * g - s; U[i][i] = f - g; for (j = l; j < n; j++) { s = (T)0.0; for (k = i; k < m; k++) s += U[k][i] * U[k][j]; f = s / h; for (k = i; k < m; k++) U[k][j] += f * U[k][i]; } for (k = i; k < m; k++) U[k][i] *= scale; } } W[i] = scale * g; g = s = scale = (T)0.0; if (i < m && i != n - 1) { for (k = l; k < n; k++) scale += fabs(U[i][k]); if (scale != (T)0.0) { for (k = l; k < n; k++) { U[i][k] /= scale; s += U[i][k] * U[i][k]; } f = U[i][l]; g = -sign(f) * sqrt(s); h = f * g - s; U[i][l] = f - g; for (k = l; k = 0; i--) { if (i < n - 1) { if (g != (T)0.0) { for (j = l; j < n; j++) V[j][i] = (U[i][j] / U[i][l]) / g; for (j = l; j < n; j++) { s = (T)0.0; for (k = l; k < n; k++) s += U[i][k] * V[k][j]; for (k = l; k < n; k++) V[k][j] += s * V[k][i]; } } for (j = l; j < n; j++) V[i][j] = V[j][i] = (T)0.0; } V[i][i] = (T)1.0; g = rv1[i]; l = i; } // Accumulation of left-hand transformations for (i = std::min(m, n) - 1; i >= 0; i--) { l = i + 1; g = W[i]; for (j = l; j < n; j++) U[i][j] = (T)0.0; if (g != (T)0.0) { g = (T)1.0 / g; for (j = l; j < n; j++) { s = (T)0.0; for (k = l; k < m; k++) s += U[k][i] * U[k][j]; f = (s / U[i][i]) * g; for (k = i; k < m; k++) U[k][j] += f * U[k][i]; } for (j = i; j < m; j++) U[j][i] *= g; } else for (j = i; j < m; j++) U[j][i] = (T)0.0; U[i][i]++; } // Diagonalization of the bidiagonal form: loop over singular values, and over allowed iterations. for (k = n - 1; k >= 0; k--) { for (unsigned int its = 0; its < max_its; its++) { flag = true; for (l = k; l >= 0; l--) // FIXME: in NR it was l >= 1 but there subscripts start from one { // Test for splitting nm = l - 1; // Note that rV[0] is always zero if ((T)(fabs(rv1[l]) + anorm) == anorm) { flag = false; break; } if ((T)(fabs(W[nm]) + anorm) == anorm) break; } if (flag) { // Cancellation of rv1[l], if l > 0 FIXME: it was l > 1 in NR c = (T)0.0; s = (T)1.0; for (i = l; i <= k; i++) { f = s * rv1[i]; rv1[i] *= c; if ((T)(fabs(f) + anorm) == anorm) break; g = W[i]; h = dist(f, g); W[i] = h; h = (T)1.0 / h; c = g * h; s = -f * h; for (j = 0; j < m; j++) { y = U[j][nm]; z = U[j][i]; U[j][nm] = y * c + z * s; U[j][i] = z * c - y * s; } } } z = W[k]; if (l == k) { // Convergence if (z < (T)0.0) { // Singular value is made nonnegative W[k] = -z; for (j = 0; j < n; j++) V[j][k] = -V[j][k]; } break; } if (its == max_its) throw std::logic_error("Error svd: no convergence in the maximum number of iterations"); x = W[l]; nm = k - 1; y = W[nm]; g = rv1[nm]; h = rv1[k]; f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y); g = dist(f, (T)1.0); f = ((x - z) * (x + z) + h * ((y / (f + sign(f)*fabs(g))) - h)) / x; c = s = (T)1.0; // Next QR transformation for (j = l; j <= nm; j++) { i = j + 1; g = rv1[i]; y = W[i]; h = s * g; g *= c; z = dist(f, h); rv1[j] = z; c = f / z; s = h / z; f = x * c + g * s; g = g * c - x * s; h = y * s; y *= c; for (jj = 0; jj < n; jj++) { x = V[jj][j]; z = V[jj][i]; V[jj][j] = x * c + z * s; V[jj][i] = z * c - x * s; } z = dist(f, h); W[j] = z; if (z != 0) // Rotation can be arbitrary if z = 0 { z = (T)1.0 / z; c = f * z; s = h * z; } f = c * g + s * y; x = c * y - s * g; for (jj = 0; jj < m; jj++) { y = U[jj][j]; z = U[jj][i]; U[jj][j] = y * c + z * s; U[jj][i] = z * c - y * s; } } rv1[l] = (T)0.0; rv1[k] = f; W[k] = x; } } } template Matrix pinv(const Matrix& A) { Matrix U, V, x, tmp(A.ncols(), A.nrows()); Vector W; CanonicalBaseVector e(0, A.nrows()); svd(A, U, W, V); for (unsigned int i = 0; i < A.nrows(); i++) { e.reset(i); tmp.setColumn(i, dot_prod(dot_prod(dot_prod(V, Matrix(DIAG, 1.0 / W, 0.0, W.size(), W.size())), t(U)), e)); } return tmp; } template int lu(const Matrix& A, Matrix& LU, Vector& index) { if (A.ncols() != A.nrows()) throw std::logic_error("Error in LU decomposition: matrix must be squared"); int i, p, j, k, n = A.ncols(), ex; T val, tmp; Vector d(n); LU = A; index.resize(n); ex = 1; for (i = 0; i < n; i++) { index[i] = i; val = (T)0.0; for (j = 0; j < n; j++) val = std::max(val, (T)fabs(LU[i][j])); if (val == (T)0.0) std::logic_error("Error in LU decomposition: matrix was singular"); d[i] = val; } for (k = 0; k < n - 1; k++) { p = k; val = fabs(LU[k][k]) / d[k]; for (i = k + 1; i < n; i++) { tmp = fabs(LU[i][k]) / d[i]; if (tmp > val) { val = tmp; p = i; } } if (val == (T)0.0) std::logic_error("Error in LU decomposition: matrix was singular"); if (p > k) { ex = -ex; std::swap(index[k], index[p]); std::swap(d[k], d[p]); for (j = 0; j < n; j++) std::swap(LU[k][j], LU[p][j]); } for (i = k + 1; i < n; i++) { LU[i][k] /= LU[k][k]; for (j = k + 1; j < n; j++) LU[i][j] -= LU[i][k] * LU[k][j]; } } if (LU[n - 1][n - 1] == (T)0.0) std::logic_error("Error in LU decomposition: matrix was singular"); return ex; } template Vector lu_solve(const Matrix& LU, const Vector& b, Vector& index) { if (LU.ncols() != LU.nrows()) throw std::logic_error("Error in LU solve: LU matrix should be squared"); unsigned int n = LU.ncols(); if (b.size() != n) throw std::logic_error("Error in LU solve: b vector must be of the same dimensions of LU matrix"); Vector x((T)0.0, n); int i, j, p; T sum; p = index[0]; x[0] = b[p]; for (i = 1; i < n; i++) { sum = (T)0.0; for (j = 0; j < i; j++) sum += LU[i][j] * x[j]; p = index[i]; x[i] = b[p] - sum; } x[n - 1] /= LU[n - 1][n - 1]; for (i = n - 2; i >= 0; i--) { sum = (T)0.0; for (j = i + 1; j < n; j++) sum += LU[i][j] * x[j]; x[i] = (x[i] - sum) / LU[i][i]; } return x; } template void lu_solve(const Matrix& LU, Vector& x, const Vector& b, Vector& index) { x = lu_solve(LU, b, index); } template Matrix lu_inverse(const Matrix& A) { if (A.ncols() != A.nrows()) throw std::logic_error("Error in LU invert: matrix must be squared"); unsigned int n = A.ncols(); Matrix A1(n, n), LU; Vector index; lu(A, LU, index); CanonicalBaseVector e(0, n); for (unsigned i = 0; i < n; i++) { e.reset(i); A1.setColumn(i, lu_solve(LU, e, index)); } return A1; } template T lu_det(const Matrix& A) { if (A.ncols() != A.nrows()) throw std::logic_error("Error in LU determinant: matrix must be squared"); unsigned int d; Matrix LU; Vector index; d = lu(A, LU, index); return d * prod(LU.extractDiag()); } template void cholesky(const Matrix A, Matrix& LL) { if (A.ncols() != A.nrows()) throw std::logic_error("Error in Cholesky decomposition: matrix must be squared"); register int i, j, k, n = A.ncols(); register double sum; LL = A; for (i = 0; i < n; i++) { for (j = i; j < n; j++) { sum = LL[i][j]; for (k = i - 1; k >= 0; k--) sum -= LL[i][k] * LL[j][k]; if (i == j) { if (sum <= 0.0) throw std::logic_error("Error in Cholesky decomposition: matrix is not postive definite"); LL[i][i] = ::std::sqrt(sum); } else LL[j][i] = sum / LL[i][i]; } for (k = i + 1; k < n; k++) LL[i][k] = LL[k][i]; } } template Matrix cholesky(const Matrix A) { Matrix LL; cholesky(A, LL); return LL; } template Vector cholesky_solve(const Matrix& LL, const Vector& b) { if (LL.ncols() != LL.nrows()) throw std::logic_error("Error in Cholesky solve: matrix must be squared"); unsigned int n = LL.ncols(); if (b.size() != n) throw std::logic_error("Error in Cholesky decomposition: b vector must be of the same dimensions of LU matrix"); Vector x, y; /* Solve L * y = b */ forward_elimination(LL, y, b); /* Solve L^T * x = y */ backward_elimination(LL, x, y); return x; } template void cholesky_solve(const Matrix& LL, Vector& x, const Vector& b) { x = cholesky_solve(LL, b); } template void forward_elimination(const Matrix& L, Vector& y, const Vector b) { if (L.ncols() != L.nrows()) throw std::logic_error("Error in Forward elimination: matrix must be squared (lower triangular)"); if (b.size() != L.nrows()) throw std::logic_error("Error in Forward elimination: b vector must be of the same dimensions of L matrix"); register int i, j, n = b.size(); y.resize(n); y[0] = b[0] / L[0][0]; for (i = 1; i < n; i++) { y[i] = b[i]; for (j = 0; j < i; j++) y[i] -= L[i][j] * y[j]; y[i] = y[i] / L[i][i]; } } template Vector forward_elimination(const Matrix& L, const Vector b) { Vector y; forward_elimination(L, y, b); return y; } template void backward_elimination(const Matrix& U, Vector& x, const Vector& y) { if (U.ncols() != U.nrows()) throw std::logic_error("Error in Backward elimination: matrix must be squared (upper triangular)"); if (y.size() != U.nrows()) throw std::logic_error("Error in Backward elimination: b vector must be of the same dimensions of U matrix"); register int i, j, n = y.size(); x.resize(n); x[n - 1] = y[n - 1] / U[n - 1][n - 1]; for (i = n - 2; i >= 0; i--) { x[i] = y[i]; for (j = i + 1; j < n; j++) x[i] -= U[i][j] * x[j]; x[i] = x[i] / U[i][i]; } } template Vector backward_elimination(const Matrix& U, const Vector y) { Vector x; forward_elimination(U, x, y); return x; } /* Setting default linear systems machinery */ #define det lu_det #define inverse lu_inverse #define solve lu_solve /* Random */ template void random(Matrix& m) { for (unsigned int i = 0; i < m.nrows(); i++) for (unsigned int j = 0; j < m.ncols(); j++) m[i][j] = (T)(rand() / double(RAND_MAX)); } /** Aggregate functions */ template Vector sum(const Matrix& m) { Vector tmp((T)0, m.ncols()); for (unsigned int j = 0; j < m.ncols(); j++) for (unsigned int i = 0; i < m.nrows(); i++) tmp[j] += m[i][j]; return tmp; } template Vector r_sum(const Matrix& m) { Vector tmp((T)0, m.nrows()); for (unsigned int i = 0; i < m.nrows(); i++) for (unsigned int j = 0; j < m.ncols(); j++) tmp[i] += m[i][j]; return tmp; } template T all_sum(const Matrix& m) { T tmp = (T)0; for (unsigned int i = 0; i < m.nrows(); i++) for (unsigned int j = 0; j < m.ncols(); j++) tmp += m[i][j]; return tmp; } template Vector prod(const Matrix& m) { Vector tmp((T)1, m.ncols()); for (unsigned int j = 0; j < m.ncols(); j++) for (unsigned int i = 0; i < m.nrows(); i++) tmp[j] *= m[i][j]; return tmp; } template Vector r_prod(const Matrix& m) { Vector tmp((T)1, m.nrows()); for (unsigned int i = 0; i < m.nrows(); i++) for (unsigned int j = 0; j < m.ncols(); j++) tmp[i] *= m[i][j]; return tmp; } template T all_prod(const Matrix& m) { T tmp = (T)1; for (unsigned int i = 0; i < m.nrows(); i++) for (unsigned int j = 0; j < m.ncols(); j++) tmp *= m[i][j]; return tmp; } template Vector mean(const Matrix& m) { Vector res((T)0, m.ncols()); for (unsigned int j = 0; j < m.ncols(); j++) { for (unsigned int i = 0; i < m.nrows(); i++) res[j] += m[i][j]; res[j] /= m.nrows(); } return res; } template Vector r_mean(const Matrix& m) { Vector res((T)0, m.rows()); for (unsigned int i = 0; i < m.nrows(); i++) { for (unsigned int j = 0; j < m.ncols(); j++) res[i] += m[i][j]; res[i] /= m.nrows(); } return res; } template T all_mean(const Matrix& m) { T tmp = (T)0; for (unsigned int i = 0; i < m.nrows(); i++) for (unsigned int j = 0; j < m.ncols(); j++) tmp += m[i][j]; return tmp / (m.nrows() * m.ncols()); } template Vector var(const Matrix& m, bool sample_correction = false) { Vector res((T)0, m.ncols()); unsigned int n = m.nrows(); double sum, ssum; for (unsigned int j = 0; j < m.ncols(); j++) { sum = (T)0.0; ssum = (T)0.0; for (unsigned int i = 0; i < m.nrows(); i++) { sum += m[i][j]; ssum += (m[i][j] * m[i][j]); } if (!sample_correction) res[j] = (ssum / n) - (sum / n) * (sum / n); else res[j] = n * ((ssum / n) - (sum / n) * (sum / n)) / (n - 1); } return res; } template Vector stdev(const Matrix& m, bool sample_correction = false) { return sqrt(var(m, sample_correction)); } template Vector r_var(const Matrix& m, bool sample_correction = false) { Vector res((T)0, m.nrows()); double sum, ssum; unsigned int n = m.ncols(); for (unsigned int i = 0; i < m.nrows(); i++) { sum = 0.0; ssum = 0.0; for (unsigned int j = 0; j < m.ncols(); j++) { sum += m[i][j]; ssum += (m[i][j] * m[i][j]); } if (!sample_correction) res[i] = (ssum / n) - (sum / n) * (sum / n); else res[i] = n * ((ssum / n) - (sum / n) * (sum / n)) / (n - 1); } return res; } template Vector r_stdev(const Matrix& m, bool sample_correction = false) { return sqrt(r_var(m, sample_correction)); } template Vector max(const Matrix& m) { Vector res(m.ncols()); double value; for (unsigned int j = 0; j < m.ncols(); j++) { value = m[0][j]; for (unsigned int i = 1; i < m.nrows(); i++) value = std::max(m[i][j], value); res[j] = value; } return res; } template Vector r_max(const Matrix& m) { Vector res(m.nrows()); double value; for (unsigned int i = 0; i < m.nrows(); i++) { value = m[i][0]; for (unsigned int j = 1; j < m.ncols(); j++) value = std::max(m[i][j], value); res[i] = value; } return res; } template Vector min(const Matrix& m) { Vector res(m.ncols()); double value; for (unsigned int j = 0; j < m.ncols(); j++) { value = m[0][j]; for (unsigned int i = 1; i < m.nrows(); i++) value = std::min(m[i][j], value); res[j] = value; } return res; } template Vector r_min(const Matrix& m) { Vector res(m.nrows()); double value; for (unsigned int i = 0; i < m.nrows(); i++) { value = m[i][0]; for (unsigned int j = 1; j < m.ncols(); j++) value = std::min(m[i][j], value); res[i] = value; } return res; } /** Single element mathematical functions */ template Matrix exp(const Matrix&m) { Matrix tmp(m.nrows(), m.ncols()); for (unsigned int i = 0; i < m.nrows(); i++) for (unsigned int j = 0; j < m.ncols(); j++) tmp[i][j] = exp(m[i][j]); return tmp; } template Matrix sqrt(const Matrix&m) { Matrix tmp(m.nrows(), m.ncols()); for (unsigned int i = 0; i < m.nrows(); i++) for (unsigned int j = 0; j < m.ncols(); j++) tmp[i][j] = sqrt(m[i][j]); return tmp; } /** Matrix operators */ template Matrix kron(const Vector& b, const Vector& a) { Matrix tmp(b.size(), a.size()); for (unsigned int i = 0; i < b.size(); i++) for (unsigned int j = 0; j < a.size(); j++) tmp[i][j] = a[j] * b[i]; return tmp; } template Matrix t(const Matrix& a) { Matrix tmp(a.ncols(), a.nrows()); for (unsigned int i = 0; i < a.nrows(); i++) for (unsigned int j = 0; j < a.ncols(); j++) tmp[j][i] = a[i][j]; return tmp; } template Matrix dot_prod(const Matrix& a, const Matrix& b) { if (a.ncols() != b.nrows()) throw std::logic_error("Error matrix dot product: dimensions of the matrices are not compatible"); Matrix tmp(a.nrows(), b.ncols()); for (unsigned int i = 0; i < tmp.nrows(); i++) for (unsigned int j = 0; j < tmp.ncols(); j++) { tmp[i][j] = (T)0; for (unsigned int k = 0; k < a.ncols(); k++) tmp[i][j] += a[i][k] * b[k][j]; } return tmp; } template Matrix dot_prod(const Matrix& a, const Vector& b) { if (a.ncols() != b.size()) throw std::logic_error("Error matrix dot product: dimensions of the matrix and the vector are not compatible"); Matrix tmp(a.nrows(), 1); for (unsigned int i = 0; i < tmp.nrows(); i++) { tmp[i][0] = (T)0; for (unsigned int k = 0; k < a.ncols(); k++) tmp[i][0] += a[i][k] * b[k]; } return tmp; } template Matrix dot_prod(const Vector& a, const Matrix& b) { if (a.size() != b.ncols()) throw std::logic_error("Error matrix dot product: dimensions of the vector and matrix are not compatible"); Matrix tmp(1, b.ncols()); for (unsigned int j = 0; j < tmp.ncols(); j++) { tmp[0][j] = (T)0; for (unsigned int k = 0; k < a.size(); k++) tmp[0][j] += a[k] * b[k][j]; } return tmp; } template inline Matrix rank(const Matrix m) { Matrix tmp(m.nrows(), m.ncols()); for (unsigned int j = 0; j < m.ncols(); j++) tmp.setColumn(j, rank(m.extractColumn(j))); return tmp; } template inline Matrix r_rank(const Matrix m) { Matrix tmp(m.nrows(), m.ncols()); for (unsigned int i = 0; i < m.nrows(); i++) tmp.setRow(i, rank(m.extractRow(i))); return tmp; } } #endif // define _ARRAY_HH_