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

#if defined(__KCC__)
using std::string ;
#endif

#include "perrors.h"
#include "nbtri.h"
#include "generalfit.h"
#include "generaldata.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);
 } CATCHALL {
   cout<<"GeneralFitData::GeneralFitData Impossible d'allouer l'espace"<<endl;
   THROW_SAME;
 } ENDTRY
 END_CONSTRUCTOR
}

//++
GeneralFitData::GeneralFitData(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));
 } CATCHALL {
   cout<<"GeneralFitData::GeneralFitData Impossible d'allouer l'espace"<<endl;
   THROW_SAME;
 } ENDTRY

 // 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(char *flnm)
//
//	Constructeur par lecture d'un fichier ppersist..
//--
  : 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)
{
  PInPersist s(flnm);
  Read(s);
  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.
//--
{
DBASSERT( 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;
} CATCHALL {
  THROW_SAME;
} ENDTRY
}

//////////////////////////////////////////////////////////////////////
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.
//--
{
 DBASSERT(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
//--
{
 DBASSERT(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
//--
{
 DBASSERT(i1 >= 0 && i1 < mNData);
 DBASSERT(i2 >= 0 && i2 < mNData);
 DBASSERT(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]).
//--
{
 DBASSERT(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.
//--
{
 DBASSERT(i1 >= 0 && i1 < mNData);
 DBASSERT(i2 >= 0 && i2 < mNData);
 DBASSERT(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}
//--
{
 DBASSERT(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}
//--
{
 DBASSERT(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)}
//--
{
 DBASSERT(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)}
//--
{
 DBASSERT(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)}
//--
{
 DBASSERT(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)}
//--
{
 DBASSERT(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).
//--
{
 DBASSERT(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.
//--
{
 DBASSERT(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
//--
{
 DBASSERT(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
//--
{
 DBASSERT(mNData>0);

 PrintData(0,mNData-1);
}

//////////////////////////////////////////////////////////////////////
//++
int GeneralFitData::GetMinMax(int var,int& imin,int& imax)
//
//	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::GetMinMax(int var,double& min,double& max)
//
//	Retourne le minimum et le maximum de la variable ``var''
//	(cf commentaires GetMinMax).
//--
{
min = 1.; max = -1.;
int imin,imax;
int ntest = GetMinMax(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 GetMinMax).
//| - 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 GetMinMax).
//| - 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) return NULL;
return gfit.DataResidus(true);
}

//++
GeneralFitData* GeneralFitData::FitFunction(GeneralFit& gfit)
//
//	Retourne une classe contenant la function du fit ``gfit''.
//--
{
if(gfit.GetNVar()!=mNVar) return NULL;
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;
}

//////////////////////////////////////////////////////////////////////
//++
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);
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::WriteSelf(POutPersist& s)  const
//
//	Ecriture ppersist
//--
{
char strg[256];

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

// Ecriture des valeurs de definitions
s.PutI4(mNVar);
s.PutI4(mNDataAlloc);
s.PutI4(mNData);
s.PutI4(mNDataGood);
s.PutU2(mOk_EXP);
if(mNVar<=0 || mNDataAlloc<=0 || mNData<=0 || mNDataGood<0) return;

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

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

return;
}

//////////////////////////////////////////////////////////////////////
//++
void GeneralFitData::ReadSelf(PInPersist& s)
//
//	Lecture ppersist
//--
{
char strg[256];

Delete();

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

// Ecriture des valeurs de definitions
int_4 nvar,ndatalloc,ndata,ndatagood;
s.GetI4(nvar);
s.GetI4(ndatalloc);
s.GetI4(ndata);
s.GetI4(ndatagood);
s.GetU2(mOk_EXP);
if(nvar<=0 || ndatalloc<=0 || ndata<=0 || ndatagood<0 || ndatalloc<ndata) return;

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

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

return;
}

//////////////////////////////////////////////////////////////////////
//++
// 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
//--