#include "machdefs.h" #include #include #include #include #include #include "pexceptions.h" #ifndef SLINPARBUF_H_SEEN #define SLINPARBUF_H_SEEN namespace SOPHYA { ////////////////////////////////////////////////////////////////////////////////// class SLinParBuff { public: friend class SLinParBuffMerger; //***************************************************************** SLinParBuff(uint_4 lenbuf,uint_4 nresynch=0 ,r_8 x0=0.,r_8 y0=0.,bool autoxy0=false) // ---- Createur { if(lenbuf==0) throw RangeCheckError("SLinParBuff: lenbuf==0 !"); mLenBuf = lenbuf; mNResynch = nresynch; mX = new r_8[mLenBuf]; mY = new r_8[mLenBuf]; mX0 = x0; mY0 = y0; mAutoXY0 = autoxy0; Reset(); } SLinParBuff(SLinParBuff& slp) // ---- Createur par copie { mLenBuf = slp.mLenBuf; mNResynch = slp.mNResynch; mX = mY = NULL; if(mLenBuf) { mX = new r_8[mLenBuf]; mY = new r_8[mLenBuf]; for(uint_4 i=0;i et // a chaque ReComputeSum { mAutoXY0 = autoxy0; ReComputeSum(); } inline void Reset(void) // ---- Reset des sommes et des compteurs { mSx=mSy=mSx2=mSy2=mSxy=mSx3=mSx2y=mSx4=0.; mNCur = mIDeb = mIResynch = mNResynchEff = mNPush = mNPop = 0; //if(mX) for(uint_4 i=0;i nouvel mX0 = mSx/N + mX0(ancien) if(mAutoXY0) {mX0 = mSx/(r_8)mNCur + mX0; mY0 = mSy/(r_8)mNCur + mY0;} mSx=mSy=mSx2=mSy2=mSxy=mSx3=mSx2y=mSx4=0.; for(uint_4 i=mIDeb;i, Var(y) // recomputeXi2=true : recalcule le xi2 (sigma) avec la courbe et les points { mean=0.; if(mNCur==0) return -1.; // Moyenne mean = mSy/(r_8) mNCur; // Sigma r_8 sigma; if(recomputeXi2) { sigma=0.; for(uint_4 i=mIDeb;i0.) return sqrt(sigma); if(sigma<0.) return -sqrt(-sigma); return sigma; } inline r_8 Compute(r_8& a0,r_8 &a1,bool recomputeXi2=false) // ---- Calcul y=a0+a1*x et slin=Var(y-(a0+a1*x))=sqrt() // recomputeXi2=true : recalcule le xi2 avec la courbe et les points { a0=a1=0.; if(mNCur==0) return -1.; // Fit lineaire r_8 w = mNCur*mSx2 - mSx*mSx; if(w==0. || mNCur==1) return -2.; a0 = (mSx2*mSy - mSx*mSxy)/w; a1 = (mNCur*mSxy - mSx*mSy )/w; // Sigma par rapport Fit lineaire // (On a XI2=mNCur*slin**2 dans notre cas ou les erreurs=1) r_8 slin; if(recomputeXi2) { slin=0.; for(uint_4 i=mIDeb;i0.) return sqrt(slin); if(slin<0.) return -sqrt(-slin); return slin; } inline r_8 Compute(r_8& a0,r_8 &a1,r_8 &a2,bool recomputeXi2=false) // ---- Calcul y=a0+a1*x+a2*x^2 et spar=Var(y-(a0+a1*x+a2*x^2))=sqrt() // recomputeXi2=true : recalcule le xi2 avec la courbe et les points { a0=a1=a2=0.; if(mNCur==0) return -1.; // Fit parabolique r_8 w = mSx2*(mSx2*mSx2-mSx3*mSx) -mSx*(mSx3*mSx2-mSx4*mSx) +mNCur*(mSx3*mSx3-mSx4*mSx2); if(w==0. || mNCur<=2) return -2.; a2 = (mSy*(mSx2*mSx2-mSx3*mSx) - mSxy*(mSx*mSx2-mSx3*mNCur) + mSx2y*(mSx*mSx-mSx2*mNCur) )/w; a1 = -(mSy*(mSx3*mSx2-mSx4*mSx) - mSxy*(mSx2*mSx2-mSx4*mNCur) + mSx2y*(mSx2*mSx-mSx3*mNCur))/w; a0 = (mSy*(mSx3*mSx3-mSx4*mSx2) - mSxy*(mSx2*mSx3-mSx4*mSx) + mSx2y*(mSx2*mSx2-mSx3*mSx) )/w; // Sigma par rapport Fit parabolique // (On a XI2=mNCur*spar**2 dans notre cas ou les erreurs=1) // Le calcul direct du Xi2 n'est pas precis r_8 spar; if(recomputeXi2) { spar=0.; for(uint_4 i=mIDeb;i0.) return sqrt(spar); if(spar<0.) return -sqrt(-spar); return spar; } //***************************************************************** void Print(int lp=0) // ---- Print de l'etat de la classe // lp = 0 : parametres // 1 : + sommes // 2 : + tableaux { cout<<"SLinParBuff(LenBuf="<=1) { cout<<"...IDeb="<