source: Sophya/trunk/SophyaLib/NTools/slinparbuff.h@ 3309

Last change on this file since 3309 was 2322, checked in by cmv, 23 years ago
  • passage xxstream.h en xxstream
  • compile avec gcc_3.2, gcc_2.96 et cxx En 3.2 le seek from ::end semble marcher (voir Eval/COS/pbseekios.cc)

rz+cmv 11/2/2003
File size: 12.4 KB
RevLine 
[1674]1#include "machdefs.h"
[2322]2#include <iostream>
[1674]3#include <stdlib.h>
4#include <stdio.h>
5#include <math.h>
6#include <algorithm>
7#include "pexceptions.h"
8
9#ifndef SLINPARBUF_H_SEEN
10#define SLINPARBUF_H_SEEN
11
12namespace SOPHYA {
13
14//////////////////////////////////////////////////////////////////////////////////
15class SLinParBuff {
16public:
17 friend class SLinParBuffMerger;
18
19 //*****************************************************************
20 SLinParBuff(uint_4 lenbuf,uint_4 nresynch=0
21 ,r_8 x0=0.,r_8 y0=0.,bool autoxy0=false)
22 // ---- Createur
23 {
24 if(lenbuf==0) throw RangeCheckError("SLinParBuff: lenbuf==0 !");
25 mLenBuf = lenbuf;
26 mNResynch = nresynch;
27 mX = new r_8[mLenBuf]; mY = new r_8[mLenBuf];
28 mX0 = x0; mY0 = y0; mAutoXY0 = autoxy0;
29 Reset();
30 }
31
32 SLinParBuff(SLinParBuff& slp)
33 // ---- Createur par copie
34 {
35 mLenBuf = slp.mLenBuf;
36 mNResynch = slp.mNResynch;
37
38 mX = mY = NULL;
39 if(mLenBuf) {
40 mX = new r_8[mLenBuf]; mY = new r_8[mLenBuf];
41 for(uint_4 i=0;i<mLenBuf;i++) {mX[i]=slp.mX[i]; mY[i]=slp.mY[i];}
42 }
43
44 mX0 = slp.mX0; mY0 = slp.mY0; mAutoXY0 = slp.mAutoXY0;
45 mNCur = slp.mNCur; mIDeb = slp.mIDeb; mIResynch = slp.mIResynch;
46 mNResynchEff = slp.mNResynchEff; mNPush = slp.mNPush; mNPop = slp.mNPop;
47
48 mSx = slp.mSx; mSy = slp.mSy;
49 mSx2 = slp.mSx2; mSy2 = slp.mSy2; mSxy = slp.mSxy;
50 mSx3 = slp.mSx3; mSx2y = slp.mSx2y;
51 mSx4 = slp.mSx4;
52 }
53
54 SLinParBuff(void)
55 // ---- Createur par defaut
56 {
57 mLenBuf = 0; mNResynch = 0;
58 mX = mY = NULL;
59 mX0 = mY0 = 0.; mAutoXY0 = false;
60 Reset();
61 }
62
63 virtual ~SLinParBuff()
64 // ---- Destructeur
65 {
66 if(mX) delete [] mX; if(mY) delete [] mY;
67 }
68
69 //*****************************************************************
70 inline void GetX0Y0(r_8& x0,r_8& y0) {x0 = mX0; y0 = mY0;}
71 inline bool GetAutoX0Y0() {return mAutoXY0;}
72
73 inline void SetX0Y0(r_8 x0=0.,r_8 y0=0.)
74 // ---- Stabilite numerique, centrage des x,y en mX0,mY0
75 {
76 mX0 = x0; mY0 = y0; mAutoXY0 = false;
77 ReComputeSum();
78 }
79
80 inline void SetAutoX0Y0(bool autoxy0=false)
81 // ---- Stabilite numerique, centrage automatique a <x> et <y>
82 // a chaque ReComputeSum
83 {
84 mAutoXY0 = autoxy0;
85 ReComputeSum();
86 }
87
88 inline void Reset(void)
89 // ---- Reset des sommes et des compteurs
90 {
91 mSx=mSy=mSx2=mSy2=mSxy=mSx3=mSx2y=mSx4=0.;
92 mNCur = mIDeb = mIResynch = mNResynchEff = mNPush = mNPop = 0;
93 //if(mX) for(uint_4 i=0;i<mLenBuf;i++) mX[i] = 0.;
94 //if(mY) for(uint_4 i=0;i<mLenBuf;i++) mY[i] = 0.;
95 }
96
97 //*****************************************************************
98 inline uint_4 Pop(void)
99 // ---- Pour enlever la donnee la plus ancienne
100 {
101 if(mNCur==0) return mNCur;
102 r_8 x=mX[mIDeb]-mX0, y=mY[mIDeb]-mY0; r_8 x2 = x*x;
103 mSx -= x; mSy -= y;
104 mSx2 -= x2; mSy2 -= y*y; mSxy -= x*y;
105 mSx3 -= x2*x; mSx2y -= x2*y;
106 mSx4 -= x2*x2;
107 mNCur--; mNPop++;
108 if(mNCur==0) {Reset(); return mNCur;}
109 mIDeb++; if(mIDeb==mLenBuf) mIDeb=0;
110 return mNCur;
111 }
112
113 inline uint_4 Push(r_8 x,r_8 y)
114 // ---- Pour ajouter une donnee, la donnee la plus ancienne
115 // est enlevee si le buffer est deja plein.
116 {
117 uint_4 ifill;
118 if(mNCur==mLenBuf) {
119 r_8 X=mX[mIDeb]-mX0, Y=mY[mIDeb]-mY0; r_8 X2 = X*X;
120 mSx -= X; mSy -=Y;
121 mSx2 -= X2; mSy2 -=Y*Y; mSxy -= X*Y;
122 mSx3 -= X2*X; mSx2y -=X*X*Y;
123 mSx4 -= X2*X2;
124 ifill = mIDeb;
125 mIDeb++; if(mIDeb==mLenBuf) mIDeb=0;
126 } else {
127 ifill = (mIDeb+mNCur)%mLenBuf;
128 mNCur++;
129 }
130 mX[ifill] = x; mY[ifill] = y;
131 x -= mX0; y -= mY0; r_8 x2 = x*x;
132 mSx += x; mSy += y;
133 mSx2 += x2; mSy2 += y*y; mSxy += x*y;
134 mSx3 += x2*x; mSx2y += x2*y;
135 mSx4 += x2*x2;
136 mNPush++;
137 // Il faut re-synchroniser pour eviter les derives numeriques
138 mIResynch++; if(mIResynch == mNResynch) ReComputeSum();
139 return mNCur;
140 }
141
142 //*****************************************************************
143 inline uint_4 ReComputeSum(void)
144 // ---- Pour re-synchroniser (eviter les derives numeriques).
145 {
146 if(mNCur==0) return 0;
147 // Re-centrage automatique a la moyenne demande:
148 // Attention, mSx = sum(x-mX0) ==> nouvel mX0 = mSx/N + mX0(ancien)
149 if(mAutoXY0) {mX0 = mSx/(r_8)mNCur + mX0; mY0 = mSy/(r_8)mNCur + mY0;}
150 mSx=mSy=mSx2=mSy2=mSxy=mSx3=mSx2y=mSx4=0.;
151 for(uint_4 i=mIDeb;i<mIDeb+mNCur;i++) {
152 uint_4 ii = i%mLenBuf;
153 r_8 x=mX[ii]-mX0, y=mY[ii]-mY0; r_8 x2 = x*x;
154 mSx += x; mSy += y;
155 mSx2 += x2; mSy2 += y*y; mSxy += x*y;
156 mSx3 += x2*x; mSx2y += x2*y;
157 mSx4 += x2*x2;
158 }
159 mIResynch=0;
160 mNResynchEff++;
161 return mNCur;
162 }
163
164 //*****************************************************************
165 // ---- Retourne le nombre de points
166 inline uint_4 NPoints(void) {return mNCur;}
167
168 inline r_8 Compute(r_8& mean,bool recomputeXi2=false)
169 // ---- Calcul <y>, Var(y)
170 // recomputeXi2=true : recalcule le xi2 (sigma) avec la courbe et les points
171 {
172 mean=0.;
173 if(mNCur==0) return -1.;
174 // Moyenne
175 mean = mSy/(r_8) mNCur;
176 // Sigma
177 r_8 sigma;
178 if(recomputeXi2) {
179 sigma=0.;
180 for(uint_4 i=mIDeb;i<mIDeb+mNCur;i++) {
181 uint_4 ii = i%mLenBuf;
182 r_8 s = mean - (mY[ii]-mY0);
183 sigma += s*s;
184 }
185 sigma /= (r_8) mNCur;
186 } else {
187 sigma = mSy2/(r_8) mNCur - mean*mean;
188 }
189 // gestion du decalage
190 mean += mY0;
191
192 if(sigma>0.) return sqrt(sigma);
193 if(sigma<0.) return -sqrt(-sigma);
194 return sigma;
195 }
196
197 inline r_8 Compute(r_8& a0,r_8 &a1,bool recomputeXi2=false)
198 // ---- Calcul y=a0+a1*x et slin=Var(y-(a0+a1*x))=sqrt(<dy^2>)
199 // recomputeXi2=true : recalcule le xi2 avec la courbe et les points
200 {
201 a0=a1=0.;
202 if(mNCur==0) return -1.;
203 // Fit lineaire
204 r_8 w = mNCur*mSx2 - mSx*mSx;
205 if(w==0. || mNCur==1) return -2.;
206 a0 = (mSx2*mSy - mSx*mSxy)/w;
207 a1 = (mNCur*mSxy - mSx*mSy )/w;
208 // Sigma par rapport Fit lineaire
209 // (On a XI2=mNCur*slin**2 dans notre cas ou les erreurs=1)
210 r_8 slin;
211 if(recomputeXi2) {
212 slin=0.;
213 for(uint_4 i=mIDeb;i<mIDeb+mNCur;i++) {
214 uint_4 ii = i%mLenBuf;
215 r_8 s = a0+a1*(mX[ii]-mX0) - (mY[ii]-mY0);
216 slin += s*s;
217 }
218 slin /= (r_8) mNCur;
219 } else {
220 slin = (mSy2 +a0*a0*mNCur +a1*a1*mSx2 -2.*a0*mSy -2.*a1*mSxy +2.*a0*a1*mSx)
221 / (r_8)mNCur;
222 }
223 // gestion du decalage y-y0 = a0 + a1*(x-x0)
224 a0 = mY0 + a0 - a1*mX0;
225
226 if(slin>0.) return sqrt(slin);
227 if(slin<0.) return -sqrt(-slin);
228 return slin;
229 }
230
231 inline r_8 Compute(r_8& a0,r_8 &a1,r_8 &a2,bool recomputeXi2=false)
232 // ---- Calcul y=a0+a1*x+a2*x^2 et spar=Var(y-(a0+a1*x+a2*x^2))=sqrt(<dy^2>)
233 // recomputeXi2=true : recalcule le xi2 avec la courbe et les points
234 {
235 a0=a1=a2=0.;
236 if(mNCur==0) return -1.;
237 // Fit parabolique
238 r_8 w = mSx2*(mSx2*mSx2-mSx3*mSx) -mSx*(mSx3*mSx2-mSx4*mSx) +mNCur*(mSx3*mSx3-mSx4*mSx2);
239 if(w==0. || mNCur<=2) return -2.;
240 a2 = (mSy*(mSx2*mSx2-mSx3*mSx) - mSxy*(mSx*mSx2-mSx3*mNCur) + mSx2y*(mSx*mSx-mSx2*mNCur) )/w;
241 a1 = -(mSy*(mSx3*mSx2-mSx4*mSx) - mSxy*(mSx2*mSx2-mSx4*mNCur) + mSx2y*(mSx2*mSx-mSx3*mNCur))/w;
242 a0 = (mSy*(mSx3*mSx3-mSx4*mSx2) - mSxy*(mSx2*mSx3-mSx4*mSx) + mSx2y*(mSx2*mSx2-mSx3*mSx) )/w;
243 // Sigma par rapport Fit parabolique
244 // (On a XI2=mNCur*spar**2 dans notre cas ou les erreurs=1)
245 // Le calcul direct du Xi2 n'est pas precis
246 r_8 spar;
247 if(recomputeXi2) {
248 spar=0.;
249 for(uint_4 i=mIDeb;i<mIDeb+mNCur;i++) {
250 uint_4 ii = i%mLenBuf;
251 r_8 s = a0+(a1+a2*(mX[ii]-mX0))*(mX[ii]-mX0) - (mY[ii]-mY0);
252 spar += s*s;
253 }
254 spar /= (r_8) mNCur;
255 } else {
256 spar = (mSy2 +a0*a0*mNCur +a1*a1*mSx2 +a2*a2*mSx4 -2.*mSy*a0 -2.*a1*mSxy
257 -2.*a2*mSx2y +2.*a0*a1*mSx +2.*a0*a2*mSx2 +2.*a1*a2*mSx3)
258 / (r_8) mNCur;
259 }
260 // gestion du decalage y-y0 = a0 + a1*(x-x0) + a2*(x-x0)^2
261 a0 = mY0 + a0 - a1*mX0 + a2*mX0*mX0;
262 a1 = a1 - 2.*a2*mX0;
263
264 if(spar>0.) return sqrt(spar);
265 if(spar<0.) return -sqrt(-spar);
266 return spar;
267 }
268
269 //*****************************************************************
270 void Print(int lp=0)
271 // ---- Print de l'etat de la classe
272 // lp = 0 : parametres
273 // 1 : + sommes
274 // 2 : + tableaux
275 {
276 cout<<"SLinParBuff(LenBuf="<<mLenBuf<<",NResynch="<<mNResynch
277 <<",auto_xy0="<<mAutoXY0
278 <<"): mX0="<<mX0<<" mY0="<<mY0<<" mNCur="<<mNCur<<endl
279 <<" NPush="<<mNPush<<" NPop="<<mNPop
280 <<" NSynchro="<<mNResynchEff<<endl;
281 if(mNCur==0) return;
282 if(lp>=2) {
283 cout<<"X:";
284 if(mX) for(uint_4 i=0;i<mNCur;i++) {
285 uint_4 ii = (mIDeb+i)%mLenBuf;
286 cout<<" "<<mX[ii]; if(i%10==9 || i==mNCur-1) cout<<endl;
287 }
288 if(mY) cout<<"Y:";
289 for(uint_4 i=0;i<mNCur;i++) {
290 uint_4 ii = (mIDeb+i)%mLenBuf;
291 cout<<" "<<mY[ii]; if(i%10==9 || i==mNCur-1) cout<<endl;
292 }
293 }
294 if(lp>=1) {
295 cout<<"...IDeb="<<mIDeb<<" IResynch="<<mIResynch<<endl;
296 cout<<" Sx="<<mSx<<" Sx2="<<mSx2<<" Sx3="<<mSx3<<" Sx4="<<mSx4<<endl
297 <<" Sy="<<mSy<<" Sy2="<<mSy2<<" Sxy="<<mSxy<<" Sx2y="<<mSx2y<<endl;
298 }
299 }
300
301 void PrintCompute(int lp=0)
302 // ---- Print des valeurs numeriques calculees
303 // lp = 0 : valeurs calculees
304 // 1 : + sommes
305 {
306 bool recompute = (mX==NULL || mY==NULL) ? false: true;
307 r_8 mean;
308 r_8 sigma = Compute(mean,recompute), sigmar = Compute(mean,false);
309 cout<<"SLinParBuff: n="<<NPoints()<<" mean="<<mean
310 <<" sigma="<<sigma<<" (raw="<<sigmar<<")"<<endl;
311 r_8 a0,a1,a2;
312 r_8 slin = Compute(a0,a1,recompute), slinr = Compute(a0,a1,false);
313 cout<<" a0="<<a0<<" a1="<<a1
314 <<" slin="<<slin<<" (raw="<<slinr<<")"<<endl;
315 r_8 spar = Compute(a0,a1,a2,recompute), sparr = Compute(a0,a1,a2,false);
316 cout<<" a0="<<a0<<" a1="<<a1<<" a2="<<a2
317 <<" spar="<<spar<<" (raw="<<sparr<<")"<<endl;
318 if(lp<1) return;
319 cout<<"...Sx="<<mSx<<" Sx2="<<mSx2<<" Sx3="<<mSx3<<" Sx4="<<mSx4<<endl
320 <<" Sy="<<mSy<<" Sy2="<<mSy2<<" Sxy="<<mSxy<<" Sx2y="<<mSx2y<<endl;
321 }
322
323protected:
324 uint_4 mLenBuf,mNResynch;
325 r_8 *mX,*mY, mX0, mY0;
326 bool mAutoXY0;
327 uint_4 mNCur,mIDeb,mIResynch;
328 r_8 mSx,mSy,mSx2,mSy2,mSxy, mSx3,mSx2y, mSx4;
329 uint_4 mNResynchEff, mNPush, mNPop;
330};
331
332///////////////////////////////////////////////////////////////////
333// That Class allows merging of SLinParBuff Class for computing //
334// parameters.
335// ex:
336// SLinParBuff s1; -> Fill s1
337// SLinParBuff s2; -> Fill s2
338// SLinParBuff s3; -> Fill s3
339// SLinParBuffMerger smerge(s1);
340// sig = smerge.Compute(mean); -> same as sig = s1.Compute(mean);
341// smerge.Add(s2);
342// sig = smerge.Compute(mean); -> sig and mean are those for
343class SLinParBuffMerger {
344public:
345 SLinParBuffMerger(void) {mFirst=true; mSlp.Reset();}
346 SLinParBuffMerger(SLinParBuff& s,bool recompute=false)
347 {mFirst=true; mSlp.Reset(); Add(s,recompute);}
348 virtual ~SLinParBuffMerger(void) {}
349 inline uint_4 NPoints(void) {return mSlp.NPoints();}
350 inline void Reset(void) {mSlp.Reset(); mFirst=true;}
351 inline void Add(SLinParBuff& s,bool recompute=false)
352 {
353 bool changex0y0=false, AutoXY0_Save; r_8 X0_Save,Y0_Save;
354 if(mFirst) {
355 // Cas ou c'est le premier SLinParBuff additionne.
356 mSlp.mX0=s.mX0; mSlp.mY0=s.mY0; mFirst=false;
357 } else if(mSlp.mX0!=s.mX0 || mSlp.mY0!=s.mY0) {
358 // Attention: pour merger il faut avoir les memes mX0,mY0
359 changex0y0=true;
360 X0_Save=s.mX0; Y0_Save=s.mY0; AutoXY0_Save=s.mAutoXY0;
361 s.mX0=mSlp.mX0; s.mY0=mSlp.mY0; s.mAutoXY0=false;
362 recompute=true;
363 }
364 if(recompute) s.ReComputeSum();
365 mSlp.mNCur += s.mNCur;
366 mSlp.mSx += s.mSx;
367 mSlp.mSy += s.mSy;
368 mSlp.mSx2 += s.mSx2;
369 mSlp.mSy2 += s.mSy2;
370 mSlp.mSxy += s.mSxy;
371 mSlp.mSx3 += s.mSx3;
372 mSlp.mSx2y += s.mSx2y;
373 mSlp.mSx4 += s.mSx4;
374 // Dans le cas ou on a change les X0,Y0, on remet en etat
375 if(changex0y0) {
376 s.mX0=X0_Save; s.mY0=Y0_Save;
377 s.ReComputeSum();
378 s.mAutoXY0=AutoXY0_Save;
379 }
380 }
381 inline r_8 Compute(r_8& mean)
382 {return mSlp.Compute(mean);}
383 inline r_8 Compute(r_8& a0,r_8 &a1)
384 {return mSlp.Compute(a0,a1);}
385 inline r_8 Compute(r_8& a0,r_8 &a1,r_8 &a2)
386 {return mSlp.Compute(a0,a1,a2);}
387 inline void Print(int lp=0) {mSlp.Print(lp);}
388 inline void PrintCompute(int lp=0) {mSlp.PrintCompute(lp);}
389
390protected:
391 bool mFirst;
392 SLinParBuff mSlp;
393};
394
395} // Fin du namespace
396
397#endif
Note: See TracBrowser for help on using the repository browser.