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

/* ---------------------------------------------------------- 
   Projet BAORadio - (C) LAL/IRFU  2008-2010

   Programme de lecture des fichiers vecteurs/matrices de 
   visibilites produits par mcrd / vismfib 
   R. Ansari, C. Magneville   -  LAL/Irfu
   V : Mai 2009
   ---------------------------------------------------------- */

// 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"


//--------------------------- Fonctions de ce fichier   ------------------- 
int Usage(bool fgshort=true);

int DecodeProc(int narg, char* arg[]);
int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, 
		vector<sa_size_t> tfrlist);
int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > vismtx_, TVector< uint_4> chanum_,
		  sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) ;

int DecodeProcVJun09(int narg, char* arg[]);
int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card);
//------------------------------------------------------------------------------------------------------------

/* --Fonction-- */
int Usage(bool fgshort)
{
  cout << " --- svv2mtx.cc : Read PPF files produced by mcrd/visfmib to make time-frequency\n"
       << "     matrices and Visibilites=f(time) datatable " << endl;
  cout << " Usage: svv2mtx -mcrd InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [card=1] \n" 
       << "   OR \n" 
       << "        svv2mtx InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq VisMtxRowList" << endl; 
  if (fgshort) {
    cout << " svv2mtx -h for detailed instructions" << endl;
    return 1;
  }
  cout << " -mcrd : Read mcrd output files -> ProcSVFilesVJun09() \n" 
       << " InOutPath : Input/Output directory name \n" 
       << " Imin,Imax,IStep: Input PPF files sequence number \n"
       << "    FileNames=InOutPath/Ch12_II.fits Imin<=II<=Imax II+=IStep \n" 
       << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple\n"
       << " VisMtxRowList : List of visibiliy matrix rows (0,2,...) -> time-freq map \n"
       << " card=1 Ch12 , card=2 Ch34    (mcrd output files) " << endl;
  return 1;
}

//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
  if (narg<4) return Usage(true);
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    bool fgvjun09=false;
    if ((strcmp(arg[1],"-mcrd")==0))  fgvjun09=true;
    ResourceUsage resu;
    if (fgvjun09)  DecodeProcVJun09(narg-2, arg+2);
    else  DecodeProc(narg-1, arg+1);
    resu.Update();
    cout << resu;
  }
  catch (PException& exc) {
    cerr << " svv2mtx.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n svv2mtx.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " svv2mtx.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

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

}

//-------------------------------------------------------------------- 
//       Traitement fichiers produits par vismfib  (V Nov09)
/* --Fonction-- */
int DecodeProc(int narg, char* arg[])
{
  // Decodage des arguments et traitement 
  string inoutpath = arg[0];    
  int imin=0;
  int imax=0;
  int istep=1;
  sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
  int jf1=0;
  int jf2=0;
  int nfreq=0;
  sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
  int card=1;
  vector<sa_size_t> rowlist;
  if (narg>3) { 
    EnumeratedSequence es;
    int nbad;
    es.Append(arg[3], nbad, ",");
    for(int j=0; j<es.Size(); j++) rowlist.push_back((int)es.Value(j));
  }
  if (rowlist.size()<1) rowlist.push_back(0);

  cout << " ----- svv2mtx/DecodeProc - Start - InOutPath= " << inoutpath << " IMin,Max,Step=" 
       << imin << "," << imax << "," << istep << " Card=" << card << endl;
  cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
  cout << " RowList: " ;
  for(int j=0; j<rowlist.size(); j++)  cout << rowlist[j] << " , " ;
  cout << endl;
  int rc=ProcSVFiles(inoutpath, imin, imax, istep, jf1, jf2, nfreq, rowlist);
  return rc;
}

// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
/* --Fonction-- */
int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, 
		vector<sa_size_t> rowlist)
{
  Timer tm("ProcSVFiles");
  
  DataTable visdt;
  visdt.AddDoubleColumn("mfc");
  visdt.AddDoubleColumn("mtt");
  visdt.AddIntegerColumn("jfreq");
  visdt.AddIntegerColumn("numch");
  visdt.AddFloatColumn("vre");
  visdt.AddFloatColumn("vim");

  vector< TMatrix< complex<r_4> > > vmtf;

  if (jf1<1)  jf1=1;
  if ((jf2<1)||(jf2<jf1))  jf2=jf1; 
  if (nfreq<1) nfreq=1;
  int djf=(jf2-jf1)/nfreq;
  if (djf<1) djf=1;

  char fname[1024];

  cout << " ---- ProcSVFiles()-Begin - Reading chanum vector" << endl;
  /*
  TVector< uint_4 > chanum(8*17);
  int ku=0; 
  for(int iu=1; iu<=16; iu++) 
    for(int ju=iu; ju<=16; ju++)  {
      chanum(ku)=iu*1000+ju;  ku++;
    }
  */
  TVector< uint_4 > chanum;      
  {
  sprintf(fname, "%s/chanum_0.ppf",inoutpath.c_str());
  PInPersist pic(fname);
  pic >> PPFNameTag("chanpairnum") >> chanum;
  }
  
  cout << " ---- ProcSVFiles()-Read chanum done " << endl;
  sa_size_t nrows = (imax-imin+1)/istep;
  sa_size_t kr=0;

  sa_size_t mtf_binfreq=25;
  sa_size_t mtf_bintime=5;

  sa_size_t ncols=1;

  for(int ifile=imin; ifile<=imax; ifile+=istep) {
    sprintf(fname, "%s/vismtx_0_%d.ppf",inoutpath.c_str(),ifile);
    cout << " ProcSVFiles[" << ifile << "] opening file " << fname << endl;
    PInPersist pin(fname);
    TMatrix< complex<r_4> > vismtx;
    pin >> vismtx;

    if (ifile==imin) {
      ncols = vismtx.NCols();
      cout << " ProcSVFilesVJun09/Info: Output Time-Frequency matrices NRows=NFiles" 
           << nrows << " NCols=NFreq=" << ncols << endl;
      for(size_t j=0; j<rowlist.size(); j++)   
	vmtf.push_back(TMatrix< complex<r_4> >(ncols/mtf_binfreq+1, nrows/mtf_bintime+1));
    }
    /*
    for(size_t j=0; j<rowlist.size(); j++) 
      vmtf[j].Row(kr) = vismtx.Row(rowlist[j]);
    */
    for(size_t j=0; j<rowlist.size(); j++) 
      for(sa_size_t icf=0; icf<ncols; icf++)  {
	vmtf[j](icf/mtf_binfreq,kr/mtf_bintime)+=vismtx(rowlist[j],icf);
    }
    //    cout << " DBG* kr=" << kr << " kr/mtf_bintime=" << kr/mtf_bintime 
    //	   << " ncols/2/mtf_binfreq=" << ncols/2/mtf_binfreq << endl;
    kr++;

     
// Calcul moyenne dans des bandes en frequence 
    FillVisDTable(visdt, vismtx, chanum, jf1, jf2, djf);
  }

  sprintf(fname, "%s/vistfreqmtx.ppf",inoutpath.c_str());
  cout << "ProcSVFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;  
  POutPersist po(fname);
  char tagb[64];
  for(size_t j=0; j<rowlist.size(); j++)  {
    sprintf(tagb,"visft%d", chanum(rowlist[j]));
    po << PPFNameTag(tagb) << vmtf[j];
  }
  cout << visdt;
  sprintf(fname, "%s/svvdt.ppf",inoutpath.c_str());
  cout << "ProcSVFile: writing visibility datatable to file " << fname << endl;
  POutPersist pod(fname);
  pod << visdt;

  return 0;
}

/* --Fonction-- */
int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > vismtx_, TVector< uint_4> chanum_,
		  sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) 
{
  double xnt_[20];
  xnt_[0]=(double)vismtx_.Info()["MeanFC"];  
  xnt_[1]=(double)vismtx_.Info()["MeanTT"]/1.25e8;

  uint_4 nmean_=vismtx_.Info()["NPAQSUM"];
  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);
	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);
	xnt_[4]=moyreal/(r_4)(nmean_*djf_);
	xnt_[5]=moyimag/(r_4)(nmean_*djf_);
	visdt_.AddRow(xnt_);
      }
    }
  }
  return 0;
}


//-------------------------------------------------------------------- 
//       Traitement fichiers produits par mcrd (V Jun09) 
/* --Fonction-- */
int DecodeProcVJun09(int narg, char* arg[])
{
  // Decodage des arguments et traitement 
  string inoutpath = arg[0];    
  int imin=0;
  int imax=0;
  int istep=1;
  sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
  int jf1=0;
  int jf2=0;
  int nfreq=0;
  sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
  int card=1;
  if (narg>3) card=atoi(arg[3]); 
  cout << " ----- svv2mtx/DecodeProcVJun09 - Start - InOutPath= " << inoutpath << " IMin,Max,Step=" 
       << imin << "," << imax << "," << istep << " Card=" << card << endl;
  cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
  int rc=ProcSVFilesVJun09(inoutpath, imin, imax, istep, jf1, jf2, nfreq, card);
  return rc;
}


// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
/* --Fonction-- */
int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card)
{
  Timer tm("ProcSVFiles");
  char fname[512];
// NTuple 
  const char* nnames[10] = {"fcsm","ttsm","jfreq","s1","s2","s12","s12re","s12im","s12phi","s12mod"};
  NTuple nt(10, nnames);
  double xnt[15];
  uint_4 nmnt = 0;
  double ms1,ms2,ms12,ms12re,ms12im,ms12phi,ms12mod;

  TMatrix<r_4> s1, s2;
  TMatrix<r_4> v12re, v12im, v12phi,v12mod;
  sa_size_t ncols = (imax-imin+1)/istep;
  sa_size_t nrows = 10;
  sa_size_t kc=0;
  for(int ifile=imin; ifile<=imax; ifile+=istep) {
    if (card==2) 
      sprintf(fname, "%s/Ch34_%d.ppf",inoutpath.c_str(),ifile);
    else 
      sprintf(fname, "%s/Ch12_%d.ppf",inoutpath.c_str(),ifile);
    cout << " ProcSVFilesVJun09[" << ifile << "] opening file " << fname << endl;
    PInPersist pin(fname);
    string tag1="specV1";
    string tag2="specV2";
    string tag12="visiV12";
    if (card==2) {
      tag1 = "specV3";
      tag2 = "specV4";
      tag12="visiV34";
    }
    TVector<r_4> sv1;
    TVector<r_4> sv2; 
    TVector< complex<r_4> > vv12;
    pin >> PPFNameTag(tag1) >> sv1; 
    pin >> PPFNameTag(tag2) >> sv2; 
    pin >> PPFNameTag(tag12) >> vv12; 
    if (ifile==imin) {
      nrows = sv1.Size();
      cout << " ProcSVFilesVJun09/Info: Output s1,s2 matrix size NRows=NFreq=" 
           << nrows << " NCols=NFiles=" << ncols << endl;
      s1.SetSize(nrows, ncols);
      s2.SetSize(nrows, ncols);
      nrows = vv12.Size();
      cout << " ProcSVFilesVJun09/Info: Output v12 matrix size NRows=NFreq=" 
           << nrows << " NCols=NFiles=" << ncols << endl;
      v12re.SetSize(nrows, ncols);
      v12im.SetSize(nrows, ncols);
      v12phi.SetSize(nrows, ncols);
      v12mod.SetSize(nrows, ncols);
    }
    s1.Column(kc) = sv1;
    s2.Column(kc) = sv2;
    v12re.Column(kc) = real(vv12);
    v12im.Column(kc) = imag(vv12);
    v12phi.Column(kc) = phase(vv12);
    v12mod.Column(kc) = module(vv12);

// Calcul moyenne dans des bandes en frequence 
    int deltajf=(jf2-jf1)/nfreq;
    if (deltajf<1) deltajf=1; 
    for(int kf=0; kf<nfreq; kf++) {
      sa_size_t jfstart=jf1+kf*deltajf;
      sa_size_t jfend=jfstart+deltajf;
      if (jfend>jf2)  break;
      nmnt=0;  ms1=ms2=ms12=ms12re=ms12im=ms12phi=ms12mod=0.;
      for(sa_size_t jf=jfstart; jf<jfend; jf++) { 
        ms1 += s1(jf,kc);
        ms2 += s2(jf,kc);
        ms12re += v12re(jf,kc);
        ms12im += v12im(jf,kc);
        ms12phi += v12phi(jf,kc);
        ms12mod += v12mod(jf,kc);
      }
      nmnt = (jfend-jfstart);
      if (nmnt>0)  {
        double fnorm = (double)nmnt; 
        xnt[0] = ((int_8)(sv1.Info()["StartFC"])+(int_8)(sv1.Info()["EndFC"]))*0.5;
        xnt[1] = ((int_8)(sv1.Info()["StartTT"])+(int_8)(sv1.Info()["EndTT"]))*0.5;
        xnt[2] = kf;
        xnt[3] = ms1/fnorm;
        xnt[4] = ms2/fnorm;
        xnt[5] = ms12/fnorm;
        xnt[6] = ms12re/fnorm;
        xnt[7] = ms12im/fnorm;
        xnt[8] = ms12phi/fnorm;
        xnt[9] = ms12mod/fnorm;
        nt.Fill(xnt);
      }
    }
    kc++;

  }
  if (card==2) 
    sprintf(fname, "%s/Ch34mtx.ppf",inoutpath.c_str());
  else 
    sprintf(fname, "%s/Ch12mtx.ppf",inoutpath.c_str());

  cout << nt;
  cout << "ProcSVFilesVJun09: Opening file " << fname << " for writing" << endl;  
  POutPersist po(fname);
  string tag1="s1";
  string tag2="s2";
  string tag12r="v12re";
  string tag12i="v12im";
  string tag12p="v12phi";
  string tagnt="nt12";
  if (card==2) {
    tag1="s3";
    tag2="s4";
    tag12r="v34re";
    tag12i="v34im";
    tag12p="v34phi";
    tagnt="nt34";
    }
  po <<  PPFNameTag(tag1) << s1;
  po <<  PPFNameTag(tag2) << s2;
  po <<  PPFNameTag(tag12r) << v12re;
  po <<  PPFNameTag(tag12i) << v12im;
  po <<  PPFNameTag(tag12p) << v12phi;
  po <<  PPFNameTag(tagnt) << nt;
  cout << "ProcSVFilesVJun09: Matrices s1, s2, v12re, v12im, v12phi, NTuple nt written to file " << fname << endl;  
  return 0;
}