// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

/* ---------------------------------------------------------- 
   Projet BAORadio/PAON4 - (C) LAL/IRFU  2008-2015 

   msvis2dt: programme de lecture des fichiers meanspecII.ppf 
   (sortie de monitoring), et vismtxII.ppf (visibilities) 
   Configuration de debug avec une carte avec firmware RAW et 
   une avec firmware FFT - Septembre 2015 
   R. Ansari, J.E. Campagne, C. Magneville   -  LAL/Irfu
   ---------------------------------------------------------- */

// include standard c/c++
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <iostream>
#include <string>

#include "pexceptions.h"
#include "tvector.h"
#include "fioarr.h"
// #include "tarrinit.h"
#include "ntuple.h" 
#include "datatable.h" 
#include "histinit.h" 
#include "matharr.h" 
#include "timestamp.h"
#include <utilarr.h>

// include sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"

// include lecteur de fichiers visibilites 
// #include "visip4reader.h"

//--------------------------- Fonctions de ce fichier   ------------------- 
int Usage(bool fgshort=true);
// int DecodeArgs(int narg, char* arg[]);

/* --Fonction-- */
int Usage(bool fgea)
{
  cout << " --- msvis2dt.cc : Read PPF files produced by mfacq (meanspecII.ppf vismtxII.ppf) \n" << endl;
  if (fgea) cout << " msvis2dt: Missing or wrong argument ! " << endl;
  cout << " Usage: msvis2dt [-tfm Nf] InPath_vismtx InPath_meanspec Imin,Imax OutPPFFile Freq1 Freq2 DeltaFreq [PrtLev=0] \n"
       << "  -tfm Nf : create time-frequncy maps of auto-correlation binning Nf frequencies \n" 
       << "  InPath_vismtx: path for input vismtxII.ppf (on bao5) \n"
       << "  InPath_meanspec: path for input meanspecII.ppf (on bao6) \n"
       << "  Imin,Imax : read files for Imin<=iii<=Imax  iii+=step \n"
       << "  OutPPFFile: Output PPF file name \n"
       << "  Freq1 Freq2: the two frequencies for which cross-correlations are computed [=1375 1410 MHz] \n"
       << "  DeltaFreq: frequency band (in MHz) [=1 MHz] \n"
       << "    Note: bande_freq=[freq1,2-deltafreq , freq1,2+deltafreq] \n"<< endl;
  /*
  if (fgshort) {
    cout << " rdvisip4 -h for detailed instructions" << endl;
    return 1;
  }
  */
  return 1;
}

//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
  if (narg<9) return Usage(true);
  int offa=0;
  bool FgTFM=false;   // true -> create time-frequency maps of cross-correlations 
  sa_size_t TFMfbin=8;
  if (strcmp(arg[1],"-tfm")==0) {
    if (narg<11) return Usage(true);
    FgTFM=true;
    TFMfbin=atoi(arg[2]);
    if (TFMfbin<1) TFMfbin=8;
    offa=2;
  }
  string path5 = arg[offa+1];
  string path6 = arg[offa+2];
  int Imin,Imax,Istep=1; 
  sscanf(arg[offa+3],"%d,%d",&Imin,&Imax);
  Istep=1;
  string outfile=arg[offa+4];
  // frequency range definition 
  double deltafreq=1.;
  double freq1=1375.;
  double freq2=1410.;
  int nargo=narg-offa;
  if ((nargo>5)&&(strcmp(arg[offa+5],"-")!=0))  freq1=atof(arg[offa+5]);
  if ((nargo>6)&&(strcmp(arg[offa+6],"-")!=0))  freq2=atof(arg[offa+6]);
  if ((nargo>7)&&(strcmp(arg[offa+7],"-")!=0))  deltafreq=atof(arg[offa+7]);
  // time averaging interval definition, in term of visibility files
  int deltaIavg=1;
  //  if ((nargo>8)&&(strcmp(arg[offa+8],"-")!=0))  deltaIavg=atoi(arg[offa+8]);
  if (deltaIavg<1)  deltaIavg=1;
  int prtlev=1;
  if (nargo>8) prtlev=atoi(arg[offa+8]);

  cout << " msvis2dt/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"  OutFile:"<<outfile<<"\n"
       << " freq1="<<freq1<<" freq2="<<freq2<<" deltafreq="<<deltafreq<<" MHz\n"
    //       << JFmin<<" <=NumFreq<= "<<JFmax<<"  "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
       <<"  DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;

  string TFMoutfile;
  if (FgTFM) {
    size_t l = outfile.length();
    size_t pp = outfile.rfind('.');
    if (pp<l) TFMoutfile = outfile.substr(0,pp) + "_tfm.ppf";
    cout << " msvis2dt/Info: Time-Frequency maps of Xcor will be created and saved to file "
	 <<TFMoutfile<<endl;
  }

  //---- calcul des index de bandes de frequences 
  // pas en frequence du firmware FFT 
  double deltanufft=250./4096;    // 250 MHz en 4096 frequences 
  sa_size_t DJF=deltafreq/deltanufft;
  double freqstart=1250.;   // Bande de 1250-1500 MHz 
  sa_size_t JFmin=(freq1-1250.-deltafreq)/deltanufft;
  sa_size_t JFmax=(freq1-1250.+deltafreq)/deltanufft;
  sa_size_t JFminB=(freq2-1250.-deltafreq)/deltanufft;
  sa_size_t JFmaxB=(freq2-1250.+deltafreq)/deltanufft;
  cout << " msvis2dt/Info:   frequency bands ... " <<endl;
  cout << " Band 1: "<<freq1<<" +/-"<<deltafreq<<" MHz -> "<<JFmin<<" <=NumFreq<= "<<JFmax<<endl;
  cout << " Band 2: "<<freq2<<" +/-"<<deltafreq<<" MHz -> "<<JFminB<<" <=NumFreq<= "<<JFmaxB<<endl;

  // Bande de frequence de 2 MHz autour du 21 cm (1420 MHz -> JF=2792) 
  sa_size_t JFmin21=2776,JFmax21=2808;   // +/- 1 MHz autour de 1420 MHz 
  sa_size_t JFmin21_5=2710,JFmax21_5=2874;  //  +/- 5 MHz autour de 1420 MHz 
  cout << " pjHI +/-1 MHz @1420 MHz " << JFmin21<<" <=NumFreq<= "<<JFmax21
       << " pjHI +/-5 MHz @1420 MHz " << JFmin21_5<<" <=NumFreq<= "<<JFmax21_5 << endl;
 
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    ResourceUsage resu;

    //  Numero des lignes des auto-correlations dans la matrice des visibilites 
    sa_size_t KVAC[4]={0,4,7,9}; 
    //  Numero des lignes des cross-correlations 1H-2H 1H-3H 1H-4H 2H-3H 2H-4H 3H-4H dans la matrice des visibilites 
    sa_size_t KVCXHH[6]={1,2,3,5,6,8}; 

    const char* dtnames[43]={"k","day","time",
			     "p1h","p2h","p3h","p4h",
			     "p1hB","p2hB","p3hB","p4hB",
			     "p1hHI","p2hHI","p3hHI","p4hHI",
			     "p1hHI5","p2hHI5","p3hHI5","p4hHI5",
                             "xh1h2","xh1h3","xh1h4","xh2h3","xh2h4","xh3h4",
                             "xh1h2B","xh1h3B","xh1h4B","xh2h3B","xh2h4B","xh3h4B",
			     "xh1h2HI","xh1h3HI","xh1h4HI","xh2h3HI","xh2h4HI","xh3h4HI",
			     "dayraw","timeraw",
			     "p3hraw","p3hrawB","p3hrawHI","p3hrawHI5"};

    DataTable dt(64);
    dt.AddIntegerColumn(dtnames[0]);
    dt.AddIntegerColumn(dtnames[1]);
    dt.AddFloatColumn(dtnames[2]);
    for(int kk=0; kk<16; kk++) 
      dt.AddFloatColumn(dtnames[3+kk]);
    for(int kk=0; kk<18; kk++) 
      dt.AddComplexColumn(dtnames[3+16+kk]);

    dt.AddIntegerColumn(dtnames[3+16+18]);    
    dt.AddFloatColumn(dtnames[3+16+18+1]);
    for(int kk=0; kk<4; kk++) 
      dt.AddFloatColumn(dtnames[3+16+18+2+kk]);

    Range freqs(JFmin,JFmax);
    Range freqs21(JFmin21,JFmax21);
    Range freqs21_5(JFmin21_5,JFmax21_5);
    Range freqsB(JFminB,JFmaxB);

    
    double acrdt[4];  // 4 bandes x 1 autocor raw/meanspec 
    double acdt[16];   // 4 bandes x [ 4 autocor vismtx ]
    complex<double> cxdt[18];   // les 3*6=18 valeurs de cross-correlations pour remplissage dans la table 

    TimeStamp dateorg(2015,1,1,12,0,0.);  // Date origine 1 jan 2015
    double mttag;
    int cnt=0;
    //----- 5 TimeFrequency maps of auto-correlations 
    vector< TArray< r_4 > > vtfm;
    
    // Getting en empty row_ptr
    DataTableRowPtr rowp = dt.EmptyRowPtr();

    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

    char filename[1024];
    for(int jf=Imin; jf<=Imax; jf+=Istep) {
      sprintf(filename,"%s/vismtx%d.ppf",path5.c_str(),jf);
      if (prtlev>1)  cout << "JFile="<<jf<<" -> Opening file " << filename << endl;
      PInPersist pin5(filename);
      TMatrix< complex<float> >  vismtx;
      time_t tm5 = pin5.CreationDate();
      struct tm * ptms5 = gmtime(&tm5);
      TimeStamp tms5(ptms5->tm_year, ptms5->tm_mon, ptms5->tm_mday, ptms5->tm_hour, ptms5->tm_min, ptms5->tm_sec);
      pin5 >> vismtx;

      sprintf(filename,"%s/meanspec%d.ppf",path6.c_str(),jf);
      if (prtlev>1)  cout << "JFile="<<jf<<" -> Opening file " << filename << endl;
      PInPersist pin6(filename);
      TMatrix<float> spectre;
      time_t tm6 = pin6.CreationDate();
      struct tm * ptms6 = gmtime(&tm6);
      TimeStamp tms6(ptms6->tm_year, ptms6->tm_mon, ptms6->tm_mday, ptms6->tm_hour, ptms6->tm_min, ptms6->tm_sec);
      pin6 >> PPFNameTag("spectre") >> spectre;

      if (vismtx.NCols() != spectre.NCols())  
	throw PError("msvis2dt/ERROR: vismtx and spectre do NOT have same number of frequencies");
      if (jf == Imin) { 
	if (FgTFM) {  //allocating time-frequency maps 
	  tfmSX=(Imax-Imin+1)/Istep/deltaIavg;
	  tfmSY=vismtx.NCols()/TFMfbin;
	  cout << " msvis2dt/Info: allocating Time-Freqncy maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
	  for(int k=0; k<5; k++) vtfm.push_back( TArray< r_4 >(tfmSX, tfmSY) ); 
	}
      }
      //---------------------------------------------------------------------------------
      // pour les 4 autocorrelations provenant des visibilites  
      for(int k=0; k<4; k++)  {
	TVector< complex<float> > acv = vismtx.Row( KVAC[k] );
	acdt[k] = acv(freqs).Sum().real();       // integration en bande 1
	acdt[k+4] = acv(freqsB).Sum().real();      // integration en bande large B
	acdt[k+8] = acv(freqs21).Sum().real();     // integration a +/- 1 MHz @ 1420 MHz 
	acdt[k+12] = acv(freqs21_5).Sum().real();   // integration a +/- 5 MHz @ 1420 MHz 

	if (FgTFM && (TFMtmidx<tfmSX)) {  //filling time-frequency maps 
	  TArray< r_4 > & tfmap = vtfm[k];
	  for(sa_size_t jy=0; jy<tfmSY; jy++) {
	    complex<r_4> zz = acv( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	    tfmap(TFMtmidx, jy) = zz.real();
	  } 
	}  //----- end of time-frequency maps filling 
      }
      {
      // Voie 3H-raw 
	TVector<float> acr = spectre.Row(0);
	acrdt[0] = acr(freqs).Sum();    // integration en bande 1
	acrdt[1] = acr(freqsB).Sum();    // integration en bande large B
	acrdt[2] = acr(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
	acrdt[3] = acr(freqs21_5).Sum();   // integration a +/- 5 MHz @ 1420 MHz 
	if (FgTFM && (TFMtmidx<tfmSX)) {  //filling time-frequency maps 
	  TArray< r_4 > & tfmap = vtfm[4];
	  for(sa_size_t jy=0; jy<tfmSY; jy++) {
	    complex<r_4> zz = acr( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	    tfmap(TFMtmidx, jy) = zz.real();
	  } 
	}  //----- end of time-frequency maps filling 
      }
      //---- pour les 6 cross-cor
      for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
	TVector< complex<float> > cxv = vismtx.Row( KVCXHH[k] );
	cxdt[k] = cxv(freqs).Sum();    // integration en bande 1
	cxdt[k+6] = cxv(freqsB).Sum();    // integration en bande large B
	cxdt[k+12] = cxv(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
      }
      //  On calcule les moyennes par bande de frequence
      // Moyennes dans la bande 1
      //       double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
      double fnorm=1./(double)(JFmax-JFmin+1);
      for(int k=0; k<4; k++)  acdt[k]*=fnorm; 
      acrdt[0]*=fnorm; 
      for(int k=0; k<6; k++)  cxdt[k]*=complex<r_8>(fnorm,0.);  
      // Moyennes dans la bande 2 
      fnorm=1./(double)(JFmaxB-JFminB+1);
      for(int k=0; k<4; k++)  acdt[k+4]*=fnorm; 
      acrdt[1]*=fnorm; 
      for(int k=0; k<6; k++)  cxdt[k+6]*=complex<r_8>(fnorm,0.);   
      // Moyennes ds +/5 MHz @1420 , exclu +/- 1 MHz centre a 1420 
      fnorm=1./(double)((JFmax21_5-JFmin21_5)-(JFmax21-JFmin21));
      for(int k=0; k<8; k++)  acdt[k+12] = (acdt[k+12]-acdt[k+8])*fnorm;
      acrdt[3] = (acrdt[3]-acrdt[2])*fnorm; 
      // Moyennes ds +/1 MHz @1420 
      fnorm=1./(double)(JFmax21-JFmin21+1);
      for(int k=0; k<4; k++)  acdt[k+8]*=fnorm; 
      acrdt[2]*=fnorm; 
      for(int k=0; k<6; k++)  cxdt[k+12]*=complex<r_8>(fnorm,0.);   
      
      // remplissage datatable 
      dt.NextRowPtr(rowp);
      rowp(0)=(int_4)cnt;   
      rowp(1)=(int_4)(tms5.DaysPart()-dateorg.DaysPart());   
      rowp(2)=(r_4)tms5.SecondsPart();
      for(int k=0; k<16; k++)   rowp(3+k) = (float)acdt[k];
      for(int k=0; k<18; k++)   rowp(19+k) = complex<float>((float)cxdt[k].real(), (float)cxdt[k].imag());
      rowp(37)=(int_4)(tms6.DaysPart()-dateorg.DaysPart());   
      rowp(38)=(r_4)tms6.SecondsPart();
      for(int k=0; k<4; k++)   rowp(39+k) = (float)acrdt[k];
      //  ... done 
      cnt++;
      
      if (prtlev>1) cout << "JFile="<<jf<<" DateFile5:"<<tms5<<" DateFile6:"<<tms6<<" p3hraw="<<acrdt[0]<<endl;
    }
 
    cout << " msvis2dt/Info: count="<<cnt<<" visimtx/meanspec read "<<endl;
    dt.SetShowMinMaxFlag(true);
    cout << dt;
    cout<<"  msvis2dt/Info: Saving auto/cross-corr=f(time) datatable to PPF file "<<outfile<<endl;
    POutPersist po(outfile);
    po<<dt; 

    if (FgTFM) {   // --- renormalizing and saving time-frequency maps 
      cout<<"  msvis2dt/Info: Saving time-frequency maps to PPF file "<<TFMoutfile<<endl;
      POutPersist potfm(TFMoutfile);
      const char* tfm_names[6]={"TFM_3Hraw", "TFM_1H", "TFM_2H", "TFM_3H", "TFM_4H"};
      for(int k=0; k<5; k++)  {   // loop over the 6 Xcor 
	TArray< r_4 > & tfmap = vtfm[k];
	tfmap /= (r_4)(1./(double)deltaIavg/(double)TFMfbin);
	potfm << PPFNameTag(tfm_names[k]) << tfmap;
      }
    }
    //    resu.Update();
    cout << resu;   // Update est fait lors du print
  }
  catch (PException& exc) {
    cerr << " msvis2dt.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n msvisi2dt.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " msvisi2dt.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

  cout << ">>>> msvisi2dt.cc ------- END ----------- RC=" << rc << endl;
  return rc;

}