source: Sophya/trunk/SophyaLib/NTools/generaldata.cc@ 2743

#include "sopnamsp.h"
#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#ifndef NO_VALUES_H
#include <values.h>
#endif
#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 SOPHYA::GeneralFitData
  \ingroup NTools
  Classe de stoquage de donnees a plusieurs variables avec erreur
  sur l'ordonnee et sur les abscisses (options) :

  \f$ {x0(i),Ex0(i), x1(i),Ex1(i), x2(i),Ex2(i) ... ; Y(i),EY(i)} \f$
  \verbatim
   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]
  \endverbatim
*/

//////////////////////////////////////////////////////////////////////
/*!
  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.
*/
GeneralFitData::GeneralFitData(unsigned int nVar, unsigned int ndatalloc, uint_2 errx)
  : 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;
}
}

/*!
  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.
*/
GeneralFitData::GeneralFitData(const GeneralFitData& data, bool clean)
  : 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++;
   }
 }

}

/*!
  Constructeur par defaut.
*/
GeneralFitData::GeneralFitData()
  : 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)
{
}

/*!
  Destructeur
*/
GeneralFitData::~GeneralFitData()
{
 Delete();
}

//////////////////////////////////////////////////////////////////////
/*!
  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.
*/
void GeneralFitData::Alloc(unsigned int nVar, unsigned int ndatalloc, int_2 errx)
{
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"));
}
}

//////////////////////////////////////////////////////////////////////
/*!
   Gestion des des-allocations
*/
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;}
}

//////////////////////////////////////////////////////////////////////
/*!
  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.
*/
void GeneralFitData::SetDataPtr(int ptr)
{
 ASSERT(ptr >= 0 && ptr < mNDataAlloc);
 mNData = ptr;
 mNDataGood = 0;
 if(ptr==0) return;
 for(int i=0;i<mNData;i++) if(mOK[i]) mNDataGood++;
}

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

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

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

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

//////////////////////////////////////////////////////////////////////
/*!
  Pour re-valider le point numero i ([0,NData]).
*/
void GeneralFitData::ValidData(int i)
{
 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++;
}

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

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

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

//////////////////////////////////////////////////////////////////////
/*!
  Pour redefinir un point a \f$ {x,[errx] ; f,err} \f$
*/
void GeneralFitData::RedefineData1(int i,double x,double f,double err,double errx)
{
 RedefineData(i,&x,f,err,&errx);
}

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

/*!
  Pour redefinir un point a
  \f$ {xp[0],[errxp[0]], xp[1],[errxp[1]], xp[2],[errxp[2]],... ; f,err} \f$ 
*/
void GeneralFitData::RedefineData(int i,double* xp,double f,double err,double* errxp)
{
 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--;}
 }
}

//////////////////////////////////////////////////////////////////////
/*!
  Pour ajouter un point \f$ {x,[errx] ; f,err} \f$
*/
void GeneralFitData::AddData1(double x, double f, double err, double errx)
{
 AddData(&x,f,err,&errx);
}

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

/*!
  Pour ajouter un point
  \f$ {xp[0],[errxp[0]], xp[1],[errxp[1]], xp[2],[errxp[2]],... ; f,err} \f$
*/
void GeneralFitData::AddData(double* xp, double f, double err,double* errxp)
{
 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++;
}

/*!
  Pour ajouter un point
 \f$ {xp[0],[errxp[0]], xp[1],[errxp[1]], xp[2],[errxp[2]],... ; f,err} \f$ 
*/
void GeneralFitData::AddData(float* xp, float f, float err, float* errxp)
{
 {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);
}

//////////////////////////////////////////////////////////////////////
/*!
  Pour remplir la structure de donnees d'un seul coup avec
  \f$ {x(i),[errx(i)] ; f(i),err(i)} \f$
*/
void GeneralFitData::SetData1(int nData
            , double* x, double* f, double *err, double *errx)
{
 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);
 }
}

/*!
  Pour remplir la structure de donnees d'un seul coup avec
  \f$ {x(i),[errx(i)] ; f(i),err(i)} \f$
*/
void GeneralFitData::SetData1(int nData
            , float* x, float* f, float* err, float *errx)
{
 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);
 }
}

/*!
  Pour remplir la structure de donnees d'un seul coup avec
  \f$ {x(i),[errx(i)], y(i),[erry(i)], ; f(i),err(i)} \f$
*/
void GeneralFitData::SetData2(int nData, double* x, double* y, double* f
                             ,double *err,double *errx,double *erry)
{
 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);
 }
}

/*!
  Pour remplir la structure de donnees d'un seul coup avec
 \f$ {x(i),[errx(i)], y(i),[erry(i)], ; f(i),err(i)} \f$
*/
void GeneralFitData::SetData2(int nData, float* x, float* y, float* f
                             ,float *err,float *errx,float *erry)
{
 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);
 }
}

/*!
  Pour remplir la structure de donnees d'un seul coup avec
  \f$ {X0(i),[EX0(i)], X1(i),[EX1(i)], X2(i),[EX2(i)], ... ; F(i),Err(i)} \f$

  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).
  \verbatim
      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).
  \endverbatim
*/
void GeneralFitData::SetData(int nData,double** xp, double *f
                            , double *err, double** 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);
 }
}

/*!
  Voir commentaire ci-dessus.
*/
void GeneralFitData::SetData(int nData,float** xp, float *f
                            , float *err, float** 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] = (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);
 }
}

//////////////////////////////////////////////////////////////////////
/*!
  Impression de l'etat de la structure de donnees
*/
void GeneralFitData::PrintStatus() const
{
  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;
}

/*!
  Impression du point i
*/
void GeneralFitData::PrintData(int i) const
{
 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";
}

/*!
  Impression des points i1 a i2
*/
void GeneralFitData::PrintData(int i1,int i2) const
{
 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;
}

/*!
  Impression de tous les points
*/
void GeneralFitData::PrintData() const
{
 ASSERT(mNData>0);

 PrintData(0,mNData-1);
}

/*!
  Impression de l'etat de la structure de donnees avec bornes sur "s"
*/
void GeneralFitData::Show(ostream& os) const
{
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;
}

//////////////////////////////////////////////////////////////////////
/*!
  Retourne les numeros des points de valeurs minimum et maximum
  de la variable ``var'':
  \verbatim
    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.
  \endverbatim
*/
int GeneralFitData::GetMnMx(int var,int& imin,int& imax) const
{
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=1., max=-1.;
int ntest = 0;
for(int i=0;i<mNData;i++) {
  if( ! IsValid(i) ) continue;
  double v;
  if(var==1)          v = EVal(i);
    else if(var==2)   v = Absc(ix,i);
      else if(var==3) v = EAbsc(ix,i);
        else          v = Val(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;
}

/*!
  Retourne le minimum et le maximum de la variable ``var''
  (cf commentaires GetMnMx).
*/
int GeneralFitData::GetMnMx(int var,double& min,double& max) const
{
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;
}

//////////////////////////////////////////////////////////////////////
/*!
//
  Retourne la moyenne et le sigma de la variable ``var''
  (cf commentaires GetMnMx).
  \verbatim
    - 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.
  \endverbatim
*/
int GeneralFitData::GetMeanSigma(int var,double& mean,double& sigma,double min,double max) const
{
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==1)          v = EVal(i);
    else if(var==2)   v = Absc(ix,i);
      else if(var==3) v = EAbsc(ix,i);
        else          v = Val(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;
}

/*!
  Retourne le mode de la variable ``var''
  (cf commentaires GetMnMx).
  \verbatim
    - 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.
  \endverbatim
*/
int GeneralFitData::GetMoMeMed(int var,double& mode,double& mean,double& median,
                               double min,double max,double coeff) const
{
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==1)          v = EVal(i);
    else if(var==2)   v = Absc(ix,i);
      else if(var==3) v = EAbsc(ix,i);
        else          v = Val(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;
}

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

//////////////////////////////////////////////////////////////////////
/*!
  Pour fiter un polynome de degre ``degre''. On fite
  Y=f(X-xc) 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''. On re-centre les abscisses X de ``xc''.
  \verbatim
  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 differents de NDataGood
  \endverbatim
*/
double GeneralFitData::PolFit(int varx,Poly& pol,int degre,bool ey,double xc) const
{
if(degre<0) return -1.;
if(varx<0 || varx>=mNVar) return -2.;
if(mNDataGood<=0) return -4.;
TVector<r_8> x(mNDataGood);
TVector<r_8> y(mNDataGood);
TVector<r_8> 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) - xc;
  y(ntest) = Val(i);
  if(ey) ey2(ntest) = EVal(i)*EVal(i);
  ntest++;
}
double res = 0.;
if(ey) {
  TVector<r_8> errcoef(1);
  res = pol.Fit(x,y,ey2,degre,errcoef);
} else {
  res = pol.Fit(x,y,degre);
}
return res;
}

/*!
  Pour fiter un polynome de degre ``degre1''. On fite
  Z=f(X-xc,Y-yc) 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''. On re-centre les abscisses X de ``xc''
  et Y de ``yc''.
  \verbatim
    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
  \endverbatim
*/
double GeneralFitData::PolFit(int varx,int vary,Poly2& pol,int degre1,int degre2,bool ez
                             ,double xc,double yc) const
{
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.;
TVector<r_8> x(mNDataGood);
TVector<r_8> y(mNDataGood);
TVector<r_8> z(mNDataGood);
TVector<r_8> 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) - xc;
  y(ntest) = Absc(vary,i) - yc;
  z(ntest) = Val(i);
  if(ez) ez2(ntest) = EVal(i)*EVal(i);
  ntest++;
}
double res = 0.;
if(ez) {
  TVector<r_8> 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;
}

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

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

//////////////////////////////////////////////////////////////////////
/*!
//
  Retourne la donnee `n' dans le vecteur de double `ret'.
  \verbatim
    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
  \endverbatim
*/
r_8* GeneralFitData::GetVec(int n, r_8* ret)   const
{
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);
}

/*!
  Retourne la donnee `n' dans le vecteur de float `ret'
  (meme commentaires que pour GetVec).
*/
r_4* GeneralFitData::GetVecR4(int n, r_4* ret)   const
{
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;
}

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

/*!
  Retourne le nombre de ligne = NData() (pour interface NTuple)
*/
sa_size_t GeneralFitData::NbLines() const
{
return(NData());
}

/*!
  Retourne le nombre de colonnes du ntuple equivalent:
  \verbatim
    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.
  \endverbatim
  (pour interface NTuple)
*/
sa_size_t GeneralFitData::NbColumns() const
{
return(2*NVar()+3);
}

//! Pour interface NTuple
r_8 * GeneralFitData::GetLineD(sa_size_t n) const
{
return(GetVec(n,NULL));
}

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

//! Pour interface NTuple
r_8 GeneralFitData::GetCell(sa_size_t n, string const & nom) const
{
sa_size_t k = ColumnIndex(nom);
return(GetCell(n,k));
}

/*!
  Retourne le minimum et le maximum de la variable `k' (pour interface NTuple).
*/
void GeneralFitData::GetMinMax(sa_size_t k, double& min, double& max)  const
{
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;
}

//! Pour interface NTuple
void GeneralFitData::GetMinMax(string const & nom, double& min, double& max)   const
{
sa_size_t k = ColumnIndex(nom);
GetMinMax(k,min,max);
}

//! Pour interface NTuple
sa_size_t 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;
}

//! Pour interface NTuple
string GeneralFitData::ColumnName(sa_size_t 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);
}

/*!
  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 (pour interface NTuple).
*/
string GeneralFitData::VarList_C(const char* nomx)  const
{
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);
}


//! Compute errors according to specifications
/*!
  \param val : value of the function
  \param err : value of the default error
  \param errtype : type of error according to enum FitErrType (def=DefaultError)
  \param errscale : scaling (or value) of the error (def=1.)
  \param errmin : minimum value of the error (def=0.)
  \param nozero : if true, do not return negative errors but
                  set them to zero (def=false)
  \return : return the error computed according to specifications
  \verbatim
  - val is the value to be fitted ex: val = f(x,y,...)
  - err is the error by default we want to set.
  - errtype = DefaultError  : errtmp = errscale*err
    errtype = ConstantError : errtmp = errscale
    errtype = SqrtError     : errtmp = errscale*sqrt(|val|)
    errtype = ProporError   : errtmp = errscale*|val|
  - errscale <=0 then errscale=1
  - errmin >=0 if errtmp>0  return max(errtmp,errmin)
               if errtmp<=0 return errtmp
    errmin <0  if errtmp>0  return max(errtmp,|errmin|)
               if errtmp<=0 return |errmin|
  \endverbatim
 */
double GeneralFitData::ComputeError(double val,double err,FitErrType errtype
                                   ,double errscale,double errmin,bool nozero)
{
 bool errminneg=false;
 if(errmin<0.) {errminneg=true; errmin*=-1.;}
 if(errscale<0.) errscale=1.;

 // Choix du type d'erreur
 if(errtype==ConstantError)    err = errscale;
 else if(errtype==SqrtError)   err = errscale*sqrt(fabs(val));
 else if(errtype==ProporError) err = errscale*fabs(val);

 // Gestion du minimum a partir de la valeur calculee precedemment "err"
 // Ex1: errmin=1.,  err=10.     ==> 10.
 //                  err=0.5     ==> 1.
 //                  err=0.      ==> 0.
 //                  err=-2.     ==> -2.
 // Ex2: errmin=-1., err=10.     ==> 10.
 //                  err=0.5     ==> 1.
 //                  err=0.      ==> 1.
 //                  err=-2.     ==> 11.
 if(err>0.) err = (err>errmin) ? err: errmin;
 else if(errminneg) err = errmin;

 // ne pas retourner d'erreurs negatives si demande
 if(nozero && err<0.) err=0.;

 return err;
}

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


DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
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;
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
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