Skip to content

Commit

Permalink
Prop cov, blocage ?
Browse files Browse the repository at this point in the history
  • Loading branch information
deseilligny committed Jul 10, 2022
1 parent 95a640d commit 3a4ea9a
Show file tree
Hide file tree
Showing 9 changed files with 317 additions and 68 deletions.
2 changes: 1 addition & 1 deletion MMVII/bin/Mk-MMVII.makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@ HEADER=$(wildcard ${MMV2DirIncl}*.h)
# Binaries
#== CFLAGS etc...
#
CXX=g++-8
CXX=g++-9
CFlags= "-fopenmp" "-std=c++17" "-Wall" "-Werror" "-O4" "-fPIC" -I${MMV2Dir} -I${MMV2Dir}/ExternalInclude -I${MMDir}/include/ -I${MMDir} -D'MMVII_INSTALL_PATH="${MMV2_INSTALL_PATH}"'
#CFlags= "-fopenmp" "-std=c++17" "-Wall" "-Werror" "-O4" "-march=native" "-fPIC" -I${MMV2Dir} -I${MMV2Dir}/ExternalInclude -I${MMDir}/include/ -I${MMDir} -D'MMVII_INSTALL_PATH="${MMV2_INSTALL_PATH}"'

Expand Down
27 changes: 24 additions & 3 deletions MMVII/include/MMVII_SysSurR.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ template <class Type> class cResolSysNonLinear
typedef cSetIORSNL_SameTmp<Type> tSetIO_ST;


typedef cLinearOverCstrSys<Type> tSysSR;
typedef cLinearOverCstrSys<Type> tLinearSysSR;
typedef cDenseVect<Type> tDVect;
typedef cSparseVect<Type> tSVect;
typedef std::vector<Type> tStdVect;
Expand All @@ -42,7 +42,7 @@ template <class Type> class cResolSysNonLinear
/// basic constructor, using a mode of matrix + a solution init
cResolSysNonLinear(eModeSSR,const tDVect & aInitSol);
/// constructor using linear system, allow finer control
cResolSysNonLinear(tSysSR *,const tDVect & aInitSol);
cResolSysNonLinear(tLinearSysSR *,const tDVect & aInitSol);
/// destructor
~cResolSysNonLinear();

Expand All @@ -51,6 +51,8 @@ template <class Type> class cResolSysNonLinear
/// Value of a given num var
const Type & CurSol(int aNumV) const;

tLinearSysSR * SysLinear() ;

/// Solve solution, update the current solution, Reset the least square system
const tDVect & SolveUpdateReset() ;

Expand Down Expand Up @@ -82,7 +84,7 @@ template <class Type> class cResolSysNonLinear

int mNbVar; ///< Number of variable, facility
tDVect mCurGlobSol; ///< Curent solution
tSysSR* mSys; ///< Sys to solve equations, equation are concerning the differences with current solution
tLinearSysSR* mSysLinear; ///< Sys to solve equations, equation are concerning the differences with current solution
};

/** Class for weighting residuals : compute the vector of weight from a
Expand Down Expand Up @@ -219,6 +221,15 @@ template <class Type> class cLinearOverCstrSys : public cMemCheck
/// Accessor
int NbVar() const;

/// Normal Matrix defined 4 now only in cLeasSqtAA, maybe later defined in other classe, else error
virtual cDenseMatrix<Type> V_tAA() const;
/// Idem "normal" vector
virtual cDenseVect<Type> V_tARhs() const;
/// Indicate if it gives acces to these normal "stuff"
virtual bool Acces2NormalEq() const;


virtual void AddCov(const cDenseMatrix<Type> &,const cDenseVect<Type>& ,const std::vector<int> &aVInd);

protected :
int mNbVar;
Expand Down Expand Up @@ -264,6 +275,8 @@ template <class Type> class cLeasSq : public cLinearOverCstrSys<Type>

/// Basic dense systems => cLeasSqtAA
static cLeasSq<Type>* AllocDenseLstSq(int aNbVar);


};

/** Implemant least by suming tA A , simple and efficient, by the way known to have
Expand Down Expand Up @@ -292,6 +305,14 @@ template <class Type> class cLeasSqtAA : public cLeasSq<Type>
cDenseMatrix<Type> & tAA () ; ///< Accessor
cDenseVect<Type> & tARhs () ; ///< Accessor

/// access to tAA via virtual interface
cDenseMatrix<Type> V_tAA() const override;
/// access to tARhs via virtual interface
cDenseVect<Type> V_tARhs() const override;
/// true because acces is given
bool Acces2NormalEq() const override;

void AddCov(const cDenseMatrix<Type> &,const cDenseVect<Type>& ,const std::vector<int> &aVInd) override;
private :
cDenseMatrix<Type> mtAA; /// Som(W tA A)
cDenseVect<Type> mtARhs; /// Som(W tA Rhs)
Expand Down
12 changes: 12 additions & 0 deletions MMVII/include/MMVII_Tpl_Images.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,10 +153,22 @@ template<class T1,class T2,class T3>
return aI1;
}

template<class T1,class T2,class T3>
cIm1D<T1> AddImage(T1* /*Type specifier*/ ,const cIm1D<T2> & aI2,const cIm1D<T3> & aI3)
{
cIm1D<T1> aI1(aI2.DIm().X0(),aI2.DIm().X1());
AddImageInPlace(aI1.DIm(),aI2.DIm(),aI3.DIm());
return aI1;
}

template<class T2,class T3> cIm2D<T2> operator + (const cIm2D<T2> & aI2,const cIm2D<T3> & aI3)
{
return AddImage((T2 *)nullptr,aI2,aI3);
}
template<class T2,class T3> cIm1D<T2> operator + (const cIm1D<T2> & aI2,const cIm1D<T3> & aI3)
{
return AddImage((T2 *)nullptr,aI2,aI3);
}

template<class T2> cDenseMatrix<T2> operator + (const cDenseMatrix<T2> & aI2,const cDenseMatrix<T2> & aI3)
{
Expand Down
62 changes: 62 additions & 0 deletions MMVII/src/Matrix/cLeasSqtAA.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,6 +262,39 @@ template<class Type> cDenseVect<Type> cLeasSqtAA<Type>::SparseSolve()
return EigenSolveCholeskyarseFromV3(aVCoeff,mtARhs);
}

template<class Type> cDenseMatrix<Type> cLeasSqtAA<Type>::V_tAA() const
{
cDenseMatrix<Type> aRes = mtAA.Dup();
aRes.SelfSymetrizeBottom();
return aRes;
}

template<class Type> cDenseVect<Type> cLeasSqtAA<Type>::V_tARhs() const
{
return mtARhs;
}
template<class Type> bool cLeasSqtAA<Type>::Acces2NormalEq() const
{
return true;
}

template<class Type> void cLeasSqtAA<Type>::AddCov
(const cDenseMatrix<Type> & aMat,const cDenseVect<Type>& aVect,const std::vector<int> &aVInd)
{
for (int aKx = 0 ; aKx<int(aVInd.size()) ; aKx++)
{
mtARhs(aVInd[aKx]) += aVect(aKx);
for (int aKy = 0 ; aKy<int(aVInd.size()) ; aKy++)
{
// Only triangular sup used
if (aVInd[aKx] >= aVInd[aKy])
mtAA.AddElem(aVInd[aKx],aVInd[aKy],aMat.GetElem(aKx,aKy));
}
}
}



/* *********************************** */
/* */
/* cLeasSq */
Expand Down Expand Up @@ -360,6 +393,35 @@ template<class Type> void cLinearOverCstrSys<Type>::AddObsWithTmpUK(const cSetIO
}


template<class Type> cDenseMatrix<Type> cLinearOverCstrSys<Type>::V_tAA() const
{
MMVII_INTERNAL_ERROR("No acces to tAA for this class");

return *((cDenseMatrix<Type> *)nullptr);
}

template<class Type> cDenseVect<Type> cLinearOverCstrSys<Type>::V_tARhs() const
{
MMVII_INTERNAL_ERROR("No acces to tARhs for this class");

return *((cDenseVect<Type> *)nullptr);
}

template <class Type> void cLinearOverCstrSys<Type>::AddCov
(const cDenseMatrix<Type> &,const cDenseVect<Type>& ,const std::vector<int> &aVInd)
{
MMVII_INTERNAL_ERROR("No AddCov for this class");
}



template<class Type> bool cLinearOverCstrSys<Type>::Acces2NormalEq() const
{
return false;
}




/* ===================================================== */
/* ===================================================== */
Expand Down
23 changes: 14 additions & 9 deletions MMVII/src/Matrix/cResolSysNonLinear.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,10 +154,10 @@ template <class Type> std::vector<Type> cResidualWeighter<Type>::WeightOfResid
/* */
/* ************************************************************ */

template <class Type> cResolSysNonLinear<Type>:: cResolSysNonLinear(tSysSR * aSys,const tDVect & aInitSol) :
template <class Type> cResolSysNonLinear<Type>::cResolSysNonLinear(tLinearSysSR * aSys,const tDVect & aInitSol) :
mNbVar (aInitSol.Sz()),
mCurGlobSol (aInitSol.Dup()),
mSys (aSys)
mSysLinear (aSys)
{
}

Expand All @@ -172,7 +172,7 @@ template <class Type> void cResolSysNonLinear<Type>::AddEqFixVar(const int & a
tSVect aSV;
aSV.AddIV(aNumV,1.0);
// Dont forget that the linear system compute the difference with current solution ...
mSys->AddObservation(aWeight,aSV,aVal-CurSol(aNumV));
mSysLinear->AddObservation(aWeight,aSV,aVal-CurSol(aNumV));
// mSys->AddObservation(aWeight,aSV,CurSol(aNumV)+aVal);
}

Expand All @@ -188,13 +188,18 @@ template <class Type> const Type & cResolSysNonLinear<Type>::CurSol(int aNumV) c

template <class Type> const cDenseVect<Type> & cResolSysNonLinear<Type>::SolveUpdateReset()
{
// mCurGlobSol += mSys->Solve();
mCurGlobSol += mSys->SparseSolve();
mSys->Reset();
mCurGlobSol += mSysLinear->Solve();
// mCurGlobSol += mSysLinear->SparseSolve();
mSysLinear->Reset();

return mCurGlobSol;
}

template <class Type> cLinearOverCstrSys<Type> * cResolSysNonLinear<Type>::SysLinear()
{
return mSysLinear;
}


template <class Type> void cResolSysNonLinear<Type>::AddEqFixCurVar(const int & aNumV,const Type& aWeight)
{
Expand All @@ -218,7 +223,7 @@ template <class Type> void cResolSysNonLinear<Type>::CalcAndAddObs

template <class Type> cResolSysNonLinear<Type>::~cResolSysNonLinear()
{
delete mSys;
delete mSysLinear;
}


Expand All @@ -243,7 +248,7 @@ template <class Type> void cResolSysNonLinear<Type>::AddObs ( const std::vector<
aSV.AddIV(aIO.mVIndUk[aKUk],aVDer[aKUk]);
}
// Note the minus sign : F(X0+dx) = F(X0) + Gx.dx => Gx.dx = -F(X0)
mSys->AddObservation(aW,aSV,-aIO.mVals[aKVal]);
mSysLinear->AddObservation(aW,aSV,-aIO.mVals[aKVal]);
}
}

Expand All @@ -265,7 +270,7 @@ template <class Type> void cResolSysNonLinear<Type>::AddEq2Subst

template <class Type> void cResolSysNonLinear<Type>::AddObsWithTmpUK (const tSetIO_ST & aSetIO)
{
mSys->AddObsWithTmpUK(aSetIO);
mSysLinear->AddObsWithTmpUK(aSetIO);
}

template <class Type> void cResolSysNonLinear<Type>::CalcVal
Expand Down
6 changes: 5 additions & 1 deletion MMVII/src/TutoBenchTrianguRSNL/BenchNetwRSNL.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,8 +82,12 @@ void BenchSSRNL(cParamExeBench & aParam)
{
if (! aParam.NewBench("SSRNL")) return;

if (DEBUG_RSNL)
StdOut() << "DEBUG_RSNLlllllllllllllllllllllllllllllllllllllll\n";

for (int aK=0 ; aK<10 ; aK++)
{
cMainNetwork <tREAL8> aNet(eModeSSR::eSSR_LsqDense,5,false);
cMainNetwork <tREAL8> aNet(eModeSSR::eSSR_LsqDense,2,false);

aNet.TestCov();
}
Expand Down
20 changes: 16 additions & 4 deletions MMVII/src/TutoBenchTrianguRSNL/TrianguRSNL.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@ namespace NS_Bench_RSNL
{
#define DEBUG_RSNL true
#if (DEBUG_RSNL)
constexpr double AMPL_Grid2Real= 0.00; // Amplitude of random differerence between real position and regular grid
constexpr double AMPL_Real2Init = 0.00; // Amplitude of random and syst differerence betwen real an init position
constexpr double AMPL_Grid2Real= 0.1; // Amplitude of random differerence between real position and regular grid
constexpr double AMPL_Real2Init = 0.1; // Amplitude of random and syst differerence betwen real an init position
#else
constexpr double AMPL_Grid2Real= 0.1; // Amplitude of random differerence between real position and regular grid
constexpr double AMPL_Real2Init = 0.1; // Amplitude of random and syst differerence betwen real an init position
Expand All @@ -46,6 +46,9 @@ constexpr double AMPL_Real2Init = 0.1; // Amplitude of random and syst differer
template <class Type> class cMainNetwork;
template <class Type> class cPNetwork;

// This class is used only in covariance propagation
template <class Type> class cElemCalcCoordInit ;

template <class Type> class cPNetwork
{
public :
Expand Down Expand Up @@ -91,11 +94,12 @@ template <class Type> class cMainNetwork
typedef tPNet * tPNetPtr;
typedef cResolSysNonLinear<Type> tSys;
typedef NS_SymbolicDerivative::cCalculator<tREAL8> tCalc;
typedef cElemCalcCoordInit<Type> tECCI;

/// initial simplify constructor, take N a parameterand construct [-N,N]x[N,N]
cMainNetwork(eModeSSR aMode,int aN,bool WithSchurr,cParamSparseNormalLstSq * = nullptr);
cMainNetwork(eModeSSR aMode,int aN,bool WithSchurr,cParamSparseNormalLstSq * = nullptr,tECCI * =nullptr);

cMainNetwork(eModeSSR aMode,cRect2,bool WithSchurr,cParamSparseNormalLstSq * = nullptr);
cMainNetwork(eModeSSR aMode,cRect2,bool WithSchurr,cParamSparseNormalLstSq * = nullptr, tECCI * =nullptr);
~cMainNetwork();

//int N() const;
Expand All @@ -106,6 +110,10 @@ template <class Type> class cMainNetwork
/// If we use this iteration for covariance calculation , we dont add constraint, and dont solve
Type OneItereCompensation(bool ForCovCalc);

Type CalcResidual();
void AddGaugeConstraint(Type aWeight);


/// Access to CurSol of mSys
const Type & CurSol(int aK) const;

Expand All @@ -123,6 +131,8 @@ template <class Type> class cMainNetwork

/// Compute the geometry of an index using internal parameters => global simi + some random value
tPt Ind2Geom(const cPt2di & anInd) const;
/// Compute the geometry in case of cov propag
tPt CovPropInd2Geom(const cPt2di & anInd) const;

/** Classically for the gauge fixing the direction by fixing some specific var, we must take precaution
i.e if P0=(0,0) is fixed and P1=(1,0), if we fix x1=Cste for the gauge, the
Expand All @@ -132,13 +142,15 @@ template <class Type> class cMainNetwork
bool AxeXIsHoriz() const;

const cSim2D<Type> & SimInd2G() const; ///<Accessor
tSys * Sys();

void TestCov();

protected :
/// Acces to reference of a adress if point from pixel value
tPNetPtr & PNetPtrOfGrid(const cPt2di & aP) {return mMatrP[aP.y()-mBoxInd.P0().y()][aP.x()-mBoxInd.P0().x()];}

tECCI * mECCI;
cRect2 mBoxInd; ///< rectangle of the network
int mX_SzM; ///< 1+2*aN = Sz of Matrix of point
int mY_SzM; ///< 1+2*aN = Sz of Matrix of point
Expand Down
Loading

0 comments on commit 3a4ea9a

Please sign in to comment.