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

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

   visi2dtacx: programme de lecture des fichiers matrices de 
   visibilites de PAON4, remplissage et ecriture d'un 
   DataTable avec les cross-correlations en fonction du temps 
   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 << " --- visi2dtacx.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
  if (fgea) cout << " visi2dtacx: Missing or wrong argument ! " << endl;
  cout << " Usage: visi2dtacx [-tfm Nf] InPathBAO5 InPathBAO6 Imin,Imax GainPPFFile OutPPFFile Freq1 Freq2 DeltaFreq DeltaIAvg [PrtLev=0] \n"
       << "  -tfm Nf : create time-frequncy maps of cross-correlation binning Nf frequencies \n" 
       << "  InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
       << "  Imin,Imax : read files for Imin<=iii<=Imax  iii+=step \n"
       << "  GainPPFFile: Input gains PPF file name \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"
       << "  DeltaIAvg: compute average power every DeltaI input vismtx files, def=10 \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<10) 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<12) 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 gainfile=arg[offa+4];
  string outfile=arg[offa+5];
  // frequency range definition 
  double deltafreq=1.;
  double freq1=1375.;
  double freq2=1410.;
  int nargo=narg-offa;
  if ((nargo>7)&&(strcmp(arg[offa+6],"-")!=0))  freq1=atof(arg[offa+6]);
  if ((nargo>8)&&(strcmp(arg[offa+7],"-")!=0))  freq2=atof(arg[offa+7]);
  if ((nargo>6)&&(strcmp(arg[offa+8],"-")!=0))  deltafreq=atof(arg[offa+8]);
  // time averaging interval definition, in term of visibility files
  int deltaIavg=10;
  if ((nargo>9)&&(strcmp(arg[offa+9],"-")!=0))  deltaIavg=atoi(arg[offa+9]);
  if (deltaIavg<1)  deltaIavg=1;
  int prtlev=1;
  if (nargo>10) prtlev=atoi(arg[offa+10]);

  cout << " visi2dtacx/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 = "tfm_"+outfile;
  if (FgTFM) {
    cout << " visi2dtacx/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 << " visi2dtacx/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[8]={0,8,15,21,26,30,33,35};
    //  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,9,10,16}; 

    cout << " visi2dtacx/Info: reading input gain file"<<gainfile<<" ..."; 
    Matrix gains;
    Vector gn;
    {
      PInPersist pin(gainfile);
      pin>>PPFNameTag("gains")>>gains;
      pin>>PPFNameTag("gn")>>gn;
      for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";  
      // on transforme gains en facteur multiplicatif
      for(sa_size_t r=0;r<gains.NRows();r++) {
	for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=1./gains(r,c);
      }
      cout<<endl;
    }
    // Gains pour les cross-correlations g_ij(nu) = sqrt(g_ii(nu) 8 g_jj(nu) )
    Matrix gaincx(6, gains.NCols());
    for(sa_size_t c=0; c<gains.NCols(); c++)   {
      gaincx(0,c)=sqrt(gains(0,c)*gains(1,c));     //  1H-2H   
      gaincx(1,c)=sqrt(gains(0,c)*gains(2,c));     //  1H-3H   
      gaincx(2,c)=sqrt(gains(0,c)*gains(3,c));     //  1H-4H   
      gaincx(3,c)=sqrt(gains(1,c)*gains(2,c));     //  2H-3H   
      gaincx(4,c)=sqrt(gains(1,c)*gains(3,c));     //  2H-4H 
      gaincx(5,c)=sqrt(gains(2,c)*gains(3,c));     //  3H-4H     
    }
    cout << " DONE"<<endl;

    const char* dtnames[53]={"k","day","time",
			     "p1","p2","p3","p4","p5","p6","p7","p8",
			     "p1B","p2B","p3B","p4B","p5B","p6B","p7B","p8B",
			     "p1HI","p2HI","p3HI","p4HI","p5HI","p6HI","p7HI","p8HI",
			     "p1HI5","p2HI5","p3HI5","p4HI5","p5HI5","p6HI5","p7HI5","p8HI5",
                             "xh1h2","xh1h3","xh1h4","xh2h3","xh2h4","xh3h4",
                             "xh1h2B","xh1h3B","xh1h4B","xh2h3B","xh2h4B","xh3h4B",
			     "xh1h2HI","xh1h3HI","xh1h4HI","xh2h3HI","xh2h4HI","xh3h4HI"};

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

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

    vector<string> paths; 
    paths.push_back(path5); 
    paths.push_back(path6);
    VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
    vreader.setPrintLevel(prtlev);
    bool fgok=true;
    TMatrix< complex<r_4> > vismtx;
    TMatrix< complex<r_4> > acsum;
    TMatrix< complex<r_4> > cxsum;
    double acdt[32];   // les 4*8=32 valeurs d'autocorrelation pour remplissage dans la table 
    complex<double> cxdt[18];   // les 3*6=18 valeurs de cross-correlations pour remplissage dans la table 

    TimeStamp dateobs, cfdate;
    TimeStamp dateorg(2015,1,1,12,0,0.);  // Date origine 1 jan 2015
    double mttag;
    int cnt=0, cntnt=0, pcntnt=0;
    int I=0; 
    //----- 6 TimeFrequency maps 
    vector< TArray< complex<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;

    while (fgok) {
      fgok=vreader.ReadNext(vismtx, cfdate, mttag);
      if (fgok)  {
	if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
	  acsum.SetSize(8, vismtx.NCols());
	  cxsum.SetSize(6, vismtx.NCols());
	  if (FgTFM) {  //allocating time-frequency maps 
	    tfmSX=(Imax-Imin)/Istep/deltaIavg;
	    tfmSY=vismtx.NCols()/TFMfbin;
	    cout << " visi2dtacx/Info: allocating Time-Freqncy maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
	    for(int k=0; k<6; k++) vtfm.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) ); 
	  }
	}
	if (I==0) {   // start filling a new DataTable row 
	  dateobs=cfdate;
	  if (prtlev>0) 
	    cout << "visi2dtacx/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
	  for(int k=0; k<32; k++)  acdt[k]=0.;
	  for(int k=0; k<18; k++)  cxdt[k]=complex<r_8>(0.,0.);
	  acsum = complex<r_4>(0.,0.);
	  cxsum = complex<r_4>(0.,0.);
	}

	//   sum (integration) along the time axis 
	for(int k=0; k<8; k++)    acsum.Row(k) += vismtx.Row(KVAC[k]);     // Les auto-correlations 
	for(int k=0; k<6; k++)    cxsum.Row(k) += vismtx.Row(KVCXHH[k]);   // les cross-correlations 
	I++;    // incrementing DeltaTime counter 

	if (I==deltaIavg) {
	  //---- On s'occupe d'abord des autocorrelations P1 ... P8 
	  Vector vac(acsum.NCols());
	  for(int k=0; k<8; k++)  {
	  // correct for gain(nu)
	    for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=acsum(k,c).real()*gains(k,c);
	    acdt[k] += vac(freqs).Sum();    // integration en bande 1
	    acdt[k+8] += vac(freqsB).Sum();    // integration en bande large B
	    acdt[k+16] += vac(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
	    acdt[k+24] += vac(freqs21_5).Sum();   // integration a +/- 5 MHz @ 1420 MHz 
	  }
	  //---- On s'occupe des cross-correlations  1H-2H ... 3H-4H 
	  TVector< complex<r_8> > vcx(acsum.NCols());
	  for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
	    for(sa_size_t c=0; c<vcx.Size(); c++) 
	      vcx(c)=complex<r_8>((r_8)cxsum(k,c).real(), (r_8)cxsum(k,c).imag())*complex<r_8>(gaincx(k,c), 0.);
	    cxdt[k] += vcx(freqs).Sum();    // integration en bande 1
	    cxdt[k+6] += vcx(freqsB).Sum();    // integration en bande large B
	    cxdt[k+12] += vcx(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 

	    if (FgTFM && (TFMtmidx<tfmSX)) {  //filling time-frequency maps 
	      TArray< complex<r_4> > & tfmap = vtfm[k];
	      for(sa_size_t jy=0; jy<tfmSY; jy++) {
		complex<r_8> zz = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
		tfmap(TFMtmidx, jy) = complex<r_4>(zz.real(), zz.imag());
	      } 
	    }  //----- end of time-frequency maps filling 
	  }  //----- end of loop over the 6 Xcor 
	  TFMtmidx++;

	  // Moyennes dans la bande 1
	  double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
	  for(int k=0; k<8; k++)  acdt[k]*=fnorm; 
	  for(int k=0; k<6; k++)  cxdt[k]*=complex<r_8>(fnorm,0.);   
	  // Moyennes dans la bande 2 
	  fnorm=(1./(double)deltaIavg)/(double)(JFmaxB-JFminB+1);
	  for(int k=0; k<8; k++)  acdt[k+8]*=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)deltaIavg)/(double)((JFmax21_5-JFmin21_5)-(JFmax21-JFmin21));
	  for(int k=0; k<8; k++)  acdt[k+24] = (acdt[k+24]-acdt[k+16])*fnorm;
	  // Moyennes ds +/1 MHz @1420 
	  fnorm=(1./(double)deltaIavg)/(double)(JFmax21-JFmin21+1);
	  for(int k=0; k<8; k++)  acdt[k+16]*=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)(dateobs.DaysPart()-dateorg.DaysPart());   
	  rowp(2)=(r_4)dateobs.SecondsPart();
	  for(int k=0; k<32; k++)   rowp(3+k) = (float)acdt[k];
	  for(int k=0; k<18; k++)   rowp(35+k) = complex<float>((float)cxdt[k].real(), (float)cxdt[k].imag());
	  //  ... done 
	  I=0;  cntnt++;
	}
	cnt++;
	if ((cntnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
	  cout << " visi2dtacx/Info: fill-count="<<cntnt<<" fileCount="<<cnt<<" /Max="<<Imax
	       <<" DateObs="<<dateobs<<endl;
	  pcntnt=cntnt;
	}
      } 
    }
    cout << " visi2dtacx/Info: count="<<cnt<<" visimtx read "<<endl;
    dt.SetShowMinMaxFlag(true);
    cout << dt;
    cout<<"  visi2dtacx/Info: Saving auto-corr=f(time) datatable to PPF file "<<outfile<<endl;
    POutPersist po(outfile);
    po<<dt; 

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

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

}