Skip to content

Commit

Permalink
Root Poly/Eigen :OK
Browse files Browse the repository at this point in the history
deseilligny committed Dec 23, 2024
1 parent 81c6749 commit d9a18a3
Showing 4 changed files with 158 additions and 161 deletions.
5 changes: 4 additions & 1 deletion MMVII/include/MMVII_nums.h
Original file line number Diff line number Diff line change
@@ -993,6 +993,7 @@ template <class Type> class cPolynom
{
public :
typedef std::vector<Type> tCoeffs;
typedef cPtxd<Type,2> tCompl;
cPolynom(const tCoeffs &);
cPolynom(const cPolynom &);
cPolynom(size_t aDegre);
@@ -1007,10 +1008,12 @@ template <class Type> class cPolynom
static cPolynom<Type> RandomPolyg(std::vector<Type> & 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<Type> operator * (const cPolynom<Type> & aP2) const;
cPolynom<Type> operator + (const cPolynom<Type> & aP2) const;
cPolynom<Type> operator - (const cPolynom<Type> & aP2) const;
8 changes: 6 additions & 2 deletions MMVII/src/Matrix/cEigenEigenDecompos.cpp
Original file line number Diff line number Diff line change
@@ -93,9 +93,13 @@ void Bench_EigenDecompos(cParamExeBench & aParam)
Tpl_Bench_EigenDecompos<tREAL8>(aParam);
}

#define INSTANTIATE_EIGEN_DECOMP(TYPE)\
template class cResulEigenDecomp<TYPE>;\
template class cDenseMatrix<TYPE>;

template class cResulEigenDecomp<tREAL8>;
template class cDenseMatrix<tREAL8>;
INSTANTIATE_EIGEN_DECOMP(tREAL4);
INSTANTIATE_EIGEN_DECOMP(tREAL8);
INSTANTIATE_EIGEN_DECOMP(tREAL16);


};
11 changes: 11 additions & 0 deletions MMVII/src/Serial/ElemStrToVal.cpp
Original file line number Diff line number Diff line change
@@ -1395,6 +1395,17 @@ std::string FixDigToStr(double aSignedVal,int aNbBef,int aNbAfter)
return aBuf;
}

// ================ double ==============================================

template <> std::string cStrIO<tREAL16>::ToStr(const tREAL16 & aD)
{
return cStrIO<tREAL8>::ToStr((tREAL8) (aD));
}

template <> std::string cStrIO<tREAL4>::ToStr(const tREAL4 & aD)
{
return cStrIO<tREAL8>::ToStr((tREAL8) (aD));
}


// ================ std::string ==============================================
295 changes: 137 additions & 158 deletions MMVII/src/UtiMaths/Polynoms.cpp
Original file line number Diff line number Diff line change
@@ -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 <Eigen/Dense>
// #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 <iostream>
//#include <vector>

using namespace Eigen;
using namespace std;

namespace MMVII
{
@@ -36,73 +21,136 @@ namespace MMVII
template <class Type> class cEigenPolynRoots
{
public :
cEigenPolynRoots(const cPolynom<Type> &) ;
const std::vector<Type> & RealRoots() const;
const VectorXcd & ComplexRoots() const;
bool RootIsReal(const std::complex<tREAL8> &,std::string * sayWhy=nullptr);
const MatrixXd & CompM() const;

typedef cDenseMatrix<Type> tMatComp;
typedef cPtxd<Type,2> tCompl;

cEigenPolynRoots(const cPolynom<Type> &,Type aEps,int aNbIterMax) ;
bool RootIsReal(const tCompl & ,std::string * sayWhy=nullptr);

const std::vector<Type> & RealRoots() const {return mRR;}
const std::vector<tCompl> & 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<Type> mPol;
cPolynom<Type> 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<Type> mRR; ///< Real roots
std::vector<tCompl> mCR; ///< Real roots

};


template <class Type> cEigenPolynRoots<Type>::cEigenPolynRoots(const cPolynom<Type> & aPol) :
template <class Type> cEigenPolynRoots<Type>::cEigenPolynRoots(const cPolynom<Type> & 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<MatrixXd> aSolver(mCompM);
mCEV = aSolver.eigenvalues() ;
cResulEigenDecomp<Type> 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 <class Type> const std::vector<Type> & cEigenPolynRoots<Type>::RealRoots() const {return mRR;}
template <class Type> cPtxd<Type,2> cEigenPolynRoots<Type>::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<aNbIter ; aKIt++)
{
tCompl aDeriv = mDPol.Value(aLastVal);
if (IsNotNull(aDeriv))
{
tCompl aNewVal = aLastVal - aLastEval/aDeriv;
tCompl aNewEval = mPol.Value(aNewVal);
Type aNewSqN2 = SqN2(aNewEval);
if (aNewSqN2 < aSqEps)
return aNewVal;
else if (aNewSqN2<aLastSqN2)
{

aLastVal = aNewVal;
aLastEval = aNewEval;
aLastSqN2 = aNewSqN2;
}
else
return aLastVal;
}
else
return aLastVal;
}
return aLastVal;
}

template <class Type> const VectorXcd & cEigenPolynRoots<Type>::ComplexRoots() const {return mCEV;}
// tCompl Refine(const tCompl & aV0,int aNbIter=5);

template <class Type> const MatrixXd & cEigenPolynRoots<Type>::CompM() const {return mCompM;}


/** Also the question seems pretty basic, it becomes more complicated due to numericall approximation */

template <class Type> bool cEigenPolynRoots<Type>::RootIsReal(const std::complex<tREAL8> & aC,std::string * sayWhy)
template <> tREAL4 cEigenPolynRoots<tREAL4>::PolRelAccuracy() {return 1e-5;}
template <> tREAL8 cEigenPolynRoots<tREAL8>::PolRelAccuracy() {return 1e-9;}
template <> tREAL16 cEigenPolynRoots<tREAL16>::PolRelAccuracy(){return 1e-11;}

template <> tREAL4 cEigenPolynRoots<tREAL4>::ComplRelAccuracy() {return 1e-8;}
template <> tREAL8 cEigenPolynRoots<tREAL8>::ComplRelAccuracy() {return 1e-8;}
template <> tREAL16 cEigenPolynRoots<tREAL16>::ComplRelAccuracy() {return 1e-8;}

template <> tREAL4 cEigenPolynRoots<tREAL4>::PolAbsAccuracy() {return 0.01;}
template <> tREAL8 cEigenPolynRoots<tREAL8>::PolAbsAccuracy() {return 1e-5;}
template <> tREAL16 cEigenPolynRoots<tREAL16>::PolAbsAccuracy() {return 1e-7;}
//static Type ComplRelAccuracy() ;
//static Type ComplAbsAccuracy() ;


template <class Type> bool cEigenPolynRoots<Type>::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,28 +159,28 @@ template <class Type> bool cEigenPolynRoots<Type>::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));
return false;
}

// [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 <class Type> bool cEigenPolynRoots<Type>::RootIsReal(const std::complex
return true;
}


template class cEigenPolynRoots<tREAL4>;
template class cEigenPolynRoots<tREAL8>;
template class cEigenPolynRoots<tREAL16>;


template<class Type> void My_Roots(const cPolynom<Type> & aPol1)
{
}
template <> void My_Roots(const cPolynom<tREAL8> & 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<tREAL8> & aPol1)
cAutoTimerSegm aTimeEigen(GlobAppTS(),"Eigen");
for (int aK=0 ; aK<aNb ; aK++)
{
cEigenPolynRoots<tREAL8> aEPR(aPol1);
cEigenPolynRoots<Type> aEPR(aPol1,1e-3,10);
}

cAutoTimerSegm aTimeV1(GlobAppTS(),"V1");
@@ -172,12 +225,13 @@ template <> void My_Roots(const cPolynom<tREAL8> & aPol1)
}
cAutoTimerSegm aTimeOthers(GlobAppTS(),"Others");

cEigenPolynRoots<tREAL8> aEPR(aPol1);
cEigenPolynRoots<Type> 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<tREAL8> & 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<double> coeffs = {-1,1,-1,1};
}


#if (0)

/*
MatrixXd createCompanionMatrix(const vector<double>& 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<double> coeffs = {1, -6, 11, -6};
vector<double> 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<MatrixXd> 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 <class Type> std::vector<Type> V2RealRoots(const cPolynom<Type> & aPol, Type aTol,int aNbMaxIter)
{
// cEigenMMVIIPoly aPoly({-1,0,3});
std::vector<tREAL8> aPoly {-1,0,3};
// bool hasArealRoot;
// tREAL8 aS = greatestRealRoot(hasArealRoot);

StdOut() << "EEE:V= " << Eigen::poly_eval (aPoly,2) << "\n";

// polynomialsolver<float,Dynamic>( internal::random<int>(9,13));
PolynomialSolver<tREAL8,10> aSolver;

aSolver.compute(aPoly);
// aSolver.roots();
// FakeUseIt(aSolver);
cEigenPolynRoots<Type> aEPR(aPol,aTol,aNbMaxIter);

return aEPR.RealRoots();
}
#endif


/* ************************************************************************ */
@@ -385,15 +348,15 @@ template <class Type>
template <class Type> std::vector<Type> cPolynom<Type>::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 =========================

template <class Type> size_t cPolynom<Type>::Degree() const { return mVCoeffs.size() - 1; }



template <class Type> Type cPolynom<Type>::Value(const Type & aVal) const
{
Type aResult = 0.0;
@@ -406,6 +369,22 @@ template <class Type> Type cPolynom<Type>::Value(const Type & aVal) const
return aResult;
}

template <class Type> cPtxd<Type,2> cPolynom<Type>::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 <class Type> Type cPolynom<Type>::AbsValue(const Type & aVal) const
{
Type aResult = 0.0;
@@ -567,7 +546,7 @@ template<class Type> void TplBenchPolynome()
Type aDif = RelativeSafeDifference(aVRootsGen[aK],aVRootsCalc[aK]);
MMVII_INTERNAL_ASSERT_bench(aDif<aEps,"roots size check");
}
My_Roots(aPol);
// My_Roots(aPol);

}
}
@@ -578,7 +557,7 @@ void BenchPolynome(cParamExeBench & aParam)

// TestPolynEigen();

TplBenchPolynome<tREAL4>();
// TplBenchPolynome<tREAL4>(); // => with eigen , impossible to have always acceptable accuracy
TplBenchPolynome<tREAL8>();
TplBenchPolynome<tREAL16>();

0 comments on commit d9a18a3

Please sign in to comment.