Changeset 3572 in Sophya for trunk/SophyaLib/NTools
- Timestamp:
- Feb 7, 2009, 10:50:34 PM (17 years ago)
- Location:
- trunk/SophyaLib/NTools
- Files:
-
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/NTools/Makefile
r3406 r3572 4 4 5 5 clean: 6 rm -f $(SOPHYAOBJP)FSAppIrrSmpl.o $(SOPHYAOBJP)cimage.o $(SOPHYAOBJP)cspline.o $(SOPHYAOBJP)datatypes.o $(SOPHYAOBJP)dates.o $(SOPHYAOBJP)difeq.o $(SOPHYAOBJP)dynccd.o $(SOPHYAOBJP)fct1dfit.o $(SOPHYAOBJP)fct2dfit.o $(SOPHYAOBJP)fftmserver.o $(SOPHYAOBJP)fftpserver.o $(SOPHYAOBJP)fftservintf.o $(SOPHYAOBJP)functab.o $(SOPHYAOBJP)generaldata.o $(SOPHYAOBJP)generalfit.o $(SOPHYAOBJP)generalfunc.o $(SOPHYAOBJP)imageop.o $(SOPHYAOBJP)integ.o $(SOPHYAOBJP)median.o $(SOPHYAOBJP)ntoolsinit.o $(SOPHYAOBJP)objfitter.o $(SOPHYAOBJP)perandom.o $(SOPHYAOBJP)poly.o $(SOPHYAOBJP)rk4cdifeq.o $(SOPHYAOBJP)simplesort.o $(SOPHYAOBJP)simplex.o $(SOPHYAOBJP) tabmath.o $(SOPHYAOBJP)toeplitzMatrix.o $(SOPHYAOBJP)datime.o $(SOPHYAOBJP)fftmayer_r4.o $(SOPHYAOBJP)fftmayer_r8.o $(SOPHYAOBJP)fftpackc.o $(SOPHYAOBJP)matxop.o $(SOPHYAOBJP)nbmath.o $(SOPHYAOBJP)nbtri.o6 rm -f $(SOPHYAOBJP)FSAppIrrSmpl.o $(SOPHYAOBJP)cimage.o $(SOPHYAOBJP)cspline.o $(SOPHYAOBJP)datatypes.o $(SOPHYAOBJP)dates.o $(SOPHYAOBJP)difeq.o $(SOPHYAOBJP)dynccd.o $(SOPHYAOBJP)fct1dfit.o $(SOPHYAOBJP)fct2dfit.o $(SOPHYAOBJP)fftmserver.o $(SOPHYAOBJP)fftpserver.o $(SOPHYAOBJP)fftservintf.o $(SOPHYAOBJP)functab.o $(SOPHYAOBJP)generaldata.o $(SOPHYAOBJP)generalfit.o $(SOPHYAOBJP)generalfunc.o $(SOPHYAOBJP)imageop.o $(SOPHYAOBJP)integ.o $(SOPHYAOBJP)median.o $(SOPHYAOBJP)ntoolsinit.o $(SOPHYAOBJP)objfitter.o $(SOPHYAOBJP)perandom.o $(SOPHYAOBJP)poly.o $(SOPHYAOBJP)rk4cdifeq.o $(SOPHYAOBJP)simplesort.o $(SOPHYAOBJP)simplex.o $(SOPHYAOBJP)slinparbuff.o $(SOPHYAOBJP)tabmath.o $(SOPHYAOBJP)toeplitzMatrix.o $(SOPHYAOBJP)datime.o $(SOPHYAOBJP)fftmayer_r4.o $(SOPHYAOBJP)fftmayer_r8.o $(SOPHYAOBJP)fftpackc.o $(SOPHYAOBJP)matxop.o $(SOPHYAOBJP)nbmath.o $(SOPHYAOBJP)nbtri.o 7 7 rm -f $(SOPHYALIBP)libNTools.a 8 8 9 $(SOPHYALIBP)libNTools.a : $(SOPHYAOBJP)cimage.o $(SOPHYAOBJP)cspline.o $(SOPHYAOBJP)datatypes.o $(SOPHYAOBJP)dates.o $(SOPHYAOBJP)datime.o $(SOPHYAOBJP)difeq.o $(SOPHYAOBJP)dynccd.o $(SOPHYAOBJP)fct1dfit.o $(SOPHYAOBJP)fct2dfit.o $(SOPHYAOBJP)fftmayer_r4.o $(SOPHYAOBJP)fftmayer_r8.o $(SOPHYAOBJP)fftmserver.o $(SOPHYAOBJP)fftpackc.o $(SOPHYAOBJP)fftpserver.o $(SOPHYAOBJP)fftservintf.o $(SOPHYAOBJP)FSAppIrrSmpl.o $(SOPHYAOBJP)functab.o $(SOPHYAOBJP)generaldata.o $(SOPHYAOBJP)generalfit.o $(SOPHYAOBJP)generalfunc.o $(SOPHYAOBJP)integ.o $(SOPHYAOBJP)matxop.o $(SOPHYAOBJP)median.o $(SOPHYAOBJP)nbmath.o $(SOPHYAOBJP)nbtri.o $(SOPHYAOBJP)ntoolsinit.o $(SOPHYAOBJP)objfitter.o $(SOPHYAOBJP)perandom.o $(SOPHYAOBJP)poly.o $(SOPHYAOBJP)rk4cdifeq.o $(SOPHYAOBJP)simplesort.o $(SOPHYAOBJP)simplex.o $(SOPHYAOBJP) tabmath.o $(SOPHYAOBJP)toeplitzMatrix.o9 $(SOPHYALIBP)libNTools.a : $(SOPHYAOBJP)cimage.o $(SOPHYAOBJP)cspline.o $(SOPHYAOBJP)datatypes.o $(SOPHYAOBJP)dates.o $(SOPHYAOBJP)datime.o $(SOPHYAOBJP)difeq.o $(SOPHYAOBJP)dynccd.o $(SOPHYAOBJP)fct1dfit.o $(SOPHYAOBJP)fct2dfit.o $(SOPHYAOBJP)fftmayer_r4.o $(SOPHYAOBJP)fftmayer_r8.o $(SOPHYAOBJP)fftmserver.o $(SOPHYAOBJP)fftpackc.o $(SOPHYAOBJP)fftpserver.o $(SOPHYAOBJP)fftservintf.o $(SOPHYAOBJP)FSAppIrrSmpl.o $(SOPHYAOBJP)functab.o $(SOPHYAOBJP)generaldata.o $(SOPHYAOBJP)generalfit.o $(SOPHYAOBJP)generalfunc.o $(SOPHYAOBJP)integ.o $(SOPHYAOBJP)matxop.o $(SOPHYAOBJP)median.o $(SOPHYAOBJP)nbmath.o $(SOPHYAOBJP)nbtri.o $(SOPHYAOBJP)ntoolsinit.o $(SOPHYAOBJP)objfitter.o $(SOPHYAOBJP)perandom.o $(SOPHYAOBJP)poly.o $(SOPHYAOBJP)rk4cdifeq.o $(SOPHYAOBJP)simplesort.o $(SOPHYAOBJP)simplex.o $(SOPHYAOBJP)slinparbuff.o $(SOPHYAOBJP)tabmath.o $(SOPHYAOBJP)toeplitzMatrix.o 10 10 $(AR) $(ARFLAGS) $@ $? 11 11 touch $(SOPHYAINCP)/SophyaConfInfo/libsophya.objlist … … 39 39 $(SOPHYAINCP)matharr.h \ 40 40 $(SOPHYAINCP)fioarr.h \ 41 $(SOPHYAINCP)sopemtx.h fftservintf.h \ 41 $(SOPHYAINCP)sopemtx.h \ 42 $(SOPHYAINCP)arrctcast.h fftservintf.h \ 42 43 $(SOPHYAINCP)pexceptions.h \ 43 44 $(SOPHYAINCP)tmatrix.h \ … … 640 641 $(CXXCOMPILE) $(CXXTEMPFLG) -o $@ simplex.cc 641 642 643 $(SOPHYAOBJP)slinparbuff.o: slinparbuff.cc \ 644 $(SOPHYAINCP)machdefs.h \ 645 $(SOPHYAINCP)pexceptions.h \ 646 $(SOPHYAINCP)machdefs.h slinparbuff.h 647 $(CXXCOMPILE) $(CXXTEMPFLG) -o $@ slinparbuff.cc 648 642 649 $(SOPHYAOBJP)tabmath.o: tabmath.cc $(SOPHYAINCP)sopnamsp.h \ 643 650 $(SOPHYAINCP)machdefs.h tabmath.h peida.h \ … … 677 684 $(SOPHYAINCP)matharr.h \ 678 685 $(SOPHYAINCP)fioarr.h \ 679 $(SOPHYAINCP)sopemtx.h fftservintf.h \ 686 $(SOPHYAINCP)sopemtx.h \ 687 $(SOPHYAINCP)arrctcast.h fftservintf.h \ 680 688 $(SOPHYAINCP)pexceptions.h \ 681 689 $(SOPHYAINCP)tmatrix.h \ -
trunk/SophyaLib/NTools/datatypes.cc
r2615 r3572 9 9 10 10 /* Nouvelle-Fonction */ 11 c har * DataName(PBaseDataTypes typ)11 const char * DataName(PBaseDataTypes typ) 12 12 { 13 13 switch (typ) … … 33 33 34 34 /* Nouvelle-Fonction */ 35 c har * DataLongName(PBaseDataTypes typ)35 const char * DataLongName(PBaseDataTypes typ) 36 36 { 37 37 switch (typ) -
trunk/SophyaLib/NTools/datatypes.h
r244 r3572 23 23 inline PBaseDataTypes DataType(r_8) {return kr_8;} 24 24 25 c har * DataName(PBaseDataTypes);26 c har * DataLongName(PBaseDataTypes);25 const char * DataName(PBaseDataTypes); 26 const char * DataLongName(PBaseDataTypes); 27 27 int DataSize(PBaseDataTypes); 28 28 -
trunk/SophyaLib/NTools/dates.cc
r2808 r3572 39 39 TimeZone::TimeZone() 40 40 { 41 c har* p = getenv("ACQ_TZ");41 const char* p = getenv("ACQ_TZ"); 42 42 #if defined(__DECCXX) || defined(__KCC__) || defined(__aCC__) 43 43 if (!p) p = const_cast<char *>("France"); -
trunk/SophyaLib/NTools/functab.cc
r2618 r3572 53 53 Dx = (r_8)(XMax-XMin)/(r_8)(NTab-1); 54 54 Tabul = new r_8[NTab+2]; 55 for( int_4 i=0;i<NTab+2;i++) Tabul[i] = 0.;55 for(uint_4 i=0;i<NTab+2;i++) Tabul[i] = 0.; 56 56 PInterpL = PInterpP[0] = PInterpP[1] = NULL; 57 57 if(!Tabul) return; 58 for( int_4 i=0;i<NTab;i++) Tabul[i+1] = func(XMin + i*Dx);58 for(uint_4 i=0;i<NTab;i++) Tabul[i+1] = func(XMin + i*Dx); 59 59 Tabul[0] = Tabul[1]; 60 60 Tabul[NTab+1] = Tabul[NTab]; … … 106 106 if(PInterpL) delete [] PInterpL; 107 107 PInterpL = new r_8[NTab+2]; 108 for( int_4 i=0;i<NTab+1;i++)108 for(uint_4 i=0;i<NTab+1;i++) 109 109 PInterpL[i] = (Tabul[i+1]-Tabul[i])/Dx; 110 110 PInterpL[NTab+1] = 0.; … … 116 116 PInterpP[0] = new r_8[NTab+2]; 117 117 PInterpP[1] = new r_8[NTab+2]; 118 for( int_4 i=1;i<=NTab;i++) {118 for(uint_4 i=1;i<=NTab;i++) { 119 119 PInterpP[0][i] = (Tabul[i+1]+Tabul[i-1]-2.*Tabul[i])/(2.*Dx*Dx); 120 120 PInterpP[1][i] = (Tabul[i+1]-Tabul[i-1])/(2.*Dx); -
trunk/SophyaLib/NTools/generalfit.cc
r3205 r3572 279 279 Pour ecrire les iterations dans le fichier filename 280 280 */ 281 void GeneralFit::WriteStep(c har *filename)281 void GeneralFit::WriteStep(const char *filename) 282 282 { 283 283 -
trunk/SophyaLib/NTools/generalfit.h
r3083 r3572 53 53 ~GeneralFit(); 54 54 55 void WriteStep(c har *filename = NULL);55 void WriteStep(const char *filename = NULL); 56 56 void SetDebug(int level = 0); 57 57 void SetMaxStep(int n = 100); -
trunk/SophyaLib/NTools/ntools.h
r1316 r3572 15 15 #include "poly.h" 16 16 #include "datatypes.h" 17 #include "slinparbuff.h" 17 18 18 19 #include "cimage.h" -
trunk/SophyaLib/NTools/objlist.list
r3083 r3572 31 31 simplesort.o 32 32 simplex.o 33 slinparbuff.o 33 34 tabmath.o 34 35 toeplitzMatrix.o -
trunk/SophyaLib/NTools/pclassids.h
r2602 r3572 30 30 ClassId_HistoErr = 0x0204, 31 31 ClassId_NTuple = 0x0210, 32 ClassId_XNTuple = 0x0211,33 32 34 33 ClassId_GeneralFitData = 0x0290, -
trunk/SophyaLib/NTools/poly.cc
r3235 r3572 710 710 Poly2& Poly2::operator *= (double a) 711 711 { 712 for ( uint_4 i=0; i<NElts(); i++) Element(i) *= a;712 for (int_4 i=0; i<NElts(); i++) Element(i) *= a; 713 713 return *this; 714 714 } -
trunk/SophyaLib/NTools/simplex.cc
r2808 r3572 366 366 int ilo, ihi, inhi; 367 367 int move = 0; 368 c har* smov[6] = { "None", "Reflection", "ReflecExpand", "ContractHigh", "ContractLow", "ExpandHigh" };368 const char* smov[6] = { "None", "Reflection", "ReflecExpand", "ContractHigh", "ContractLow", "ExpandHigh" }; 369 369 int movcnt[6] = {0,0,0,0,0,0}; 370 370 … … 511 511 int MinZSimplex::StopReason(string& s) 512 512 { 513 c har* sr[5] = { "NoConverg, MaxIterReached", "OK, fm<Tol0", "OK, Df/f<Tol1",513 const char* sr[5] = { "NoConverg, MaxIterReached", "OK, fm<Tol0", "OK, Df/f<Tol1", 514 514 "OK, [Df/f max]Iter<Tol2" "Error - Wrong StopReason" }; 515 515 int stop = mStop; -
trunk/SophyaLib/NTools/slinparbuff.h
r2322 r3572 1 #include "machdefs.h"2 #include <iostream>3 #include <stdlib.h>4 #include <stdio.h>5 #include <math.h>6 #include <algorithm>7 1 #include "pexceptions.h" 8 2 … … 17 11 friend class SLinParBuffMerger; 18 12 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 } 13 SLinParBuff(uint_4 lenbuf,uint_4 nresynch=0,r_8 x0=0.,r_8 y0=0.,bool autoxy0=false); 14 SLinParBuff(SLinParBuff& slp); 15 SLinParBuff(void); 16 virtual ~SLinParBuff(); 31 17 32 SLinParBuff(SLinParBuff& slp) 33 // ---- Createur par copie 34 { 35 mLenBuf = slp.mLenBuf; 36 mNResynch = slp.mNResynch; 18 void GetX0Y0(r_8& x0,r_8& y0); 19 bool GetAutoX0Y0(); 20 void SetX0Y0(r_8 x0=0.,r_8 y0=0.); 21 void SetAutoX0Y0(bool autoxy0=false); 37 22 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 23 void Reset(void); 24 uint_4 Pop(void); 25 uint_4 Push(r_8 x,r_8 y); 166 26 inline uint_4 NPoints(void) {return mNCur;} 167 27 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; 28 r_8 SumX(void) {return mSx+mNCur*mX0;} 29 r_8 SumY(void) {return mSy+mNCur*mY0;} 191 30 192 if(sigma>0.) return sqrt(sigma);193 if(sigma<0.) return -sqrt(-sigma);194 return sigma;195 }31 uint_4 ReComputeSum(void); 32 r_8 Compute(r_8& mean,bool recomputeXi2=false); 33 r_8 Compute(r_8& a0,r_8 &a1,bool recomputeXi2=false); 34 r_8 Compute(r_8& a0,r_8 &a1,r_8 &a2,bool recomputeXi2=false); 196 35 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 } 36 void Print(int lp=0); 37 void PrintCompute(int lp=0); 322 38 323 39 protected: … … 329 45 uint_4 mNResynchEff, mNPush, mNPop; 330 46 }; 47 331 48 332 49 /////////////////////////////////////////////////////////////////// … … 343 60 class SLinParBuffMerger { 344 61 public: 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 62 SLinParBuffMerger(void); 63 SLinParBuffMerger(SLinParBuff& s,bool recompute=false); 64 virtual ~SLinParBuffMerger(void); 65 66 inline uint_4 NPoints(void) {return mSlp.NPoints();} 350 67 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);} 68 void Add(SLinParBuff& s,bool recompute=false); 69 inline r_8 Compute(r_8& mean) {return mSlp.Compute(mean);} 70 inline r_8 Compute(r_8& a0,r_8 &a1) {return mSlp.Compute(a0,a1);} 71 inline r_8 Compute(r_8& a0,r_8 &a1,r_8 &a2) {return mSlp.Compute(a0,a1,a2);} 387 72 inline void Print(int lp=0) {mSlp.Print(lp);} 388 73 inline void PrintCompute(int lp=0) {mSlp.PrintCompute(lp);} -
trunk/SophyaLib/NTools/smakefile
r3503 r3572 4 4 5 5 clean: 6 rm -f $(SOPHYAOBJP)FSAppIrrSmpl.o $(SOPHYAOBJP)cimage.o $(SOPHYAOBJP)cspline.o $(SOPHYAOBJP)datatypes.o $(SOPHYAOBJP)dates.o $(SOPHYAOBJP)difeq.o $(SOPHYAOBJP)dynccd.o $(SOPHYAOBJP)fct1dfit.o $(SOPHYAOBJP)fct2dfit.o $(SOPHYAOBJP)fftmserver.o $(SOPHYAOBJP)fftpserver.o $(SOPHYAOBJP)fftservintf.o $(SOPHYAOBJP)functab.o $(SOPHYAOBJP)generaldata.o $(SOPHYAOBJP)generalfit.o $(SOPHYAOBJP)generalfunc.o $(SOPHYAOBJP)imageop.o $(SOPHYAOBJP)integ.o $(SOPHYAOBJP)median.o $(SOPHYAOBJP)ntoolsinit.o $(SOPHYAOBJP)objfitter.o $(SOPHYAOBJP)perandom.o $(SOPHYAOBJP)poly.o $(SOPHYAOBJP)rk4cdifeq.o $(SOPHYAOBJP)simplesort.o $(SOPHYAOBJP)simplex.o $(SOPHYAOBJP) tabmath.o $(SOPHYAOBJP)toeplitzMatrix.o $(SOPHYAOBJP)datime.o $(SOPHYAOBJP)fftmayer_r4.o $(SOPHYAOBJP)fftmayer_r8.o $(SOPHYAOBJP)fftpackc.o $(SOPHYAOBJP)matxop.o $(SOPHYAOBJP)nbmath.o $(SOPHYAOBJP)nbtri.o6 rm -f $(SOPHYAOBJP)FSAppIrrSmpl.o $(SOPHYAOBJP)cimage.o $(SOPHYAOBJP)cspline.o $(SOPHYAOBJP)datatypes.o $(SOPHYAOBJP)dates.o $(SOPHYAOBJP)difeq.o $(SOPHYAOBJP)dynccd.o $(SOPHYAOBJP)fct1dfit.o $(SOPHYAOBJP)fct2dfit.o $(SOPHYAOBJP)fftmserver.o $(SOPHYAOBJP)fftpserver.o $(SOPHYAOBJP)fftservintf.o $(SOPHYAOBJP)functab.o $(SOPHYAOBJP)generaldata.o $(SOPHYAOBJP)generalfit.o $(SOPHYAOBJP)generalfunc.o $(SOPHYAOBJP)imageop.o $(SOPHYAOBJP)integ.o $(SOPHYAOBJP)median.o $(SOPHYAOBJP)ntoolsinit.o $(SOPHYAOBJP)objfitter.o $(SOPHYAOBJP)perandom.o $(SOPHYAOBJP)poly.o $(SOPHYAOBJP)rk4cdifeq.o $(SOPHYAOBJP)simplesort.o $(SOPHYAOBJP)simplex.o $(SOPHYAOBJP)slinparbuff.o $(SOPHYAOBJP)tabmath.o $(SOPHYAOBJP)toeplitzMatrix.o $(SOPHYAOBJP)datime.o $(SOPHYAOBJP)fftmayer_r4.o $(SOPHYAOBJP)fftmayer_r8.o $(SOPHYAOBJP)fftpackc.o $(SOPHYAOBJP)matxop.o $(SOPHYAOBJP)nbmath.o $(SOPHYAOBJP)nbtri.o 7 7 rm -f $(SOPHYALIBP)libNTools.a 8 8 9 $(SOPHYALIBP)libNTools.a : $(SOPHYAOBJP)cimage.o $(SOPHYAOBJP)cspline.o $(SOPHYAOBJP)datatypes.o $(SOPHYAOBJP)dates.o $(SOPHYAOBJP)datime.o $(SOPHYAOBJP)difeq.o $(SOPHYAOBJP)dynccd.o $(SOPHYAOBJP)fct1dfit.o $(SOPHYAOBJP)fct2dfit.o $(SOPHYAOBJP)fftmayer_r4.o $(SOPHYAOBJP)fftmayer_r8.o $(SOPHYAOBJP)fftmserver.o $(SOPHYAOBJP)fftpackc.o $(SOPHYAOBJP)fftpserver.o $(SOPHYAOBJP)fftservintf.o $(SOPHYAOBJP)FSAppIrrSmpl.o $(SOPHYAOBJP)functab.o $(SOPHYAOBJP)generaldata.o $(SOPHYAOBJP)generalfit.o $(SOPHYAOBJP)generalfunc.o $(SOPHYAOBJP)integ.o $(SOPHYAOBJP)matxop.o $(SOPHYAOBJP)median.o $(SOPHYAOBJP)nbmath.o $(SOPHYAOBJP)nbtri.o $(SOPHYAOBJP)ntoolsinit.o $(SOPHYAOBJP)objfitter.o $(SOPHYAOBJP)perandom.o $(SOPHYAOBJP)poly.o $(SOPHYAOBJP)rk4cdifeq.o $(SOPHYAOBJP)simplesort.o $(SOPHYAOBJP)simplex.o $(SOPHYAOBJP) tabmath.o $(SOPHYAOBJP)toeplitzMatrix.o9 $(SOPHYALIBP)libNTools.a : $(SOPHYAOBJP)cimage.o $(SOPHYAOBJP)cspline.o $(SOPHYAOBJP)datatypes.o $(SOPHYAOBJP)dates.o $(SOPHYAOBJP)datime.o $(SOPHYAOBJP)difeq.o $(SOPHYAOBJP)dynccd.o $(SOPHYAOBJP)fct1dfit.o $(SOPHYAOBJP)fct2dfit.o $(SOPHYAOBJP)fftmayer_r4.o $(SOPHYAOBJP)fftmayer_r8.o $(SOPHYAOBJP)fftmserver.o $(SOPHYAOBJP)fftpackc.o $(SOPHYAOBJP)fftpserver.o $(SOPHYAOBJP)fftservintf.o $(SOPHYAOBJP)FSAppIrrSmpl.o $(SOPHYAOBJP)functab.o $(SOPHYAOBJP)generaldata.o $(SOPHYAOBJP)generalfit.o $(SOPHYAOBJP)generalfunc.o $(SOPHYAOBJP)integ.o $(SOPHYAOBJP)matxop.o $(SOPHYAOBJP)median.o $(SOPHYAOBJP)nbmath.o $(SOPHYAOBJP)nbtri.o $(SOPHYAOBJP)ntoolsinit.o $(SOPHYAOBJP)objfitter.o $(SOPHYAOBJP)perandom.o $(SOPHYAOBJP)poly.o $(SOPHYAOBJP)rk4cdifeq.o $(SOPHYAOBJP)simplesort.o $(SOPHYAOBJP)simplex.o $(SOPHYAOBJP)slinparbuff.o $(SOPHYAOBJP)tabmath.o $(SOPHYAOBJP)toeplitzMatrix.o 10 10 $(AR) $(ARFLAGS) $@ $? 11 11 touch $(SOPHYAINCP)/SophyaConfInfo/libsophya.objlist … … 641 641 $(CXXCOMPILE) $(CXXTEMPFLG) -o $@ simplex.cc 642 642 643 $(SOPHYAOBJP)slinparbuff.o: slinparbuff.cc \ 644 $(SOPHYAINCP)machdefs.h \ 645 $(SOPHYAINCP)pexceptions.h \ 646 $(SOPHYAINCP)machdefs.h slinparbuff.h 647 $(CXXCOMPILE) $(CXXTEMPFLG) -o $@ slinparbuff.cc 648 643 649 $(SOPHYAOBJP)tabmath.o: tabmath.cc $(SOPHYAINCP)sopnamsp.h \ 644 650 $(SOPHYAINCP)machdefs.h tabmath.h peida.h \
Note:
See TracChangeset
for help on using the changeset viewer.