source: Sophya/trunk/Poubelle/DPC:FitsIOServer/NTools/generaldata.cc@ 658

#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <values.h>
#include <math.h>
#include <string.h>
#include <string>

#include "strutil.h"
#include "nbtri.h"
#include "generalfit.h"
#include "generaldata.h"
#include "pexceptions.h"
#include "objfio.h"

//================================================================
// GeneralFitData
//================================================================

//++
// Class        GeneralFitData
// Lib          Outils++ 
// include      generaldata.h
//
//      Classe de stoquage de donnees a plusieurs variables avec erreur
//      sur l'ordonnee et sur les abscisses (options).
//|   {x0(i),Ex0(i), x1(i),Ex1(i), x2(i),Ex2(i) ... ; Y(i),EY(i)}
//--

// Pour memoire, structure du rangement (n=mNVar):
// - Valeur des abscisses mXP (idem pour mErrXP):
//   x0,x1,x2,...,xn   x0,x1,x2,...,xn  ....  x0,x1,x2,....,xn
//   |  1er point  |   |  2sd point  |  ....  | point mNData |
//   Donc abscisse J=[0,mNVar[ du point numero I=[0,mNData[: mXP[I*mNVar+J]
// - Valeur de l'ordonnee mF (idem pour mErr et mOK):
//          f                f                      f
//   |  1er point  |  |  2sd point  |  .... | point mNData |
//   Donc point numero I [0,mNData[ : mF[i]

//////////////////////////////////////////////////////////////////////
//++
GeneralFitData::GeneralFitData(unsigned int nVar, unsigned int ndatalloc, uint_2 errx)
//
//      Constructeur. ``nVar'' represente la dimension de l'espace des abscisses,
//      ``ndatalloc'' le nombre maximum de points et ``errx'' si non nul
//      indique que l'on fournit des erreurs sur les ``nVar'' variables en abscisse.
//--
  : mNVar(0), mNDataAlloc(0), mNData(0), mNDataGood(0), mOk_EXP(0)
  , mXP(NULL), mErrXP(NULL), mF(NULL), mErr(NULL), mOK(NULL)
  , BuffVar(NULL), BuffVarR4(NULL)
{
try {
  Alloc(nVar,ndatalloc,errx);
} catch(PException e) {
  cout << "Exception : " << typeid(e).name() << " " << e.Msg() << endl;
  throw;
}
}

//++
GeneralFitData::GeneralFitData(const GeneralFitData& data, bool clean)
//
//      Constructeur par copie. Si ``clean'' est ``true''
//      seules les donnees valides de ``data'' sont copiees.
//      Si ``clean'' est ``false'' (defaut) toutes les donnees
//      sont copiees et la taille totale de ``data'' est allouee
//      meme si elle est plus grande que la taille des donnees stoquees.
//--
  : mNVar(0), mNDataAlloc(0), mNData(0), mNDataGood(0), mOk_EXP(0)
  , mXP(NULL), mErrXP(NULL), mF(NULL), mErr(NULL), mOK(NULL)
  , BuffVar(NULL), BuffVarR4(NULL)
{
try {
  Alloc(data.mNVar,((clean)?data.mNDataGood:data.mNDataAlloc),((data.mErrXP)?1:0));
} catch(PException e) {
  cout << "Exception : " << typeid(e).name() << " " << e.Msg() << endl;
  throw;
}

 // Remplissage
 if(data.mNData>0) {
   r_8* ret;
   for(int i=0;i<data.mNData;i++) {
     if( clean && data.mOK[i]==0 ) continue;
       ret = data.GetVec(i,NULL);
       memcpy(mXP+mNData*mNVar,ret,mNVar*sizeof(r_8));
       if(mErrXP) memcpy(mErrXP+mNData*mNVar,ret+mNVar,mNVar*sizeof(r_8));
       mF[mNData]   = ret[2*mNVar];
       mErr[mNData] = ret[2*mNVar+1];
       mOK[mNData]  = (uint_2) (ret[2*mNVar+2]+0.001);
       if(mOK[mNData]!=0) mNDataGood++;
       mNData++;
   }
 }

 END_CONSTRUCTOR
}

//++
GeneralFitData::GeneralFitData()
//
//      Constructeur par defaut.
//--
  : mNVar(0), mNDataAlloc(0), mNData(0), mNDataGood(0), mOk_EXP(0)
  , mXP(NULL), mErrXP(NULL), mF(NULL), mErr(NULL), mOK(NULL)
  , BuffVar(NULL), BuffVarR4(NULL)
{
 END_CONSTRUCTOR
}

//++
GeneralFitData::~GeneralFitData()
//
//      Destructeur
//--
{
 Delete();
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::Alloc(unsigned int nVar, unsigned int ndatalloc, int_2 errx)
//
//      Pour redefinir la structure de donnees (ou la creer si on a utilise
//      le createur par defaut). Voir les explications des arguments
//      dans les commentaires du constructeur. Si ``errx''<0 alors
//      la valeur prise est celle definie auparavent.
//--
{
ASSERT( nVar>0 && ndatalloc>0 );

Delete();

if(errx>=0) mOk_EXP = (uint_2) errx;
mNVar = nVar;
mNDataAlloc = ndatalloc;

try {
  mXP  = new r_8[nVar*ndatalloc];
  if(mOk_EXP) mErrXP  = new r_8[nVar*ndatalloc];
  mF   = new r_8[ndatalloc];
  mErr = new r_8[ndatalloc];
  mOK  = new uint_2[ndatalloc];
  BuffVar   = new r_8[2*nVar+3];
  BuffVarR4 = (r_4 *) BuffVar;
} catch(PException e) {
  throw(AllocationError("GeneralFitData::Alloc allocation error\n"));
}
}

//////////////////////////////////////////////////////////////////////
void GeneralFitData::Delete()
{
mNVar = mNDataAlloc = mNData = mNDataGood = 0;
if( mXP  != NULL ) {delete [] mXP;  mXP  = NULL;}
if( mErrXP  != NULL ) {delete [] mErrXP;  mErrXP  = NULL;}
if( mF   != NULL ) {delete [] mF;   mF   = NULL;}
if( mErr != NULL ) {delete [] mErr; mErr = NULL;}
if( mOK  != NULL ) {delete [] mOK;  mOK  = NULL;}
if( BuffVar  != NULL ) {delete [] BuffVar; BuffVar = NULL; BuffVarR4 = NULL;}
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::SetDataPtr(int ptr)
//
//      Remise a zero de la structure pour nouveau remplissage (pas d'arg)
//      ou remise a la position ``ptr'' (si arg). Les donnees apres ``ptr''
//      sont sur-ecrites.
//--
{
 ASSERT(ptr >= 0 && ptr < mNDataAlloc);
 mNData = ptr;
 mNDataGood = 0;
 if(ptr==0) return;
 for(int i=0;i<mNData;i++) if(mOK[i]) mNDataGood++;
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::KillData(int i)
//
//      Pour tuer un point
//--
{
 ASSERT(i >= 0 && i < mNData);

 if( ! mOK[i] ) return;
 mOK[i] = 0;
 mNDataGood--;
}

//++
void GeneralFitData::KillData(int i1,int i2)
//
//      Pour tuer une serie de points
//--
{
 ASSERT(i1 >= 0 && i1 < mNData);
 ASSERT(i2 >= 0 && i2 < mNData);
 ASSERT(i1 <= i2 );

 for(int i=i1;i<=i2;i++) KillData(i);
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::ValidData(int i)
//
//      Pour re-valider le point numero i ([0,NData]).
//--
{
 ASSERT(i >= 0 && i < mNData);

 if( mOK[i] ) return;
 if( mErr[i]<=0. ) return;
 if(mOk_EXP) {
   for(int j=0;j<mNVar;j++) if(mErrXP[i*mNVar+j]<=0.) return;
 }
 mOK[i] = 1;
 mNDataGood++;
}

//++
void GeneralFitData::ValidData(int i1,int i2)
//
//      Pour re-valider les points numeros i1 a i2.
//--
{
 ASSERT(i1 >= 0 && i1 < mNData);
 ASSERT(i2 >= 0 && i2 < mNData);
 ASSERT(i1 <= i2 );

 for(int i=i1;i<=i2;i++) ValidData(i);
}

//++
void GeneralFitData::ValidData()
//
//      Pour re-valider tous les points.
//--
{
 for(int i=0;i<mNData;i++) ValidData(i);
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::RedefineData1(int i,double x,double f,double err,double errx)
//
//      Pour redefinir un point a
//|     {x,[errx] ; f,err}
//--
{
 RedefineData(i,&x,f,err,&errx);
}

//++
void GeneralFitData::RedefineData2(int i,double x,double y,double f
                                  ,double err,double errx,double erry)
//
//      Pour redefinir un point a
//|     {x,[errx], y,[erry] ; f,err}
//--
{
 double xp[2] = {x,y};
 double errxp[2] = {errx,erry};
 RedefineData(i,xp,f,err,errxp);
}

//++
void GeneralFitData::RedefineData(int i,double* xp,double f,double err,double* errxp)
//
//      Pour redefinir un point a
//| {xp[0],[errxp[0]], xp[1],[errxp[1]], xp[2],[errxp[2]],... ; f,err}
//--
{
 ASSERT(i>=0 && i<mNData);
 bool ok = true;

 int ip = mNVar*i;
 for(int j=0;j<mNVar;j++) mXP[ip+j] = xp[j];
 if(mOk_EXP) {
   if(errxp) {
     for(int j=0;j<mNVar;j++)
       {mErrXP[ip+j] = errxp[j]; if(errxp[j]<=0.) ok=false;}
   } else {
     for(int j=0;j<mNVar;j++) mErrXP[ip+j] = Def_ErrX;
     ok=false;
   }
 }
 mF[i] = f;
 mErr[i] = err;  if(err<=0.) ok = false;
 if(ok) {
   if(! mOK[i]) {mOK[i]=1; mNDataGood++;}
 } else {
   if(  mOK[i]) {mOK[i]=0; mNDataGood--;}
 }
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::AddData1(double x, double f, double err, double errx)
//
//      Pour ajouter un point
//|     {x,[errx] ; f,err}
//--
{
 AddData(&x,f,err,&errx);
}

//++
void GeneralFitData::AddData2(double x, double y, double f
                             , double err, double errx, double erry)
//
//      Pour ajouter un point
//|     {x,[errx], y,[erry] ; f,err}
//--
{
 double xp[2] = {x,y};
 double errxp[2] = {errx,erry};
 AddData(xp,f,err,errxp);
}

//++
void GeneralFitData::AddData(double* xp, double f, double err,double* errxp)
//
//      Pour ajouter un point
//| {xp[0],[errxp[0]], xp[1],[errxp[1]], xp[2],[errxp[2]],... ; f,err}
//--
{
 ASSERT(mNData < mNDataAlloc);
 bool ok = true;

 int ip = mNVar*mNData;
 for(int i=0;i<mNVar;i++) mXP[ip+i] = xp[i];
 if(mOk_EXP) {
   if(errxp) {
     for(int j=0;j<mNVar;j++)
       {mErrXP[ip+j] = errxp[j]; if(errxp[j]<=0.) ok=false;}
   } else {
     for(int j=0;j<mNVar;j++) mErrXP[ip+j] = Def_ErrX;
     ok=false;
   }
 }
 mF[mNData] = f;
 mErr[mNData] = err;
 if(err<=0.) ok = false;
 if(ok) { mOK[mNData]=1; mNDataGood++; } else mOK[mNData]=0;
 mNData++;
}

//++
void GeneralFitData::AddData(float* xp, float f, float err, float* errxp)
//
//      Pour ajouter un point
//| {xp[0],[errxp[0]], xp[1],[errxp[1]], xp[2],[errxp[2]],... ; f,err}
//--
{
 {for(int i=0;i<mNVar;i++) BuffVar[i] = (double) xp[i];}
 if(errxp) for(int i=0;i<mNVar;i++) BuffVar[mNVar+i] = (double) errxp[i];
 AddData(BuffVar,(double) f,(double)  err,BuffVar+mNVar);
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::SetData1(int nData
            , double* x, double* f, double *err, double *errx)
//
//      Pour remplir la structure de donnees d'un seul coup avec
//|     {x(i),[errx(i)] ; f(i),err(i)}
//--
{
 ASSERT(nData>0 && mNData+nData<=mNDataAlloc);

 for(int i=0;i<nData;i++) {
   double ex = (errx) ? errx[i]: Def_ErrX;
   double ef = (err ) ? err[i]:  Def_ErrF;
   AddData1(x[i],f[i],ef,ex);
 }
}

//++
void GeneralFitData::SetData1(int nData
            , float* x, float* f, float* err, float *errx)
//
//      Pour remplir la structure de donnees d'un seul coup avec
//|     {x(i),[errx(i)] ; f(i),err(i)}
//--
{
 ASSERT(nData>0 && mNData+nData<=mNDataAlloc);

 for(int i=0;i<nData;i++) {
   double ex = (errx) ? (double) errx[i]: Def_ErrX;
   double ef = (err ) ? (double) err[i]:  Def_ErrF;
   AddData1((double) x[i],(double) f[i],ef,ex);
 }
}

//++
void GeneralFitData::SetData2(int nData, double* x, double* y, double* f
                             ,double *err,double *errx,double *erry)
//
//      Pour remplir la structure de donnees d'un seul coup avec
//|     {x(i),[errx(i)], y(i),[erry(i)], ; f(i),err(i)}
//--
{
 ASSERT(nData>0 && mNData+nData<=mNDataAlloc);

 for(int i=0;i<nData;i++) {
   double ex = (errx) ? (double) errx[i]: Def_ErrX;
   double ey = (erry) ? (double) erry[i]: Def_ErrX;
   double ef = (err ) ? (double) err[i]:  Def_ErrF;
   AddData2(x[i],y[i],f[i],ef,ex,ey);
 }
}

//++
void GeneralFitData::SetData2(int nData, float* x, float* y, float* f
                             ,float *err,float *errx,float *erry)
//
//      Pour remplir la structure de donnees d'un seul coup avec
//|     {x(i),[errx(i)], y(i),[erry(i)], ; f(i),err(i)}
//--
{
 ASSERT(nData>0 && mNData+nData<=mNDataAlloc);

 for(int i=0;i<nData;i++) {
   double ex = (errx) ? (double) errx[i]: Def_ErrX;
   double ey = (erry) ? (double) erry[i]: Def_ErrX;
   double ef = (err ) ? (double) err[i]:  Def_ErrF;
   AddData2((double) x[i],(double) y[i],(double) f[i],ef,ex,ey);
 }
}

//++
void GeneralFitData::SetData(int nData,double** xp, double *f
                            , double *err, double** errxp)
//
//      Pour remplir la structure de donnees d'un seul coup avec
//|   {X0(i),[EX0(i)], X1(i),[EX1(i)], X2(i),[EX2(i)], ... ; F(i),Err(i)}
//      Attention: si la structure est n'est pas vide, les tableaux sont copies
//      apres les donnees pre-existantes (qui ne sont donc pas detruites). Pour
//      effacer les donnees pre-existantes utiliser SetDataPtr(0).
//|   Ici **xp est un pointeur sur un tableau de pointeurs tq
//|     xp[0] = &X0[0], xp[1] = &X1[0], xp[2] = &X2[0] ...
//|   ou X0,X1,X2,... sont les tableaux X0[nData] X1[nData] X2[nData] ...
//|   des variables (meme commentaire pour errxp).
//--
{
 ASSERT(nData>0 && mNData+nData<=mNDataAlloc);
 if(mOk_EXP && !errxp) {for(int j=0;j<mNVar;j++) BuffVar[mNVar+j] = Def_ErrX;}

 for(int i=0;i<nData;i++) {
   {for(int j=0;j<mNVar;j++) BuffVar[j] = *(xp[j]+i);}
   if(mOk_EXP && errxp)
     {for(int j=0;j<mNVar;j++) BuffVar[mNVar+j] = *(errxp[j]+i);}
   double ef = (err) ? err[i]:  Def_ErrF;
   AddData(BuffVar,f[i],ef,BuffVar+mNVar);
 }
}

//++
void GeneralFitData::SetData(int nData,float** xp, float *f
                            , float *err, float** errxp)
//
//      Voir commentaire ci-dessus.
//--
{
 ASSERT(nData>0 && mNData+nData<=mNDataAlloc);

 if(mOk_EXP && !errxp) {for(int j=0;j<mNVar;j++) BuffVar[mNVar+j] = Def_ErrX;}

 for(int i=0;i<nData;i++) {
   {for(int j=0;j<mNVar;j++) BuffVar[j] = (double) *(xp[j]+i);}
   if(mOk_EXP && errxp)
     {for(int j=0;j<mNVar;j++) BuffVar[mNVar+j] = (double) *(errxp[j]+i);}
   double ef = (err) ? err[i]:  Def_ErrF;
   AddData(BuffVar,(double) f[i],ef,BuffVar+mNVar);
 }
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::PrintStatus()
//
//      Impression de l'etat de la structure de donnees
//--
{
  cout<<"GeneralFitData:: "<<endl
      <<"NVar="<<mNVar<<" NDataAlloc="<<mNDataAlloc<<" Ok_EXP="<<mOk_EXP
      <<" ,NData="<<mNData<<" NDataGood="<<mNDataGood<<endl
      <<"  mXP="<<mXP<<"  [mErrXP="<<mErrXP<<"] mF="<<mF<<" mErr="<<mErr
      <<" mOK="<<mOK<<endl;
}

//++
void GeneralFitData::PrintData(int i)
//
//      Impression du point i
//--
{
 ASSERT(i>=0 && i<mNData);

 cout<<" "<<i<<" F( ";
 {for(int j=0;j<mNVar;j++) cout<<" "<<Absc(j,i);}
 if(mOk_EXP) {
   cout<<"  ;  ";
   for(int j=0;j<mNVar;j++) cout<<" "<<EAbsc(j,i);
 }
 cout<<")= "<<Val(i)<<" "<<EVal(i)<<" ("<<IsValid(i)<<")\n";
}

//++
void GeneralFitData::PrintData(int i1,int i2)
//
//      Impression des points i1 a i2
//--
{
 if(i1<0) i1=0;
 if(i1>=mNData) i1 = mNData-1;
 if(i2>=mNData) i2 = mNData-1;
 if(i1>i2) i2 = mNData-1;

 cout<<"GeneralFitData::PrintData[NData="
     <<mNData<<"/ NDataGood="<<mNDataGood<<"]"<<endl;
 for(int i=i1;i<=i2;i++) PrintData(i);
 cout<<flush;
}

//++
void GeneralFitData::PrintData()
//
//      Impression de tous les points
//--
{
 ASSERT(mNData>0);

 PrintData(0,mNData-1);
}

//++
void GeneralFitData::Show(ostream& os) const
//
//      Impression de l'etat de la structure de donnees avec bornes sur "s"
//--
{
double min,max;
os<<"GeneralFitData:: NVar,ErrX="<<mNVar<<","<<mOk_EXP
  <<" Data: "<<mNData<<" Good,Alloc="<<mNDataGood<<","<<mNDataAlloc<<endl;
for(int k=0;k<2*NVar()+3;k++) {
  GetMinMax(k,min,max);
  os<<" - "<<k<<" "<<ColumnName(k)<<"  ,  "<<min<<","<<max<<endl;
}
return;
}

//////////////////////////////////////////////////////////////////////
//++
int GeneralFitData::GetMnMx(int var,int& imin,int& imax) const
//
//      Retourne les numeros des points de valeurs minimum et maximum
//      de la variable ``var'':
//| La variable "var" est de la forme : var = AB avec
//|  B = 0 : variable d'ordonnee Y (valeur de A indifferente)
//|  B = 1 : erreur variable d'ordonnee EY (valeur de A indifferente)
//|  B = 2 : variable d'abscisse X numero A #[0,NVar[
//|  B = 3 : erreur variable d'abscisse EX numero A #[0,NVar[
//| - Return NData checked si ok, -1 si probleme.
//--
{
imin = imax = -1;
int ix = var/10;
var = var%10;
if(var<0 || var>3) return -1;
if(var>=2 && (ix<0 || ix>=mNVar) ) return -1;
double min, max;
int ntest = 0;
for(int i=0;i<mNData;i++) {
  if( ! IsValid(i) ) continue;
  double v;
  if(var==0)      v = Val(i);
  else if(var==1) v = EVal(i);
  else if(var==2) v = Absc(ix,i);
  else if(var==3) v = EAbsc(ix,i);
  if(ntest==0) {min = max = v; imin = imax = i;}
  if(v<min) {min = v; imin = i;}
  if(v>max) {max = v; imax = i;}
  ntest++;
}
return ntest;
}

//++
int GeneralFitData::GetMnMx(int var,double& min,double& max) const
//
//      Retourne le minimum et le maximum de la variable ``var''
//      (cf commentaires GetMnMx).
//--
{
min = 1.; max = -1.;
int imin,imax;
int ntest = GetMnMx(var,imin,imax);
if(ntest<=0) return ntest;
int ix = var/10;
var = var%10;
if(var==0) {
  if(imin>=0) min = Val(imin);
  if(imax>=0) max = Val(imax);
} else if(var==1) {
  if(imin>=0) min = EVal(imin);
  if(imax>=0) max = EVal(imax);
} else if(var==2) {
  if(imin>=0) min = Absc(ix,imin);
  if(imax>=0) max = Absc(ix,imax);
} else if(var==3) {
  if(imin>=0) min = EAbsc(ix,imin);
  if(imax>=0) max = EAbsc(ix,imax);
}
return ntest;
}

//////////////////////////////////////////////////////////////////////
//++
int GeneralFitData::GetMeanSigma(int var,double& mean,double& sigma,double min,double max)
//
//      Retourne la moyenne et le sigma de la variable ``var''
//      (cf commentaires GetMnMx).
//| - Return : nombre de donnees utilisees, -1 si pb, -2 si sigma<0.
//| - Seuls les points valides de valeur entre min,max sont utilises.
//|   Si min>=max pas de coupures sur les valeurs.
//--
{
mean = sigma = 0.;
int ix = var/10;
var = var%10;
if(var<0 || var>3) return -1;
if(var>=2 && (ix<0 || ix>=mNVar) ) return -1;
int ntest = 0;
for(int i=0;i<mNData;i++) {
  if( ! IsValid(i) ) continue;
  double v;
  if(var==0)      v = Val(i);
  else if(var==1) v = EVal(i);
  else if(var==2) v = Absc(ix,i);
  else if(var==3) v = EAbsc(ix,i);
  if(min<max && (v<min || max<v)) continue;
  mean += v;
  sigma += v*v;
  ntest++;
}
if(ntest==0) {
  mean = sigma = 0.;
} else {
  mean /= (double)ntest;
  sigma = sigma/(double)ntest - mean*mean;
  if(sigma<0.) ntest = -2;
  else if(sigma>0.) sigma = sqrt(sigma);
}
return ntest;
}

//++
int GeneralFitData::GetMoMeMed(int var,double& mode,double& mean,double& median,
                               double min,double max,double coeff)
//
//      Retourne le mode de la variable ``var''
//      (cf commentaires GetMnMx).
//| - Return : nombre de donnees utilisees, -1 si pb.
//| - Seuls les points valides de valeur entre min,max sont utilises.
//|   Si min>=max pas de coupures sur les valeurs.
//| - Le calcul du mode est approximee par la formule:
//|     Mode = Median - coeff*(Mean-Median)   (def: coeff=0.8)
//| - Kendall and Stuart donne coeff=2., mais coeff peut etre regle.
//--
{
mode = mean = median = 0.;
if(mNData<=0) return -1;
int ix = var/10;
var = var%10;
if(var<0 || var>3) return -1;
if(var>=2 && (ix<0 || ix>=mNVar) ) return -1;
double* buff = new double[mNData];
int ntest = 0;
for(int i=0;i<mNData;i++) {
  if( ! IsValid(i) ) continue;
  double v;
  if(var==0)      v = Val(i);
  else if(var==1) v = EVal(i);
  else if(var==2) v = Absc(ix,i);
  else if(var==3) v = EAbsc(ix,i);
  if(min<max && (v<min || max<v)) continue;
  buff[ntest] = v;
  mean += v;
  ntest++;
}
if(ntest==0) {
  mean = 0.;
} else {
  mean /= (double)ntest;
  qsort(buff,(size_t) ntest,sizeof(double),qSort_Dble);
  int im;
  if(ntest%2==1) {
    // nombre impair de points
    im = ntest/2;
    median = buff[im];
  } else {
    // nombre pair de points
    im = (ntest-1)/2;
    median = (buff[im]+buff[im+1])/2.;
  }
  mode = median - coeff*(mean-median);
}
delete [] buff;
return ntest;
}

//++
int GeneralFitData::GetMode(int var,double& mode,double min,double max,double coeff)
//
//      Cf description ci-dessus ``GetMoMeMed''.
//--
{
double mean,median;
return GetMoMeMed(var,mode,mean,median,min,max,coeff);
}

//////////////////////////////////////////////////////////////////////
//++
double GeneralFitData::PolFit(int varx,Poly& pol,int degre,bool ey)
//
//      Pour fiter un polynome de degre ``degre''. On fite
//      Y=f(X) ou Y=Val et X=Absc(varx). Si ``ey'' est ``true''
//      le fit prend en compte les erreurs stoquees dans EVal,
//      sinon fit sans erreurs. Le resultat du fit est retourne
//      dans le polynome ``pol''.
//| Return:
//|  -   Res = le residu du fit
//|  -   -1 si degre<0
//|  -   -2 si probleme sur numero de variable X
//|  -   -4 si NDataGood<0
//|  -   -5 si nombre de data trouves different de NDataGood
//--
{
if(degre<0) return -1.;
if(varx<0 || varx>=mNVar) return -2.;
if(mNDataGood<=0) return -4.;
Vector x(mNDataGood);
Vector y(mNDataGood);
Vector ey2(1);
if(ey) ey2.Realloc(mNDataGood,true);
int ntest = 0;
for(int i=0;i<mNData;i++) {
  if( ! IsValid(i) ) continue;
  if(ntest>=mNDataGood) return -5.;
  x(ntest) = Absc(varx,i);
  y(ntest) = Val(i);
  if(ey) ey2(ntest) = EVal(i)*EVal(i);
  ntest++;
}
double res = 0.;
if(ey) {
  Vector errcoef(1);
  res = pol.Fit(x,y,ey2,degre,errcoef);
} else {
  res = pol.Fit(x,y,degre);
}
return res;
}

//++
double GeneralFitData::PolFit(int varx,int vary,Poly2& pol,int degre1,int degre2,bool ez)
//
//
//      Pour fiter un polynome de degre ``degre1''. On fite
//      Z=f(X,Y) ou Z=Val et X=Absc(varx) et Y=Absc(vary).
//      Si ``ey'' est ``true'' le fit prend en compte les erreurs
//      stoquees dans EVal, sinon fit sans erreurs. Si ``degre2''
//      negatif, le fit determine un polynome en X,Y de degre
//      total ``degre`''. Si ``degre2'' positif ou nul, le fit
//      demande un fit de ``degre1'' pour la variable X et de degre
//      ``degre2'' sur la variable Y. Le resultat du fit est retourne
//      dans le polynome ``pol''.
//| Return:
//|  -   Res = le residu du fit
//|  -   -1 si degre<0
//|  -   -2 si probleme sur numero de variable X
//|  -   -3 si probleme sur numero de variable Y
//|  -   -4 si NDataGood<0
//|  -   -5 si nombre de data trouves different de NDataGood
//--
{
if(degre1<0) return -1.;
if(varx<0 || varx>=mNVar) return -2.;
if(vary<0 || vary>=mNVar || vary==varx) return -3.;
if(mNDataGood<=0) return -4.;
Vector x(mNDataGood);
Vector y(mNDataGood);
Vector z(mNDataGood);
Vector ez2(1);
if(ez) ez2.Realloc(mNDataGood,true);
int ntest = 0;
for(int i=0;i<mNData;i++) {
  if( ! IsValid(i) ) continue;
  if(ntest>=mNDataGood) return -5.;
  x(ntest) = Absc(varx,i);
  y(ntest) = Absc(vary,i);
  z(ntest) = Val(i);
  if(ez) ez2(ntest) = EVal(i)*EVal(i);
  ntest++;
}
double res = 0.;
if(ez) {
  Vector errcoef(1);
  if(degre2>0) res = pol.Fit(x,y,z,ez2,degre1,degre2,errcoef);
  else         res = pol.Fit(x,y,z,ez2,degre1,errcoef);
} else {
  if(degre2>0) res = pol.Fit(x,y,z,degre1,degre2);
  else         res = pol.Fit(x,y,z,degre1);
}
return res;
}

//////////////////////////////////////////////////////////////////////
//++
GeneralFitData GeneralFitData::FitResidus(GeneralFit& gfit)
//
//      Retourne une classe contenant les residus du fit ``gfit''.
//--
{
if(gfit.GetNVar()!=mNVar)
  throw(SzMismatchError("GeneralFitData::FitResidus: size mismatch\n"));
return gfit.DataResidus(true);
}

//++
GeneralFitData GeneralFitData::FitFunction(GeneralFit& gfit)
//
//      Retourne une classe contenant la function du fit ``gfit''.
//--
{
if(gfit.GetNVar()!=mNVar)
  throw(SzMismatchError("GeneralFitData::FitFunction: size mismatch\n"));
return gfit.DataFunction(true);
}

//////////////////////////////////////////////////////////////////////
//++
r_8* GeneralFitData::GetVec(int n, r_8* ret)   const
//
//      Retourne la donnee `n' dans le vecteur de double `ret'.
//| Par defaut, ret=NULL et le buffer interne de la classe est retourne
//| - Les donnees sont rangees dans l'ordre:
//|    x0,x1,x2,... ; ex0,ex1,ex2,...  ; y ; ey ; ok(0/1)
//|    |<- NVar ->| + |<-   NVar   ->| + 1 +  1 +  1
//|   Le vecteur ret a la taille 2*NVar+2+1
//--
{
int i;
if (ret == NULL) ret = BuffVar;
for(i=0; i<2*mNVar+3; i++)  ret[i] = 0.;
if (n >= mNData) return(ret);

memcpy(ret, mXP+n*mNVar, mNVar*sizeof(r_8));
if(mErrXP) memcpy(ret+mNVar, mErrXP+n*mNVar, mNVar*sizeof(r_8));
ret[2*mNVar] = mF[n];
ret[2*mNVar+1] = mErr[n];
ret[2*mNVar+2] = (double) mOK[n];
return(ret);
}

//++
r_4* GeneralFitData::GetVecR4(int n, r_4* ret)   const
//
//      Retourne la donnee `n' dans le vecteur de float `ret'
//      (meme commentaires que pour GetVec).
//--
{
if (ret == NULL) ret = BuffVarR4;
double *buff = new double[2*mNVar+3];
GetVec(n,buff);
for(int i=0;i<2*mNVar+3;i++) ret[i] = (float) buff[i];
delete [] buff;
return ret;
}

//////////////////////////////////////////////////////////////////////
//++
// int inline int GetSpaceFree() const
//      Retourne la place restante dans la structure (nombre de
//      donnees que l'on peut encore stoquer).
//--
//++
// inline int  NVar()       const
//      Retourne le nombre de variables Xi
//--
//++
// inline int  NData()
//      Retourne le nombre de donnees
//--
//++
// inline int  NDataGood()  const
//      Retourne le nombre de bonnes donnees (utilisees pour le fit)
//--
//++
// inline int  NDataAlloc() const
//      Retourne la place maximale allouee pour les donnees
//--
//++
// inline unsigned short int IsValid(int i) const
//      Retourne 1 si point valide, sinon 0
//--
//++
// inline bool HasXErrors()
//      Retourne ``true'' si il y a des erreurs sur les variables
//      d'abscisse, ``false'' sinon.
//--
//++
// inline double X1(int i) const
//      Retourne l'abscisse pour 1 dimension (y=f(x)) donnee I
//--
//++
// inline double X(int i) const
//      Retourne la 1er abscisse (X) pour (v=f(x,y,z,...)) donnee I
//--
//++
// inline double Y(int i) const
//      Retourne la 2sd abscisse (Y) pour (v=f(x,y,z,...)) donnee I
//--
//++
// inline double Z(int i) const
//      Retourne la 3ieme abscisse (Z) pour (v=f(x,y,z,...)) donnee I
//--
//++
// inline double Absc(int j,int i) const
//      Retourne la Jieme abscisse (Xj) pour (v=f(x0,x1,x2,...)) donnee I
//--
//++
// inline double Val(int i) const
//      Retourne la valeur de la Ieme donnee
//--
//++
// inline double EX1(int i) const
//      Retourne l'erreur (dx) sur l'abscisse pour 1 dimension (y=f(x)) donnee I
//--
//++
// inline double EX(int i) const
//      Retourne l'erreur (dx) sur la 1er abscisse (X) pour (v=f(x,y,z,...)) donnee I
//--
//++
// inline double EY(int i) const
//      Retourne l'erreur (dy) sur la 2sd abscisse (Y) pour (v=f(x,y,z,...)) donnee I
//--
//++
// inline double EZ(int i) const
//      Retourne l'erreur (dz) sur la 3ieme abscisse (Z) pour (v=f(x,y,z,...)) donnee I
//--
//++
// inline double EAbsc(int j,int i) const
//      Retourne l'erreur (dxj) sur la Jieme abscisse (Xj) pour (v=f(x0,x1,x2,...)) donnee I
//--
//++
// inline double EVal(int i) const {return mErr[i];}
//      Retourne l'erreur de la Ieme donnee
//--


//////////////////////////////////////////////////////////////////////
// ------- Implementation de  l interface NTuple  ---------

uint_4 GeneralFitData::NbLines() const
{
return(NData());
}

//++
uint_4 GeneralFitData::NbColumns() const
//
//      Retourne le nombre de colonnes du ntuple equivalent:
//| Exemple: on a une fonction sur un espace a 4 dimensions:
//| "x0,x1,x2,x3    , ex0,ex1,ex2,ex3    , y,   ey ,    ok"
//|   0  1  2  3        4   5   6   7      8     9      10
//|   |        |        |           |      |     |       |        
//|   0       nv-1     nv         2*nv-1  2*nv  2*nv+1  2*nv+2
//| soit 2*nvar+3 variables/colonnes.
//--
{
return(2*NVar()+3);
}

r_8 * GeneralFitData::GetLineD(int n) const
{
return(GetVec(n,NULL));
}

r_8 GeneralFitData::GetCell(int n, int k) const
{
if(k<0 || k>=2*NVar()+3) return 0.;
r_8 * val = GetVec(n,NULL);
return val[k];
}

r_8 GeneralFitData::GetCell(int n, string const & nom) const
{
int k = ColumnIndex(nom);
return(GetCell(n,k));
}

//++
void GeneralFitData::GetMinMax(int k, double& min, double& max)  const
//
//      Retourne le minimum et le maximum de la variable `k'.
//--
{
int var;
if(k<0 || k>=2*NVar()+3) return;
else if(k<NVar())      var = 10*k+2;          // Variable Xi
else if(k<2*NVar())    var = 10*(k-NVar())+3; // Variable EXi
else if(k==2*NVar())   var = 0;               // Variable Y
else if(k==2*NVar()+1) var = 1;               // Variable EY
else {min=0.; max=1.; return;}                // Variable Ok
GetMnMx(var,min,max);
return;
}

void GeneralFitData::GetMinMax(string const & nom, double& min, double& max)   const
{
int k = ColumnIndex(nom);
GetMinMax(k,min,max);
}

int GeneralFitData::ColumnIndex(string const & nom)  const
{
char str[64]; int k = -1;
strcpy(str,nom.c_str()); strip(str,'L',' ');
if(str[0]=='y') return 2*NVar();
if(str[0]=='o') return 2*NVar()+2;
if(str[0]=='x') {sscanf(str,"x%d",&k); return k;}
if(str[0]=='e')
  if(str[1]=='y') return 2*NVar()+1;
  else if(str[1]=='x') {sscanf(str,"ex%d",&k); return NVar()+k;}
return -1;
}

string GeneralFitData::ColumnName(int k) const
{
if(k==2*NVar())                return string("y");
else if(k==2*NVar()+1)         return string("ey");
else if(k==2*NVar()+2)         return string("ok");
else if(k<0 || k>=2*NVar()+3)  return string("");

char str[64] = "";
if(k<NVar()) sprintf(str,"x%d",k);
else if(k<2*NVar()) sprintf(str,"ex%d",k-NVar());
return string(str);
}

//++
string GeneralFitData::VarList_C(const char* nomx)  const
//
//      Retourne une chaine de caracteres avec la declaration des noms de 
//      variables. si "nomx!=NULL" , des instructions d'affectation
//      a partir d'un tableau "nomx[i]" sont ajoutees. 
//--
{
char buff[256];
string rets;
int i;
rets = "\ndouble";
for(i=0; i<mNVar; i++) {
  sprintf(buff," x%d, ex%d",i,i);
  rets += buff;
  if(i!=mNVar-1) rets += ","; else rets += ";\n";
}
sprintf(buff,"\ndouble y, ey, ok;\n");
rets += buff;
if (nomx) {
  for(i=0; i<mNVar; i++) {
    sprintf(buff,"x%d=%s[%d];\n", i, nomx, i);
    rets += buff;
  }
  for(i=0; i<mNVar; i++) {
    sprintf(buff,"ex%d=%s[%d];\n", i, nomx, mNVar+i);
    rets += buff;
  }
}
sprintf(buff,"y=%s[%d];\ney=%s[%d];\nok=%s[%d];\n"
       ,nomx,2*mNVar,nomx,2*mNVar+1,nomx,2*mNVar+2);
rets += buff;

return(rets);
}

///////////////////////////////////////////////////////////
// --------------------------------------------------------
//   Les objets delegues pour la gestion de persistance 
// --------------------------------------------------------
///////////////////////////////////////////////////////////


void ObjFileIO<GeneralFitData>::ReadSelf(PInPersist& is)
{
char strg[256];

if(dobj==NULL) dobj=new GeneralFitData;
  else         dobj->Delete();

// Lecture entete
is.GetLine(strg, 255);

// Ecriture des valeurs de definitions
int_4 nvar,ndatalloc,ndata,ndatagood;
is.Get(nvar);
is.Get(ndatalloc);
is.Get(ndata);
is.Get(ndatagood);
is.Get(dobj->mOk_EXP);
if(nvar<=0 || ndatalloc<=0 || ndata<=0 || ndatagood<0 || ndatalloc<ndata) return;

// Allocation de la place (attention Alloc efface mNData,mNDataGood);
dobj->Alloc(nvar,ndatalloc,-1);
dobj->mNData = ndata;
dobj->mNDataGood = ndatagood;

// Lecture des datas
is.GetLine(strg, 255);
int blen = dobj->mNVar + 3;
if(dobj->mOk_EXP) blen += dobj->mNVar;
double *buff = new double[blen];
for(int i=0;i<dobj->mNData;i++) {
  is.Get(buff, blen);
  int ip = i*dobj->mNVar;
  {for(int j=0;j<dobj->mNVar;j++)  dobj->mXP[ip+j] = buff[j];}
  dobj->mF[i] = buff[dobj->mNVar];
  dobj->mErr[i] = buff[dobj->mNVar+1];
  dobj->mOK[i] = (uint_2)(buff[dobj->mNVar+2]+0.01);
  if(dobj->mOk_EXP) {for(int j=0;j<dobj->mNVar;j++)
                        dobj->mErrXP[ip+j] = buff[dobj->mNVar+3+j];}
}
delete [] buff;

return;
}

void ObjFileIO<GeneralFitData>::WriteSelf(POutPersist& os) const
{
if (dobj == NULL)   return;
char strg[256];

// Ecriture entete pour identifier facilement
sprintf(strg,"GeneralFitData: NVar=%d NDataAlloc=%d NData=%d NDataGood=%d Ok_EXP=%d"
       ,dobj->mNVar,dobj->mNDataAlloc,dobj->mNData,dobj->mNDataGood,dobj->mOk_EXP);
os.PutLine(strg);

// Ecriture des valeurs de definitions
os.Put(dobj->mNVar);
os.Put(dobj->mNDataAlloc);
os.Put(dobj->mNData);
os.Put(dobj->mNDataGood);
os.Put(dobj->mOk_EXP);
if(dobj->mNVar<=0 || dobj->mNDataAlloc<=0 || dobj->mNData<=0 || dobj->mNDataGood<0) return;

// Ecriture des datas (on n'ecrit que mNData / mNDataAlloc)
sprintf(strg
       ,"GeneralFitData: Abscisses, Ordonnee, Erreur Ordonnee, Flag, Erreur Abscisses");
os.PutLine(strg);

int blen = dobj->mNVar + 3;
if(dobj->mOk_EXP) blen += dobj->mNVar;
double *buff = new double[blen];
for(int i=0;i<dobj->mNData;i++) {
  {for(int j=0;j<dobj->mNVar;j++) buff[j] = dobj->Absc(j,i);}
  buff[dobj->mNVar] = dobj->Val(i);
  buff[dobj->mNVar+1] = dobj->EVal(i);
  buff[dobj->mNVar+2] = (double) dobj->IsValid(i);
  if(dobj->mOk_EXP) {for(int j=0;j<dobj->mNVar;j++) buff[dobj->mNVar+3+j] = dobj->EAbsc(j,i);}
  os.Put(buff, blen);
}
delete [] buff;

return;
}


#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ObjFileIO<GeneralFitData>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class ObjFileIO<GeneralFitData>;
#endif