Index: trunk/SophyaLib/NTools/slinparbuff.h
===================================================================
--- trunk/SophyaLib/NTools/slinparbuff.h	(revision 2322)
+++ trunk/SophyaLib/NTools/slinparbuff.h	(revision 3572)
@@ -1,8 +1,2 @@
-#include "machdefs.h"
-#include <iostream>
-#include <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-#include <algorithm>
 #include "pexceptions.h"
 
@@ -17,307 +11,29 @@
   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(uint_4 lenbuf,uint_4 nresynch=0,r_8 x0=0.,r_8 y0=0.,bool autoxy0=false);
+  SLinParBuff(SLinParBuff& slp);
+  SLinParBuff(void);
+  virtual ~SLinParBuff();
 
-  SLinParBuff(SLinParBuff& slp)
-  // ---- Createur par copie
-  {
-    mLenBuf = slp.mLenBuf;
-    mNResynch = slp.mNResynch;
+  void GetX0Y0(r_8& x0,r_8& y0);
+  bool GetAutoX0Y0();
+  void SetX0Y0(r_8 x0=0.,r_8 y0=0.);
+  void SetAutoX0Y0(bool autoxy0=false);
 
-    mX = mY = NULL;
-    if(mLenBuf) {
-      mX = new r_8[mLenBuf]; mY = new r_8[mLenBuf];
-      for(uint_4 i=0;i<mLenBuf;i++) {mX[i]=slp.mX[i]; mY[i]=slp.mY[i];}
-    }
-
-    mX0 = slp.mX0; mY0 = slp.mY0; mAutoXY0 = slp.mAutoXY0;
-    mNCur = slp.mNCur; mIDeb = slp.mIDeb; mIResynch = slp.mIResynch;
-    mNResynchEff = slp.mNResynchEff; mNPush = slp.mNPush; mNPop = slp.mNPop;
-
-    mSx = slp.mSx; mSy = slp.mSy;
-    mSx2 = slp.mSx2; mSy2 = slp.mSy2; mSxy = slp.mSxy;
-    mSx3 = slp.mSx3; mSx2y = slp.mSx2y;
-    mSx4 = slp.mSx4;
-  }
-
-  SLinParBuff(void)
-  // ---- Createur par defaut
-  {
-    mLenBuf = 0; mNResynch = 0;
-    mX = mY = NULL;
-    mX0 = mY0 = 0.; mAutoXY0 = false;
-    Reset();
-  }
-
-  virtual ~SLinParBuff()
-  // ---- Destructeur
-  {
-    if(mX) delete [] mX; if(mY) delete [] mY;
-  }
-
-  //***************************************************************** 
-  inline void GetX0Y0(r_8& x0,r_8& y0) {x0 = mX0; y0 = mY0;}
-  inline bool GetAutoX0Y0() {return mAutoXY0;}
-
-  inline void SetX0Y0(r_8 x0=0.,r_8 y0=0.)
-  // ---- Stabilite numerique, centrage des x,y en mX0,mY0
-  {
-    mX0 = x0; mY0 = y0; mAutoXY0 = false;
-    ReComputeSum();
-  }
-
-  inline void SetAutoX0Y0(bool autoxy0=false)
-  // ---- Stabilite numerique, centrage automatique a <x> et <y>
-  //      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<mLenBuf;i++) mX[i] = 0.;
-    //if(mY) for(uint_4 i=0;i<mLenBuf;i++) mY[i] = 0.;
-  }
-
-  //***************************************************************** 
-  inline uint_4 Pop(void)
-  // ---- Pour enlever la donnee la plus ancienne
-  {
-    if(mNCur==0) return mNCur;
-    r_8 x=mX[mIDeb]-mX0, y=mY[mIDeb]-mY0; r_8 x2 = x*x;
-    mSx  -= x;     mSy   -= y;
-    mSx2 -= x2;    mSy2  -= y*y;   mSxy -= x*y;
-    mSx3 -= x2*x;  mSx2y -= x2*y;
-    mSx4 -= x2*x2;
-    mNCur--; mNPop++;
-    if(mNCur==0) {Reset(); return mNCur;}
-    mIDeb++; if(mIDeb==mLenBuf) mIDeb=0;
-    return mNCur;
-  }
-
-  inline uint_4 Push(r_8 x,r_8 y)
-  // ---- Pour ajouter une donnee, la donnee la plus ancienne
-  //      est enlevee si le buffer est deja plein.
-  {
-    uint_4 ifill;
-    if(mNCur==mLenBuf) {
-      r_8 X=mX[mIDeb]-mX0, Y=mY[mIDeb]-mY0; r_8 X2 = X*X;
-      mSx  -= X;     mSy   -=Y;
-      mSx2 -= X2;    mSy2  -=Y*Y;   mSxy -= X*Y;
-      mSx3 -= X2*X;  mSx2y -=X*X*Y;
-      mSx4 -= X2*X2;
-      ifill = mIDeb;
-      mIDeb++; if(mIDeb==mLenBuf) mIDeb=0;
-    } else {
-      ifill = (mIDeb+mNCur)%mLenBuf;
-      mNCur++;
-    }
-    mX[ifill] = x; mY[ifill] = y;
-    x -= mX0; y -= mY0; r_8 x2 = x*x;
-    mSx += x; mSy += y;
-    mSx2 += x2; mSy2 += y*y; mSxy += x*y;
-    mSx3 += x2*x; mSx2y += x2*y;
-    mSx4 += x2*x2;
-    mNPush++;
-    // Il faut re-synchroniser pour eviter les derives numeriques
-    mIResynch++; if(mIResynch == mNResynch) ReComputeSum();
-    return mNCur;
-  }
-
-  //***************************************************************** 
-  inline uint_4 ReComputeSum(void)
-  // ---- Pour re-synchroniser (eviter les derives numeriques).
-  {
-    if(mNCur==0) return 0;
-    // Re-centrage automatique a la moyenne demande:
-    // Attention, mSx = sum(x-mX0) ==> 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<mIDeb+mNCur;i++) {
-      uint_4 ii = i%mLenBuf;
-      r_8 x=mX[ii]-mX0, y=mY[ii]-mY0; r_8 x2 = x*x;
-      mSx  += x;     mSy   += y;
-      mSx2 += x2;    mSy2  += y*y;   mSxy += x*y;
-      mSx3 += x2*x;  mSx2y += x2*y;
-      mSx4 += x2*x2;
-    }
-    mIResynch=0;
-    mNResynchEff++;
-    return mNCur;
-  }
-
-  //***************************************************************** 
-  // ---- Retourne le nombre de points
+  void Reset(void);
+  uint_4 Pop(void);
+  uint_4 Push(r_8 x,r_8 y);
   inline uint_4 NPoints(void) {return mNCur;}
 
-  inline r_8 Compute(r_8& mean,bool recomputeXi2=false)
-  // ---- Calcul <y>, 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;i<mIDeb+mNCur;i++) {
-        uint_4 ii = i%mLenBuf;
-        r_8 s = mean - (mY[ii]-mY0);
-	sigma += s*s;
-      }
-      sigma /= (r_8) mNCur;
-    } else {
-      sigma = mSy2/(r_8) mNCur - mean*mean;
-    }
-    // gestion du decalage
-    mean += mY0;
+  r_8 SumX(void) {return mSx+mNCur*mX0;}
+  r_8 SumY(void) {return mSy+mNCur*mY0;}
 
-    if(sigma>0.) return sqrt(sigma);
-    if(sigma<0.) return -sqrt(-sigma);
-    return sigma;
-  }
+  uint_4 ReComputeSum(void);
+  r_8 Compute(r_8& mean,bool recomputeXi2=false);
+  r_8 Compute(r_8& a0,r_8 &a1,bool recomputeXi2=false);
+  r_8 Compute(r_8& a0,r_8 &a1,r_8 &a2,bool recomputeXi2=false);
 
-  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(<dy^2>)
-  //   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;i<mIDeb+mNCur;i++) {
-        uint_4 ii = i%mLenBuf;
-        r_8 s = a0+a1*(mX[ii]-mX0) - (mY[ii]-mY0);
-	slin += s*s;
-      }
-      slin /= (r_8) mNCur;
-    } else {
-      slin = (mSy2 +a0*a0*mNCur +a1*a1*mSx2 -2.*a0*mSy -2.*a1*mSxy +2.*a0*a1*mSx)
-             / (r_8)mNCur;
-    }
-    // gestion du decalage y-y0 = a0 + a1*(x-x0)
-    a0 = mY0 + a0 - a1*mX0;
-
-    if(slin>0.) 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(<dy^2>)
-  //   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;i<mIDeb+mNCur;i++) {
-        uint_4 ii = i%mLenBuf;
-        r_8 s = a0+(a1+a2*(mX[ii]-mX0))*(mX[ii]-mX0) - (mY[ii]-mY0);
-	spar += s*s;
-      }
-      spar /= (r_8) mNCur;
-    } else {
-      spar = (mSy2 +a0*a0*mNCur +a1*a1*mSx2 +a2*a2*mSx4 -2.*mSy*a0 -2.*a1*mSxy
-              -2.*a2*mSx2y +2.*a0*a1*mSx +2.*a0*a2*mSx2 +2.*a1*a2*mSx3)
-             / (r_8) mNCur;
-    }
-    // gestion du decalage y-y0 = a0 + a1*(x-x0) + a2*(x-x0)^2
-    a0 = mY0 + a0 - a1*mX0 + a2*mX0*mX0;
-    a1 = a1 - 2.*a2*mX0;
-
-    if(spar>0.) 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="<<mLenBuf<<",NResynch="<<mNResynch
-        <<",auto_xy0="<<mAutoXY0
-        <<"): mX0="<<mX0<<" mY0="<<mY0<<" mNCur="<<mNCur<<endl
-        <<"           NPush="<<mNPush<<" NPop="<<mNPop
-        <<" NSynchro="<<mNResynchEff<<endl;
-    if(mNCur==0) return;
-    if(lp>=2) {
-      cout<<"X:";
-      if(mX) for(uint_4 i=0;i<mNCur;i++) {
-        uint_4 ii = (mIDeb+i)%mLenBuf;
-        cout<<" "<<mX[ii]; if(i%10==9 || i==mNCur-1) cout<<endl;
-      }
-      if(mY) cout<<"Y:";
-      for(uint_4 i=0;i<mNCur;i++) {
-        uint_4 ii = (mIDeb+i)%mLenBuf;
-        cout<<" "<<mY[ii]; if(i%10==9 || i==mNCur-1) cout<<endl;
-      }
-    }
-    if(lp>=1) {
-      cout<<"...IDeb="<<mIDeb<<" IResynch="<<mIResynch<<endl;
-      cout<<"   Sx="<<mSx<<" Sx2="<<mSx2<<" Sx3="<<mSx3<<" Sx4="<<mSx4<<endl
-          <<"   Sy="<<mSy<<" Sy2="<<mSy2<<" Sxy="<<mSxy<<" Sx2y="<<mSx2y<<endl;
-    }
-  }
-
-  void PrintCompute(int lp=0)
-  // ---- Print des valeurs numeriques calculees
-  //      lp = 0 : valeurs calculees
-  //           1 :             + sommes
-  {
-    bool recompute = (mX==NULL || mY==NULL) ? false: true;
-    r_8 mean;
-    r_8 sigma = Compute(mean,recompute), sigmar = Compute(mean,false);
-    cout<<"SLinParBuff: n="<<NPoints()<<" mean="<<mean
-        <<" sigma="<<sigma<<" (raw="<<sigmar<<")"<<endl;
-    r_8 a0,a1,a2;
-    r_8 slin  = Compute(a0,a1,recompute), slinr  = Compute(a0,a1,false);
-    cout<<"             a0="<<a0<<" a1="<<a1
-        <<" slin="<<slin<<" (raw="<<slinr<<")"<<endl;
-    r_8 spar = Compute(a0,a1,a2,recompute), sparr = Compute(a0,a1,a2,false);
-    cout<<"             a0="<<a0<<" a1="<<a1<<" a2="<<a2
-        <<" spar="<<spar<<" (raw="<<sparr<<")"<<endl;
-    if(lp<1) return;
-    cout<<"...Sx="<<mSx<<" Sx2="<<mSx2<<" Sx3="<<mSx3<<" Sx4="<<mSx4<<endl
-        <<"   Sy="<<mSy<<" Sy2="<<mSy2<<" Sxy="<<mSxy<<" Sx2y="<<mSx2y<<endl;
-  }
+  void Print(int lp=0);
+  void PrintCompute(int lp=0);
 
 protected:
@@ -329,4 +45,5 @@
   uint_4 mNResynchEff, mNPush, mNPop;
 };
+
 
 ///////////////////////////////////////////////////////////////////
@@ -343,46 +60,14 @@
 class SLinParBuffMerger {
 public:
-  SLinParBuffMerger(void) {mFirst=true; mSlp.Reset();}
-  SLinParBuffMerger(SLinParBuff& s,bool recompute=false)
-                   {mFirst=true; mSlp.Reset(); Add(s,recompute);}
-  virtual ~SLinParBuffMerger(void) {}
-  inline uint_4  NPoints(void) {return mSlp.NPoints();}
+  SLinParBuffMerger(void);
+  SLinParBuffMerger(SLinParBuff& s,bool recompute=false);
+  virtual ~SLinParBuffMerger(void);
+
+  inline uint_4 NPoints(void) {return mSlp.NPoints();}
   inline void Reset(void) {mSlp.Reset(); mFirst=true;}
-  inline void Add(SLinParBuff& s,bool recompute=false)
-  {
-    bool changex0y0=false, AutoXY0_Save; r_8 X0_Save,Y0_Save;
-    if(mFirst) {
-      // Cas ou c'est le premier SLinParBuff additionne.
-      mSlp.mX0=s.mX0; mSlp.mY0=s.mY0; mFirst=false;
-    } else if(mSlp.mX0!=s.mX0 || mSlp.mY0!=s.mY0) {
-      // Attention: pour merger il faut avoir les memes mX0,mY0
-      changex0y0=true;
-      X0_Save=s.mX0; Y0_Save=s.mY0; AutoXY0_Save=s.mAutoXY0;
-      s.mX0=mSlp.mX0; s.mY0=mSlp.mY0; s.mAutoXY0=false;
-      recompute=true;
-    }
-    if(recompute) s.ReComputeSum();
-    mSlp.mNCur += s.mNCur;
-    mSlp.mSx   += s.mSx;
-    mSlp.mSy   += s.mSy;
-    mSlp.mSx2  += s.mSx2;
-    mSlp.mSy2  += s.mSy2;
-    mSlp.mSxy  += s.mSxy;
-    mSlp.mSx3  += s.mSx3;
-    mSlp.mSx2y += s.mSx2y;
-    mSlp.mSx4  += s.mSx4;
-    // Dans le cas ou on a change les X0,Y0, on remet en etat
-    if(changex0y0) {
-      s.mX0=X0_Save; s.mY0=Y0_Save;
-      s.ReComputeSum();
-      s.mAutoXY0=AutoXY0_Save;
-    }
-  }
-  inline r_8 Compute(r_8& mean)
-              {return mSlp.Compute(mean);}
-  inline r_8 Compute(r_8& a0,r_8 &a1)
-              {return mSlp.Compute(a0,a1);}
-  inline r_8 Compute(r_8& a0,r_8 &a1,r_8 &a2)
-              {return mSlp.Compute(a0,a1,a2);}
+  void Add(SLinParBuff& s,bool recompute=false);
+  inline r_8 Compute(r_8& mean) {return mSlp.Compute(mean);}
+  inline r_8 Compute(r_8& a0,r_8 &a1) {return mSlp.Compute(a0,a1);}
+  inline r_8 Compute(r_8& a0,r_8 &a1,r_8 &a2) {return mSlp.Compute(a0,a1,a2);}
   inline void Print(int lp=0) {mSlp.Print(lp);}
   inline void PrintCompute(int lp=0) {mSlp.PrintCompute(lp);}
