source: Sophya/trunk/AddOn/TAcq/brviscalc.cc@ 3876

//----------------------------------------------------------------
// Projet BAORadio - (C) LAL/IRFU  2008-2010
// Classes de threads de calcul de visibilites pour BAORadio 
//----------------------------------------------------------------

#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <fstream>
#include <signal.h>

#include "pexceptions.h"
#include "fioarr.h"
#include "matharr.h"
#include "timestamp.h"
#include "ctimer.h"

#include "brviscalc.h"


//---------------------------------------------------------------------
// Classe de traitement - calcul de visibilite pour n fibres 
//---------------------------------------------------------------------

/* --Methode-- */
BRVisibilityCalculator::BRVisibilityCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, size_t nthr)
  : BRBaseProcessor(memgr), paralex_(*this, nthr), nparthr_(nthr), 
    outpath_(outpath), nmean_(nmean), nbcalc_(1), calcid_(0), vpdata_(2*memgr.NbFibres())
    // , dtfos_(outpath+"visdt.fits", Fits_Create), visdt_(dtfos_, 1024, true);
{  
  DefineRank(1,0);

  uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();
  chanum_.SetSize(maxnpairs);
  sa_size_t k=0;
  for(size_t i=0; i<2*memgr_.NbFibres(); i++)  vpdata_[i]=NULL;
  for(size_t i=0; i<2*memgr_.NbFibres(); i++) {
    for(size_t j=i; j<2*memgr_.NbFibres(); j++) {
      chanum_(k) = (i+1)*100+(j+1);  k++;
    }
  }
  SelectPairs();

  //  visdt_.AddFloatColumn("mfc");
  visdt_.AddDoubleColumn("mfc");
  visdt_.AddDoubleColumn("mtt");
  visdt_.AddIntegerColumn("jfreq");
  visdt_.AddIntegerColumn("numch");
  visdt_.AddFloatColumn("vre");
  visdt_.AddFloatColumn("vim");

  if (nmean_ < 1)  nmean_=memgr_.NbPaquets();
  if (nmean_ < 1)  nmean_=1;

  cout << " BRVisibilityCalculator::/Info  nmean=" << nmean_ << endl;

  totnbpaq_=0;
  numfile_=0;
  nb_flop_=0.;
  moyfc_=moytt_=0.;

  fgallfibok=NULL;
  fgcktt_=false;
  setNameId("viscalc", 0); 
}

/* --Methode-- */
BRVisibilityCalculator::~BRVisibilityCalculator()
{
  if (totnbpaq_<1)   return; 
  cout << " BRVisibilityCalculator - Visibility Datatable : " << endl;
  cout << visdt_;
  
  string filename;
  filename = outpath_+"visdt.ppf";
  if (nbcalc_>1) {
    char sbuff[32];
    sprintf(sbuff,"visdt_%d",(int)calcid_);
    filename = outpath_+sbuff;
  }
  POutPersist po(filename);
  po << visdt_;
  if (calcid_ == 0) {
    POutPersist poc(outpath_+"chanum.ppf");
    poc << chanum_;
    
    if (fgcktt_) {
      cout << " BRVisibilityCalculator -  Check TimeTag done: TotNPaqProc= " << totnbpaq_ << endl;
      for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)  {
        cout << " BRTTCheck-Fiber[" << fib << "] NBadTT=" << vbadtt_[fib] << " NDiffTT>5="
             << vndiff5tt_[fib] << "  NotSameTT=" << vnsamett_[fib] << endl; 
      }
      POutPersist pott(outpath_+"ttfcmtx.ppf");
      pott << PPFNameTag("FC") << fcmtx_; 
      pott << PPFNameTag("TT") << ttmtx_;     
    }
  }
}

/* --Methode-- */
void BRVisibilityCalculator::DefineRank(uint_4 nbc, uint_4 cid)
{
  if ((nbc>6)||(cid>=nbc)) 
    throw ParmError("BRVisibilityCalculator::DefineRank() NbCalc > 6 !");
  nbcalc_=nbc;
  calcid_=cid;
  if (nbcalc_>1) {
    uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();
    uint_4 npairs=maxnpairs/nbcalc_;
    if (calcid_==(nbcalc_-1)) 
      SelectPairs(calcid_*npairs, maxnpairs-calcid_*npairs);
    else 
      SelectPairs(calcid_*npairs, npairs);
    MemZaction mmzas[6]={MemZA_ProcA,MemZA_ProcB,MemZA_ProcC,MemZA_ProcD,MemZA_ProcE,MemZA_ProcF};
    SetMemZAction(mmzas[calcid_]);
    setNameId("viscalc_grp", calcid_); 
  }
  return ;
}

/* --Methode-- */
uint_4 BRVisibilityCalculator::SelectPairs(uint_4 pair1, uint_4 nbpairs)
{
  BRPaquet paq(memgr_.PaqSize());
  uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();

  if (pair1 >= maxnpairs)  pair1=maxnpairs-1;
  if (nbpairs > maxnpairs-pair1)  nbpairs=maxnpairs-pair1;
  pairst_=pair1;
  nbpairs_=nbpairs;
  vismtx_.SetSize(nbpairs_, paq.DataSize()/4); 
  return nbpairs_;
}

/* --Methode-- */
int BRVisibilityCalculator::SelectFreqBinning(uint_4 freq1, uint_4 freq2, uint_4 nbfreq)
{
  jf1_=freq1; jf2_=freq2;
  if ((jf1_<1)||(jf1_>=vismtx_.NCols()))  jf1_=1;
  if ((jf2_<1)||(jf2_>=vismtx_.NCols())||(jf2_<jf1_))  jf2_=vismtx_.NCols()-1;
  if (nbfreq<1) nbfreq=(jf2_-jf1_);
  djf_=(jf2_-jf1_)/nbfreq;
  if (djf_<1) djf_=1;
  cout << " BRVisibilityCalculator::SelectFreqBinning/Info  JF1=" << jf1_ 
       << " JF2=" << jf2_ << " DJF=" << djf_ << endl;

}


/* --Methode-- */
int BRVisibilityCalculator::ActivateTimeTagCheck(uint_8 maxnpaq)
{
  mindeltatt_=memgr_.PaqSize()/2;
  if (mindeltatt_<1) mindeltatt_=1;
  fcmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
  ttmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
  vlasttt_.resize(memgr_.NbFibres(), 0);
  vbadtt_.resize(memgr_.NbFibres(), 0);
  vnsamett_.resize(memgr_.NbFibres(), 0);
  vndiff5tt_.resize(memgr_.NbFibres(), 0);

  fgcktt_=true;
  cout << " BRVisibilityCalculator::ActivateTimeTagCheck() - TT/Fc matrix NCols=" << maxnpaq 
       << " MinDeltaTT=" << mindeltatt_ << endl;

  return 0;
}


/* --Methode-- */
void BRVisibilityCalculator::run()
{
  cout << " BRVisibilityCalculator[" << calcid_ << "/" << nbcalc_ 
       << "]::run() - Starting " << " NFibers=" << memgr_.NbFibres() 
       << " NChan=" << 2*memgr_.NbFibres() << " NPairs=" << nbpairs_ << " First:" << pairst_ << endl;   
  
  if (nparthr_ < 2)  return BRBaseProcessor::run();
  // Execution multithread parallele
  setRC(1);     
  int rc=0;
  try {
    if ((nmean_%memgr_.NbPaquets())!=0) {
      uint_4 mnmean = (nmean_/memgr_.NbPaquets()+1)*memgr_.NbPaquets();
      cout << " BRVisibilityCalculator::run()/Info changing nmean=" << nmean_  << " to multiple of"
           << " memgr_.NbPaquets() -> " << mnmean << endl;
      nmean_=mnmean;
    }
    paralex_.SetParallelTask(*this);
    cout << " BRVisibilityCalculator::run()/Info : starting ParallelExecutor with nThreads=" 
         << paralex_.nThreads() << " ... " << endl;
    paralex_.start();

    fgallfibok = new bool[memgr_.NbPaquets()];

    size_t paqsz=memgr_.PaqSize();
    bool fgrun=true;
    while (fgrun) {
      if (stop_) break;
      if (memgr_.GetRunState() == MemZR_Stopped) break;
      int mid = memgr_.FindMemZoneId(mmact_);  // (MemZA_ProcA);
      //      Byte* buffg = memgr_.GetMemZone(mid);
      //      if (buffg == NULL) {
      if (mid < 0) {
        cout << "BRVisibilityCalculator[" << calcid_ << "]::run()/ERROR FindMemZoneId(" 
             << (int)mmact_ << ") ->" << mid << ") -> NULL" << endl;
        setRC(7);      fgrun=false;             
        break;  
      }
      for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
        fbuff_[fib] = memgr_.GetMemZone(mid,fib);
        if (fbuff_[fib] == NULL) { // cela ne devrait pas arriver 
          cout << "BRVisibilityCalculator[" << calcid_ << "]::run()/ERROR memgr.GetMemZone(" << mid << "," << fib << ") -> NULL" << endl;
          setRC(9);       fgrun=false;          
          break;
        }
      }

      if (totnbpaq_%nmean_ == 0) {
        if (totnbpaq_ > 0) {
          moyfc_/=nmean_;
          moytt_/=nmean_;
          vismtx_.Info()["MeanFC"] = moyfc_;
          vismtx_.Info()["MeanTT"] = moytt_;
          vismtx_.Info()["NPAQSUM"] = nmean_;
          
          // ATTENTION : Matrice visibilites non moyennee 
          char nfile[48];
          if (nbcalc_==1) 
            sprintf(nfile,"vismtx%d.ppf",numfile_);
          else 
            sprintf(nfile,"vismtx_%d_%d.ppf",(int)calcid_,numfile_);
          string flnm=outpath_+nfile;
          POutPersist po(flnm);
          po << vismtx_;
          cout << numfile_ << "-BRVisibilityCalculator[" << calcid_ << "]::run() NPaqProc="
               << totnbpaq_ << "  -> Visibility Matrix in " << flnm << endl;
          FillVisibTable(moyfc_, moytt_);
          numfile_++;
        }
        vismtx_ = complex<r_4>((r_4)0.,(r_4)0.);
        moyfc_=moytt_=0.;
      }

      for(size_t jp=0; jp<memgr_.NbPaquets(); jp++) {   // boucle sur les paquets d'une zone  
        fgallfibok[jp]=fgokallfibers_=true;
        for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
          vpaq_[fib].Set(fbuff_[fib]+jp*paqsz);
          vfgok_[fib] = vpchk_[fib].Check(vpaq_[fib],curfc_[fib]);
          if (!vfgok_[fib])  fgallfibok[jp]=fgokallfibers_=false;
        }
        if (fgokallfibers_)  {
          if (totprocnpaq_==0) {
            for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)  {
              fcfirst_[fib]=curfc_[fib];
              ttfirst_[fib]=vpaq_[fib].TimeTag();
            }
          }
          totprocnpaq_++;
          moyfc_ += curfc_[0];
          moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
          if ((fgcktt_)&&(calcid_==0))  CheckTimeTag();
          totnbpaq_++;
        }
      }  // Fin de boucle sur les paquets 
      
      // Execution parallele  pour calcul des visibilites par bandes de frequence
      int rcpex=paralex_.execute();
      if (rcpex!=0)  cout << " BRVisibilityCalculator[" << calcid_ << "]::run() / Error Rc[paralex_.execute()]=" << rcpex << endl;

      memgr_.FreeMemZone(mid, mmsta_);  // (MemZS_ProcA);
    } // Fin de boucle sur les zones a traiter 
    //------------------------------------
    cout << " --------- END BRVisibilityCalculator[" << calcid_ << "]::run() , TotNbProcPaq=" << totprocnpaq_ << endl;
    /*
    for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)  vpchk_[fib].Print();
    cout << " ------------------------------------ " << endl;
    */
    delete[] fgallfibok;
  }
  catch (std::exception& exc) {
    cout << " BRVisibilityCalculator[" << calcid_ << "]::run()/catched std::exception " << exc.what() << endl;
    setRC(98);  
    return; 
  }
  catch(...) {
    cout << " BRVisibilityCalculator[" << calcid_ << "]::run()/catched unknown ... exception " << endl;
    setRC(99);  
    return; 
  }
  
}

/* --Methode-- */
int BRVisibilityCalculator::Process()
{
    
  for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
    vpdata_[2*fib] = vpaq_[fib].Data1C();
    vpdata_[2*fib+1] = vpaq_[fib].Data2C();
  }
  
  if (totnbpaq_%nmean_ == 0) {
    if (totnbpaq_ > 0) {
      moyfc_/=nmean_;
      moytt_/=nmean_;
      vismtx_.Info()["MeanFC"] = moyfc_;
      vismtx_.Info()["MeanTT"] = moytt_;
      vismtx_.Info()["NPAQSUM"] = nmean_;

      // ATTENTION : Matrice visibilites non moyennee 
      char nfile[48];
      if (nbcalc_==1) 
        sprintf(nfile,"vismtx%d.ppf",numfile_);
      else 
        sprintf(nfile,"vismtx_%d_%d.ppf",(int)calcid_,numfile_);
      string flnm=outpath_+nfile;
      POutPersist po(flnm);
      po << vismtx_;
      cout << numfile_ << "-BRVisibilityCalculator::Process() NPaqProc="
           << totnbpaq_ << "  -> Visibility Matrix in " << flnm << endl;
      FillVisibTable(moyfc_, moytt_);
      numfile_++;
    }
    vismtx_ = complex<r_4>((r_4)0.,(r_4)0.);
    moyfc_=moytt_=0.;
  }

  sa_size_t k=0;
  for(size_t i=0; i<vpdata_.size(); i++) {
    for(size_t j=i; j<vpdata_.size(); j++) {
      size_t kpair=i*vpdata_.size()+j;
      if (kpair<pairst_)  continue;
      if (kpair>=(pairst_+nbpairs_))  break;
      TVector< complex<r_4> > vis = vismtx_.Row(k);   k++;
      for(sa_size_t f=1; f<vis.Size(); f++) {
        vis(f) += complex<r_4>((r_4)vpdata_[i][f].realB(), (r_4)vpdata_[i][f].imagB()) * 
          complex<r_4>((r_4)vpdata_[j][f].realB(), -(r_4)vpdata_[j][f].imagB());
      }
      nb_flop_ += (7.*(r_8)vis.Size());
    }
  }

  moyfc_ += curfc_[0];
  moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
  if ((fgcktt_)&&(calcid_==0))  CheckTimeTag();
  totnbpaq_++;
  return 0;
}

/* --Methode-- */
int BRVisibilityCalculator::execute(int tid)
{
  vector<TwoByteComplex*>  pvpdata(2*memgr_.NbFibres());
  size_t paqsz=memgr_.PaqSize();
  BRPaquet ppaq(paqsz);

  sa_size_t fdelt = vismtx_.NCols()/nparthr_;
  sa_size_t fdeb = tid*fdelt;
  sa_size_t ffin = (tid+1)*fdelt;

  if (fdeb<1) fdeb=1;
  if ((ffin>vismtx_.NCols())||(tid==(nparthr_-1))) ffin=vismtx_.NCols();

  for(size_t jp=0; jp<memgr_.NbPaquets(); jp++) {   // boucle sur les paquets d'une zone  
    if (!fgallfibok[jp])  continue;
    for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
      ppaq.Set(fbuff_[fib]+jp*paqsz);
      pvpdata[2*fib] = ppaq.Data1C();
      pvpdata[2*fib+1] = ppaq.Data2C();
    }
    sa_size_t k=0;
    for(size_t i=0; i<vpdata_.size(); i++) {
      for(size_t j=i; j<vpdata_.size(); j++) {
        size_t kpair=i*vpdata_.size()+j;
        if (kpair<pairst_)  continue;
        if (kpair>=(pairst_+nbpairs_))  break;
        TVector< complex<r_4> > vis = vismtx_.Row(k);   k++;
        for(sa_size_t f=fdeb; f<ffin; f++) {
          vis(f) += complex<r_4>((r_4)pvpdata[i][f].realB(), (r_4)pvpdata[i][f].imagB()) * 
            complex<r_4>((r_4)pvpdata[j][f].realB(), -(r_4)pvpdata[j][f].imagB());
        }
        nb_flop_ += (7.*(r_8)(ffin-fdeb));
      }
    }

  }  // Fin de boucle sur les paquets 
  
  return 0;
}

/* --Methode-- */
int BRVisibilityCalculator::FillVisibTable(double fcm, double ttm)
{
  double xnt[10];
  xnt[0]=fcm;  xnt[1]=ttm/1.25e8;

  if (djf_<2) {
    for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
      for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
        xnt[2]=jf;
        xnt[3]=chanum_(rv+pairst_);
        xnt[4]=vismtx_(rv,jf).real()/(r_4)(nmean_);
        xnt[5]=vismtx_(rv,jf).imag()/(r_4)(nmean_);
        visdt_.AddRow(xnt);
      }
    }
  }
  else {
    for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
      for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
        r_4 moyreal=0.;
        r_4 moyimag=0.;
        sa_size_t jjfmx=jf+djf_;
        if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
        for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
          moyreal+=vismtx_(rv,jjf).real();
          moyimag+=vismtx_(rv,jjf).imag();
        }
        xnt[2]=jf+djf_/2;
        xnt[3]=chanum_(rv+pairst_);
        xnt[4]=moyreal/(r_4)(nmean_*djf_);
        xnt[5]=moyimag/(r_4)(nmean_*djf_);
        visdt_.AddRow(xnt);
      }
    }
  }
  return 0;
}

/* --Methode-- */
int BRVisibilityCalculator::CheckTimeTag()
{
  if (totnbpaq_==0) {
    for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)  {
      vlasttt_[fib]=ttfirst_[fib];
      if (ttmtx_.NCols()>0) {
        fcmtx_(fib,totnbpaq_) = curfc_[fib]; 
        ttmtx_(fib,totnbpaq_) = vlasttt_[fib]; 
      }
    }
    return 0;
  }
  for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
    int_8 ld = (int_8)vpaq_[fib].TimeTag()-(int_8)vlasttt_[fib];
    int_8 fd = (int_8)vpaq_[fib].TimeTag()-(int_8)ttfirst_[fib]-(int_8)vpaq_[0].TimeTag()+(int_8)ttfirst_[0];
    /*    if ( (ld < mindeltatt_) || (fd<-5) || (fd>5)) { vbadtt_[fib]++; vnsamett_[fib]++; }
    else {
      if (fd!=0)  vnsamett_[fib]++;
    }
    */
    if  (ld < mindeltatt_) vbadtt_[fib]++;
    else {
      if (fd != 0)  vnsamett_[fib]++;
      if ((fd<-5)||(fd>5))  vndiff5tt_[fib]++;
    }
    vlasttt_[fib]=vpaq_[fib].TimeTag();
    if (totnbpaq_<ttmtx_.NCols()) {
      fcmtx_(fib,totnbpaq_) = curfc_[fib]; 
      ttmtx_(fib,totnbpaq_) = vlasttt_[fib]; 
    }
  }
  return 0;
}

//-------------------------------------------------------------------------------
// Classe Groupe (ensemble) de Calculateur de Visibilites, tournant en parallele
//-------------------------------------------------------------------------------

/* --Methode-- */
BRVisCalcGroup::BRVisCalcGroup(size_t nbcalc, RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, size_t nthr)
  : tm_(false)
{
  if ((nbcalc<1)||(nbcalc>6)) 
    throw ParmError("BRVisCalcGroup::BRVisCalcGroup NbCalc > 6 !");
  for(size_t i=0; i<nbcalc; i++) {
    BRVisibilityCalculator * viscp=new BRVisibilityCalculator(memgr, outpath, nmean, nthr);
    viscp->DefineRank(nbcalc, i);
    viscalcp_.push_back(viscp); 
  }
}
/* --Methode-- */
BRVisCalcGroup::~BRVisCalcGroup()
{
  for(size_t i=0; i<viscalcp_.size(); i++) 
    delete viscalcp_[i];
}
/* --Methode-- */
int BRVisCalcGroup::SelectFreqBinning(uint_4 freq1, uint_4 freq2, uint_4 nbfreq)
{
  int rc=0;
  for(size_t i=0; i<viscalcp_.size(); i++) 
    rc=viscalcp_[i]->SelectFreqBinning(freq1, freq2, nbfreq);
  return rc;
}
/* --Methode-- */
void BRVisCalcGroup::start()
{
  for(size_t i=0; i<viscalcp_.size(); i++) 
    viscalcp_[i]->start();
  tm_.SplitQ();
}
/* --Methode-- */
void BRVisCalcGroup::join()
{
  r_8 totflop=0.;
  for(size_t i=0; i<viscalcp_.size(); i++) {
    viscalcp_[i]->join();
    totflop += viscalcp_[i]->TotNbFLOP();
  }
  tm_.SplitQ();
  cout << "      ----------------------------------------------------------" << endl;
  cout << " BRVisCalcGroup::join() : Finished , Elaspsed time: " << tm_.PartialElapsedTimems() 
       << " ms (total:" << tm_.TotalElapsedTimems() << ")" << endl;
  cout << " ... TotalMegaFLOP= " << totflop/(1024.e3) << " @ " 
       << totflop/(r_8)tm_.PartialElapsedTimems()/(1024) << " MFLOP/s" << endl; 
  cout << "      ----------------------------------------------------------" << endl;
  return;
}