source: Sophya/trunk/SophyaLib/HiStats/ntuple.cc@ 975

#include <stdio.h>
#include <string.h>

#include "strutil.h"
#include "perrors.h"
#include "ntuple.h"


#define BADVAL -1.e19
#define LENNAME 8
#define LENNAME1  (LENNAME+1)

//++
// Class        NTuple
// Lib  Outils++ 
// include      ntuple.h
//
//      Classe de ntuples
//--
//++
// Links        Parents
// PPersist
// NTupleInterface
//--

/* --Methode-- */
//++
NTuple::NTuple()
//
//      Createur par defaut
//--
{
mNVar = mNEnt = mBlk = mNBlk = 0;
mVar = NULL;
mVarD = NULL;
mNames = NULL;
mInfo = NULL;
}


//++
NTuple::NTuple(int nvar, char** noms, int blk)
//
//      Createur d'un ntuple de `nvar' variables dont les
//      noms sont dans le tableau de cahines de caracteres `noms'
//      avec `blk' d'evenements par blocks.
//--
{
mNVar = mNEnt = mBlk = mNBlk = 0;
mVar = NULL;
mVarD = NULL;
mNames = NULL;
mInfo = NULL;
if (nvar <= 0)  THROW(sizeMismatchErr); 
mNVar = nvar;
mVar = new r_4[nvar];
mVarD = new r_8[nvar];
if (blk < 10) blk = 10;
mBlk = blk;
// On prend des noms de LENNAME char pour le moment
mNames = new char[nvar*LENNAME1];
r_4* pt = new r_4[nvar*blk];
mNBlk = 1;
mPtr.push_back(pt);
int i;
for(i=0; i<nvar; i++)
  { strncpy(mNames+i*LENNAME1, noms[i], LENNAME);  
  mNames[i*LENNAME1+LENNAME] = '\0'; }
return;
}

//                                       cmv 8/10/99
//++
NTuple::NTuple(const NTuple& NT)
//
//      Createur par copie (clone).
//--
: mNVar(0), mNEnt(0), mBlk(0), mNBlk(0)
, mVar(NULL), mVarD(NULL), mNames(NULL), mInfo(NULL)
{
if(NT.mNVar<=0) return; // cas ou NT est cree par defaut
mNVar = NT.mNVar;
mBlk = NT.mBlk;
mVar = new r_4[NT.mNVar];
mVarD = new r_8[NT.mNVar];
mNames = new char[NT.mNVar*LENNAME1];

int i;
r_4* pt = new r_4[mNVar*mBlk];
mNBlk = 1; mPtr.push_back(pt);

for(i=0;i<mNVar;i++) strcpy(mNames+i*LENNAME1,NT.NomIndex(i)); 

if(NT.mInfo!=NULL) {mInfo = new DVList; *mInfo = *(NT.mInfo);}

if(NT.mNEnt<=0) return;
for(i=0;i<NT.mNEnt;i++) {r_4* x=NT.GetVec(i,NULL); Fill(x);}

return;
}

/* --Methode-- */
//++
NTuple::NTuple(char *flnm)
//
//      Createur lecture fichier ppersist.
//--
{
mNVar = mNEnt = mBlk = mNBlk = 0;
mVar = NULL;
mVarD = NULL;
mNames = NULL;
mInfo = NULL;
PInPersist s(flnm);
ObjFileIO<NTuple> fiont(this);
fiont.Read(s);
}

/* --Methode-- */
NTuple::~NTuple()
{
Clean();
}

/* --Methode-- */
void NTuple::Clean()
{
if (mVar) delete[] mVar;
if (mVarD) delete[] mVarD;
if (mNames)  delete[] mNames;
if (mInfo) delete mInfo;
int i;
if(mNBlk>0) for(i=0; i<mNBlk; i++)  delete[] mPtr[i];
mPtr.erase(mPtr.begin(), mPtr.end());
mNVar = mNEnt = mBlk = mNBlk = 0;
mVar = NULL;
mVarD = NULL;
mNames = NULL;
mInfo = NULL;
return;
}

/* --Methode--        cmv 08/10/99 */
//++
NTuple& NTuple::operator = (const NTuple& NT)
//
//      Operateur egal (clone).
//--
{
if(this == &NT) return *this;
Clean();
if(NT.mNVar<=0) return *this; // cas ou NT est cree par defaut
mNVar = NT.mNVar;
mBlk = NT.mBlk;
mVar = new r_4[NT.mNVar];
mVarD = new r_8[NT.mNVar];
mNames = new char[NT.mNVar*LENNAME1];

int i;
r_4* pt = new r_4[mNVar*mBlk];
mNBlk = 1; mPtr.push_back(pt);

for(i=0;i<mNVar;i++) strcpy(mNames+i*LENNAME1,NT.NomIndex(i)); 

if(NT.mInfo!=NULL) {mInfo = new DVList; *mInfo = *(NT.mInfo);}

if(NT.mNEnt<=0) return *this;
for(i=0;i<NT.mNEnt;i++) {r_4* x=NT.GetVec(i,NULL); Fill(x);}

// En fait il faudrait un createur par copie qui partage les donnees
// quand l'objet est temporaire... trop complique A FAIRE !  cmv.
return *this;
}

/* --Methode-- */
//++
void  NTuple::Fill(r_4* x)
//
//      Remplit le ntuple avec le tableau cd reels `x'.
//--
{
int numb = mNEnt/mBlk;
if (numb >= mNBlk) {
  r_4* pt = new r_4[mNVar*mBlk];
  mNBlk++;
  mPtr.push_back(pt);
}
int offb = mNEnt-numb*mBlk;
memcpy((mPtr[numb]+offb*mNVar), x, mNVar*sizeof(r_4));
mNEnt++;
return;
}


/* --Methode-- */
//++
float NTuple::GetVal(int n, int k)  const
//
//      Retourne la valeur de la variable `k' de l'evenement `n'.
//--
{
if (n >= mNEnt)   return(BADVAL);
if ( (k < 0) || (k >= mNVar) )    return(BADVAL);
int numb = n/mBlk;
int offb = n-numb*mBlk;
return(*(mPtr[numb]+offb*mNVar+k));
}

/* --Methode-- */
//++
void NTuple::SetVal(int n, int k, float value)  
//
//      initialise la valeur de la variable `k' de l'evenement `n'.
//--
{
if (n >= mNEnt)   return;
if ( (k < 0) || (k >= mNVar) )    return;
int numb = n/mBlk;
int offb = n-numb*mBlk;
*(mPtr[numb]+offb*mNVar+k) = value;
}

/* --Methode-- */
//++
int NTuple::IndexNom(const char* nom)  const
//
//      Retourne le numero de la variable de nom `nom'.
//--
{
int i;
for(i=0; i<mNVar; i++)  
  if ( strcmp(nom, mNames+i*LENNAME1) == 0)  return(i);
return(-1);
}


static char nomretour[2*LENNAME1];
/* --Methode-- */
//++
char* NTuple::NomIndex(int k)  const
//
//      Retourne le nom de la variable numero 'k'
//--
{
nomretour[0] = '\0';
if ((k >= 0) && (k < mNVar))  strcpy(nomretour, mNames+k*LENNAME1);
return(nomretour);
}

  
/* --Methode-- */
//++
r_4* NTuple::GetVec(int n, r_4* ret)   const
//
//      Retourne l'evenement `n' dans le vecteur `ret'.
//--
{
int i;
if (ret == NULL)   ret = mVar;
if (n >= mNEnt) {
  for(i=0; i<mNVar; i++)   ret[i] = BADVAL;
  return(ret);
}
  
int numb = n/mBlk;
int offb = n-numb*mBlk;
memcpy(ret, (mPtr[numb]+offb*mNVar), mNVar*sizeof(r_4));
return(ret);
}

/* --Methode-- */
//++
r_8* NTuple::GetVecD(int n, r_8* ret)   const
//
//      Retourne l'evenement `n' dans le vecteur `ret'.
//--
{
int i;
if (ret == NULL)   ret = mVarD;
float *  fr = GetVec(n);
for(i=0; i<mNVar; i++)   ret[i] = fr[i];
return(ret);
}



/* --Methode-- */
//++
DVList&  NTuple::Info()
//
//      Renvoie une référence sur l'objet DVList Associé
//--
{
if (mInfo == NULL)  mInfo = new DVList;
return(*mInfo);
}

/* --Methode-- */
//++
void  NTuple::Print(int num, int nmax)  const
//
//      Imprime `nmax' evenements a partir du numero `num'.
//--
{
int i,j;

printf("Num     ");
for(i=0; i<mNVar; i++)  printf("%8s ", mNames+i*LENNAME1);
putchar('\n');

if (nmax <= 0)  nmax = 1;
if (num < 0)  num = 0;
nmax += num;
if (nmax > mNEnt) nmax = mNEnt;
for(i=num; i<nmax; i++) {
  GetVec(i, NULL); 
  printf("%6d  ", i);  
  for(j=0; j<mNVar; j++)  printf("%8g ", (float)mVar[j]);
  putchar('\n');
}
return;
}

/* --Methode-- */
//++
void  NTuple::Show(ostream& os)  const
//
//      Imprime l'information generale sur le ntuple.
//--
{
os << "NTuple: NVar= " << mNVar << " NEnt=" << mNEnt  
   << " (Blk Sz,Nb= " << mBlk << " ," << mNBlk << ")\n";
os << "            Variables       Min      Max       \n";
int i;
double min, max;
char buff[128];
for(i=0; i<mNVar; i++) {
  GetMinMax(i, min, max);
  sprintf(buff, "%3d  %16s  %10lg  %10lg \n", i, mNames+i*LENNAME1, min, max);
  os << (string)buff ;
  }
os << endl;
}


/* --Methode-- */
//++
int  NTuple::FillFromASCIIFile(string const& fn, float defval)
//
//      Remplit le ntuple a partir d'un fichier ASCII.
//      Renvoie le nombre de lignes ajoutees.
//--
{
if (NbColumns() < 1)  { 
  cout << "NTuple::FillFromASCIIFile() Ntuple has " << NbColumns() << " columns" << endl;
  return(-1);
  }
FILE * fip = fopen(fn.c_str(), "r");
if (fip == NULL) {
  cout << "NTuple::FillFromASCIIFile() Error opening file " << fn << endl;
  return(-2);
  }

char lineb[4096];
char *line;
char* ccp;
int j,kk;
int postab, posb;
float* xv = new float[NbColumns()];

int nlread = 0;
int nvar = NbColumns();
// On boucle sur toutes les lignes
while (fgets(lineb,4096,fip) != NULL) {
  lineb[4095] = '\0';
  j = 0; line = lineb;
//  On enleve les espaces et tab de debut 
  while ( (line[j] != '\0') && ((line[j] == ' ') || (line[j] == '\t')) )  j++;
  line = lineb+j; 
// Il faut que le premier caractere non-espace soit un digit, ou + ou - ou .
  if (!( isdigit(line[0]) || (line[0] == '+') || (line[0] == '-') || (line[0] == '.') ))  continue;
  ccp = line;
  for(kk=0; kk<nvar; kk++)  xv[kk] = defval;
  for(kk=0; kk<nvar; kk++) {
// Les mots sont separes par des espaces ou des tab
    postab = posc(ccp, '\t' );
    posb = posc(ccp, ' ' );
    if (postab >= 0) { 
       if (posb < 0) posb = postab;
       else if (postab < posb)  posb = postab; 
       } 
    if (posb >= 0)  ccp[posb] = '\0';
    if ( isdigit(line[0]) || (line[0] == '+') || (line[0] == '-') || (line[0] == '.') ) 
      xv[kk] = atof(ccp);
    if (posb < 0)  break;
    ccp += posb+1;   j = 0;
    while ( (ccp[j] != '\0') && ((ccp[j] == ' ') || (ccp[j] == '\t')) )  j++;
    ccp += j;
    }
  Fill(xv);
  nlread++;
  }

delete[] xv;
cout << "NTuple::FillFromASCIIFile( " << fn << ") " << nlread << " fill from file " << endl;
return(nlread);
}


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

/* --Methode-- */
uint_4 NTuple::NbLines() const
{
return(NEntry());
}
/* --Methode-- */
uint_4 NTuple::NbColumns() const
{
return(NVar());
}

/* --Methode-- */
r_8 * NTuple::GetLineD(int n) const
{
return(GetVecD(n));
}

/* --Methode-- */
r_8 NTuple::GetCell(int n, int k) const
{
return(GetVal(n, k));
}

/* --Methode-- */
r_8 NTuple::GetCell(int n, string const & nom) const
{
return(GetVal(n, nom.c_str()));
}

/* --Methode-- */
//++
void  NTuple::GetMinMax(int k, double& min, double& max)  const
//
//      Retourne le minimum et le maximum de la variable `k'.
//--
{
min = 9.e19; max = -9.e19;
if ( (k < 0) || (k >= mNVar) )    return;
int jb,ib,i;
double x;
i=0;
for(jb=0; jb< mNBlk; jb++)
  for(ib=0; ib< mBlk; ib++) {
    if (i >= mNEnt)  break;
    i++;
    x = *(mPtr[jb]+ib*mNVar+k);
    if(i==1) {min = x; max = x;}
    if (x < min)  min = x;
    if (x > max)  max = x;
  }
return;
}

/* --Methode-- */
void NTuple::GetMinMax(string const & nom, double& min, double& max)   const
{
GetMinMax(IndexNom(nom.c_str()), min, max);
}

/* --Methode-- */
int NTuple::ColumnIndex(string const & nom)  const
{
return(IndexNom(nom.c_str()));
}

/* --Methode-- */
string NTuple::ColumnName(int k) const
{
return(NomIndex(k));
}

/* --Methode-- */
//++
string NTuple::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. 
//--
{
string rets="";
int i;
for(i=0; i<mNVar; i++) {
  if ( (i%5 == 0) && (i > 0) )  rets += ";";  
  if (i%5 == 0)   rets += "\ndouble "; 
  else rets += ",";
  rets += mNames+i*LENNAME1;
  }
rets += "; \n";
if (nomx) { 
  char buff[256];
  for(i=0; i<mNVar; i++) {
    sprintf(buff,"%s=%s[%d]; ",  mNames+i*LENNAME1, nomx, i);
    rets += buff;
    if ( (i%3 == 0) && (i > 0) )  rets += "\n"; 
    }
  }

return(rets);
}


/* --Methode-- */
//++
string NTuple::LineHeaderToString() const 
//      Retourne une chaine de caracteres avec la liste des noms de  
//      variables, utilisables pour une impression
//--
{
char buff[32];
string rets=" Num    ";
for(int i=0; i<mNVar; i++) {
  sprintf(buff, "%8s ", mNames+i*LENNAME1);
  rets += buff;
  }
rets += '\n';
return(rets);
}

/* --Methode-- */
//++
string NTuple::LineToString(int n) const
//      Retourne une chaine de caracteres avec le contenu de la ligne "n"
//      utilisable pour une impression
//--
{
char buff[32];
double* val;
val = GetLineD(n);
sprintf(buff,"%6d: ",n);  
string rets=buff;
int i;
for(i=0; i<mNVar; i++) {
  sprintf(buff, "%8.3g ", val[i]);
  rets += buff;
  }
rets += '\n';
return(rets);
}


/* --Methode-- */ 
//++
void   ObjFileIO<NTuple>::WriteSelf(POutPersist& s)  const
//
//      Ecriture ppersist du ntuple.
//--
{
char strg[256];
if (dobj->mInfo)  sprintf(strg, "NVar=%6d  NEnt=%9d  BlkSz=%6d NBlk=%6d  HasInfo", 
                          (int)dobj->mNVar, (int)dobj->mNEnt, (int)dobj->mBlk, (int)dobj->mNBlk);
else sprintf(strg, "NVar=%6d  NEnt=%9d  BlkSz=%6d NBlk=%6d ", 
                   (int)dobj->mNVar, (int)dobj->mNEnt, (int)dobj->mBlk, (int)dobj->mNBlk);
s.PutLine(strg);
s.PutI4(dobj->mNVar);
s.PutBytes(dobj->mNames, dobj->mNVar*LENNAME1);
s.PutI4(dobj->mNEnt);
s.PutI4(dobj->mBlk);
s.PutI4(dobj->mNBlk);
if (dobj->mInfo)  s << (*(dobj->mInfo));
int jb;
for(jb=0; jb<dobj->mNBlk; jb++)
  s.PutR4s(dobj->mPtr[jb], dobj->mNVar*dobj->mBlk); 
return;
}

/* --Methode-- */
//++
void  ObjFileIO<NTuple>::ReadSelf(PInPersist& s)
//
//      Lecture ppersist du ntuple.
//--
{

dobj->Clean();

char strg[256];
s.GetLine(strg, 255);
// Pour savoir s'il y avait un DVList Info associe
bool hadinfo = false;
if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;

s.GetI4(dobj->mNVar);
dobj->mNames = new char[dobj->mNVar*LENNAME1];
dobj->mVar = new r_4[dobj->mNVar];
dobj->mVarD = new r_8[dobj->mNVar];
s.GetBytes(dobj->mNames, dobj->mNVar*LENNAME1);
s.GetI4(dobj->mNEnt);
s.GetI4(dobj->mBlk);
s.GetI4(dobj->mNBlk);

if (hadinfo) {    // Lecture eventuelle du DVList Info
  if (dobj->mInfo == NULL)  dobj->mInfo = new DVList;
  s >> (*(dobj->mInfo));
  }

int jb; 
for(jb=0; jb<dobj->mNBlk; jb++) {
  r_4* pt = new r_4[dobj->mNVar*dobj->mBlk];
  dobj->mPtr.push_back(pt);
  s.GetR4s(dobj->mPtr[jb], dobj->mNVar*dobj->mBlk); 
}

}


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

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