source: Sophya/trunk/ArchTOIPipe/Kernel/fitstoirdr.cc@ 2998

// ArchTOIPipe           (C)     CEA/DAPNIA/SPP IN2P3/LAL
//                               Eric Aubourg
//                               Christophe Magneville
//                               Reza Ansari
// $Id: fitstoirdr.cc,v 1.43 2003-10-13 20:48:37 aubourg Exp $

#include "fitstoirdr.h"
#include "toimanager.h"
#include <sched.h>


FITSTOIReader::FITSTOIReader(string fn,int buff_sz) {
  fname = fn;
  allfn.push_back(fn);
  Buff_Sz = (buff_sz>0) ? buff_sz: 1000;
  cout<<"FITSTOIReader::FITSTOIReader"<<endl;
  cout<<"FITSTOIReader::inited "<<inited<<" bsz="<<Buff_Sz<<endl;
  name = "rdr";
  fptr = NULL;
  totnscount = 0;

  implicitSN = false;
  implicitSNStart = 0;

  // Variables rajoutee pour gerer les fichiers de flag de LevelS Reza 18/6/2002
  sepFlagfile = false;
  fptrflg = NULL;
}

FITSTOIReader::~FITSTOIReader() {
}

void FITSTOIReader::setImplicitSN(long snStart) {
  implicitSN = true;
  implicitSNStart = snStart;
}

void FITSTOIReader::setBufferSize(int buffsz)
{
  Buff_Sz = (buffsz>0) ? buffsz: 1024;
  return;
}

pthread_mutex_t fits_mutex = PTHREAD_MUTEX_INITIALIZER;
void fits_lock();
void fits_unlock();

void fits_lock() {
  pthread_mutex_lock(&fits_mutex);
}
void fits_unlock() {
  pthread_mutex_unlock(&fits_mutex);
}

void FITSTOIReader::openFile(string fn) {
  fits_lock();
  if (fptr) {
    fits_close_file(fptr,&fstatus);
    fptr = NULL;
  }
  fname = fn;
  cout << "FITSTOIReader::open FileName=" << fname << endl;
  fstatus = 0;
  // Open file
  fits_open_file(&fptr,fname.c_str(),READONLY,&fstatus);
  if (fstatus != 0) {
    fits_report_error(stderr, fstatus);
    fits_unlock();
    throw IOExc("FITSTOIReader::openFile() fitsio error");
  }
  
  // Go to first extension, which should be a BINTABLE
  
  int simple, bitpix, naxis;
  long naxes;
  long pcount, gcount;
  int extend;
  fits_read_imghdr(fptr, 1, &simple, &bitpix, 
                   &naxis, &naxes, &pcount, &gcount, &extend, &fstatus);
  
  fits_movabs_hdu(fptr, 2, NULL, &fstatus);
  
  fits_get_num_cols(fptr,&ncols,&fstatus);
  fits_get_num_rows(fptr,&nrows,&fstatus);

  cout << "FITSTOIReader cols = " << ncols << " rows=" << nrows << endl;
  if (implicitSN) {
    firstSn = implicitSNStart;
  } else {
    int anyNul;
    double y;
    fits_read_col_dbl(fptr,1,1,1,1,0,&y,&anyNul,&fstatus);
    firstSn = (int) (y+.1);
  }

  // Ouverture fichier de flag separe  de LevelS (Reza 18/6/2002)
  if (sepFlagfile) {
    cout << " FITSTOIReader::openFile() - Opening separate flag file " << sepFlagFileName 
         << " NColFlags= " << sepFlagCols.size() << endl;
    fits_open_file(&fptrflg,sepFlagFileName.c_str(),READONLY,&fstatus);
    if (fstatus != 0) {
      fits_report_error(stderr, fstatus);
      fits_unlock();
      throw IOExc("FITSTOIReader::openFile() - sepFlagfile open fitsio error");
    }
    fits_movabs_hdu(fptrflg, 2, NULL, &fstatus);
    if (fstatus != 0) {
      fits_report_error(stderr, fstatus);
      fits_unlock();
      throw IOExc("FITSTOIReader::openFile() - sepFlagfile fits_movabs_hdu(2) fitsio error");
    }
    long nrowsflg = 0;
    fits_get_num_rows(fptrflg,&nrowsflg,&fstatus);
    if (nrows != nrowsflg) {
      cerr << " FITSTOIReader::openFile()/Error: Different NRows in flag and data files!" << endl;
      fits_unlock();
      throw ParmError("FITSTOIReader::openFile() Different NRows in flag and data files");
    }
    int ncolsflg = 0;
    fits_get_num_cols(fptrflg,&ncolsflg,&fstatus);    
    if (ncolsflg < sepFlagCols.size()) {
      cerr << " FITSTOIReader::openFile()/Error: ncolsflg (=" << ncolsflg
           << ") < sepFlagCols.size() (=" << sepFlagCols.size() << ")" << endl;
      fits_unlock();
      throw ParmError("FITSTOIReader::openFile() Different NRows in flag and data files");
    }
  }
  // 
  fits_unlock();  
}

void FITSTOIReader::init() {

  // Modif pour fichiers de flag separe de LevelS (Reza 18/6/2002)
  if (sepFlagfile && (allfn.size() > 1) ) {
    cerr << "FITSTOIReader::init()/Error- Multiple files and separate flag file not allowed !"<<endl;
    throw ParmError("FITSTOIReader::init() Multiple files and separate flag file not allowed");
  }

  openFile(allfn.front());

  fits_lock();
  // si pas implicitSN, la premiere colonne est le sampleNum. 
  // Sinon, le samplenum est la fitsline + offset.
  int itoi=-1;
  int col1 = implicitSN ? 0 : 1;
  for (int i=col1; i<ncols; i++) {
    char templt[10];
    sprintf(templt, "%d", i+1);
    char colname[200];
    int colnum;
    fits_get_colname(fptr, CASESEN, templt, colname, &colnum, &fstatus);
    cout << "FITSTOIReader col " << colname << endl;
    // On verifie si c'est une colonne de flag 
    if (!strncmp(colname, "fg_", 3) || !strncmp(colname, "FLAG_", 5) ) {
      colsinput[itoi].second=true;
    } else {
      declareOutput(colname);
      itoi++;
      colsinput[itoi] = pair<int,bool>(i,false);
    }
  }
  fits_unlock();

  snBegin = firstSn;
  if (forcedMinIn > 0 && forcedMinIn > snBegin) {
    snBegin = forcedMinIn;
  }
  openFile(allfn.back());
  snEnd = firstSn+nrows-1;
  if (forcedMaxIn > 0 && forcedMaxIn < snEnd) {
    snEnd = forcedMaxIn;
  }
  /* init */
  //  snMin=snBegin;     NNNNNOOOOOONNNNNN !!!!
  //  snMax=snEnd;
  cout << "FITSTOIReader range " << snBegin << " -> " << snEnd << endl;
}

long FITSTOIReader::calcMinOut() {
  chkinit();
  TOIManager* mgr = TOIManager::getManager();
  long firstReq = mgr->getRequestedBegin();
  return snBegin > firstReq ? snBegin : firstReq;
}

long FITSTOIReader::calcMaxOut() {
  chkinit();
  TOIManager* mgr = TOIManager::getManager();
  long lastReq = mgr->getRequestedEnd();
  return snEnd < lastReq ? snEnd : lastReq;
}


// ajout vf 31/07/2002
bool FITSTOIReader::checkSampleLimits(long& min, long& max, int pass)
{
  bool sample_ok=true;
  chkinit();
  return TOIProcessor::checkSampleLimits(min, max, pass);
  /*  
  cout << "check " << name << " in  " << min << " - " << max << " ; " << snMin << " - " << snMax << endl;

  // on verifie qu'on peut effectivement produire

  if (min < snBegin) {
    min = snBegin;
  }

  if (max > snEnd) {
    max = snEnd;
  }


  // si intersection des intervals vide on leve une exception
  //if (min >= snEnd || max <= snBegin) {
  //cerr << " FITSTOIReader::checkSampleLimits(long& min long& max)/Error Sample out of bounds" << endl;
  //throw ParmError("  FITSTOIReader::checkSampleLimits(long& min long& max) min >= snEnd || max <= snBegin");
  //}
  
  // ajustements si necessaires
  if (min < snMin) {
    snMin = min;
  }
  if (max > snMax) {
    snMax = max;
  }


  min=min<snMin?snMin:min;
  max=max>snMax?snMax:max;

  // cas sans contraintes, on retourne nos bornes
  if (min>max) {
    min = snBegin;
    max = snEnd;
  }

  cout << "check " << name << " out " << min << " - " << max << " ; " << snMin << " - " << snMax << endl;
  //cout << "fitstoirdr : limites verifiees : " << snBegin << " , " << snEnd << " : " << sample_ok << endl; 
  return sample_ok;
*/
}

void FITSTOIReader::calcSampleLimits(long& min, long& max)
{
  chkinit(); 
  cout << "calc  " << name << " in  " << min << " - " << max << " ; " << snBegin << " - " << snEnd << endl;
  
  min=min<snMin?snMin:min;
  max=max>snMax?snMax:max;
  snBegin=snMin;
  snEnd=snMax;
  cout << "calc  " << name << " out " << min << " - " << max << " ; " << snBegin << " - " << snEnd << endl;
  //cout << "fitstoirdr : limites temporaires calculees : " << snBegin << " , " << snEnd << endl;

}
// fin ajout vf



void FITSTOIReader::addFile(string fn) {
  allfn.push_back(fn);
}

void FITSTOIReader::setFlagFile(string fn, vector<FlagToiDef> flags)
{
  if (flags.size() < 1) {
    cerr << " FITSTOIReader::setFlagFile()/Error flag.size() = 0 ! " << endl;
    throw ParmError("FITSTOIReader::setFlagFile() flag.size() = 0");
  }

  sepFlagfile = true;
  sepFlagFileName = fn;
  sepFlagCols =  flags;
}

void FITSTOIReader::run() {
  for (vector<string>::iterator i=allfn.begin(); i!=allfn.end(); i++) {
    openFile(*i);
    if (Buff_Sz > 1) run2();  // Lecture bufferise
    else run1();              // Lecture un echantillon a la fois
  }
  fits_lock();
  if (fptr) {
    fits_close_file(fptr,&fstatus);
    fptr = NULL;
  }
  if (sepFlagfile && fptrflg) {
    fits_close_file(fptrflg,&fstatus);
    fptrflg = NULL;
  }
  fits_unlock();
}

// run 1 : deprecated. NON MAINTENU. Incompatible avec implicit SN.
//                     ^^^^^^^^^^^^  Reza , 13/5/2002 , Je viens de rajouter 
//                                   implicitSNStart et switvh run1/run2 
//                                   suivant buffersize
  
void FITSTOIReader::run1() {
  // Il faudrait optimiser en fonction de ce qui a ete demande comme samplenum,
  // mais cela implique de gerer aussi bien echant uniforme que non.
  // On pourrait aussi lire plusieurs elements d'un coup.
  int ncols = outIx.size();
  cout << "FITSTOIReader::run1() reader reading... NRows=" << nrows << " firstSn= " 
       << firstSn << endl;

  double* tabdata = new double[getNOut()]; 
  uint_8* tabflag = new uint_8[getNOut()]; 

  for (long i=0; i<nrows; i++) {  
    int anyNul;
    double y;
    fits_lock();
    long sn;
    if (implicitSN) {
      fits_read_col_dbl(fptr,1,i+1,1,1,0,&y,&anyNul,&fstatus);
      sn = (long) (y+.1); 
    }
    else sn = implicitSNStart+i;
    //fits_unlock();
    TOIManager* mgr = TOIManager::getManager();
    if (sn > mgr->getRequestedEnd()) {fits_unlock(); break;}
    if (sn < mgr->getRequestedBegin()) {fits_unlock(); continue;}
    // if (sn < mgr->getRequestedBegin()+10) cout << "rdr out " << sn << endl;
    int k;
    for (k=0; k<getNOut(); k++) {
      int j = colsinput[k].first;
      if ( !checkOutputTOIIndex(k) )  continue;  // Reza - Si TOI non connecte
      //fits_lock();
      fits_read_col_dbl(fptr,j+1,i+1,1,1,0,&y,&anyNul,&fstatus);
      tabdata[k] = y;
      long flg = 0;
      if (colsinput[k].second) {
        fits_read_col_lng(fptr,j+2,i+1,1,1,0,&flg,&anyNul,&fstatus);
      }

      if (sepFlagfile) { // Ajout Reza (18/6/2002) pour fichier de flags separe
        unsigned char sflg;
        flg = 0;
        for(int skf=0; skf<sepFlagCols.size(); skf++) {
          fits_read_col_byt(fptrflg, skf+1, i+1,1,1,0,&sflg,&anyNul,&fstatus);
          if (sflg) flg |= sepFlagCols[skf];
        }
      }  // Fin modif pour fichier de flags separe (18/6/2002)

      tabflag[k] = flg;
      //      fits_unlock();
      //putData(k, sn, y, flg);
    }
    fits_unlock();
    for (k=0; k<getNOut(); k++) {
      putData(k, sn, tabdata[k], tabflag[k]);
    }
    /*     if (i%50==0)  {
      // fits_unlock();
      sched_yield();
      //fits_lock();
      }*/
    totnscount++;
  }
  //fits_unlock();
  delete[] tabflag;
  delete[] tabdata;
  cout << "reader done reading... " << pthread_self() << endl;
}

void FITSTOIReader::run2()
// ---- Version bufferisee (Christophe 20/10/2001)
{
 int ncols = outIx.size();
 cout<<"FITSTOIReader::run2() reader reading... NRows="<<nrows<<" firstSn= " <<firstSn<<endl;

 //////// Prepare buffer, allocate memory
 double *samplenum = new double[Buff_Sz];
 vector<double*> colval;
 vector<long*> colflg;
 for(int k=0; k<getNOut(); k++) {
   colval.push_back(NULL);
   colflg.push_back(NULL);
   if(!checkOutputTOIIndex(k)) continue;
   colval[k] = new double[Buff_Sz];
   // modif || sepFlagfile pour fichier de flag separe (Reza 18/6/2002)
   if((colsinput[k].second) || sepFlagfile) colflg[k] = new long[Buff_Sz];
 }
 uint_8 * tmpflg = new uint_8[Buff_Sz];

 // Ajout Reza (18/6/2002) pour fichier de flags separe
 unsigned char* stmpflg = NULL; 
 if (sepFlagfile) stmpflg = new unsigned char[Buff_Sz];
 // Fin modif pour fichier de flags separe (18/6/2002)
 

 TOIManager* mgr = TOIManager::getManager();

 //////// Read data and put into pipe
 long ideb=-1,ifin=-1;
 for(long i=0; i<nrows; i+=Buff_Sz) {
   fits_lock();   // lock en debut de boucle de lecture - Reza 13/05/2002
   // faut-il lire dans le fichier fits ?
   if(i<ideb || i>ifin) { // Toujours vrai avec le += Buff_Sz
     ideb = i;
     ifin = ideb+Buff_Sz-1;
     if(ifin>=nrows) ifin=nrows-1;
     long n = ifin-ideb+1;
     int anyNul;
     if (implicitSN) {
       if (ideb+implicitSNStart > snEnd
           || (forcedMaxIn > 0 && ideb+implicitSNStart > forcedMaxIn)) {
         fits_unlock();  // unlock avant de casser la boucle - Reza 16/05/2002
         break;
       }
       if (ifin+implicitSNStart < snBegin
           || (forcedMinIn > 0 && ifin+implicitSNStart < forcedMinIn)) {
         fits_unlock();  // unlock avant de continuer - Reza 16/05/2002
         continue;
       }
       for (long j=0; j<Buff_Sz; j++) {
         samplenum[j] = j+ideb+implicitSNStart;
       }
       //---- Pas la peine, fait en debut de boucle - Reza 13/05/2002  fits_lock();
     } else {
       //---- Pas la peine, fait en debut de boucle - Reza 13/05/2002 fits_lock();
       fits_read_col_dbl(fptr,1,ideb+1,1,n,0,samplenum,&anyNul,&fstatus);
       if (samplenum[0] > snEnd
           || (forcedMaxIn > 0 && samplenum[0] > forcedMaxIn)) {
         fits_unlock();  // unlock avant de casser la boucle - Reza 13/05/2002
         break;
       }
       if (samplenum[n-1] < snBegin
           || (forcedMinIn > 0 && samplenum[n-1] < forcedMinIn)) {
         fits_unlock();
         continue;
       }
     }
     for(int k=0; k<getNOut(); k++) {
       if(colval[k]==NULL)  continue;
       int j = colsinput[k].first;
       fits_read_col_dbl(fptr,j+1,ideb+1,1,n,0,colval[k],&anyNul,&fstatus);
       if (fstatus !=0) {
         cerr << " RzBug - j+1" << j+1 << " ideb= " << ideb << endl;  
         fits_report_error(stderr,fstatus);
       }

       //--- Ajout Reza (18/6/2002) pour fichier de flags separe
       if (sepFlagfile) { 
         int sjj;
         for(sjj=0; sjj<n; sjj++) colflg[k][sjj] = 0;
         for(int skf=0; skf<sepFlagCols.size(); skf++) {
           fits_read_col_byt(fptrflg, skf+1, ideb+1,1,n,0,stmpflg,&anyNul,&fstatus);
           if (fstatus !=0) {
             cerr << " RzBug2 - skf=" << skf << " ideb= " << ideb << endl;  
             fits_report_error(stderr,fstatus);
           }
           for(sjj=0; sjj<n; sjj++) {
             if (stmpflg[sjj]) colflg[k][sjj] |= sepFlagCols[skf];
           }
         }
       }  
       else if(colflg[k]!=NULL)         
         fits_read_col_lng(fptr,j+2,ideb+1,1,n,0,colflg[k],&anyNul,&fstatus);
       //--- Fin modif pour fichier de flags separe (18/6/2002)
     }
     if(fstatus!=0) {
       fits_report_error(stderr,fstatus);
       fits_unlock(); // On supprime le lock avant le throw - Reza 13/05/2002
       throw IOExc("FITSTOIReader::run2: Error Reading Fits file\n");
     }
   }
   fits_unlock();   // unlock en fin de boucle de lecture - Reza 13/05/2002

   long ip = i-ideb; // pointeurs dans les buffers
   long sn = (long) (samplenum[ip]+.1);
   if (ip == 0) { 
     for(int k=0; k<getNOut(); k++) {
       if(colval[k]==NULL)  continue;
       if (colflg[k] != NULL) {
         for (int ii=0; ii<Buff_Sz; ii++) {
           tmpflg[ii] = colflg[k][ii];
         }
         putData(k, sn, ifin-ideb+1, colval[k], tmpflg);
       } else {
         putData(k, sn, ifin-ideb+1, colval[k]);
       }
     }
   }
   totnscount += (ifin-ideb+1);
 }

 //////// des-allocate buffers
 delete [] samplenum;
 delete [] tmpflg;
 if (sepFlagfile) delete [] stmpflg;
 {for(int k=0; k<getNOut(); k++) {
   if(colval[k]!=NULL) delete [] colval[k];
   if(colflg[k]!=NULL) delete [] colflg[k];
 }}
 colval.resize(0); colflg.resize(0);

 cout << "reader (buffered) done reading... " << pthread_self() << endl;

}


// affichage des limites
void FITSTOIReader::printLimits()
{
  cout << "fitstoirdr " << name <<" : limites calculees : " << snBegin << " , " << snEnd << endl;
}