From 6e53dd0034ba7c59b54fd6977543b9d75a7601dd Mon Sep 17 00:00:00 2001 From: Abdullah Tahiri Date: Sun, 13 Dec 2020 08:43:02 +0100 Subject: [PATCH] GCS: Refactor dependent parameters identification which is only available for DenseQR --- src/Mod/Sketcher/App/planegcs/GCS.cpp | 206 ++++++++++++-------------- src/Mod/Sketcher/App/planegcs/GCS.h | 9 ++ 2 files changed, 106 insertions(+), 109 deletions(-) diff --git a/src/Mod/Sketcher/App/planegcs/GCS.cpp b/src/Mod/Sketcher/App/planegcs/GCS.cpp index ad77b95526..d76b568be8 100644 --- a/src/Mod/Sketcher/App/planegcs/GCS.cpp +++ b/src/Mod/Sketcher/App/planegcs/GCS.cpp @@ -3923,122 +3923,22 @@ SolverReportingManager::Manager().LogToFile("GCS::System::diagnose()\n"); if(qrAlgorithm==EigenDenseQR){ makeDenseQRDecomposition( J, jacobianconstraintmap, qrJT, paramsNum, constrNum, rank, R); + + if (J.rows() > 0) { + identifyDependentGeometryParametersInTransposedJacobianDenseQRDecomposition( + qrJT, + pdiagnoselist, + paramsNum, rank); + } } else if(qrAlgorithm==EigenSparseQR){ makeSparseQRDecomposition( J, jacobianconstraintmap, SqrJT, paramsNum, constrNum, rank, R); } + + if (J.rows() > 0) { - // DETECTING CONSTRAINT SOLVER PARAMETERS - // - // NOTE: This is only true for dense QR with full pivoting, because solve parameters get reordered. - // I am unable to adapt it to Sparse QR. (abdullah). See: - // - // https://stackoverflow.com/questions/49009771/getting-rows-transpositions-with-sparse-qr - // https://forum.kde.org/viewtopic.php?f=74&t=151239 - // - // R (original version, not R here which is trimmed to not have empty rows) - // has paramsNum rows, the first "rank" rows correspond to parameters that are constraint - - // Calculate the Permutation matrix from the Transposition matrix - Eigen::PermutationMatrix rowPermutations; - - rowPermutations.setIdentity(paramsNum); - - if(qrAlgorithm==EigenDenseQR){ // P.J.P' = Q.R see https://eigen.tuxfamily.org/dox/classEigen_1_1FullPivHouseholderQR.html - const MatrixIndexType rowTranspositions = qrJT.rowsTranspositions(); - - for(int k = 0; k < rank; ++k) - rowPermutations.applyTranspositionOnTheRight(k, rowTranspositions.coeff(k)); - } -#ifdef EIGEN_SPARSEQR_COMPATIBLE - else if(qrAlgorithm==EigenSparseQR){ - // J.P = Q.R, see https://eigen.tuxfamily.org/dox/classEigen_1_1SparseQR.html - // There is no rowsTransposition in this QR decomposition. - // TODO: This detection method won't work for SparseQR - } - -#endif - -#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX - std::stringstream stream; -#endif - -// params (in the order of J) shown as independent from QR - std::set indepParamCols; - std::set depParamCols; - - for (int j=0; j < rank; j++) { - - int origRow = rowPermutations.indices()[j]; - - indepParamCols.insert(origRow); - - // NOTE: Q*R = transpose(J), so the row of R corresponds to the col of J (the rows of transpose(J)). - // The cols of J are the parameters, the rows are the constraints. -#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX - stream << "R row " << j << " = J col " << origRow << std::endl; -#endif - } - -#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX - std::string tmp = stream.str(); - - SolverReportingManager::Manager().LogString(tmp); -#endif - - // If not independent, must be dependent - for(int j=0; j < paramsNum; j++) { - auto result = indepParamCols.find(j); - if(result == indepParamCols.end()) { - depParamCols.insert(j); - } - } - - // Last (NumParams-rank) rows of Q construct the dependent part of J - // in conjunction with the R matrix - // Last (NumParams-rank) cols of Q never contribute as R is zero after the rank - /*std::set associatedParamCols; - for(int i = rank; i < paramsNum; i++) { - for(int j = 0; j < rank; j++) { - if(fabs(Q(i,j))>1e-10) { - for(int k = j; k < constrNum; k++) { - if(fabs(R(j,k))>1e-10) { - associatedParamCols.insert(k); - } - } - } - } - } - - for( auto param : associatedParamCols) { - auto pos = indepParamCols.find(rowPermutations.indices()[param]); - - if(pos!=indepParamCols.end()) { - indepParamCols.erase(pos); - } - }*/ -#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX - stream.flush(); - - stream << "Indep params: ["; - for(auto indep :indepParamCols) - stream << indep; - stream << "]" << std::endl; - - stream << "Dep params: ["; - for(auto dep :depParamCols) - stream << dep; - stream << "]" << std::endl; - - tmp = stream.str(); - SolverReportingManager::Manager().LogString(tmp); -#endif - for( auto param : depParamCols) { - pdependentparameters.push_back(pdiagnoselist[param]); - } - // Detecting conflicting or redundant constraints if (constrNum > rank) { // conflicting or redundant constraints for (int i=1; i < rank; i++) { @@ -4365,6 +4265,94 @@ void System::makeSparseQRDecomposition( Eigen::MatrixXd &J, std::map &j } } +void System::identifyDependentGeometryParametersInTransposedJacobianDenseQRDecomposition( + const Eigen::FullPivHouseholderQR& qrJT, + const GCS::VEC_pD &pdiagnoselist, + int paramsNum, int rank) +{ + // DETECTING CONSTRAINT SOLVER PARAMETERS + // + // NOTE: This is only true for dense QR with full pivoting, because solve parameters get reordered. + // I am unable to adapt it to Sparse QR. (abdullah). See: + // + // https://stackoverflow.com/questions/49009771/getting-rows-transpositions-with-sparse-qr + // https://forum.kde.org/viewtopic.php?f=74&t=151239 + // + // R (original version, not R here which is trimmed to not have empty rows) + // has paramsNum rows, the first "rank" rows correspond to parameters that are constraint + + // Calculate the Permutation matrix from the Transposition matrix + Eigen::PermutationMatrix rowPermutations; + + rowPermutations.setIdentity(paramsNum); + + // P.J.P' = Q.R see https://eigen.tuxfamily.org/dox/classEigen_1_1FullPivHouseholderQR.html + const MatrixIndexType rowTranspositions = qrJT.rowsTranspositions(); + + for(int k = 0; k < rank; ++k) + rowPermutations.applyTranspositionOnTheRight(k, rowTranspositions.coeff(k)); + +#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX + std::stringstream stream; +#endif + + // params (in the order of J) shown as independent from QR + std::set indepParamCols; + std::set depParamCols; + + for (int j=0; j < rank; j++) { + + int origRow = rowPermutations.indices()[j]; + + indepParamCols.insert(origRow); + + // NOTE: Q*R = transpose(J), so the row of R corresponds to the col of J (the rows of transpose(J)). + // The cols of J are the parameters, the rows are the constraints. +#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX + stream << "R row " << j << " = J col " << origRow << std::endl; +#endif + } + +#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX + std::string tmp = stream.str(); + + SolverReportingManager::Manager().LogString(tmp); +#endif + + // If not independent, must be dependent + for(int j=0; j < paramsNum; j++) { + auto result = indepParamCols.find(j); + if(result == indepParamCols.end()) { + depParamCols.insert(j); + } + } + +#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX + stream.flush(); + + stream << "Indep params: ["; + for(auto indep :indepParamCols) + stream << indep; + stream << "]" << std::endl; + + stream << "Dep params: ["; + for(auto dep :depParamCols) + stream << dep; + stream << "]" << std::endl; + + tmp = stream.str(); + SolverReportingManager::Manager().LogString(tmp); +#endif + + + for( auto param : depParamCols) { + pdependentparameters.push_back(pdiagnoselist[param]); + } + +} + + + void System::clearSubSystems() { isInit = false; diff --git a/src/Mod/Sketcher/App/planegcs/GCS.h b/src/Mod/Sketcher/App/planegcs/GCS.h index cee752fa36..5f1fc73fdd 100644 --- a/src/Mod/Sketcher/App/planegcs/GCS.h +++ b/src/Mod/Sketcher/App/planegcs/GCS.h @@ -141,6 +141,15 @@ namespace GCS Eigen::SparseQR, Eigen::COLAMDOrdering > &SqrJT, int ¶msNum, int &constrNum, int &rank, Eigen::MatrixXd &R); + // This function name is long for a reason: + // - Only for DenseQR + // - Only for Transposed Jacobian QR decomposition + void identifyDependentGeometryParametersInTransposedJacobianDenseQRDecomposition( + const Eigen::FullPivHouseholderQR& qrJT, + const GCS::VEC_pD &pdiagnoselist, + int paramsNum, int rank + ); + #ifdef _GCS_EXTRACT_SOLVER_SUBSYSTEM_ void extractSubsystem(SubSystem *subsys, bool isRedundantsolving); #endif