diff --git a/mlpp/lin_alg/lin_alg.cpp b/mlpp/lin_alg/lin_alg.cpp index 9c14435..824f649 100644 --- a/mlpp/lin_alg/lin_alg.cpp +++ b/mlpp/lin_alg/lin_alg.cpp @@ -932,34 +932,44 @@ MLPPLinAlg::QRDResult MLPPLinAlg::qrd(const Ref &A) { return res; } -/* -MLPPLinAlg::CholeskyResult MLPPLinAlg::cholesky(std::vector> A) { - std::vector> L = zeromat(A.size(), A[0].size()); - for (uint32_t j = 0; j < L.size(); j++) { // Matrices entered must be square. No problem here. - for (uint32_t i = j; i < L.size(); i++) { +MLPPLinAlg::CholeskyResult MLPPLinAlg::cholesky(const Ref &A) { + Size2i a_size = A->size(); + + CholeskyResult res; + + ERR_FAIL_COND_V(a_size.x != a_size.y, res); + + Ref L = zeromatnm(a_size.y, a_size.x); + + for (int j = 0; j < a_size.y; ++j) { // Matrices entered must be square. No problem here. + for (int i = j; i < a_size.y; ++i) { if (i == j) { real_t sum = 0; - for (uint32_t k = 0; k < j; k++) { - sum += L[i][k] * L[i][k]; + + for (int k = 0; k < j; k++) { + real_t lik = L->element_get(i, k); + + sum += lik * lik; } - L[i][j] = Math::sqrt(A[i][j] - sum); + + L->element_set(i, j, Math::sqrt(A->element_get(i, j) - sum)); } else { // That is, i!=j real_t sum = 0; - for (uint32_t k = 0; k < j; k++) { - sum += L[i][k] * L[j][k]; + + for (int k = 0; k < j; k++) { + sum += L->element_get(i, k) * L->element_get(j, k); } - L[i][j] = (A[i][j] - sum) / L[j][j]; + + L->element_set(i, j, (A->element_get(i, j) - sum) / L->element_get(j, j)); } } } - CholeskyResult res; res.L = L; - res.Lt = transpose(L); // Indeed, L.T is our upper triangular matrix. + res.Lt = L->transposen(); // Indeed, L.T is our upper triangular matrix. return res; } -*/ /* real_t MLPPLinAlg::sum_elements(std::vector> A) { @@ -990,62 +1000,54 @@ Ref MLPPLinAlg::flattenvvnv(const Ref &A) { return res; } -/* -std::vector MLPPLinAlg::solve(std::vector> A, std::vector b) { - return mat_vec_mult(inverse(A), b); +Ref MLPPLinAlg::solve(const Ref &A, const Ref &b) { + return A->inverse()->mult_vec(b); } -bool MLPPLinAlg::positiveDefiniteChecker(std::vector> A) { - auto eig_result = eig(A); - auto eigenvectors = std::get<0>(eig_result); - auto eigenvals = std::get<1>(eig_result); +bool MLPPLinAlg::positive_definite_checker(const Ref &A) { + EigenResult eig_result = eigen(A); - std::vector eigenvals_vec; - for (uint32_t i = 0; i < eigenvals.size(); i++) { - eigenvals_vec.push_back(eigenvals[i][i]); - } - for (uint32_t i = 0; i < eigenvals_vec.size(); i++) { - if (eigenvals_vec[i] <= 0) { // Simply check to ensure all eigenvalues are positive. + Ref eigenvals = eig_result.eigen_values; + Size2i eigenvals_size = eigenvals->size(); + + for (int i = 0; i < eigenvals_size.y; ++i) { + if (eigenvals->element_get(i, i) <= 0) { // Simply check to ensure all eigenvalues are positive. return false; } } + return true; } -bool MLPPLinAlg::negativeDefiniteChecker(std::vector> A) { - auto eig_result = eig(A); - auto eigenvectors = std::get<0>(eig_result); - auto eigenvals = std::get<1>(eig_result); +bool MLPPLinAlg::negative_definite_checker(const Ref &A) { + EigenResult eig_result = eigen(A); - std::vector eigenvals_vec; - for (uint32_t i = 0; i < eigenvals.size(); i++) { - eigenvals_vec.push_back(eigenvals[i][i]); - } - for (uint32_t i = 0; i < eigenvals_vec.size(); i++) { - if (eigenvals_vec[i] >= 0) { // Simply check to ensure all eigenvalues are negative. + Ref eigenvals = eig_result.eigen_values; + Size2i eigenvals_size = eigenvals->size(); + + for (int i = 0; i < eigenvals_size.y; ++i) { + if (eigenvals->element_get(i, i) >= 0) { // Simply check to ensure all eigenvalues are negative. return false; } } + return true; } -bool MLPPLinAlg::zeroEigenvalue(std::vector> A) { - auto eig_result = eig(A); - auto eigenvectors = std::get<0>(eig_result); - auto eigenvals = std::get<1>(eig_result); +bool MLPPLinAlg::zero_eigenvalue(const Ref &A) { + EigenResult eig_result = eigen(A); - std::vector eigenvals_vec; - for (uint32_t i = 0; i < eigenvals.size(); i++) { - eigenvals_vec.push_back(eigenvals[i][i]); - } - for (uint32_t i = 0; i < eigenvals_vec.size(); i++) { - if (eigenvals_vec[i] == 0) { - return true; + Ref eigenvals = eig_result.eigen_values; + Size2i eigenvals_size = eigenvals->size(); + + for (int i = 0; i < eigenvals_size.y; ++i) { + if (eigenvals->element_get(i, i) == 0) { // TODO should it use is_equal_approx? + return false; } } - return false; + + return true; } -*/ Ref MLPPLinAlg::flattenmnv(const Vector> &A) { Ref a; diff --git a/mlpp/lin_alg/lin_alg.h b/mlpp/lin_alg/lin_alg.h index 3a5d245..ada3813 100644 --- a/mlpp/lin_alg/lin_alg.h +++ b/mlpp/lin_alg/lin_alg.h @@ -121,28 +121,22 @@ public: QRDResult qrd(const Ref &A); - /* struct CholeskyResult { - std::vector> L; - std::vector> Lt; + Ref L; + Ref Lt; }; - CholeskyResult cholesky(std::vector> A); - */ + CholeskyResult cholesky(const Ref &A); //real_t sum_elements(std::vector> A); Ref flattenvvnv(const Ref &A); + Ref solve(const Ref &A, const Ref &b); - /* - std::vector solve(std::vector> A, std::vector b); + bool positive_definite_checker(const Ref &A); + bool negative_definite_checker(const Ref &A); - bool positiveDefiniteChecker(std::vector> A); - - bool negativeDefiniteChecker(std::vector> A); - - bool zeroEigenvalue(std::vector> A); - */ + bool zero_eigenvalue(const Ref &A); // VECTOR FUNCTIONS diff --git a/test/mlpp_tests.cpp b/test/mlpp_tests.cpp index 949839b..0e05314 100644 --- a/test/mlpp_tests.cpp +++ b/test/mlpp_tests.cpp @@ -1157,7 +1157,7 @@ void MLPPTests::test_new_math_functions() { PLOG_MSG(alg.gram_schmidt_process(P)->to_string()); MLPPLinAlg::QRDResult qrd_result = alg.qrd(P); // It works! - + //[MLPPMatrix: // [ 0.857143 -0.394286 -0.331429 ] // [ 0.428571 0.902857 0.034286 ] @@ -1172,27 +1172,26 @@ void MLPPTests::test_new_math_functions() { PLOG_MSG(qrd_result.R->to_string()); } void MLPPTests::test_positive_definiteness_checker() { - /* - //MLPPStat stat; MLPPLinAlg alg; - //MLPPActivation avn; - //MLPPCost cost; - //MLPPData data; - //MLPPConvolutions conv; // Checking positive-definiteness checker. For Cholesky Decomp. - std::vector> A = { + std::vector> A_arr = { { 1, -1, -1, -1 }, { -1, 2, 2, 2 }, { -1, 2, 3, 1 }, { -1, 2, 1, 4 } }; - std::cout << std::boolalpha << alg.positiveDefiniteChecker(A) << std::endl; + Ref A(memnew(MLPPMatrix(A_arr))); + + PLOG_MSG("positive_definite_checker Example:"); + PLOG_MSG(String::bool_str(alg.positive_definite_checker(A))); + + PLOG_MSG("Cholesky Example:"); MLPPLinAlg::CholeskyResult chres = alg.cholesky(A); // works. - alg.printMatrix(chres.L); - alg.printMatrix(chres.Lt); - */ + + PLOG_MSG(chres.L->to_string()); + PLOG_MSG(chres.Lt->to_string()); } // real_t f(real_t x){