From 4668474cebfbead40522d719833a348a7d12aa93 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Mon, 16 Dec 2024 22:30:45 +0100 Subject: [PATCH 1/5] In supr dep to mmv1 --- MMVII/src/Appli/cMMVII_CalcSet.cpp | 13 +- MMVII/src/Bench/BenchGlob.cpp | 77 --------- MMVII/src/MMV1/ExportSensorMMV1.cpp | 23 ++- MMVII/src/MMV1/FluxPtsMMV1.cpp | 16 +- MMVII/src/MMV1/GenCodeTestCam.cpp | 197 +---------------------- MMVII/src/MMV1/ImageFilterMMV1.cpp | 55 +++++-- MMVII/src/MMV1/XmlMMV1.cpp | 63 +++++++- MMVII/src/Mesh/Ply.cpp | 4 + MMVII/src/Sensors/cConvCalib.cpp | 7 + MMVII/src/TestLibsExtern/TestJetsCam.cpp | 28 +--- 10 files changed, 148 insertions(+), 335 deletions(-) diff --git a/MMVII/src/Appli/cMMVII_CalcSet.cpp b/MMVII/src/Appli/cMMVII_CalcSet.cpp index 6ed470e4ae..4c3da5f22b 100755 --- a/MMVII/src/Appli/cMMVII_CalcSet.cpp +++ b/MMVII/src/Appli/cMMVII_CalcSet.cpp @@ -1,3 +1,5 @@ +#define WITH_MMV1_FUNCTION false + #include "MMVII_2Include_Serial_Tpl.h" #include "MMVII_DeclareAllCmd.h" #include "MMVII_Sensor.h" @@ -90,7 +92,10 @@ static void OneBenchEditSet const std::string & ExpSet // Expect set ) { - // StdOut() << "OneBenchEditSet " << anOp << std::endl; +#if (!WITH_MMV1_FUNCTION) + aNumAskedOut = 0; + aRealNumOut = 2; +#endif cMMVII_Appli & anAp = cMMVII_Appli::CurrentAppli(); std::string aDirI = anAp.InputDirTestMMVII() + "Files/" ; @@ -196,15 +201,15 @@ int cAppli_EditSet::ExecuteBench(cParamExeBench & aParam) aK, // Change test condition : Use or Not Dir Project "+=", // operator for mpdi false, // If true, Previous OutPut is moved on input, else Input is purged at end of process - ".*txt", // Pattern of used files + "F.*txt", // Pattern of used files 0, // Required num version 2, // Real Num Version 10, // Number of element expected (become obsolet with expected set) "", // Interval modifying the pattern if != "" C09 // Ground truth, what the string should be ); // - OneBenchEditSet(aK,"+=",false,".*txt" ,1,1,10,"",C09); - OneBenchEditSet(aK,"+=",false,".*txt" ,2,2,10,"",C09); + OneBenchEditSet(aK,"+=",false,"F.*txt" ,1,1,10,"",C09); + OneBenchEditSet(aK,"+=",false,"F.*txt" ,2,2,10,"",C09); OneBenchEditSet(aK,"+=",false,"F[02468].txt",2,2,5,"","02468"); // here we init from previous OneBenchEditSet(aK,"+=",true ,"F[3-5].txt" ,2,2,7,"","0234568"); // 0234568 diff --git a/MMVII/src/Bench/BenchGlob.cpp b/MMVII/src/Bench/BenchGlob.cpp index 3868e6ca60..0db0872455 100755 --- a/MMVII/src/Bench/BenchGlob.cpp +++ b/MMVII/src/Bench/BenchGlob.cpp @@ -1056,88 +1056,11 @@ int cAppli_MPDTest::Exe() return EXIT_SUCCESS; } TTT (); -#if 1 { StdOut() << "T0:" << cName2Calc::CalcFromName("toto",10,true) << std::endl; StdOut() << "T1:" << cName2Calc::CalcFromName("EqDist_Dist_Rad3_Dec1_XY1",10) << std::endl; return EXIT_SUCCESS; } -#else - if (mMMV1_GenCodeTestCam) - { - //StdOut() << "kkk=[" << mTopDirMMVII <<"]" << std::endl; - MMV1_GenerateCodeTestCam(); - return EXIT_SUCCESS; - } - { - // Si on le met a 10h => reveil a 6h20 - double t = 8.0; - //sleep(3600.0 * t); - sleepcp(3600.0 * t * 1000); - std::string aName= "/home/mpd/Bureau/Perso1/Musik/Bach/bach-goldberg-variations-bwv-988-glenn-gould-1981.mp3"; - aName = "cvlc " + aName; - StdOut() << system(aName.c_str()) << std::endl;; - } - { - cPt3dr * anAdr = nullptr; - StdOut () << "ADDDDDr " << anAdr << "\n"; - StdOut () << "ADDDDDrx " << &(anAdr->x()) << "\n"; - StdOut () << "ADDDDDry " << &(anAdr->y()) << "\n"; - StdOut () << "ADDDDDrz " << &(anAdr->z()) << "\n"; - ShowAdr(anAdr->y()); - } - { - double aV= 3.3333; - printf("VVVVV=%05.2f\n",aV); - - } - if ((UN>DEUX) && PrintAndTrue("aaaa")) - { - PrintAndTrue("bbbb"); - } - PrintAndTrue("ccccc"); - cRotation3D::RandomRot(); - -/* - cSparseVect aSV; - for (const auto & aP : aSV) - { - StdOut() << aP.mI << std::endl; - } -*/ - -/* - cIm2D aIm(cPt2di(3,3)); - aIm.DIm().SetV(cPt2di(0,0),13); - // aIm.DIm().SetV(cPt2di(0,0),1000); - // aIm.DIm().SetV(cPt2di(-1,0),1); - // new cIm2D(cPt2di(3,3)); - cDataIm2D & aDIm = aIm.DIm(); - tU_INT1* aPtr = aDIm.RawDataLin(); - StdOut() << "aIm=" << int(aPtr[0]) << std::endl; - aPtr[0] = 14; - StdOut() << "aIm=" << (int)aDIm.GetV(cPt2di(0,0)) << std::endl; - // aPtr[-1] = 0; -*/ - -/* - TestVectBool(); - cMMVII_Ofs aOs1("toto1.txt"); - cMMVII_Ofs aOs2("toto2.txt"); - - - cMultipleOfs amOs; // (aOs1.Ofs(),aOs2.Ofs()); - amOs.Add(aOs1.Ofs()); - amOs.Add(aOs2.Ofs()); - amOs << "1+1=" << 1+1 << "\n"; - cMMVII_Ofs aFile("toto.txt"); - std::ostream & anOFS = aFile.Ofs(); - anOFS << "TEST OFFFSSSSSSSSSSSS\n"; -*/ - - return EXIT_SUCCESS; - -#endif } tMMVII_UnikPApli Alloc_MPDTest(const std::vector & aVArgs,const cSpecMMVII_Appli & aSpec) diff --git a/MMVII/src/MMV1/ExportSensorMMV1.cpp b/MMVII/src/MMV1/ExportSensorMMV1.cpp index d50f9a6646..f751bb73c0 100755 --- a/MMVII/src/MMV1/ExportSensorMMV1.cpp +++ b/MMVII/src/MMV1/ExportSensorMMV1.cpp @@ -1,9 +1,15 @@ +#define WITH_MMV1_FUNCTION false + +#if (WITH_MMV1_FUNCTION) #include "V1VII.h" -#include "MMVII_util.h" +#endif + +#include "MMVII_MMV1Compat.h" namespace MMVII { +#if (WITH_MMV1_FUNCTION) /* ************************************************* */ /* */ /* cExportV1StenopeCalInterne */ @@ -154,5 +160,20 @@ cExportV1StenopeCalInterne::cExportV1StenopeCalInterne if ((aNbLayer>0) && (aNbPointPerDim>0)) DoCorresp(mCorresp,aCamV1,aNbPointPerDim,aNbLayer,aDS); } +#else +cExportV1StenopeCalInterne::cExportV1StenopeCalInterne +( + bool isForCalib, + const std::string& aFile, + int aNbPointPerDim, + int aNbLayer, + tREAL8 aDS, + const std::string& aFileInterneCalib +) : + mPose (cIsometry3D::Identity()) +{ + MMVII_INTERNAL_ERROR("Reading of MMV1 Camera is deprecated"); +} +#endif }; diff --git a/MMVII/src/MMV1/FluxPtsMMV1.cpp b/MMVII/src/MMV1/FluxPtsMMV1.cpp index d0a309d409..5dee36e857 100755 --- a/MMVII/src/MMV1/FluxPtsMMV1.cpp +++ b/MMVII/src/MMV1/FluxPtsMMV1.cpp @@ -1,3 +1,5 @@ +#define WITH_MMV1_FUNCTION false + #include "V1VII.h" #include "MMVII_Geom2D.h" @@ -35,6 +37,10 @@ Neighbourhood DiscNeich(const tREAL8 & aRay) void FluxToV2Points(tResFlux & aRes,Flux_Pts aFlux) { +//StdOut() << "FluxToV2PointsFluxToV2PointsFluxToV2PointsuuuuuuuuuuuuFluxToV2PointsFluxToV2PointsFluxToV2PointsFluxToV2PointsFluxToV2Points\n"; +//getchar(); + + aRes.clear(); Liste_Pts aLV1(2); @@ -80,7 +86,6 @@ void GetPts_Ellipse(tResFlux & aRes,const cPt2dr & aC,double aRayA,double aRayB } - void GetPts_Line(tResFlux & aRes,const cPt2dr & aP1,const cPt2dr &aP2,tREAL8 aDil) { Flux_Pts aFlux = line(ToMMV1(ToI(aP1)),ToMMV1(ToI(aP2))); @@ -95,15 +100,6 @@ void GetPts_Line(tResFlux & aRes,const cPt2dr & aP1,const cPt2dr &aP2) } -/* -void Txxxxxxxxx() -{ - Flux_Pts disc(Pt2dr,REAL,bool front = true); - - Neighbourhood -} -*/ - }; diff --git a/MMVII/src/MMV1/GenCodeTestCam.cpp b/MMVII/src/MMV1/GenCodeTestCam.cpp index c86bbc4e26..8e1782da63 100755 --- a/MMVII/src/MMV1/GenCodeTestCam.cpp +++ b/MMVII/src/MMV1/GenCodeTestCam.cpp @@ -1,197 +1,2 @@ -#include "V1VII.h" -#include "cGeneratedCodeTestCam.h" -#include "MMVII_Derivatives.h" +// Empty file - -namespace MMVII -{ -/** \file GenCodeTestCam.cpp - \brief Make benchmark on generated formal code - - This file contain test to generate formal class, and use them - to make benchmark between jets and computed analyticall - -*/ - - -class cGenCodeTestCam -{ - public : - cGenCodeTestCam(); - void Generate(); - void Init(cGeneratedCodeTestCam &aGCTC); - - private : - cIncListInterv mLInterv; - cSetEqFormelles mSet; - AllocateurDInconnues & mAlloc; - cMatr_Etat_PhgrF mRot0; - cP2d_Etat_PhgrF mMesIm; - - cP3dFormel mPTer; - cP3dFormel mCCam; - cP3dFormel mOmega; - cP2dFormel mCDist; - cValFormel mK2; - cValFormel mK4; - cValFormel mK6; - cP2dFormel mPP; - cValFormel mFoc; -}; - -class cMMV1_TestCam : public cInterfaceTestCam -{ - public : - void InitFromParams(const std::vector &) override; - void Compute(std::vector & Vals,std::vector > & ) override; - void Compute(int aNb) override; - cMMV1_TestCam (); - private : - cGenCodeTestCam mGenerator; - cGeneratedCodeTestCam mCamGen; -}; - -/* *********************************************** */ -/* */ -/* cMMV1_TestCam */ -/* */ -/* *********************************************** */ - -void cMMV1_TestCam::InitFromParams(const std::vector & aVals) -{ - mCamGen.SetCoordCur(aVals.data()); -// MesIm_x - *(mCamGen.AdrVarLocFromString("MesIm_x")) = 0.0; - *(mCamGen.AdrVarLocFromString("MesIm_y")) = 0.0; - - *(mCamGen.AdrVarLocFromString("R0_0_0")) = 1.0; - *(mCamGen.AdrVarLocFromString("R0_1_0")) = 0.0; - *(mCamGen.AdrVarLocFromString("R0_2_0")) = 0.0; - - *(mCamGen.AdrVarLocFromString("R0_0_1")) = 0.0; - *(mCamGen.AdrVarLocFromString("R0_1_1")) = 1.0; - *(mCamGen.AdrVarLocFromString("R0_2_1")) = 0.0; - - *(mCamGen.AdrVarLocFromString("R0_0_2")) = 0.0; - *(mCamGen.AdrVarLocFromString("R0_1_2")) = 0.0; - *(mCamGen.AdrVarLocFromString("R0_2_2")) = 1.0; - -} - -cMMV1_TestCam::cMMV1_TestCam () -{ - mGenerator.Init(mCamGen); -} - -void cMMV1_TestCam::Compute(std::vector & aVals,std::vector > & aDeriv) -{ - mCamGen.ComputeValDeriv(); - aVals = mCamGen.ValSsVerif(); - aDeriv = mCamGen.CompDerSsVerif(); -} - -void cMMV1_TestCam::Compute(int aNb) -{ - for (int aK=0 ; aK aPCam = mRot0.Mat() * (mPTer.FPt()-mCCam.FPt()); - aPCam = aPCam + (mOmega.FPt() ^ aPCam); - - Pt2d aPi (aPCam.x/aPCam.z,aPCam.y/aPCam.z); - Pt2d aCdPi = aPi-mCDist.FPt(); - Fonc_Num aRho2 = Square(aCdPi.x) + Square(aCdPi.y); - Fonc_Num aRho4 = Square(aRho2); - Fonc_Num aRho6 = Cube(aRho2); - - Fonc_Num aCoeffDist = mK2.FVal()*aRho2 + mK4.FVal()*aRho4 + mK6.FVal()*aRho6; - Pt2d aPDist = aPi + Pt2d(aCdPi.x*aCoeffDist,aCdPi.y*aCoeffDist); - Pt2d aPProj = mPP.FPt() + Pt2d(aPDist.x*mFoc.FVal(),aPDist.y*mFoc.FVal()); - Pt2d aResidu = aPProj-mMesIm.PtF(); - - std::string aDirApp = cMMVII_Appli::CurrentAppli().TopDirMMVII() ; - cElCompileFN::DoEverything - ( - aDirApp + "src/MMV1/", // Directory ou est localise le code genere - "cGeneratedCodeTestCam", // donne les noms de .cpp et .h de classe - // std::vector ({aPi.x,aPi.y}), // expressions formelles - std::vector ({aResidu.x,aResidu.y}), // expressions formelles - mLInterv // intervalle de reference - ); - - -} - - -void MMV1_GenerateCodeTestCam() -{ - if (1) - { - cGenCodeTestCam aGCTC; - aGCTC.Generate(); - } - if (0) - { - double aVVals[30]; - cGenCodeTestCam aGen; - cGeneratedCodeTestCam aCamGen; - aGen.Init(aCamGen); - aCamGen.SetCoordCur(aVVals); - aCamGen.ComputeValDeriv(); - - } -} - - - -}; diff --git a/MMVII/src/MMV1/ImageFilterMMV1.cpp b/MMVII/src/MMV1/ImageFilterMMV1.cpp index 94302bf1a7..5ba74495f6 100755 --- a/MMVII/src/MMV1/ImageFilterMMV1.cpp +++ b/MMVII/src/MMV1/ImageFilterMMV1.cpp @@ -1,6 +1,6 @@ -#define WITH_MMV1_BENCH false +#define WITH_MMV1_FUNCTION false -#if (WITH_MMV1_BENCH) +#if (WITH_MMV1_FUNCTION) #include "V1VII.h" #endif @@ -61,10 +61,17 @@ void Bench_LabMaj() } -#if (WITH_MMV1_BENCH) +#if (WITH_MMV1_FUNCTION) + +/* ************************************************* */ +/* */ +/* UNUSED */ +/* */ +/* ************************************************* */ + + // Stuffs not implemanted because not used in MMVII 4 now -// Not implemanted/used in MMVII 4 now void MMV1_MakeImageDist(cIm2D aImIn,const std::string & aNameChamfer) { const Chamfer & aChamf = Chamfer::ChamferFromName(aNameChamfer); @@ -97,8 +104,16 @@ void ExportHomMMV1(const std::string & aIm1,const std::string & aIm2,const std:: ExportHomMMV1(aIm1,aIm2,SH,aVR); } +/* ************************************************* */ +/* */ +/* DEPRECATED V1 SOLUTION */ +/* */ +/* ************************************************* */ - //================= V1 Solution ================== + // implemantation of methods that used to do be done unsing MMV1 + // Activated with WITH_MMV1_FUNCTION , to check correctness + + //================= Curvature of tangent line -> Courb Tgt ================== template cIm2D MMV1_CourbTgt(cIm2D aImIn) { @@ -117,6 +132,8 @@ template void MMV1_SelfCourbTgt(cIm2D aImIn) ELISE_COPY(aV1In.all_pts(),courb_tgt(aV1In.in_proj(),0.5),aV1In.out()); } + //================= Laplacian ================== + cIm2D MMV1_Lapl(cIm2D aImIn) { cIm2D aRes(aImIn.DIm().Sz()); @@ -130,19 +147,9 @@ cIm2D MMV1_Lapl(cIm2D aImIn) } - //================= Global Solution ================== -template void BenchImFilterV1V2(cIm2D anI1,cIm2D anI2,tREAL8 aEps) -{ - - tREAL8 aRD = anI1.DIm().SafeMaxRelDif(anI2.DIm(),1e-2); - if ( aRD > aEps ) - { - MMVII_INTERNAL_ASSERT_bench( false,"BenchImFilterV1V2 RD="+ToStr(aRD)); - } - -} + //================= Extremal value of image ================== template std::pair MMV1_ValExtre(cIm2D aImIn) { @@ -151,6 +158,7 @@ template std::pair MMV1_ValExtre(cIm2D aImIn) ELISE_COPY(aV1In.all_pts(),aV1In.in(),VMin(aVMin)|VMax(aVMax)); return std::pair((Type)aVMin,(Type)aVMax); } + //================= Average of absolute value ================== template double MMV1_MoyAbs(cIm2D aImIn) { @@ -160,6 +168,8 @@ template double MMV1_MoyAbs(cIm2D aImIn) return aSom[0] / aSom[1]; } + //================= Deriche gradient ================== + template void MMV1_ComputeDeriche(cImGrad & aResGrad,const cDataIm2D & aImIn,double aAlpha) { auto aV1In = cMMV1_Conv::ImToMMV1(aImIn); @@ -180,6 +190,8 @@ template cImGrad MMV1_Deriche(const cDataIm2D & aImIn,do return aResGrad; } + //================= Make a visualisable image with adapative dyn ================== + void MMV1_MakeStdIm8BIts(cIm2D aImIn,const std::string& aName) { Im2D aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); @@ -204,6 +216,17 @@ void MMV1_MakeStdIm8BIts(cIm2D aImIn,const std::string& aName) } + //============== BENCH with old MMV1 implementation and new one =============== + +template void BenchImFilterV1V2(cIm2D anI1,cIm2D anI2,tREAL8 aEps) +{ + tREAL8 aRD = anI1.DIm().SafeMaxRelDif(anI2.DIm(),1e-2); + if ( aRD > aEps ) + { + MMVII_INTERNAL_ASSERT_bench( false,"BenchImFilterV1V2 RD="+ToStr(aRD)); + } +} + template void Tpl_BenchImFilterV1V2(const cPt2di& aSz) { diff --git a/MMVII/src/MMV1/XmlMMV1.cpp b/MMVII/src/MMV1/XmlMMV1.cpp index 909398f624..74e2ce7fbb 100755 --- a/MMVII/src/MMV1/XmlMMV1.cpp +++ b/MMVII/src/MMV1/XmlMMV1.cpp @@ -1,16 +1,21 @@ +#define WITH_MMV1_FUNCTION false + +#if (WITH_MMV1_FUNCTION) #include "V1VII.h" +#endif #include "MMVII_Mappings.h" #include "MMVII_AimeTieP.h" +#if (WITH_MMV1_FUNCTION) #include "../src/uti_image/NewRechPH/cParamNewRechPH.h" -// #include "../CalcDescriptPCar/AimeTieP.h" - #include "im_tpl/cPtOfCorrel.h" #include "algo_geom/qdt.h" #include "algo_geom/qdt_implem.h" #include "ext_stl/heap.h" +#endif + namespace MMVII @@ -22,7 +27,11 @@ namespace MMVII /* */ /* *********************************** */ - +#if (WITH_MMV1_FUNCTION) +/* + This class is an interface to use in MMV2, the 3d masq produce by SaisieMasqQT of MMV1, + using the MMV1 library +*/ class cMasq_MMV1asBoundeSet : public cDataBoundedSet { public : @@ -54,10 +63,18 @@ cDataBoundedSet * MMV1_Masq(const cBox3dr &aBox,const std::string & a { return new cMasq_MMV1asBoundeSet(aBox,aNameFile); } +#else +cDataBoundedSet * MMV1_Masq(const cBox3dr &aBox,const std::string & aNameFile) +{ + MMVII_INTERNAL_ERROR("Handling MMV1 3d-Masq is deprecated"); + return nullptr; +} +#endif //============= tNameRel ==================== +#if (WITH_MMV1_FUNCTION) void TestTimeV1V2() { @@ -595,13 +612,47 @@ template void cImplem_ExportAimeTiep::FiltrageSpatialPts() } } +template class cImplem_ExportAimeTiep; +template class cImplem_ExportAimeTiep; +#else +template cInterf_ExportAimeTiep::~cInterf_ExportAimeTiep() +{ +} + +template cInterf_ExportAimeTiep * cInterf_ExportAimeTiep::Alloc(const cPt2di& aSzIm0,bool IsMax,eTyPyrTieP ATypePt,const std::string & aName,bool ForInspect,const cGP_Params & aParam ) +{ + + MMVII_INTERNAL_ERROR("Creating cInterf_ExportAimeTiep is deprecated"); + return nullptr; +} + +template<> void MMv1_SaveInFile(const tNameRel & aSet,const std::string & aName) +{ + MMVII_INTERNAL_ERROR("Handling writing of MMV1 Relations is deprecated"); +} +template<> void MMv1_SaveInFile(const tNameSet & aSet,const std::string & aName) +{ + MMVII_INTERNAL_ERROR("Handling writing MMV1 Sets is deprecated"); +} +tNameSet MMV1InitSet(const std::string & aName) +{ + MMVII_INTERNAL_ERROR("Handling reaing of MMV1 Sets is deprecated"); + return tNameSet(); +} + +tNameRel MMV1InitRel(const std::string & aName) +{ + MMVII_INTERNAL_ERROR("Handling reaing of MMV1 Rels is deprecated"); + return tNameRel(); +} +#endif // ============ INSTANTIATION ====================== template class cInterf_ExportAimeTiep; template class cInterf_ExportAimeTiep; -template class cImplem_ExportAimeTiep; -template class cImplem_ExportAimeTiep; - }; + + + diff --git a/MMVII/src/Mesh/Ply.cpp b/MMVII/src/Mesh/Ply.cpp index e03e9b60c2..4be24d3899 100644 --- a/MMVII/src/Mesh/Ply.cpp +++ b/MMVII/src/Mesh/Ply.cpp @@ -2,6 +2,8 @@ #include "MMVII_Geom3D.h" #include "MMVII_Mappings.h" +#define WITH_MMV1_FUNCTION false + namespace MMVII { @@ -2114,12 +2116,14 @@ template void cTriangulation3D::Bench() MMVII_INTERNAL_ASSERT_bench(aTriBin3D.HeuristikAlmostEqual(aNewTriBin3D,1e-5,1e-5),"cTriangulation3D::Bench"); MMVII_INTERNAL_ASSERT_bench(aTriTxt3D.HeuristikAlmostEqual(aNewTriTxt3D,1e-5,1e-5),"cTriangulation3D::Bench"); +#if (WITH_MMV1_FUNCTION) cDataBoundedSet * aMasq= MMV1_Masq(aNewTriBin3D.BoxEngl().ToR(),aDirI+"AperiCloud_Basc_selectionInfo.xml"); aNewTriBin3D.Filter(*aMasq); aNewTriBin3D.WriteFile(aDirTmp+"MeshFilteredBin.ply",true); // BREAK_POINT("BENCHTRI"); delete aMasq; +#endif } void BenchPly(cParamExeBench & aParam) diff --git a/MMVII/src/Sensors/cConvCalib.cpp b/MMVII/src/Sensors/cConvCalib.cpp index 547a05cdcb..52dd2f5cb1 100644 --- a/MMVII/src/Sensors/cConvCalib.cpp +++ b/MMVII/src/Sensors/cConvCalib.cpp @@ -1,8 +1,12 @@ +#define WITH_MMV1_FUNCTION false + #include "MMVII_PCSens.h" #include "MMVII_MMV1Compat.h" #include "MMVII_DeclareCste.h" #include "MMVII_BundleAdj.h" + + /** \file cConvCalib.cpp testgit @@ -298,6 +302,9 @@ void BenchPoseImportV1(const std::string & aNameOriV1,double anAccuracy) void BenchCentralePerspective_ImportV1(cParamExeBench & aParam) { +#if (!WITH_MMV1_FUNCTION) + return; +#endif BenchPoseImportV1("Orientation-_DSC8385.tif.xml" ,1e-5); BenchPoseImportV1("Orientation-Img0937.tif.xml" ,1e-5); diff --git a/MMVII/src/TestLibsExtern/TestJetsCam.cpp b/MMVII/src/TestLibsExtern/TestJetsCam.cpp index 9a1c4b7a7a..90d29c2a77 100755 --- a/MMVII/src/TestLibsExtern/TestJetsCam.cpp +++ b/MMVII/src/TestLibsExtern/TestJetsCam.cpp @@ -284,8 +284,8 @@ void BenchJetsCam(cParamExeBench & aParam) for (int aK=0 ; aK aVParam = StdParamTestCam(1.0); - std::vector aValJet,aValV1,aVarValJet; - std::vector > aDerJet,aDerV1,aVarDerJet; + std::vector aValJet,aVarValJet; + std::vector > aDerJet,aVarDerJet; cJetsTestCam aCamJet; aCamJet.InitFromParams(aVParam); @@ -295,33 +295,11 @@ void BenchJetsCam(cParamExeBench & aParam) aVarCamJet.InitFromParams(aVParam); aVarCamJet.Compute(aVarValJet,aVarDerJet); - - std::shared_ptr aCamV1 (cInterfaceTestCam::AllocMMV1()); - aCamV1->InitFromParams(aVParam); - aCamV1->Compute(aValV1,aDerV1); - - for (int aKXY=0 ; aKXY<2 ; aKXY++) - { - double aDif = aValJet[aKXY]-aValV1[aKXY]; - MMVII_INTERNAL_ASSERT_bench(std::abs(aDif)<1e-5,"Jets/V1 vals"); - - aDif = aVarValJet[aKXY]-aValV1[aKXY]; - - // StdOut() << "DDD=" << aDif<< " " << aValJet[aKXY] << " " << aValV1[aKXY] << std::endl; - for (int aKD=0 ; aKDCompute(aNbTest); + // aCamV1->Compute(aNbTest); double aT2 = cMMVII_Appli::CurrentAppli().SecFromT0(); // Now we now that sparse jets are slow save some time double aT3 = aT2; From c9bb90c7ce79873e30f76bc8c80ea1926f2b56d0 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Mon, 16 Dec 2024 22:37:24 +0100 Subject: [PATCH 2/5] In rm V1 ... --- MMVII/src/MMV1/GenCodeTestCam.cpp | 2 - MMVII/src/MMV1/cGeneratedCodeTestCam.cpp | 542 ----------------------- 2 files changed, 544 deletions(-) delete mode 100755 MMVII/src/MMV1/GenCodeTestCam.cpp diff --git a/MMVII/src/MMV1/GenCodeTestCam.cpp b/MMVII/src/MMV1/GenCodeTestCam.cpp deleted file mode 100755 index 8e1782da63..0000000000 --- a/MMVII/src/MMV1/GenCodeTestCam.cpp +++ /dev/null @@ -1,2 +0,0 @@ -// Empty file - diff --git a/MMVII/src/MMV1/cGeneratedCodeTestCam.cpp b/MMVII/src/MMV1/cGeneratedCodeTestCam.cpp index 43478da79f..e69de29bb2 100755 --- a/MMVII/src/MMV1/cGeneratedCodeTestCam.cpp +++ b/MMVII/src/MMV1/cGeneratedCodeTestCam.cpp @@ -1,542 +0,0 @@ -// File Automatically generated by eLiSe -#include "StdAfx.h" -#include "cGeneratedCodeTestCam.h" - - -cGeneratedCodeTestCam::cGeneratedCodeTestCam(): - cElCompiledFonc(2) -{ - AddIntRef (cIncIntervale("CCam",3,6)); - AddIntRef (cIncIntervale("CD",9,11)); - AddIntRef (cIncIntervale("F",16,17)); - AddIntRef (cIncIntervale("K2",11,12)); - AddIntRef (cIncIntervale("K4",12,13)); - AddIntRef (cIncIntervale("K6",13,14)); - AddIntRef (cIncIntervale("PGr",0,3)); - AddIntRef (cIncIntervale("PP",14,16)); - AddIntRef (cIncIntervale("W",6,9)); - Close(false); -} - - - -void cGeneratedCodeTestCam::ComputeVal() -{ - double tmp0_ = mCompCoord[0]; - double tmp1_ = mCompCoord[3]; - double tmp2_ = tmp0_ - tmp1_; - double tmp3_ = mCompCoord[1]; - double tmp4_ = mCompCoord[4]; - double tmp5_ = tmp3_ - tmp4_; - double tmp6_ = mCompCoord[2]; - double tmp7_ = mCompCoord[5]; - double tmp8_ = tmp6_ - tmp7_; - double tmp9_ = mLocR0_0_2 * (tmp2_); - double tmp10_ = mLocR0_1_2 * (tmp5_); - double tmp11_ = tmp9_ + tmp10_; - double tmp12_ = mLocR0_2_2 * (tmp8_); - double tmp13_ = tmp11_ + tmp12_; - double tmp14_ = mLocR0_0_1 * (tmp2_); - double tmp15_ = mLocR0_1_1 * (tmp5_); - double tmp16_ = tmp14_ + tmp15_; - double tmp17_ = mLocR0_2_1 * (tmp8_); - double tmp18_ = tmp16_ + tmp17_; - double tmp19_ = mCompCoord[7]; - double tmp20_ = mLocR0_0_0 * (tmp2_); - double tmp21_ = mLocR0_1_0 * (tmp5_); - double tmp22_ = tmp20_ + tmp21_; - double tmp23_ = mLocR0_2_0 * (tmp8_); - double tmp24_ = tmp22_ + tmp23_; - double tmp25_ = tmp19_ * (tmp13_); - double tmp26_ = mCompCoord[8]; - double tmp27_ = tmp26_ * (tmp18_); - double tmp28_ = tmp25_ - tmp27_; - double tmp29_ = tmp24_ + tmp28_; - double tmp30_ = mCompCoord[6]; - double tmp31_ = tmp30_ * (tmp18_); - double tmp32_ = tmp19_ * (tmp24_); - double tmp33_ = tmp31_ - tmp32_; - double tmp34_ = tmp13_ + tmp33_; - double tmp35_ = (tmp29_) / (tmp34_); - double tmp36_ = mCompCoord[9]; - double tmp37_ = tmp35_ - tmp36_; - double tmp38_ = ElSquare(tmp37_); - double tmp39_ = tmp26_ * (tmp24_); - double tmp40_ = tmp30_ * (tmp13_); - double tmp41_ = tmp39_ - tmp40_; - double tmp42_ = tmp18_ + tmp41_; - double tmp43_ = (tmp42_) / (tmp34_); - double tmp44_ = mCompCoord[10]; - double tmp45_ = tmp43_ - tmp44_; - double tmp46_ = ElSquare(tmp45_); - double tmp47_ = tmp38_ + tmp46_; - double tmp48_ = mCompCoord[11]; - double tmp49_ = tmp48_ * (tmp47_); - double tmp50_ = mCompCoord[12]; - double tmp51_ = ElSquare(tmp47_); - double tmp52_ = tmp50_ * tmp51_; - double tmp53_ = tmp49_ + tmp52_; - double tmp54_ = mCompCoord[13]; - double tmp55_ = VCube(tmp47_); - double tmp56_ = tmp54_ * tmp55_; - double tmp57_ = tmp53_ + tmp56_; - double tmp58_ = mCompCoord[16]; - - mVal[0] = (mCompCoord[14] + (tmp35_ + (tmp37_) * (tmp57_)) * tmp58_) - mLocMesIm_x; - - mVal[1] = (mCompCoord[15] + (tmp43_ + (tmp45_) * (tmp57_)) * tmp58_) - mLocMesIm_y; - -} - - -void cGeneratedCodeTestCam::ComputeValDeriv() -{ - double tmp0_ = mCompCoord[0]; - double tmp1_ = mCompCoord[3]; - double tmp2_ = tmp0_ - tmp1_; - double tmp3_ = mCompCoord[1]; - double tmp4_ = mCompCoord[4]; - double tmp5_ = tmp3_ - tmp4_; - double tmp6_ = mCompCoord[2]; - double tmp7_ = mCompCoord[5]; - double tmp8_ = tmp6_ - tmp7_; - double tmp9_ = mLocR0_0_2 * (tmp2_); - double tmp10_ = mLocR0_1_2 * (tmp5_); - double tmp11_ = tmp9_ + tmp10_; - double tmp12_ = mLocR0_2_2 * (tmp8_); - double tmp13_ = tmp11_ + tmp12_; - double tmp14_ = mLocR0_0_1 * (tmp2_); - double tmp15_ = mLocR0_1_1 * (tmp5_); - double tmp16_ = tmp14_ + tmp15_; - double tmp17_ = mLocR0_2_1 * (tmp8_); - double tmp18_ = tmp16_ + tmp17_; - double tmp19_ = mCompCoord[7]; - double tmp20_ = mLocR0_0_0 * (tmp2_); - double tmp21_ = mLocR0_1_0 * (tmp5_); - double tmp22_ = tmp20_ + tmp21_; - double tmp23_ = mLocR0_2_0 * (tmp8_); - double tmp24_ = tmp22_ + tmp23_; - double tmp25_ = tmp19_ * (tmp13_); - double tmp26_ = mCompCoord[8]; - double tmp27_ = tmp26_ * (tmp18_); - double tmp28_ = tmp25_ - tmp27_; - double tmp29_ = tmp24_ + tmp28_; - double tmp30_ = mCompCoord[6]; - double tmp31_ = tmp30_ * (tmp18_); - double tmp32_ = tmp19_ * (tmp24_); - double tmp33_ = tmp31_ - tmp32_; - double tmp34_ = tmp13_ + tmp33_; - double tmp35_ = (tmp29_) / (tmp34_); - double tmp36_ = mCompCoord[9]; - double tmp37_ = tmp35_ - tmp36_; - double tmp38_ = ElSquare(tmp37_); - double tmp39_ = tmp26_ * (tmp24_); - double tmp40_ = tmp30_ * (tmp13_); - double tmp41_ = tmp39_ - tmp40_; - double tmp42_ = tmp18_ + tmp41_; - double tmp43_ = (tmp42_) / (tmp34_); - double tmp44_ = mCompCoord[10]; - double tmp45_ = tmp43_ - tmp44_; - double tmp46_ = ElSquare(tmp45_); - double tmp47_ = tmp38_ + tmp46_; - double tmp48_ = mLocR0_0_2 * tmp19_; - double tmp49_ = mLocR0_0_1 * tmp26_; - double tmp50_ = tmp48_ - tmp49_; - double tmp51_ = mLocR0_0_0 + tmp50_; - double tmp52_ = (tmp51_) * (tmp34_); - double tmp53_ = mLocR0_0_1 * tmp30_; - double tmp54_ = mLocR0_0_0 * tmp19_; - double tmp55_ = tmp53_ - tmp54_; - double tmp56_ = mLocR0_0_2 + tmp55_; - double tmp57_ = (tmp29_) * (tmp56_); - double tmp58_ = tmp52_ - tmp57_; - double tmp59_ = ElSquare(tmp34_); - double tmp60_ = (tmp58_) / tmp59_; - double tmp61_ = mCompCoord[11]; - double tmp62_ = tmp61_ * (tmp47_); - double tmp63_ = mCompCoord[12]; - double tmp64_ = ElSquare(tmp47_); - double tmp65_ = tmp63_ * tmp64_; - double tmp66_ = tmp62_ + tmp65_; - double tmp67_ = mCompCoord[13]; - double tmp68_ = VCube(tmp47_); - double tmp69_ = tmp67_ * tmp68_; - double tmp70_ = tmp66_ + tmp69_; - double tmp71_ = 2 * (tmp60_); - double tmp72_ = tmp71_ * (tmp37_); - double tmp73_ = mLocR0_0_0 * tmp26_; - double tmp74_ = mLocR0_0_2 * tmp30_; - double tmp75_ = tmp73_ - tmp74_; - double tmp76_ = mLocR0_0_1 + tmp75_; - double tmp77_ = (tmp76_) * (tmp34_); - double tmp78_ = (tmp42_) * (tmp56_); - double tmp79_ = tmp77_ - tmp78_; - double tmp80_ = (tmp79_) / tmp59_; - double tmp81_ = 2 * (tmp80_); - double tmp82_ = tmp81_ * (tmp45_); - double tmp83_ = tmp72_ + tmp82_; - double tmp84_ = mCompCoord[16]; - double tmp85_ = mLocR0_1_2 * tmp19_; - double tmp86_ = mLocR0_1_1 * tmp26_; - double tmp87_ = tmp85_ - tmp86_; - double tmp88_ = mLocR0_1_0 + tmp87_; - double tmp89_ = (tmp88_) * (tmp34_); - double tmp90_ = mLocR0_1_1 * tmp30_; - double tmp91_ = mLocR0_1_0 * tmp19_; - double tmp92_ = tmp90_ - tmp91_; - double tmp93_ = mLocR0_1_2 + tmp92_; - double tmp94_ = (tmp29_) * (tmp93_); - double tmp95_ = tmp89_ - tmp94_; - double tmp96_ = (tmp95_) / tmp59_; - double tmp97_ = 2 * (tmp96_); - double tmp98_ = tmp97_ * (tmp37_); - double tmp99_ = mLocR0_1_0 * tmp26_; - double tmp100_ = mLocR0_1_2 * tmp30_; - double tmp101_ = tmp99_ - tmp100_; - double tmp102_ = mLocR0_1_1 + tmp101_; - double tmp103_ = (tmp102_) * (tmp34_); - double tmp104_ = (tmp42_) * (tmp93_); - double tmp105_ = tmp103_ - tmp104_; - double tmp106_ = (tmp105_) / tmp59_; - double tmp107_ = 2 * (tmp106_); - double tmp108_ = tmp107_ * (tmp45_); - double tmp109_ = tmp98_ + tmp108_; - double tmp110_ = mLocR0_2_2 * tmp19_; - double tmp111_ = mLocR0_2_1 * tmp26_; - double tmp112_ = tmp110_ - tmp111_; - double tmp113_ = mLocR0_2_0 + tmp112_; - double tmp114_ = (tmp113_) * (tmp34_); - double tmp115_ = mLocR0_2_1 * tmp30_; - double tmp116_ = mLocR0_2_0 * tmp19_; - double tmp117_ = tmp115_ - tmp116_; - double tmp118_ = mLocR0_2_2 + tmp117_; - double tmp119_ = (tmp29_) * (tmp118_); - double tmp120_ = tmp114_ - tmp119_; - double tmp121_ = (tmp120_) / tmp59_; - double tmp122_ = 2 * (tmp121_); - double tmp123_ = tmp122_ * (tmp37_); - double tmp124_ = mLocR0_2_0 * tmp26_; - double tmp125_ = mLocR0_2_2 * tmp30_; - double tmp126_ = tmp124_ - tmp125_; - double tmp127_ = mLocR0_2_1 + tmp126_; - double tmp128_ = (tmp127_) * (tmp34_); - double tmp129_ = (tmp42_) * (tmp118_); - double tmp130_ = tmp128_ - tmp129_; - double tmp131_ = (tmp130_) / tmp59_; - double tmp132_ = 2 * (tmp131_); - double tmp133_ = tmp132_ * (tmp45_); - double tmp134_ = tmp123_ + tmp133_; - double tmp135_ = -(1); - double tmp136_ = tmp135_ * mLocR0_0_2; - double tmp137_ = tmp135_ * mLocR0_0_1; - double tmp138_ = tmp135_ * mLocR0_0_0; - double tmp139_ = tmp136_ * tmp19_; - double tmp140_ = tmp137_ * tmp26_; - double tmp141_ = tmp139_ - tmp140_; - double tmp142_ = tmp138_ + tmp141_; - double tmp143_ = (tmp142_) * (tmp34_); - double tmp144_ = tmp137_ * tmp30_; - double tmp145_ = tmp138_ * tmp19_; - double tmp146_ = tmp144_ - tmp145_; - double tmp147_ = tmp136_ + tmp146_; - double tmp148_ = (tmp29_) * (tmp147_); - double tmp149_ = tmp143_ - tmp148_; - double tmp150_ = (tmp149_) / tmp59_; - double tmp151_ = 2 * (tmp150_); - double tmp152_ = tmp151_ * (tmp37_); - double tmp153_ = tmp138_ * tmp26_; - double tmp154_ = tmp136_ * tmp30_; - double tmp155_ = tmp153_ - tmp154_; - double tmp156_ = tmp137_ + tmp155_; - double tmp157_ = (tmp156_) * (tmp34_); - double tmp158_ = (tmp42_) * (tmp147_); - double tmp159_ = tmp157_ - tmp158_; - double tmp160_ = (tmp159_) / tmp59_; - double tmp161_ = 2 * (tmp160_); - double tmp162_ = tmp161_ * (tmp45_); - double tmp163_ = tmp152_ + tmp162_; - double tmp164_ = tmp135_ * mLocR0_1_2; - double tmp165_ = tmp135_ * mLocR0_1_1; - double tmp166_ = tmp135_ * mLocR0_1_0; - double tmp167_ = tmp164_ * tmp19_; - double tmp168_ = tmp165_ * tmp26_; - double tmp169_ = tmp167_ - tmp168_; - double tmp170_ = tmp166_ + tmp169_; - double tmp171_ = (tmp170_) * (tmp34_); - double tmp172_ = tmp165_ * tmp30_; - double tmp173_ = tmp166_ * tmp19_; - double tmp174_ = tmp172_ - tmp173_; - double tmp175_ = tmp164_ + tmp174_; - double tmp176_ = (tmp29_) * (tmp175_); - double tmp177_ = tmp171_ - tmp176_; - double tmp178_ = (tmp177_) / tmp59_; - double tmp179_ = 2 * (tmp178_); - double tmp180_ = tmp179_ * (tmp37_); - double tmp181_ = tmp166_ * tmp26_; - double tmp182_ = tmp164_ * tmp30_; - double tmp183_ = tmp181_ - tmp182_; - double tmp184_ = tmp165_ + tmp183_; - double tmp185_ = (tmp184_) * (tmp34_); - double tmp186_ = (tmp42_) * (tmp175_); - double tmp187_ = tmp185_ - tmp186_; - double tmp188_ = (tmp187_) / tmp59_; - double tmp189_ = 2 * (tmp188_); - double tmp190_ = tmp189_ * (tmp45_); - double tmp191_ = tmp180_ + tmp190_; - double tmp192_ = tmp135_ * mLocR0_2_2; - double tmp193_ = tmp135_ * mLocR0_2_1; - double tmp194_ = tmp135_ * mLocR0_2_0; - double tmp195_ = tmp192_ * tmp19_; - double tmp196_ = tmp193_ * tmp26_; - double tmp197_ = tmp195_ - tmp196_; - double tmp198_ = tmp194_ + tmp197_; - double tmp199_ = (tmp198_) * (tmp34_); - double tmp200_ = tmp193_ * tmp30_; - double tmp201_ = tmp194_ * tmp19_; - double tmp202_ = tmp200_ - tmp201_; - double tmp203_ = tmp192_ + tmp202_; - double tmp204_ = (tmp29_) * (tmp203_); - double tmp205_ = tmp199_ - tmp204_; - double tmp206_ = (tmp205_) / tmp59_; - double tmp207_ = 2 * (tmp206_); - double tmp208_ = tmp207_ * (tmp37_); - double tmp209_ = tmp194_ * tmp26_; - double tmp210_ = tmp192_ * tmp30_; - double tmp211_ = tmp209_ - tmp210_; - double tmp212_ = tmp193_ + tmp211_; - double tmp213_ = (tmp212_) * (tmp34_); - double tmp214_ = (tmp42_) * (tmp203_); - double tmp215_ = tmp213_ - tmp214_; - double tmp216_ = (tmp215_) / tmp59_; - double tmp217_ = 2 * (tmp216_); - double tmp218_ = tmp217_ * (tmp45_); - double tmp219_ = tmp208_ + tmp218_; - double tmp220_ = (tmp29_) * (tmp18_); - double tmp221_ = -(tmp220_); - double tmp222_ = tmp221_ / tmp59_; - double tmp223_ = 2 * (tmp222_); - double tmp224_ = tmp223_ * (tmp37_); - double tmp225_ = -(tmp13_); - double tmp226_ = tmp225_ * (tmp34_); - double tmp227_ = (tmp42_) * (tmp18_); - double tmp228_ = tmp226_ - tmp227_; - double tmp229_ = (tmp228_) / tmp59_; - double tmp230_ = 2 * (tmp229_); - double tmp231_ = tmp230_ * (tmp45_); - double tmp232_ = tmp224_ + tmp231_; - double tmp233_ = (tmp13_) * (tmp34_); - double tmp234_ = -(tmp24_); - double tmp235_ = (tmp29_) * tmp234_; - double tmp236_ = tmp233_ - tmp235_; - double tmp237_ = (tmp236_) / tmp59_; - double tmp238_ = 2 * (tmp237_); - double tmp239_ = tmp238_ * (tmp37_); - double tmp240_ = (tmp42_) * tmp234_; - double tmp241_ = -(tmp240_); - double tmp242_ = tmp241_ / tmp59_; - double tmp243_ = 2 * (tmp242_); - double tmp244_ = tmp243_ * (tmp45_); - double tmp245_ = tmp239_ + tmp244_; - double tmp246_ = -(tmp18_); - double tmp247_ = tmp246_ * (tmp34_); - double tmp248_ = (tmp247_) / tmp59_; - double tmp249_ = 2 * (tmp248_); - double tmp250_ = tmp249_ * (tmp37_); - double tmp251_ = (tmp24_) * (tmp34_); - double tmp252_ = (tmp251_) / tmp59_; - double tmp253_ = 2 * (tmp252_); - double tmp254_ = tmp253_ * (tmp45_); - double tmp255_ = tmp250_ + tmp254_; - double tmp256_ = 2 * tmp135_; - double tmp257_ = tmp256_ * (tmp37_); - double tmp258_ = tmp256_ * (tmp45_); - double tmp259_ = (tmp37_) * (tmp70_); - double tmp260_ = tmp35_ + tmp259_; - double tmp261_ = (tmp83_) * tmp61_; - double tmp262_ = 2 * (tmp83_); - double tmp263_ = tmp262_ * (tmp47_); - double tmp264_ = tmp263_ * tmp63_; - double tmp265_ = tmp261_ + tmp264_; - double tmp266_ = 3 * (tmp83_); - double tmp267_ = tmp266_ * tmp64_; - double tmp268_ = tmp267_ * tmp67_; - double tmp269_ = tmp265_ + tmp268_; - double tmp270_ = (tmp109_) * tmp61_; - double tmp271_ = 2 * (tmp109_); - double tmp272_ = tmp271_ * (tmp47_); - double tmp273_ = tmp272_ * tmp63_; - double tmp274_ = tmp270_ + tmp273_; - double tmp275_ = 3 * (tmp109_); - double tmp276_ = tmp275_ * tmp64_; - double tmp277_ = tmp276_ * tmp67_; - double tmp278_ = tmp274_ + tmp277_; - double tmp279_ = (tmp134_) * tmp61_; - double tmp280_ = 2 * (tmp134_); - double tmp281_ = tmp280_ * (tmp47_); - double tmp282_ = tmp281_ * tmp63_; - double tmp283_ = tmp279_ + tmp282_; - double tmp284_ = 3 * (tmp134_); - double tmp285_ = tmp284_ * tmp64_; - double tmp286_ = tmp285_ * tmp67_; - double tmp287_ = tmp283_ + tmp286_; - double tmp288_ = (tmp163_) * tmp61_; - double tmp289_ = 2 * (tmp163_); - double tmp290_ = tmp289_ * (tmp47_); - double tmp291_ = tmp290_ * tmp63_; - double tmp292_ = tmp288_ + tmp291_; - double tmp293_ = 3 * (tmp163_); - double tmp294_ = tmp293_ * tmp64_; - double tmp295_ = tmp294_ * tmp67_; - double tmp296_ = tmp292_ + tmp295_; - double tmp297_ = (tmp191_) * tmp61_; - double tmp298_ = 2 * (tmp191_); - double tmp299_ = tmp298_ * (tmp47_); - double tmp300_ = tmp299_ * tmp63_; - double tmp301_ = tmp297_ + tmp300_; - double tmp302_ = 3 * (tmp191_); - double tmp303_ = tmp302_ * tmp64_; - double tmp304_ = tmp303_ * tmp67_; - double tmp305_ = tmp301_ + tmp304_; - double tmp306_ = (tmp219_) * tmp61_; - double tmp307_ = 2 * (tmp219_); - double tmp308_ = tmp307_ * (tmp47_); - double tmp309_ = tmp308_ * tmp63_; - double tmp310_ = tmp306_ + tmp309_; - double tmp311_ = 3 * (tmp219_); - double tmp312_ = tmp311_ * tmp64_; - double tmp313_ = tmp312_ * tmp67_; - double tmp314_ = tmp310_ + tmp313_; - double tmp315_ = (tmp232_) * tmp61_; - double tmp316_ = 2 * (tmp232_); - double tmp317_ = tmp316_ * (tmp47_); - double tmp318_ = tmp317_ * tmp63_; - double tmp319_ = tmp315_ + tmp318_; - double tmp320_ = 3 * (tmp232_); - double tmp321_ = tmp320_ * tmp64_; - double tmp322_ = tmp321_ * tmp67_; - double tmp323_ = tmp319_ + tmp322_; - double tmp324_ = (tmp245_) * tmp61_; - double tmp325_ = 2 * (tmp245_); - double tmp326_ = tmp325_ * (tmp47_); - double tmp327_ = tmp326_ * tmp63_; - double tmp328_ = tmp324_ + tmp327_; - double tmp329_ = 3 * (tmp245_); - double tmp330_ = tmp329_ * tmp64_; - double tmp331_ = tmp330_ * tmp67_; - double tmp332_ = tmp328_ + tmp331_; - double tmp333_ = (tmp255_) * tmp61_; - double tmp334_ = 2 * (tmp255_); - double tmp335_ = tmp334_ * (tmp47_); - double tmp336_ = tmp335_ * tmp63_; - double tmp337_ = tmp333_ + tmp336_; - double tmp338_ = 3 * (tmp255_); - double tmp339_ = tmp338_ * tmp64_; - double tmp340_ = tmp339_ * tmp67_; - double tmp341_ = tmp337_ + tmp340_; - double tmp342_ = tmp257_ * tmp61_; - double tmp343_ = 2 * tmp257_; - double tmp344_ = tmp343_ * (tmp47_); - double tmp345_ = tmp344_ * tmp63_; - double tmp346_ = tmp342_ + tmp345_; - double tmp347_ = 3 * tmp257_; - double tmp348_ = tmp347_ * tmp64_; - double tmp349_ = tmp348_ * tmp67_; - double tmp350_ = tmp346_ + tmp349_; - double tmp351_ = tmp135_ * (tmp70_); - double tmp352_ = tmp258_ * tmp61_; - double tmp353_ = 2 * tmp258_; - double tmp354_ = tmp353_ * (tmp47_); - double tmp355_ = tmp354_ * tmp63_; - double tmp356_ = tmp352_ + tmp355_; - double tmp357_ = 3 * tmp258_; - double tmp358_ = tmp357_ * tmp64_; - double tmp359_ = tmp358_ * tmp67_; - double tmp360_ = tmp356_ + tmp359_; - double tmp361_ = (tmp45_) * (tmp70_); - double tmp362_ = tmp43_ + tmp361_; - - mVal[0] = (mCompCoord[14] + (tmp260_) * tmp84_) - mLocMesIm_x; - - mCompDer[0][0] = (tmp60_ + (tmp60_) * (tmp70_) + (tmp269_) * (tmp37_)) * tmp84_; - mCompDer[0][1] = (tmp96_ + (tmp96_) * (tmp70_) + (tmp278_) * (tmp37_)) * tmp84_; - mCompDer[0][2] = (tmp121_ + (tmp121_) * (tmp70_) + (tmp287_) * (tmp37_)) * tmp84_; - mCompDer[0][3] = (tmp150_ + (tmp150_) * (tmp70_) + (tmp296_) * (tmp37_)) * tmp84_; - mCompDer[0][4] = (tmp178_ + (tmp178_) * (tmp70_) + (tmp305_) * (tmp37_)) * tmp84_; - mCompDer[0][5] = (tmp206_ + (tmp206_) * (tmp70_) + (tmp314_) * (tmp37_)) * tmp84_; - mCompDer[0][6] = (tmp222_ + (tmp222_) * (tmp70_) + (tmp323_) * (tmp37_)) * tmp84_; - mCompDer[0][7] = (tmp237_ + (tmp237_) * (tmp70_) + (tmp332_) * (tmp37_)) * tmp84_; - mCompDer[0][8] = (tmp248_ + (tmp248_) * (tmp70_) + (tmp341_) * (tmp37_)) * tmp84_; - mCompDer[0][9] = (tmp351_ + (tmp350_) * (tmp37_)) * tmp84_; - mCompDer[0][10] = (tmp360_) * (tmp37_) * tmp84_; - mCompDer[0][11] = (tmp47_) * (tmp37_) * tmp84_; - mCompDer[0][12] = tmp64_ * (tmp37_) * tmp84_; - mCompDer[0][13] = tmp68_ * (tmp37_) * tmp84_; - mCompDer[0][14] = 1; - mCompDer[0][15] = 0; - mCompDer[0][16] = tmp260_; - mVal[1] = (mCompCoord[15] + (tmp362_) * tmp84_) - mLocMesIm_y; - - mCompDer[1][0] = (tmp80_ + (tmp80_) * (tmp70_) + (tmp269_) * (tmp45_)) * tmp84_; - mCompDer[1][1] = (tmp106_ + (tmp106_) * (tmp70_) + (tmp278_) * (tmp45_)) * tmp84_; - mCompDer[1][2] = (tmp131_ + (tmp131_) * (tmp70_) + (tmp287_) * (tmp45_)) * tmp84_; - mCompDer[1][3] = (tmp160_ + (tmp160_) * (tmp70_) + (tmp296_) * (tmp45_)) * tmp84_; - mCompDer[1][4] = (tmp188_ + (tmp188_) * (tmp70_) + (tmp305_) * (tmp45_)) * tmp84_; - mCompDer[1][5] = (tmp216_ + (tmp216_) * (tmp70_) + (tmp314_) * (tmp45_)) * tmp84_; - mCompDer[1][6] = (tmp229_ + (tmp229_) * (tmp70_) + (tmp323_) * (tmp45_)) * tmp84_; - mCompDer[1][7] = (tmp242_ + (tmp242_) * (tmp70_) + (tmp332_) * (tmp45_)) * tmp84_; - mCompDer[1][8] = (tmp252_ + (tmp252_) * (tmp70_) + (tmp341_) * (tmp45_)) * tmp84_; - mCompDer[1][9] = (tmp350_) * (tmp45_) * tmp84_; - mCompDer[1][10] = (tmp351_ + (tmp360_) * (tmp45_)) * tmp84_; - mCompDer[1][11] = (tmp47_) * (tmp45_) * tmp84_; - mCompDer[1][12] = tmp64_ * (tmp45_) * tmp84_; - mCompDer[1][13] = tmp68_ * (tmp45_) * tmp84_; - mCompDer[1][14] = 0; - mCompDer[1][15] = 1; - mCompDer[1][16] = tmp362_; -} - - -void cGeneratedCodeTestCam::ComputeValDerivHessian() -{ - ELISE_ASSERT(false,"Foncteur cGeneratedCodeTestCam Has no Der Sec"); -} - -void cGeneratedCodeTestCam::SetMesIm_x(double aVal){ mLocMesIm_x = aVal;} -void cGeneratedCodeTestCam::SetMesIm_y(double aVal){ mLocMesIm_y = aVal;} -void cGeneratedCodeTestCam::SetR0_0_0(double aVal){ mLocR0_0_0 = aVal;} -void cGeneratedCodeTestCam::SetR0_0_1(double aVal){ mLocR0_0_1 = aVal;} -void cGeneratedCodeTestCam::SetR0_0_2(double aVal){ mLocR0_0_2 = aVal;} -void cGeneratedCodeTestCam::SetR0_1_0(double aVal){ mLocR0_1_0 = aVal;} -void cGeneratedCodeTestCam::SetR0_1_1(double aVal){ mLocR0_1_1 = aVal;} -void cGeneratedCodeTestCam::SetR0_1_2(double aVal){ mLocR0_1_2 = aVal;} -void cGeneratedCodeTestCam::SetR0_2_0(double aVal){ mLocR0_2_0 = aVal;} -void cGeneratedCodeTestCam::SetR0_2_1(double aVal){ mLocR0_2_1 = aVal;} -void cGeneratedCodeTestCam::SetR0_2_2(double aVal){ mLocR0_2_2 = aVal;} - - - -double * cGeneratedCodeTestCam::AdrVarLocFromString(const std::string & aName) -{ - if (aName == "MesIm_x") return & mLocMesIm_x; - if (aName == "MesIm_y") return & mLocMesIm_y; - if (aName == "R0_0_0") return & mLocR0_0_0; - if (aName == "R0_0_1") return & mLocR0_0_1; - if (aName == "R0_0_2") return & mLocR0_0_2; - if (aName == "R0_1_0") return & mLocR0_1_0; - if (aName == "R0_1_1") return & mLocR0_1_1; - if (aName == "R0_1_2") return & mLocR0_1_2; - if (aName == "R0_2_0") return & mLocR0_2_0; - if (aName == "R0_2_1") return & mLocR0_2_1; - if (aName == "R0_2_2") return & mLocR0_2_2; - return 0; -} - - -cElCompiledFonc::cAutoAddEntry cGeneratedCodeTestCam::mTheAuto("cGeneratedCodeTestCam",cGeneratedCodeTestCam::Alloc); - - -cElCompiledFonc * cGeneratedCodeTestCam::Alloc() -{ return new cGeneratedCodeTestCam(); -} - - From 88628cf863b44a523d395e88e3cd7359705c2426 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Mon, 16 Dec 2024 22:54:04 +0100 Subject: [PATCH 3/5] In supr V1 ... --- MMVII/src/MMV1/DenseMapCorresp.cpp | 8 ----- MMVII/src/MMV1/ExportHomMMV1.cpp | 22 ++++++++++-- MMVII/src/MMV1/ExportMeasuresImGCPMMV1.cpp | 25 +++++++++---- MMVII/src/MMV1/cGeneratedCodeTestCam.cpp | 0 MMVII/src/MMV1/cGeneratedCodeTestCam.h | 42 ---------------------- 5 files changed, 38 insertions(+), 59 deletions(-) delete mode 100755 MMVII/src/MMV1/DenseMapCorresp.cpp delete mode 100755 MMVII/src/MMV1/cGeneratedCodeTestCam.cpp delete mode 100755 MMVII/src/MMV1/cGeneratedCodeTestCam.h diff --git a/MMVII/src/MMV1/DenseMapCorresp.cpp b/MMVII/src/MMV1/DenseMapCorresp.cpp deleted file mode 100755 index 231d57c2c9..0000000000 --- a/MMVII/src/MMV1/DenseMapCorresp.cpp +++ /dev/null @@ -1,8 +0,0 @@ -#include "V1VII.h" - - -//extern std::string MM3DFixeByMMVII; - -namespace MMVII -{ -}; diff --git a/MMVII/src/MMV1/ExportHomMMV1.cpp b/MMVII/src/MMV1/ExportHomMMV1.cpp index ad1f232ed5..81dcb8a96c 100755 --- a/MMVII/src/MMV1/ExportHomMMV1.cpp +++ b/MMVII/src/MMV1/ExportHomMMV1.cpp @@ -1,9 +1,16 @@ +#define WITH_MMV1_FUNCTION false + +#if (WITH_MMV1_FUNCTION) #include "V1VII.h" +#endif + #include "MMVII_util.h" #include "MMVII_MeasuresIm.h" namespace MMVII { + +#if (WITH_MMV1_FUNCTION) cHomogCpleIm ToMMVII(const cNupletPtsHomologues & aNUp) { @@ -75,13 +82,24 @@ void cImportHomV1::GetHom(cSetHomogCpleIm & aPackV2,const std::string & aNameIm /* */ /* ************************************************* */ -cInterfImportHom::~cInterfImportHom() +cInterfImportHom * cInterfImportHom::CreateImportV1(const std::string&aDir,const std::string&aSubDir,const std::string&aExt) { + return new cImportHomV1(aDir,aSubDir,aExt); } +#else cInterfImportHom * cInterfImportHom::CreateImportV1(const std::string&aDir,const std::string&aSubDir,const std::string&aExt) { - return new cImportHomV1(aDir,aSubDir,aExt); + MMVII_INTERNAL_ERROR("No CreateImportV1 "); + return nullptr; } +cInterfImportHom::~cInterfImportHom() +{ +} + +#endif + + + }; diff --git a/MMVII/src/MMV1/ExportMeasuresImGCPMMV1.cpp b/MMVII/src/MMV1/ExportMeasuresImGCPMMV1.cpp index 4106976e55..349140c569 100755 --- a/MMVII/src/MMV1/ExportMeasuresImGCPMMV1.cpp +++ b/MMVII/src/MMV1/ExportMeasuresImGCPMMV1.cpp @@ -1,4 +1,9 @@ +#define WITH_MMV1_FUNCTION false + +#if (WITH_MMV1_FUNCTION) #include "V1VII.h" +#endif + #include "MMVII_util.h" #include "MMVII_MeasuresIm.h" @@ -10,6 +15,8 @@ namespace MMVII /* */ /* ************************************************* */ +#if (WITH_MMV1_FUNCTION) + void ImportMesImV1(std::list & aResult,const std::string & aNameFileMesImV1) { aResult.clear(); @@ -36,13 +43,6 @@ cSetMesGCP ImportMesGCPV1(const std::string & aNameFileMesGCPV1,const std::strin cMes1GCP aMesV2(ToMMVII(aMesV1.Pt()),aMesV1.NamePt(),1.0); aMesV2.SetSigma2(ToMMVII(aMesV1.Incertitude())); - /* - aMesV2.Sigma2() = {0,0,0,0,0,0}; - - (aMesV2.Sigma2())[cMes1GCP::IndXX] = Square(aMesV1.Incertitude().x); - (aMesV2.Sigma2())[cMes1GCP::IndYY] = Square(aMesV1.Incertitude().y); - (aMesV2.Sigma2())[cMes1GCP::IndZZ] = Square(aMesV1.Incertitude().z); - */ aResult.AddMeasure(aMesV2); } @@ -50,6 +50,17 @@ cSetMesGCP ImportMesGCPV1(const std::string & aNameFileMesGCPV1,const std::strin return aResult; } +#else +void ImportMesImV1(std::list & aResult,const std::string & aNameFileMesImV1) +{ + MMVII_INTERNAL_ERROR("No ImportMesImV1 "); +} +cSetMesGCP ImportMesGCPV1(const std::string & aNameFileMesGCPV1,const std::string & aNameSet) +{ + MMVII_INTERNAL_ERROR("No ImportMesGCPV1"); + return cSetMesGCP (aNameSet); +} +#endif diff --git a/MMVII/src/MMV1/cGeneratedCodeTestCam.cpp b/MMVII/src/MMV1/cGeneratedCodeTestCam.cpp deleted file mode 100755 index e69de29bb2..0000000000 diff --git a/MMVII/src/MMV1/cGeneratedCodeTestCam.h b/MMVII/src/MMV1/cGeneratedCodeTestCam.h deleted file mode 100755 index 1f803d250e..0000000000 --- a/MMVII/src/MMV1/cGeneratedCodeTestCam.h +++ /dev/null @@ -1,42 +0,0 @@ -// File Automatically generated by eLiSe -#include "StdAfx.h" - - -class cGeneratedCodeTestCam: public cElCompiledFonc -{ - public : - - cGeneratedCodeTestCam(); - void ComputeVal(); - void ComputeValDeriv(); - void ComputeValDerivHessian(); - double * AdrVarLocFromString(const std::string &); - void SetMesIm_x(double); - void SetMesIm_y(double); - void SetR0_0_0(double); - void SetR0_0_1(double); - void SetR0_0_2(double); - void SetR0_1_0(double); - void SetR0_1_1(double); - void SetR0_1_2(double); - void SetR0_2_0(double); - void SetR0_2_1(double); - void SetR0_2_2(double); - - - static cAutoAddEntry mTheAuto; - static cElCompiledFonc * Alloc(); - private : - - double mLocMesIm_x; - double mLocMesIm_y; - double mLocR0_0_0; - double mLocR0_0_1; - double mLocR0_0_2; - double mLocR0_1_0; - double mLocR0_1_1; - double mLocR0_1_2; - double mLocR0_2_0; - double mLocR0_2_1; - double mLocR0_2_2; -}; From 81c6749cbcea3c7ec97e87a0e0592476e2316f9a Mon Sep 17 00:00:00 2001 From: deseilligny Date: Sat, 21 Dec 2024 18:44:22 +0100 Subject: [PATCH 4/5] In supress V1 for root poly, Algo OK --- .../Eigen/src/Polynomials/PolynomialUtils.h | 3 +- MMVII/include/MMVII_nums.h | 6 +- MMVII/src/MMV1/Numerics.cpp | 10 - MMVII/src/UtiMaths/Polynoms.cpp | 335 +++++++++++++++++- 4 files changed, 336 insertions(+), 18 deletions(-) diff --git a/MMVII/ExternalInclude/eigen-3.4.0/unsupported/Eigen/src/Polynomials/PolynomialUtils.h b/MMVII/ExternalInclude/eigen-3.4.0/unsupported/Eigen/src/Polynomials/PolynomialUtils.h index 394e857acf..b490a03588 100644 --- a/MMVII/ExternalInclude/eigen-3.4.0/unsupported/Eigen/src/Polynomials/PolynomialUtils.h +++ b/MMVII/ExternalInclude/eigen-3.4.0/unsupported/Eigen/src/Polynomials/PolynomialUtils.h @@ -53,7 +53,8 @@ T poly_eval( const Polynomials& poly, const T& x ) { T val=poly[0]; T inv_x = T(1)/x; - for( DenseIndex i=1; i class cPolynom Type Value(const Type & 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; cPolynom operator * (const Type & aVal) const; - std::vector RealRoots(const Type & aTol,int ItMax); + cPolynom Deriv() const; + + std::vector RealRoots(const Type & aTol,int ItMax) const; Type& operator [] (size_t aK) {return mVCoeffs[aK];} diff --git a/MMVII/src/MMV1/Numerics.cpp b/MMVII/src/MMV1/Numerics.cpp index 072cfdf883..9546bf2a0f 100755 --- a/MMVII/src/MMV1/Numerics.cpp +++ b/MMVII/src/MMV1/Numerics.cpp @@ -54,16 +54,6 @@ INST_V1ROOTS(tREAL8) INST_V1ROOTS(tREAL16) -/* -template void - RealRootsOfRealPolynome - ( - ElSTDNS vector & Sols, - const ElPolynome &aPol, - Type tol, - INT ItMax - ) - */ }; diff --git a/MMVII/src/UtiMaths/Polynoms.cpp b/MMVII/src/UtiMaths/Polynoms.cpp index 3cd3ddd8e3..3064ccc416 100755 --- a/MMVII/src/UtiMaths/Polynoms.cpp +++ b/MMVII/src/UtiMaths/Polynoms.cpp @@ -1,17 +1,303 @@ + +#if defined(__GNUC__) && !defined(__clang__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wmaybe-uninitialized" +#endif + #include "cMMVII_Appli.h" #include "V1VII.h" +#include +// #include "unsupported/Eigen/Polynomials" + +//#if defined(__GNUC__) && !defined(__clang__) +// # pragma GCC diagnostic pop +// #endif -/* // #include "include/MMVII_nums.h" // #include "include/MMVII_Bench.h" - //#include "Eigen/unsupported/Eigen/Polynomials" - */ +//#include +//#include + +using namespace Eigen; +using namespace std; namespace MMVII { +/** + Class for extraction of roots of polynoms using "Companion matrix method" + + See https://en.wikipedia.org/wiki/Companion_matrix +*/ + + +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; + private : + + cPolynom mPol; + cPolynom mDPol; + size_t mDeg; + MatrixXd mCompM; ///< companion matrix + VectorXcd mCEV; ///< Complex eigen values + std::vector mRR; ///< Real roots + +}; + + +template cEigenPolynRoots::cEigenPolynRoots(const cPolynom & aPol) : + mPol (aPol), + mDPol (aPol.Deriv()), + mDeg (mPol.Degree()), + mCompM (mDeg,mDeg) +{ + + 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; + } + + 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; + } + + EigenSolver aSolver(mCompM); + mCEV = aSolver.eigenvalues() ; + + for (const auto & aC : mCEV) + if (RootIsReal(aC)) + mRR.push_back(aC.real()); + std::sort(mRR.begin(),mRR.end()); +} + +template const std::vector & cEigenPolynRoots::RealRoots() const {return mRR;} + +template const VectorXcd & cEigenPolynRoots::ComplexRoots() const {return mCEV;} + +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) +{ + // [1] Test is "aC" is a real number + + 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) + { + if (sayWhy) + *sayWhy = "ABS REAL COMPLEX=" + ToStr(std::abs(C_i)); + return false; + } + + tREAL8 C_r = aC.real(); + // [1.1] if relative imaginary part is "big" + if (std::abs(C_i) > 1e-8 * (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)); + // [2.1] if absolute value of polynom is big + if (aAbsVP > 1e-5) + { + if (sayWhy) + *sayWhy = "ABS VALUE POL " + ToStr(aAbsVP); + return false; + } + + tREAL8 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 (sayWhy) + *sayWhy = "RELATIVE VALUE POL " + ToStr(aAbsVP/(aAVA+1e-5)); + return false; + } + + if (sayWhy) + *sayWhy = "is real"; + + return true; +} + + +template class cEigenPolynRoots; + +template void My_Roots(const cPolynom & aPol1) +{ +} +template <> void My_Roots(const cPolynom & aPol1) +{ +// StdOut() << "DDDD " << aPol1.Degree() << "\n"; + // (X2+1)(X-1) = X3-X2+X-1 + int aNb=300; + // vector aCoeffs1 = {-1,1,-1,1,5,-2,0.12}; + // cPolynom aPol1(aCoeffs1); + + cAutoTimerSegm aTimeEigen(GlobAppTS(),"Eigen"); + for (int aK=0 ; aK aEPR(aPol1); + } + + cAutoTimerSegm aTimeV1(GlobAppTS(),"V1"); + for (int aK=0 ; aK aEPR(aPol1); + auto aV2 = aEPR.RealRoots(); + auto aV1 = aPol1.RealRoots(1e-20,60); + StdOut() << " SZzzzZ= " << aV1.size() << " " << aV2.size() << "\n"; + if (aV1.size() != aV2.size()) + { + StdOut() << "Coeffs=" << aPol1.VCoeffs() << "\n"; + StdOut() << "V1=" << aV1 << "\n"; + StdOut() << "V2=" << aV2 << "\n"; + for (const auto & aC : aEPR.ComplexRoots()) + { + std::string strWhy; + bool isR = aEPR.RootIsReal(aC,&strWhy); + 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); + +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() +{ + // 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); + +} +#endif + /* ************************************************************************ */ /* */ @@ -24,7 +310,7 @@ namespace MMVII template cPolynom::cPolynom(const tCoeffs & aVCoeffs) : mVCoeffs (aVCoeffs) { - MMVII_INTERNAL_ASSERT_tiny(!mVCoeffs.empty(),"Empty polynom not handled"); + // MMVII_INTERNAL_ASSERT_tiny(!mVCoeffs.empty(),"Empty polynom not handled"); } template cPolynom::cPolynom(const cPolynom & aPol) : cPolynom(aPol.mVCoeffs) @@ -96,8 +382,9 @@ template return aRes; } -template std::vector cPolynom::RealRoots(const Type & aTol,int ItMax) +template std::vector cPolynom::RealRoots(const Type & aTol,int ItMax) const { +// StdOut() << "RealRootsRealRootsRealRootsRealRootsRealRootsRealRootsRealRootsRealRoots \n"; getchar(); return V1RealRoots(mVCoeffs,aTol,ItMax); } @@ -119,6 +406,22 @@ template Type cPolynom::Value(const Type & aVal) const return aResult; } +template Type cPolynom::AbsValue(const Type & aVal) const +{ + Type aResult = 0.0; + Type aPowV = 1.0; + for (const auto & aCoef : mVCoeffs) + { + aResult += std::abs(aCoef * aPowV); + aPowV *= aVal; + } + return aResult; +} + + + + + template cPolynom cPolynom::operator * (const cPolynom & aP2) const { cPolynom aRes(Degree() + aP2.Degree()); @@ -159,6 +462,15 @@ template cPolynom cPolynom::operator - (const cPolynom return aRes; } +template cPolynom cPolynom::Deriv() const +{ + std::vector aVCD; // Vector Coeff Derivates + for (size_t aDeg=1 ; aDeg(aVCD); +} + template cPolynom cPolynom::operator * (const Type & aMul) const { @@ -201,6 +513,10 @@ template void TplBenchPolynome() cPolynom aP1min2 = aPol1 - aPol2; cPolynom aP2min1 = aPol2 - aPol1; + cPolynom aDerP1P2_A = aP1mul2.Deriv(); + cPolynom aDerP1P2_B = aPol1.Deriv() * aPol2 + aPol1*aPol2.Deriv(); + + Type aEps = tElemNumTrait ::Accuracy(); for (int aK=0 ; aK< 20 ; aK++) @@ -210,6 +526,10 @@ template void TplBenchPolynome() Type aChekP = aPol1.Value(aV) + aPol2.Value(aV); Type aChekMin = aPol1.Value(aV) - aPol2.Value(aV); + Type aDerA = aDerP1P2_A.Value(aV) ; + Type aDerB = aDerP1P2_B.Value(aV) ; + MMVII_INTERNAL_ASSERT_bench(std::abs(aDerA-aDerB) void TplBenchPolynome() MMVII_INTERNAL_ASSERT_bench(std::abs(RelativeSafeDifference(-aChekMin,aP2min1.Value(aV))) aVRootsGen; Type aAmpl = 10*RandUnif_NotNull(1e-2); @@ -247,6 +567,7 @@ template void TplBenchPolynome() Type aDif = RelativeSafeDifference(aVRootsGen[aK],aVRootsCalc[aK]); MMVII_INTERNAL_ASSERT_bench(aDif(); TplBenchPolynome(); TplBenchPolynome(); From d9a18a3d63680cd56815a54c4a96dc8150315120 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Mon, 23 Dec 2024 15:59:16 +0100 Subject: [PATCH 5/5] 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();