// 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 <sys/stat.h>
#include <iostream>
#include <string>

#include "pexceptions.h"
#include "tvector.h"
#include "fioarr.h"

//--------------------------- Fonctions de ce fichier   ------------------- 
int Usage(void);
int Usage(void)
{
cout<<" --- svv2mtx2.cc : Read PPF files produced by mcrd/visfmib"<<endl
    <<"     and make time(row) - mean_frequency(col) matrices"<<endl
    <<"Usage: svv2mtx2 InOutPath Imin,Imax NumFreq1,NumFreq2,NBinFreq"<<endl 
    << " svv2mtx2 -h for detailed instructions"<<endl;
return 1;
}

//----------------------------------------------------
int main(int narg, char* arg[])
{
  // --- Decodage des arguments et traitement 
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage();
  if (narg<=1) return Usage();
  string inoutpath = arg[1];
  int ifilmin=0, ifilmax=-1;
  if(narg>2) sscanf(arg[2],"%d,%d",&ifilmin,&ifilmax);
  if(ifilmin<0) ifilmin=0;
  int jfr1=0, jfr2=0, nbinfreq=0;
  if(narg>3) sscanf(arg[3],"%d,%d,%d",&jfr1,&jfr2,&nbinfreq);
  char str[1024];

  // --- recherche des fichiers de visibilites
  int nfile = 0;
  {
  struct stat buffer;
  int i1=ifilmin, i2=(ifilmax<ifilmin) ? 999999999: ifilmax;
  for(int ifile=i1; ifile<=i2; ifile++) {
    sprintf(str, "%s/vismtx%d.ppf",inoutpath.c_str(),ifile);
    int status = stat(str,&buffer);
    if(status) break;
    ifilmax = ifile;
  }
  nfile = ifilmax-ifilmin+1;
  cout<<"Found "<<nfile<<" files from "<<ifilmin<<" to "<<ifilmax<<endl;
  if(nfile==0 || ifilmax<ifilmin) return -3;
  }


  //--------------------------------------------------------------------
  int rc = 0;
  try {

  // --- read visibilities coding vector
  cout<<"reading chanum vector"<<endl;
  TVector<uint_4> chanum;
  {
  sprintf(str,"%s/chanum.ppf",inoutpath.c_str());
  PInPersist pic(str);
  pic >> chanum;
  cout<<"number of visib is "<<chanum.Size()<<endl;
  }
  vector< TMatrix< complex<r_4> >* > vMVis;
  for(int i=0;i<chanum.Size();i++) vMVis.push_back(new TMatrix< complex<r_4> >);

  // --- read visibility files
  int nfreq0, nvisi, nmoyfreq = 0;
  TVector<r_8> Ttag(nfile), Freq;
  TVector<int_4> Npaqsum(nfile), NFreq;
  for(int ifile=ifilmin; ifile<=ifilmax; ifile++) {
    sprintf(str, "%s/vismtx%d.ppf",inoutpath.c_str(),ifile);
    cout<<ifile<<" opening: "<<str<<endl;
    PInPersist pin(str);
    TMatrix< complex<r_4> > vismtx;
    pin >> vismtx;

    Ttag(ifile-ifilmin) = (double)vismtx.Info()["MeanTT"]/125.e6;
    if(ifile==ifilmin) printf("Reference TimeTag is %.2f = %.5f sec\n"
                             ,(double)vismtx.Info()["MeanTT"],Ttag(0));
    uint_4 nmean = vismtx.Info()["NPAQSUM"];
    Npaqsum(ifile-ifilmin) = (int_4)nmean;
    nfreq0 = vismtx.NCols();
    nvisi = vismtx.NRows();

    // --- For initialisation purposes
    if (ifile==ifilmin) {
      cout<<"vismtx: number of frequencies = "<<nfreq0
          <<" , visib = "<<nvisi<<" , nmean = "<<nmean<<endl;
      if(chanum.Size()!=nvisi)
        throw ParmError("ERROR: nvisi != chanum.Size() !!!");
      // how many average frequencies ?
      if(jfr1<=0) jfr1=1;
      if(jfr1>=nfreq0) jfr1=nfreq0-1;
      if(jfr2<jfr1 || jfr2>=nfreq0) jfr2=nfreq0-1;
      if(nbinfreq<=0) nbinfreq=1;
      if(nbinfreq>jfr2-jfr1+1) nbinfreq=jfr2-jfr1+1;
      nmoyfreq = (jfr2-jfr1+1)/nbinfreq;
      if(nmoyfreq*nbinfreq<jfr2-jfr1+1) nmoyfreq++;
      cout<<"frequency averaging: jfr1="<<jfr1<<" jfr2="<<jfr2
          <<" nbinfreq="<<nbinfreq<<" (nmoyfreq="<<nmoyfreq<<")"<<endl;
      // average frequencies value
      Freq.ReSize(nbinfreq); Freq = 0.;
      NFreq.ReSize(nbinfreq); NFreq = 0;
      cout<<"Frequency average:"<<endl;
      for(int i=0;i<nbinfreq;i++) {
        Freq(i) = 0.;
        int nj = 0;
        for(int j=0;j<nmoyfreq;j++) {
          int ip = jfr1 + i*nmoyfreq + j;
          if(ip>=nfreq0) break;
          Freq(i) += ip;
          nj++;
	}
        if(nj>0) Freq(i) /= (double)nj;
        NFreq(i) = nj;
        cout<<"  F("<<i<<") = "<<Freq(i)<<" (nj="<<NFreq(i)
            <<" ["<<jfr1+i*nmoyfreq<<","<<jfr1+i*nmoyfreq+nj-1<<"])"<<endl;
      }
      // allocate visib matrice <f> vs t
      cout<<"allocating "<<nvisi<<" visibility matrices ("<<nfile<<","<<nbinfreq<<") "
          <<nvisi*nfile*nbinfreq*sizeof(complex<r_4>)/1.e6<<" Mo"<<endl;
      for(int ivi=0;ivi<nvisi;ivi++) {
        vMVis[ivi]->ReSize(nfile,nbinfreq);
        (*vMVis[ivi]) = complex<r_4>(0.,0.);
      }
    }

    // --- Fill visibility matrix
    for(int ivi=0;ivi<nvisi;ivi++) {
      for(int ifr=0;ifr<nbinfreq;ifr++) {
        int nj = 0;
        complex<r_8> v(0.,0.);
        for(int j=0;j<nmoyfreq;j++) {
          int ifrp = jfr1 + ifr*nmoyfreq + j;
          if(ifrp>=nfreq0) break;
          v += vismtx(ivi,ifrp);
          nj++;
	}
        if(nj>0) v /= (double)nj;
        (*vMVis[ivi])(ifile-ifilmin,ifr) = v / (double)nmean;
      }
    }
  }

  sprintf(str, "%s/svvdt2.ppf",inoutpath.c_str());
  cout<<"writing visibility matrix to file "<<str<<endl;
  POutPersist pos(str);
  pos.PutObject(chanum,"chanum");
  pos.PutObject(Ttag,"ttsec");
  pos.PutObject(Npaqsum,"npaqsum");
  pos.PutObject(Freq,"ifreq");
  pos.PutObject(NFreq,"nifreq");
  for(unsigned int i=0;i<vMVis.size();i++) {
    sprintf(str,"visft%d",chanum(i));
    pos.PutObject(*vMVis[i],str);
    delete vMVis[i];
  }
  cout<<" -------- End of Job --------"<<endl;
  }
  //--------------------------------------------------------------------
  catch(PException& exc) {
    cerr<<" svv2mtx2.cc catched PException "<<exc.Msg()<<endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr<<"\n svv2mtx2.cc std::exception :" 
        <<(string)typeid(sex).name() << "\n msg= "<<sex.what()<<endl;
    rc = 78;
  }
  catch(...) {
    cerr<<" svv2mtx2.cc catched unknown (...) exception  "<<endl; 
    rc = 79; 
  } 

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