From d9a18a3d63680cd56815a54c4a96dc8150315120 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Mon, 23 Dec 2024 15:59:16 +0100 Subject: [PATCH] Root Poly/Eigen :OK --- MMVII/include/MMVII_nums.h | 5 +- MMVII/src/Matrix/cEigenEigenDecompos.cpp | 8 +- MMVII/src/Serial/ElemStrToVal.cpp | 11 + MMVII/src/UtiMaths/Polynoms.cpp | 295 +++++++++++------------ 4 files changed, 158 insertions(+), 161 deletions(-) diff --git a/MMVII/include/MMVII_nums.h b/MMVII/include/MMVII_nums.h index faeb68c3fd..418c7e9de3 100755 --- a/MMVII/include/MMVII_nums.h +++ b/MMVII/include/MMVII_nums.h @@ -993,6 +993,7 @@ template class cPolynom { public : typedef std::vector tCoeffs; + typedef cPtxd tCompl; cPolynom(const tCoeffs &); cPolynom(const cPolynom &); cPolynom(size_t aDegre); @@ -1007,10 +1008,12 @@ template class cPolynom static cPolynom RandomPolyg(std::vector & aVRoots,int aNbRoot,int aNbNoRoot,Type Interv,Type MinDist); - Type Value(const Type & aVal) const; + Type Value(const Type & aVal) const; + tCompl Value(const tCompl & aVal) const; /// return som(|a_k x^k|) , used for some bounding stuffs Type AbsValue(const Type & aVal) const; + cPolynom operator * (const cPolynom & aP2) const; cPolynom operator + (const cPolynom & aP2) const; cPolynom operator - (const cPolynom & aP2) const; diff --git a/MMVII/src/Matrix/cEigenEigenDecompos.cpp b/MMVII/src/Matrix/cEigenEigenDecompos.cpp index 084ba7c645..ad451b538d 100755 --- a/MMVII/src/Matrix/cEigenEigenDecompos.cpp +++ b/MMVII/src/Matrix/cEigenEigenDecompos.cpp @@ -93,9 +93,13 @@ void Bench_EigenDecompos(cParamExeBench & aParam) Tpl_Bench_EigenDecompos(aParam); } +#define INSTANTIATE_EIGEN_DECOMP(TYPE)\ +template class cResulEigenDecomp;\ +template class cDenseMatrix; -template class cResulEigenDecomp; -template class cDenseMatrix; +INSTANTIATE_EIGEN_DECOMP(tREAL4); +INSTANTIATE_EIGEN_DECOMP(tREAL8); +INSTANTIATE_EIGEN_DECOMP(tREAL16); }; diff --git a/MMVII/src/Serial/ElemStrToVal.cpp b/MMVII/src/Serial/ElemStrToVal.cpp index 67716b10bf..1f4a5d03e3 100644 --- a/MMVII/src/Serial/ElemStrToVal.cpp +++ b/MMVII/src/Serial/ElemStrToVal.cpp @@ -1395,6 +1395,17 @@ std::string FixDigToStr(double aSignedVal,int aNbBef,int aNbAfter) return aBuf; } + // ================ double ============================================== + +template <> std::string cStrIO::ToStr(const tREAL16 & aD) +{ + return cStrIO::ToStr((tREAL8) (aD)); +} + +template <> std::string cStrIO::ToStr(const tREAL4 & aD) +{ + return cStrIO::ToStr((tREAL8) (aD)); +} // ================ std::string ============================================== diff --git a/MMVII/src/UtiMaths/Polynoms.cpp b/MMVII/src/UtiMaths/Polynoms.cpp index 3064ccc416..977f6bf95b 100755 --- a/MMVII/src/UtiMaths/Polynoms.cpp +++ b/MMVII/src/UtiMaths/Polynoms.cpp @@ -1,27 +1,12 @@ +#define WITH_MMV1_FUNCTION true -#if defined(__GNUC__) && !defined(__clang__) -# pragma GCC diagnostic push -# pragma GCC diagnostic ignored "-Wmaybe-uninitialized" -#endif +// #include "cMMVII_Appli.h" +#include "MMVII_Matrix.h" -#include "cMMVII_Appli.h" +#if (WITH_MMV1_FUNCTION) #include "V1VII.h" -#include -// #include "unsupported/Eigen/Polynomials" +#endif -//#if defined(__GNUC__) && !defined(__clang__) -// # pragma GCC diagnostic pop -// #endif - - - // #include "include/MMVII_nums.h" - // #include "include/MMVII_Bench.h" - -//#include -//#include - -using namespace Eigen; -using namespace std; namespace MMVII { @@ -36,73 +21,136 @@ namespace MMVII template class cEigenPolynRoots { public : - cEigenPolynRoots(const cPolynom &) ; - const std::vector & RealRoots() const; - const VectorXcd & ComplexRoots() const; - bool RootIsReal(const std::complex &,std::string * sayWhy=nullptr); - const MatrixXd & CompM() const; + + typedef cDenseMatrix tMatComp; + typedef cPtxd tCompl; + + cEigenPolynRoots(const cPolynom &,Type aEps,int aNbIterMax) ; + bool RootIsReal(const tCompl & ,std::string * sayWhy=nullptr); + + const std::vector & RealRoots() const {return mRR;} + const std::vector & ComplexRoots() const {return mCR;} + const tMatComp & CompM() const {return mCompM;} + + tCompl Refine(const tCompl & aV0,Type aEps,int aNbIter) const; + private : + static Type PolRelAccuracy() ; + static Type PolAbsAccuracy() ; + static Type ComplRelAccuracy() ; + static Type ComplAbsAccuracy() ; cPolynom mPol; cPolynom mDPol; size_t mDeg; - MatrixXd mCompM; ///< companion matrix - VectorXcd mCEV; ///< Complex eigen values + size_t mSzMat; /// to avoid size 0 + tMatComp mCompM; ///< companion matrix std::vector mRR; ///< Real roots + std::vector mCR; ///< Real roots }; -template cEigenPolynRoots::cEigenPolynRoots(const cPolynom & aPol) : +template cEigenPolynRoots::cEigenPolynRoots(const cPolynom & aPol,Type aEps,int aNbIter) : mPol (aPol), mDPol (aPol.Deriv()), mDeg (mPol.Degree()), - mCompM (mDeg,mDeg) + mSzMat (std::max((size_t)1,mDeg)), + mCompM (mSzMat,mSzMat,eModeInitImage::eMIA_Null) { if (mDeg==0) return; - for (size_t aI = 0; aI < mDeg; ++aI) - for (size_t aJ = 0; aJ < mDeg; ++aJ) - mCompM(aI,aJ) = 0; - // fill the diagonal up the principal diag for (size_t aK = 0; aK < mDeg-1; ++aK) { - mCompM(aK+1, aK ) = 1; + mCompM.SetElem(aK+1, aK,1); } const Type & aHighCoeff = mPol[mDeg]; // Fill last line with normalized coeff for (size_t aK = 0; aK < mDeg; ++aK) { - mCompM(aK, mDeg - 1) = -mPol[aK] / aHighCoeff; + mCompM.SetElem(aK, mDeg - 1, -mPol[aK] / aHighCoeff); } - EigenSolver aSolver(mCompM); - mCEV = aSolver.eigenvalues() ; + cResulEigenDecomp aRED = mCompM.Eigen_Decomposition() ; + - for (const auto & aC : mCEV) - if (RootIsReal(aC)) - mRR.push_back(aC.real()); - std::sort(mRR.begin(),mRR.end()); + for (size_t aK = 0; aK < mDeg; ++aK) + { + tCompl aCR(aRED.mEigenVal_R(aK),aRED.mEigenVal_I(aK)); + + aCR = Refine(aCR,aEps*PolAbsAccuracy(),aNbIter); + mCR.push_back(aCR); + + if (RootIsReal(aCR)) + mRR.push_back(aCR.x()); + } + std::sort(mRR.begin(),mRR.end()); } -template const std::vector & cEigenPolynRoots::RealRoots() const {return mRR;} +template cPtxd cEigenPolynRoots::Refine(const tCompl & aVal0,Type aEps,int aNbIter) const +{ + Type aSqEps = Square(aEps); + tCompl aLastVal = aVal0; + tCompl aLastEval = mPol.Value(aLastVal); + Type aLastSqN2 = SqN2(aLastEval); + + for (int aKIt=0 ; aKIt const VectorXcd & cEigenPolynRoots::ComplexRoots() const {return mCEV;} + // tCompl Refine(const tCompl & aV0,int aNbIter=5); -template const MatrixXd & cEigenPolynRoots::CompM() const {return mCompM;} /** Also the question seems pretty basic, it becomes more complicated due to numericall approximation */ -template bool cEigenPolynRoots::RootIsReal(const std::complex & aC,std::string * sayWhy) +template <> tREAL4 cEigenPolynRoots::PolRelAccuracy() {return 1e-5;} +template <> tREAL8 cEigenPolynRoots::PolRelAccuracy() {return 1e-9;} +template <> tREAL16 cEigenPolynRoots::PolRelAccuracy(){return 1e-11;} + +template <> tREAL4 cEigenPolynRoots::ComplRelAccuracy() {return 1e-8;} +template <> tREAL8 cEigenPolynRoots::ComplRelAccuracy() {return 1e-8;} +template <> tREAL16 cEigenPolynRoots::ComplRelAccuracy() {return 1e-8;} + +template <> tREAL4 cEigenPolynRoots::PolAbsAccuracy() {return 0.01;} +template <> tREAL8 cEigenPolynRoots::PolAbsAccuracy() {return 1e-5;} +template <> tREAL16 cEigenPolynRoots::PolAbsAccuracy() {return 1e-7;} + //static Type ComplRelAccuracy() ; + //static Type ComplAbsAccuracy() ; + + +template bool cEigenPolynRoots::RootIsReal(const tCompl & aC,std::string * sayWhy) { + // [1] Test is "aC" is a real number + Type C_i =aC.y(); - tREAL8 C_i = aC.imag(); // [1.1] if absolute value of imaginary part is "big" it's not if (std::abs(C_i) > 1e-5) { @@ -111,9 +159,9 @@ template bool cEigenPolynRoots::RootIsReal(const std::complex return false; } - tREAL8 C_r = aC.real(); + Type C_r =aC.x(); // [1.1] if relative imaginary part is "big" - if (std::abs(C_i) > 1e-8 * (std::abs(C_r)+1e-5)) + if (std::abs(C_i) > ComplRelAccuracy() * (std::abs(C_r)+1e-5)) { if (sayWhy) *sayWhy = "RELAT REAL COMPLEX=" + ToStr(std::abs(C_i)/(std::abs(C_r)+1e-5)); @@ -121,18 +169,18 @@ template bool cEigenPolynRoots::RootIsReal(const std::complex } // [2] Test - tREAL8 aAbsVP = std::abs(mPol.Value(C_r)); + Type aAbsVP = std::abs(mPol.Value(C_r)); // [2.1] if absolute value of polynom is big - if (aAbsVP > 1e-5) + if (aAbsVP > PolAbsAccuracy()) { if (sayWhy) *sayWhy = "ABS VALUE POL " + ToStr(aAbsVP); return false; } - tREAL8 aAVA = mPol.AbsValue(C_r); + Type aAVA = mPol.AbsValue(C_r); // [2.1] if absolute value of polynom is big relatively to norm - if (aAbsVP > 1e-8 * (aAVA+1e-5)) + if (aAbsVP > PolRelAccuracy() * (aAVA+1e-5)) { if (sayWhy) *sayWhy = "RELATIVE VALUE POL " + ToStr(aAbsVP/(aAVA+1e-5)); @@ -145,14 +193,19 @@ template bool cEigenPolynRoots::RootIsReal(const std::complex return true; } - +template class cEigenPolynRoots; template class cEigenPolynRoots; +template class cEigenPolynRoots; + template void My_Roots(const cPolynom & aPol1) { -} -template <> void My_Roots(const cPolynom & aPol1) -{ + // 4 => low accuracy + // 16 => low timing + // As we just want to see that there is no regression for standard double + if (sizeof(Type) != 8) return; + + // StdOut() << "DDDD " << aPol1.Degree() << "\n"; // (X2+1)(X-1) = X3-X2+X-1 int aNb=300; @@ -162,7 +215,7 @@ template <> void My_Roots(const cPolynom & aPol1) cAutoTimerSegm aTimeEigen(GlobAppTS(),"Eigen"); for (int aK=0 ; aK aEPR(aPol1); + cEigenPolynRoots aEPR(aPol1,1e-3,10); } cAutoTimerSegm aTimeV1(GlobAppTS(),"V1"); @@ -172,12 +225,13 @@ template <> void My_Roots(const cPolynom & aPol1) } cAutoTimerSegm aTimeOthers(GlobAppTS(),"Others"); - cEigenPolynRoots aEPR(aPol1); + cEigenPolynRoots aEPR(aPol1,1e-3,10); auto aV2 = aEPR.RealRoots(); auto aV1 = aPol1.RealRoots(1e-20,60); - StdOut() << " SZzzzZ= " << aV1.size() << " " << aV2.size() << "\n"; if (aV1.size() != aV2.size()) { + StdOut() << " SZzzzZ= " << aV1.size() << " " << aV2.size() << " SIZOFTYPE=" << sizeof(Type) << "\n"; + StdOut() << aV1 << aV2 << "\n"; StdOut() << "Coeffs=" << aPol1.VCoeffs() << "\n"; StdOut() << "V1=" << aV1 << "\n"; StdOut() << "V2=" << aV2 << "\n"; @@ -188,115 +242,24 @@ template <> void My_Roots(const cPolynom & aPol1) StdOut() << "R=" << isR << " C=" << aC << " W=" << strWhy << "\n"; } - StdOut() << "DET=" << aEPR.CompM().determinant() << "\n"; StdOut() << " ------------------ MAT ---------------------\n"; StdOut() << aEPR.CompM() << "\n"; getchar(); } - else - StdOut() << aV1 << aV2 << "\n"; // (X2+1)(X-1) = X3-X2+X-1 // vector coeffs = {-1,1,-1,1}; } -#if (0) - -/* -MatrixXd createCompanionMatrix(const vector& coeffs) { -StdOut() << "BEGIN createCompanionMatrixcreateCompanionMatrixcreateCompanionMatrix\n"; - int n = coeffs.size(); - MatrixXd companionMatrix(n - 1, n - 1); + // return V1RealRoots(mVCoeffs,aTol,ItMax); -StdOut() << "LLLLL " << __LINE__ << "\n"; - - // Remplir la matrice compagnon - // for (int i = 0; i < n - 1; ++i) { - for (int i = 0; i < n - 2; ++i) { - companionMatrix(i+1, i ) = 1; // Remplir la diagonale au-dessus de la diagonale principale - } -StdOut() << "LLLLL " << __LINE__ << "\n"; - - // Remplir la dernière colonne de la matrice avec les coefficients du polynôme - for (int i = 0; i < n - 1; ++i) { - // companionMatrix(i, n - 1) = -coeffs[i] / coeffs[n - 1]; - companionMatrix(i, n - 2) = -coeffs[i] / coeffs[n - 1]; - } - -StdOut() << "END createCompanionMatrixcreateCompanionMatrixcreateCompanionMatrix\n"; - return companionMatrix; -} -// Fonction principale -int My_Roots() { - // (X2+1)(X-1) = X3-X2+X-1 - // Exemple de coefficients d'un polynôme : x^3 - 6x^2 + 11x - 6 - // (x-1) (x2+ 5x +6) - // vector coeffs = {1, -6, 11, -6}; - vector coeffs = {1,-1,1,-1}; - - - - // Créer la matrice compagnon - MatrixXd companionMatrix = createCompanionMatrix(coeffs); - - std::cout << "mmmm " << companionMatrix << "\n"; - - // Calculer les valeurs propres (racines du polynôme) - EigenSolver solver(companionMatrix); - std::cout << "eeeee " << solver.eigenvalues() << "\n"; - - // const auto& theEigenV = solver.eigenvalues(); - VectorXcd eivals = solver.eigenvalues(); - std::cout << "EEEEEE " << eivals << "\n"; - - VectorXd R_roots = solver.eigenvalues().real(); - VectorXd I_roots = solver.eigenvalues().imag(); - - for (int i = 0; i < R_roots.size(); ++i) { - cout << "R=" << R_roots[i] << " I=" << I_roots[i] << endl; - } -*/ - -/* - VectorXd roots = solver.eigenvalues().real(); - - // Afficher les racines - cout << "Les racines du polynôme sont : " << endl; - for (int i = 0; i < roots.size(); ++i) { - cout << roots[i] << endl; - } -*/ - - return 0; -} -#endif - -/* - */ - - - - -#if (0) -void TestPolynEigen() +template std::vector V2RealRoots(const cPolynom & aPol, Type aTol,int aNbMaxIter) { - // cEigenMMVIIPoly aPoly({-1,0,3}); - std::vector aPoly {-1,0,3}; - // bool hasArealRoot; - // tREAL8 aS = greatestRealRoot(hasArealRoot); - - StdOut() << "EEE:V= " << Eigen::poly_eval (aPoly,2) << "\n"; - - // polynomialsolver( internal::random(9,13)); - PolynomialSolver aSolver; - - aSolver.compute(aPoly); - // aSolver.roots(); - // FakeUseIt(aSolver); + cEigenPolynRoots aEPR(aPol,aTol,aNbMaxIter); + return aEPR.RealRoots(); } -#endif /* ************************************************************************ */ @@ -385,7 +348,8 @@ template template std::vector cPolynom::RealRoots(const Type & aTol,int ItMax) const { // StdOut() << "RealRootsRealRootsRealRootsRealRootsRealRootsRealRootsRealRootsRealRoots \n"; getchar(); - return V1RealRoots(mVCoeffs,aTol,ItMax); + // return V1RealRoots(mVCoeffs,aTol,ItMax); + return V2RealRoots(*this,aTol,ItMax); } // =========== others ========================= @@ -393,7 +357,6 @@ template std::vector cPolynom::RealRoots(const Type & a template size_t cPolynom::Degree() const { return mVCoeffs.size() - 1; } - template Type cPolynom::Value(const Type & aVal) const { Type aResult = 0.0; @@ -406,6 +369,22 @@ template Type cPolynom::Value(const Type & aVal) const return aResult; } +template cPtxd cPolynom::Value(const tCompl & aVal) const +{ + tCompl aResult (0.0,0.0); + tCompl aPowV (1.0,0.0); + for (const auto & aCoef : mVCoeffs) + { + aResult += aCoef * aPowV; + aPowV = aVal * aPowV; + } + return aResult; +} + + + + + template Type cPolynom::AbsValue(const Type & aVal) const { Type aResult = 0.0; @@ -567,7 +546,7 @@ template void TplBenchPolynome() Type aDif = RelativeSafeDifference(aVRootsGen[aK],aVRootsCalc[aK]); MMVII_INTERNAL_ASSERT_bench(aDif(); + // TplBenchPolynome(); // => with eigen , impossible to have always acceptable accuracy TplBenchPolynome(); TplBenchPolynome();