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

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

   Programme de lecture des fichiers matrices de 
   visibilites de PAON4 produits mfacq 
   R. Ansari, J.E. Campagne, C. Magneville   -  LAL/Irfu
   V0 : Fev 2015 
   ---------------------------------------------------------- */

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

// ---- compute gain for each channel from visibilies and save it to file 
void computeSaveGain(TMatrix< complex<r_4> > & vismtx_mean, string const & gainfile);
// ----- adaptation de la fonction de calcul de gain, ecrit par cmv pour l'analyse Amas (Aout 2013) 
Matrix computeGain2(Matrix & acs, int SZW=96,float frackeep=0.75);
// ----- fonction pour normaliser les gains 
Vector normalizeGain2(Matrix& gain);
 
//--------------------------- Fonctions de ce fichier   ------------------- 
int Usage(bool fgshort=true);
// int DecodeArgs(int narg, char* arg[]);

/* --Fonction-- */
int Usage(bool fgea)
{
  cout << " --- rdvisip4.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
  if (fgea) cout << " rdvisip4: Missing or wrong argument ! " << endl;
  cout << " Usage: rdvisip4 [-mr] InPathBAO5 InPathBAO6 Imin,Imax,step OutPPFFile GainPPFFile [PrintLev=1] \n"
       << "  -m: read also the mean Fourier coeff matrix \n"
       << "  -r: reorder frequence (for FFT firmeware)  \n"
       << "  -mr: do both  \n"
       << "  InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
       << "  Imin,Imax,step : read files for Imin<=iii<=Imax  iii+=step \n"
       << "  OutPPFFile: Output PPF file name (with the average visibility matrix) \n"
       << "  GainPPFFile: PPF file name containing the computed gain g(nu) \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<6) return Usage(true);
  if ((narg>1)&&((strcmp(arg[1],"-m")==0)||(strcmp(arg[1],"-r")==0)||(strcmp(arg[1],"-mr")==0)) && (narg<7))  return Usage(false);

  int off = 0; 
  bool fgrdmean = false;
  bool fgreorderfreq = false;
  if (strcmp(arg[1],"-m")==0 ){ fgrdmean = true; off++; }
  else if (strcmp(arg[1],"-r")==0 ){ fgreorderfreq = true; off++; }
  else if (strcmp(arg[1],"-mr")==0 ){ fgrdmean = fgreorderfreq = true; off++; }

  string path5 = arg[1+off];
  string path6 = arg[2+off];
  int Imin=0,Imax=10,Istep=1; 
  sscanf(arg[3+off],"%d,%d,%d",&Imin,&Imax,&Istep);
  string outfile=arg[4+off];
  string gainfile=arg[5+off];
  int prtlev=1;
  if (narg>6+off) prtlev=atoi(arg[6+off]);

  cout << " rdvisip4/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"  OutFile:"<<outfile
       << " PrtLev="<<prtlev
       <<endl;
  if(fgrdmean) cout << "Trying to read also the FFT mean coef matrix" <<endl;
  
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    ResourceUsage resu;
    vector<string> paths; 
    paths.push_back(path5); 
    paths.push_back(path6);
 
    int nfiles=(Imax-Imin+1)/Istep;
    int nmod=nfiles/20;
    if (nmod<1) nmod=1;
    VisiP4Reader vreader(paths, Imin,Imax,Istep,fgreorderfreq);
    vreader.setPrintLevel(prtlev);
    bool fgok=true;
    //vismtx: visib matrice, meanmtx: mean FFT coef, vismtx_center: visib matrix corrigee de l'offset (normalisee avec NPAQSUM)
    TMatrix< complex<r_4> > vismtx, meanmtx, vismtx_center;
    TMatrix< complex<r_4> > vismtx_mean, meanmtx_mean, vismtx_center_mean;
    TimeStamp dateobs;
    double mttag;
    int cnt=0;
    while (fgok) {
      if(fgrdmean) {
	fgok=vreader.ReadNext(vismtx, meanmtx, dateobs, mttag);
	vismtx_center = vreader.SubtractOffset(vismtx, meanmtx);
      } else {
	fgok=vreader.ReadNext(vismtx, dateobs, mttag);
      }
      if (fgok)  {
	if (cnt==0)  vismtx_mean=vismtx;
	else vismtx_mean+=vismtx;
	if(fgrdmean) {
	  if (cnt==0)  meanmtx_mean = meanmtx;
	  else meanmtx_mean += meanmtx;
	  if (cnt==0)  vismtx_center_mean = vismtx_center;
	  else  vismtx_center_mean += vismtx_center;
	}
	cnt++;
	if (cnt%nmod==0)  cout<<"rdvisip4/Info file read count="<<cnt<<" / nfiles="<<nfiles<<endl;
      } 
    }
    cout << " rdvisip4/Info: count="<<cnt<<" visimtx read "<<endl;
    if (cnt>0)  { 
      vismtx_mean /= complex<r_4>((r_4)cnt, (r_4)0.);
      cout << " rdvisip4/Info: computing gains ..."<<endl;
      computeSaveGain(vismtx_mean, gainfile);
      cout<<"  rdvisip4/Info: Saving vismtx_mean to PPF file "<<outfile<<endl;
      POutPersist po(outfile);
      po<<vismtx_mean;
      if(fgrdmean){
	meanmtx_mean /=  complex<r_4>((r_4)cnt, (r_4)0.);
	cout<<"  rdvisip4/Info: Saving meanmtx_mean to PPF file "<<outfile<<endl;
	po<<meanmtx_mean;

	vismtx_center_mean /=  complex<r_4>((r_4)cnt, (r_4)0.);
	cout<<"  rdvisip4/Info: Saving vismtx_center_mean to PPF file "<<outfile<<endl;
	po<<vismtx_center_mean;
      }

    } 
    //    resu.Update();
    cout << resu;   // Update est fait lors du print
  }
  catch (PException& exc) {
    cerr << " rdvisip4.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n rdvisip4.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " rdvisip4.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

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

}

//---------------------  Fonctions pour le calcul de gain -----------------------
//----- include specifique pour computeGain2()
#include "nbtrixx.h"
#include "cspline.h"

/*--Fonction--*/
void computeSaveGain(TMatrix< complex<r_4> > & vismtx_mean, string const & gainfile)
{
  sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
  Matrix acs(8, vismtx_mean.NCols());
  for(sa_size_t r=0; r<8; r++) 
    for(sa_size_t c=0; c<acs.NCols(); c++)  acs(r,c)=vismtx_mean(KVAC[r],c).real();
  int SZW=96;  float frackeep=0.75; 
  cout << " computeSaveGain/Info: computing gains with SZW="<<SZW<<" frackeep="<<frackeep<<" ..."<<endl;
  Matrix gains = computeGain2(acs, SZW, frackeep);
  
  Vector gn = normalizeGain2(gains);
  cout << " computeSaveGain/Info: saving gains, gn to "<<gainfile<<endl;
  POutPersist po(gainfile);
  po<<PPFNameTag("gains")<<gains;
  po<<PPFNameTag("gn")<<gn;
  return;
}

// ----- adaptation de la fonction de calcul de gain, ecrit par cmv pour l'analyse Amas (Aout 2013) 
/* --- fonction de calcul de gain avec fit de spline qui flag tout seul 
   les mauvais points du spectre - calcule la moyenne par fenetre (en 
   utilisant une mediane) et fit de spline sur les un point / fenetre 
   base sur la fonction rawgain2 de cmv (en date du gain (fonction de normalisation) 
   des spectres ON et OFF separes 
*/ 
/*--Fonction--*/
Matrix computeGain2(Matrix& acs, int SZW,float frackeep)
{
if(SZW%2) SZW++;
int nrows = acs.NRows(), ncols = acs.NCols();
double *val = new double[2*SZW+2];
Matrix gain(nrows,ncols); gain = 0.;
for(int r=0;r<nrows;r++) {
  TVector<r_8> DVon(ncols); DVon = 0.;
  for(int c=1; c<ncols-1; c++) DVon(c) = fabs(acs(r,c) - 0.5*(acs(r,c-1)+acs(r,c+1)));
  vector<double> vsplx, vsply;
  for(int c=0; c<ncols; c+=SZW/4) {
    // on enleve d'abord les RFI tres fortes
    int ng = 0;
    for(int cc=c;cc<c+SZW&&cc<ncols;cc++) {val[ng] = DVon(cc); ng++;}
    if(ng<2) continue;
    TabSort(ng,val);
    double dvcut = val[int(ng*frackeep)];
    //
    ng = 0;
    double x = 0.;
    for(int cc=c;cc<c+SZW&&cc<ncols;cc++) if(DVon(cc)<dvcut) {val[ng] = acs(r,cc); x+=cc; ng++;}
    if(ng!=0) x /= (double)ng;
    ////cout<<"rawgain2: r="<<r<<" c="<<c<<" dvcut="<<dvcut<<" ng="<<ng<<" x="<<x<<" v="<<val[ng/2]<<endl;
    if(ng<2) continue;
    TabSort(ng,val);
    vsplx.push_back(x);
    vsply.push_back(val[ng/2]);
  }
  TVector<r_8> vvx(vsplx), vvy(vsply);
  CSpline csp(vvx.Size(),vvx.Data(),vvy.Data());
  csp.ComputeCSpline();
  for(int c=0;c<ncols;c++) gain(r,c) = csp.CSplineInt(c);
 }
 delete [] val;
 return gain;
}

/*--Fonction--*/
Vector normalizeGain2(Matrix& gain)
{
  char buff[16];
  Vector gnorm(gain.NRows());
  for(sa_size_t r=0;r<gain.NRows();r++) {
    double gn=gain.Row(r).Sum();
    double fcg=1000./gn;    // we decide arbitray to renormalize the gain to have the sum[gains(i)]=1000 
    gain.Row(r) *= fcg;
    sprintf(buff,"GNORM%d",(int)r);
    gain.Info()[buff]=gn;
    gnorm(r)=gn/1000.;
    size_t badcnt=0;
#define LOWGTHR 1.e-9 
    for(sa_size_t c=0; c<gain.NCols(); c++) 
      if (gain(r,c)<LOWGTHR)  { 
	badcnt++;   gain(r,c)=LOWGTHR; 
      }
    cout << " normalizeGain2()/Info GNorm["<<r<<"]= "<<gn<<"  (badcount="<<badcnt<<")"<<endl;
  }
  return gnorm;
}