54 throw ExceptionDimension(
StringUTF8(
"SquareMatrixDouble::SquareMatrixDouble(const std::vector<std::vector<double>> &m): ") +
_(
"the matrix is not square."));
61 throw ExceptionDimension(
StringUTF8(
"SquareMatrixDouble::SquareMatrixDouble(const std::vector<std::vector<double>> &m): ") +
_(
"the matrix is not square."));
65 #define MAX_GAUSS_W 40
83 if (
Gauss(
double(tmp), sigma) < 0.1)
91 for (
size_t r = 0; r < d; r++)
93 for (
size_t r = d; r < mat.
rows; r++)
96 for (
size_t c = 0; c < d; c++)
98 for (
size_t c = d; c < mat.
cols; c++)
101 for (
size_t r = 0; r < d; r++)
102 for (
size_t c = d + 1; c < mat.
cols; c++)
103 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
105 for (
size_t r = 0; r < d; r++)
106 for (
size_t c = 0; c < d; c++)
107 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
109 for (
size_t r = d + 1; r < mat.
rows; r++)
110 for (
size_t c = 0; c < d; c++)
111 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
113 for (
size_t r = d + 1; r < mat.
rows; r++)
114 for (
size_t c = d + 1; c < mat.
cols; c++)
115 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
140 if (
Gauss(
double(tmp), sigma) < 0.1)
149 for (
size_t r = 0; r < d; r++)
151 for (
size_t r = d; r < mat.
rows; r++)
154 for (
size_t c = 0; c < d; c++)
156 for (
size_t c = d; c < mat.
cols; c++)
159 for (
size_t r = 0; r < d; r++)
160 for (
size_t c = d + 1; c < mat.
cols; c++)
161 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
163 for (
size_t r = 0; r < d; r++)
164 for (
size_t c = 0; c < d; c++)
165 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
167 for (
size_t r = d + 1; r < mat.
rows; r++)
168 for (
size_t c = 0; c < d; c++)
169 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
171 for (
size_t r = d + 1; r < mat.
rows; r++)
172 for (
size_t c = d + 1; c < mat.
cols; c++)
173 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
176 for (
size_t r = 0; r < mat.
rows; r++)
177 for (
size_t c = 0; c < mat.
cols; c++)
179 mat[r][c] *= double(d) - double(c);
180 mat[r][c] /= double(d);
206 if (
Gauss(
double(tmp), sigma) < 0.1)
214 for (
size_t r = 0; r < d; r++)
216 for (
size_t r = d; r < mat.
rows; r++)
219 for (
size_t c = 0; c < d; c++)
221 for (
size_t c = d; c < mat.
cols; c++)
224 for (
size_t r = 0; r < d; r++)
225 for (
size_t c = d + 1; c < mat.
cols; c++)
226 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
228 for (
size_t r = 0; r < d; r++)
229 for (
size_t c = 0; c < d; c++)
230 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
232 for (
size_t r = d + 1; r < mat.
rows; r++)
233 for (
size_t c = 0; c < d; c++)
234 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
236 for (
size_t r = d + 1; r < mat.
rows; r++)
237 for (
size_t c = d + 1; c < mat.
cols; c++)
238 mat[r][c] = mat[r][d] * mat[d][c] /
MULT;
241 for (
size_t r = 0; r < mat.
rows; r++)
242 for (
size_t c = 0; c < mat.
cols; c++)
244 mat[r][c] *= double(d) - double(r);
245 mat[r][c] /= double(d);
261 while ((Status) && (r <
rows))
265 while ((Status) && (c < r))
267 Status = (
At(r, c) == 0.0);
287 while ((Status) && (r <
rows - 1))
291 while ((Status) && (c <
cols))
293 Status = (
At(r, c) == 0.0);
313 while (status && (c <
int(
cols) - 2))
317 while (status && (r <
int(
rows)))
319 status = status && !
At(r, c);
339 while ((Status) && (r <
rows - 1))
343 while ((Status) && (c <
cols))
345 Status = (
At(r, c) == 0.0) && (
At(c, r) == 0.0);
364 for (
size_t r = 0; r <
rows; ++r)
379 for (
size_t r = 0; r <
rows; ++r)
390 for (
size_t l = 0; l <
rows - 1; ++l)
391 for (
size_t c = l + 1; c <
cols; ++c)
392 std::swap(
At(l, c),
At(c, l));
410 _(
"row or column index out of range"));
416 for (
size_t rindex = 0; rindex <
rows; ++rindex)
422 for (
size_t cindex = 0; cindex <
cols; ++cindex)
426 M[i][j] =
At(rindex, cindex);
452 _(
"row or column index out of range"));
457 if ((r + c) % 2 == 1)
476 return At(0, 0) *
At(1, 1) -
At(0, 1) *
At(1, 0);
483 det +=
At(0, 0) *
At(1, 1) *
At(2, 2);
484 det +=
At(1, 0) *
At(2, 1) *
At(0, 2);
485 det +=
At(2, 0) *
At(0, 1) *
At(1, 2);
486 det -=
At(0, 2) *
At(1, 1) *
At(2, 0);
487 det -=
At(1, 2) *
At(2, 1) *
At(0, 0);
488 det -=
At(2, 2) *
At(0, 1) *
At(1, 0);
505 size_t MaxNullInRow = 0;
506 size_t MaxNullInColumn = 0;
510 for (
size_t k = 0; k <
rows; ++k)
515 if (nr > MaxNullInRow)
521 if (nc > MaxNullInColumn)
523 MaxNullInColumn = nc;
530 if (MaxNullInRow > MaxNullInColumn)
532 for (
size_t k = 0; k <
cols; ++k)
534 double Erk =
At(r, k);
542 for (
size_t k = 0; k <
rows; k++)
544 double Ekc =
At(k, c);
567 for (
size_t r = 0; r <
rows; ++r)
568 for (
size_t c = 0; c <
cols; ++c)
591 for (
size_t r = 0; r < n; ++r)
592 for (
size_t c = 0; c < n; ++c)
618 for (
size_t c = 0; c < n - 1; ++c)
622 double Pivot = M[c][c];
623 double AbsMaxPivot = fabs(Pivot);
626 for (
size_t r = c + 1 ; r < n; ++r)
628 double Candidate = M[r][c];
630 if (fabs(Candidate) > AbsMaxPivot)
633 AbsMaxPivot = fabs(Pivot);
651 for (
size_t r = c + 1; r < n; ++r)
653 double Coeff = M[r][c];
657 double Scale = - Coeff / Pivot;
659 for (
size_t k = 0; k < n; ++k)
661 M[r][k] += M[c][k] *
Scale;
662 Id[r][k] += Id[c][k] *
Scale;
670 for (
int c =
int(n) - 1; c > 0; --c)
672 double DCoeff = M[c][c];
676 for (
int r = c - 1; r >= 0; --r)
678 double Coeff = M[r][c];
682 double Scale = - Coeff / DCoeff;
684 for (
size_t k = 0; k < n; ++k)
686 Id[r][k] += Id[c][k] *
Scale;
687 M[r][k] += M[c][k] *
Scale;
696 for (
size_t r = 0; r < n; r++)
718 L[0][0] = sqrt(
At(0, 0));
722 for (
int j = 0; j < int(n); ++j)
726 for (
int k = 0; k <= j - 1; ++k)
729 L[j][j] = sqrt(
At(j, j) - Sum);
731 for (
int i = j + 1 ; i < int(n); ++i)
735 for (
int k = 0; k <= j - 1; ++k)
736 Sum += L[i][k] * L[j][k];
738 L[i][j] = (
At(i, j) - Sum) / L[j][j];
744 L[1][0] =
At(0, 1) / L[0][0];
745 L[1][1] = sqrt(
At(1, 1) -
Sqr(
At(0, 1)) /
At(0, 0));
759 auto hess_mat(*
this);
761 std::vector<std::vector<double>> reflector_seeds;
763 reflector_seeds.reserve(n - 2);
765 for (
auto k = 0; k <= n - 3; ++k)
768 std::vector<double> u_k;
772 for (
auto i = k + 1; i < n; ++i)
773 u_k.push_back(hess_mat.At(i, k));
776 Scale(u_k.begin(), u_k.end(), 1.0 /
Pythagoras(u_k.begin(), u_k.end()));
777 reflector_seeds.push_back(u_k);
780 for (
auto j = k; j < n; ++j)
782 auto scale_current_column = 0.0;
783 auto it = u_k.begin();
785 for (
auto i = k + 1; i < n; ++i)
787 scale_current_column += (*it) * hess_mat.At(i, j);
793 for (
auto i = k + 1; i < n; ++i)
795 hess_mat.At(i, j) -= 2 * (*it) * scale_current_column;
801 for (
auto i = 0; i < n; ++i)
803 auto scale_current_row = 0.0;
804 auto it = u_k.begin();
806 for (
auto j = k + 1; j < n; ++j)
808 scale_current_row += hess_mat.At(i, j) * (*it);
814 for (
auto j = k + 1; j < n; ++j)
816 hess_mat.At(i, j) -= 2 * scale_current_row * (*it);
836 std::multimap<double, MatrixDouble> eigen_pairs;
847 eigen_vector_1[0][0] = 1.0;
848 eigen_vector_1[1][0] = 0;
849 eigen_vector_2[0][0] = 0.0;
850 eigen_vector_2[1][0] = 1.0;
852 eigen_pairs.insert(std::make_pair(g_1, eigen_vector_1));
853 eigen_pairs.insert(std::make_pair(g_3, eigen_vector_2));
857 auto delta =
Sqr(g_1 - g_3) + 4 *
Sqr(g_2);
858 auto sqrt_delta = sqrt(delta);
860 auto eigen_value_1 = (g_1 + g_3 + sqrt_delta) / 2.0;
861 auto eigen_value_2 = (g_1 + g_3 - sqrt_delta) / 2.0;
863 eigen_vector_1[0][0] = 1.0;
864 eigen_vector_1[1][0] = (eigen_value_1 - g_1) / g_2;
868 eigen_vector_2[0][0] = - eigen_vector_1[1][0];
869 eigen_vector_2[1][0] = eigen_vector_1[0][0];
871 eigen_pairs.insert(std::make_pair(eigen_value_1, eigen_vector_1));
872 eigen_pairs.insert(std::make_pair(eigen_value_2, eigen_vector_2));
895 size_t MaxIter = 100;
902 for (k = 0; k < n; k++)
903 Eigenvectors[k][k] = 1.0;
918 while ((k < MaxIter) && Loop)
926 for (i = 0; i < n; ++i)
927 for (j = 0; j < n; ++j)
930 double a = fabs(A[i][j]);
944 for (i = 0; i < n; i++)
947 double Xi = (A[q][q] - A[p][p]) / (2.0 * A[p][q]);
954 t =
Min(*Roots.begin(), *Roots.rbegin());
957 double c = 1.0 / sqrt(1 + t * t);
958 double s = t / sqrt(1 + t * t);
967 Eigenvectors *= Omega;
986 std::multimap<double, MatrixDouble> EigenPairs;
987 std::vector<bool> Flags(n);
989 for (k = 0; k < n; ++k)
992 for (i = 0; i < n; ++i)
994 double MaxEigenvalue = -1.0;
997 for (j = 0; j < n; ++j)
1000 double Value = A[j][j];
1002 if (Value > MaxEigenvalue)
1004 MaxEigenvalue = Value;
1011 EigenPairs.insert(std::make_pair(A[Argmax][Argmax], Eigenvectors.
MakeColumn(Argmax)));
1026 std::vector<double> eigenvalues, offdiag;
1027 tred2(z, eigenvalues, offdiag);
1030 for (
int i = 1; i < int(
rows); ++i)
1031 offdiag[i - 1] = offdiag[i];
1032 offdiag[
rows - 1] = 0.0;
1033 for (
int l = 0; l < int(
rows); ++l)
1039 for (m = l; m < int(
rows) - 1; ++m)
1041 double dd =
Abs(eigenvalues[m]) +
Abs(eigenvalues[m + 1]);
1042 if (
Abs(offdiag[m]) <= std::numeric_limits<double>::epsilon() * dd)
1047 if (iter++ >= maxiter)
1048 throw ExceptionRuntime(
StringUTF8(
"std::multimap<double, SMatrixDouble> SquareMatrixDouble::MakeTQLIEigensystem(size_t maxiter) const: ") +
_(
"Too many iterations."));
1049 double g = (eigenvalues[l + 1] - eigenvalues[l]) / (2.0 * offdiag[l]);
1051 g = eigenvalues[m] - eigenvalues[l] + offdiag[l] / (g + ((g < 0) ? -
Abs(r) :
Abs(r)));
1052 double s = 1.0, c = 1.0, p = 0.0;
1054 for (i = m - 1; i >= l; --i)
1056 double f = s * offdiag[i];
1057 double b = c * offdiag[i];
1062 eigenvalues[i + 1] -= p;
1068 g = eigenvalues[i + 1] - p;
1069 r = (eigenvalues[i] - g) * s + 2.0 * c * b;
1071 eigenvalues[i + 1] = g + p;
1073 for (
int k = 0; k < int(
rows); ++k)
1076 z[k][i + 1] = s * z.
At(k, i) + c * f;
1077 z.
At(k, i) = c * z.
At(k, i) - s * f;
1080 if ((r == 0.0) && (i >= l))
1082 eigenvalues[l] -= p;
1090 std::multimap<double, MatrixDouble> out;
1091 for (
size_t tmp = 0; tmp < eigenvalues.size(); ++tmp)
1094 for (
size_t r = 0; r <
rows; ++r)
1095 ev[r][0] = z[r][tmp];
1096 out.insert(std::make_pair(eigenvalues[tmp], ev));
1106 void SquareMatrixDouble::tred2(
SquareMatrixDouble &z, std::vector<double> &diag, std::vector<double> &offdiag)
const
1110 std::vector<double>(
rows).swap(diag);
1111 std::vector<double>(
rows).swap(offdiag);
1113 for (
int i =
int(
rows) - 1; i > 0; --i)
1120 for (
int k = 0; k < i; ++k)
1121 scale +=
Abs(z.
At(i, k));
1123 offdiag[i] = z.
At(i, l);
1126 for (
int k = 0; k < i; ++k)
1128 z.
At(i, k) /= scale;
1129 h +=
Sqr(z.
At(i, k));
1131 double f = z.
At(i, l);
1132 double g = (f > 0.0) ? -sqrt(h) : sqrt(h);
1133 offdiag[i] = scale * g;
1137 for (
int j = 0; j < i; ++j)
1139 z.
At(j, i) = z.
At(i, j) / h;
1141 for (
int k = 0; k < j + 1; ++k)
1142 g += z.
At(j, k) * z.
At(i, k);
1143 for (
int k = j + 1; k < i; ++k)
1144 g += z.
At(k, j) * z.
At(i, k);
1146 f += offdiag[j] * z.
At(i, j);
1148 double hh = f / (h + h);
1149 for (
int j = 0; j < i; ++j)
1152 offdiag[j] = g = offdiag[j] - hh * f;
1153 for (
int k = 0; k < j + 1; ++k)
1154 z.
At(j, k) -= f * offdiag[k] + g * z.
At(i, k);
1159 offdiag[i] = z.
At(i, l);
1164 for (
size_t i = 0; i <
rows; ++i)
1168 for (
size_t j = 0; j < i; ++j)
1171 for (
size_t k = 0; k < i; ++k)
1172 g += z.
At(i, k) * z.
At(k, j);
1173 for (
size_t k = 0; k < i; ++k)
1174 z.
At(k, j) -= g * z.
At(k, i);
1177 diag[i] = z.
At(i, i);
1179 for (
size_t j = 0; j < i; ++j)
1180 z.
At(j, i) = z.
At(i, j) = 0.0;
1197 const auto n = int(
GetCols());
1201 std::vector<std::complex<double>> eigenvalues;
1203 for (
auto k = 0; k < n; ++k)
1204 eigenvalues.push_back(std::complex<double>(
At(k, k)));
1218 for (
auto r = 0; r < n; ++r)
1219 for (
auto c = 0; c < n; ++c)
1220 a[r + 1][c + 1] =
At(r, c);
1225 for (
auto r = 0; r < n; ++r)
1226 for (
auto c = 0; c < n; ++c)
1227 a[r + 1][c + 1] = hess_mat.
At(r, c);
1230 std::vector<double> wr(n + 1);
1231 std::vector<double> wi(n + 1);
1234 double z, y, x, w, v, u, t = 0.0, s, r, q, p;
1240 for (
auto i = 1; i <= n; ++i)
1241 for (
auto j =
Max(i - 1, 1); j <= n; ++j)
1242 anorm +=
Abs(a[i][j]);
1252 for (l = nn; l >= 2; --l)
1254 s =
Abs(a[l - 1][l - 1]) +
Abs(a[l][l]);
1258 if (
Abs(a[l][l - 1]) + s == s)
1271 y = a[nn - 1][nn - 1];
1272 w = a[nn][nn - 1] * a[nn - 1][nn];
1285 wr[nn - 1] = wr[nn] = x + z;
1290 wi[nn - 1] = wi[nn] = 0.0;
1294 wr[nn - 1] = wr[nn] = x + p;
1295 wi[nn - 1] = -(wi[nn] = z);
1303 std::cout <<
"ARGH !!!" << std::endl;
1309 for (
auto i = 1; i <= nn; ++i)
1312 s =
Abs(a[nn][nn - 1]) +
Abs(a[nn - 1][nn - 2]);
1314 w = -0.4375 *
Sqr(s);
1319 for (m = (nn - 2); m >= l; --m)
1325 p = (r * s - w) / a[m + 1][m] + a[m][m + 1];
1326 q = a[m + 1][m + 1] - z - r - s;
1327 r = a[m + 2][m + 1];
1338 u =
Abs(a[m][m - 1]) * (
Abs(q) +
Abs(r));
1339 v =
Abs(p) * (
Abs(a[m - 1][m - 1]) +
Abs(z) +
Abs(a[m + 1][m + 1]));
1345 for (
auto i = m + 2; i <= nn; ++i)
1353 for (
auto k = m; k <= nn - 1; ++k)
1361 q = a[k + 1][k - 1];
1365 r = a[k + 2][k - 1];
1368 if ((x =
Abs(p) +
Abs(q) +
Abs(r)) != 0.0)
1382 a[k][k - 1] = -a[k][k - 1];
1385 a[k][k - 1] = -s * x;
1394 for (
auto j = k; j <= nn; ++j)
1396 p = a[k][j] + q * a[k + 1][j];
1400 p += r * a[k + 2][j];
1401 a[k + 2][j] -= p * z;
1404 a[k + 1][j] -= p * y;
1408 auto mmin = nn < k + 3 ? nn : k + 3;
1410 for (
auto i = l; i <= mmin; ++i)
1412 p = x * a[i][k] + y * a[i][k + 1];
1417 a[i][k + 2] -= p * r;
1420 a[i][k + 1] -= p * q;
1431 std::vector<std::complex<double>> eigenvalues;
1433 for (
auto k = 1; k <=n; ++k)
1434 eigenvalues.push_back(std::complex<double>(wr[k], wi[k]));
1438 std::sort(eigenvalues.begin(), eigenvalues.end(), [](
const std::complex<double> &c1,
const std::complex<double> &c2)
1440 return norm(c1) < norm(c2);
1448 Cloner::Register<SquareMatrixDouble>();
std::vector< std::complex< double > > Eigenvalues(size_t max_iter=30) const
Extract eigenvalues for matrix of real (eigenvalues may be complex)
static SquareMatrixDouble NewGaussianSobelX(double sigma)
Create Gaussian Sobel X derivation mask.
bool IsDiagonal() const
Check if matrix is diagonal.
bool IsUpperTriangular() const
Check if matrix is upper triangular.
size_t GetRows() const noexcept
Returns the number of rows.
size_t GetCols() const noexcept
Returns the number of columns.
SquareMatrixDouble MakeGaussJordanInverse() const
Invert using Gauss-Jordan elimination.
virtual SquareMatrixDouble & Transpose() override
Transposition.
SquareMatrixDouble(size_t size, double value=0.0)
Constructor.
double MeanGauss(double x, double sigma)
Computes Gauss function at x for a given standard deviation (centered in 0) – to use with matrices...
const T & Max(const T &a, const T &b)
Returns the max of two values.
Matrix MakeColumn(size_t k) const
Extracts one column from matrix.
size_t CountNullCellsInColumn(size_t c) const
Counts null cells in column.
size_t CountNullCellsInRow(size_t r) const
Counts null cells in row.
SquareMatrixDouble MakeUpperHessenberg() const
Get the upper Hessenberg form of matrix.
#define CRN_END_CLASS_CONSTRUCTOR(classname)
Defines a class constructor.
double Trace() const
Trace.
std::multimap< double, MatrixDouble > MakeJacobiEigensystem(size_t MaxIteration=100) const
Perform diagonalization for symmetric matrix.
static SquareMatrixDouble NewGaussian(double sigma)
Create Gaussian matrix.
void Scale(ITER it_begin, ITER it_end, const double s)
Scale a collection of numbers.
static SquareMatrixDouble NewIdentity(size_t n)
Create identity matrix.
void MultRow(size_t r, double v)
Scale one row from matrix.
bool IsUpperHessenberg() const
Check if matrix is upper Hessenberg.
SquareMatrixDouble MakeInverse() const
Invert using determinants.
std::multimap< double, MatrixDouble > MakeTQLIEigensystem(size_t maxiter=30) const
Perform diagonalization for positive symmetric matrix.
std::set< double > RealRoots(double a, double b, double c) noexcept
Real roots of trinom .
#define CRN_DATA_FACTORY_REGISTER(elemname, classname)
Registers a class to the data factory.
constexpr SumType< T > Sqr(const T &v) noexcept(noexcept(v *v))
Returns the square of a value.
void SetAll(const T &v)
Set all elements.
double Gauss(double x, double sigma)
Computes Gauss function at x for a given standard deviation (centered in 0)
constexpr int SignOf(const T &x) noexcept(noexcept(x< 0))
Returns the sign (-1 or 1) of a value.
void Abs(Image< T > &img, typename std::enable_if< std::is_arithmetic< T >::value >::type *dummy=nullptr) noexcept
Replaces each pixel by its absolute value.
const T & Min(const T &a, const T &b)
Returns the min of two values.
const double & At(size_t pos) const noexcept
static SquareMatrixDouble NewGaussianSobelY(double sigma)
Create Gaussian Sobel Y derivation mask.
Square double matrix class.
bool IsLowerTriangular() const
Check if matrix is lower triangular.
std::multimap< double, MatrixDouble > MakeSpectralEigensystem() const
Perform diagonalization for 2x2 symmetric matrix.
SquareMatrixDouble MakeMinor(size_t r, size_t c) const
Minor matrix.
A character string class.
std::pair< size_t, size_t > Argmax() const
double Cofactor(size_t r, size_t c) const
Cofactor.
double Determinant() const
Determinant.
SquareMatrixDouble MakeCholesky() const
Get the lower triangular factor in Cholesky decomposition.
double DiagonalProduct() const
Diagonal product.
#define CRN_BEGIN_CLASS_CONSTRUCTOR(classname)
Defines a class constructor.
double Pythagoras(double a, double b) noexcept
Computes sqrt(a²+b²) without destructive underflow or overflow.
void SwapRows(size_t r1, size_t r2)
Swap two rows in matrix.