source for visiavg (OP+RA)

this computes averages for time and visi (auto and HH cross correlations)

visiavg.cc

+// Utilisation de SOPHYA pour faciliter les tests ...
+#include "sopnamsp.h"
+#include "machdefs.h"
+
+/* ---------------------------------------------------------- 
+   Projet BAORadio/PAON4 - (C) LAL/IRFU  2017
+
+   visiavg: programme de lecture des fichiers matrices de 
+   visibilites de PAON4, calcul de visibilities moyennes 
+    en bin de temps et de frequence
+   O. Perdereau, R.Ansari   -  LAL
+   ---------------------------------------------------------- */
+
+// 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 "p4autils.h"
+#include "visip4reader.h"
+#include "p4gnugain.h"
+
+int Usage(void);
+int Usage(void)
+{
+  cout << " --- visiavg.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
+  cout << " Usage: visiavg [-arguments] \n" << endl;
+  P4AnaParams::UsageOptions();
+  cout<< endl;
+  return 1;
+}
+
+//----------------------------------------------------
+int main(int narg, const char* arg[])
+{
+  // --- Decoding parameters 
+  if( (narg<2) || ((narg>1)&&(strcmp(arg[1],"-h")==0) ) )  return Usage();
+  P4AnaParams params;
+  params.DecodeArgs(narg, arg);
+  string outfile = params.outfile_;
+  if (outfile.length()<1)  outfile = "visavg.ppf";
+  int deltaIavg = params.TFMtimebin_;
+  sa_size_t TFMfbin = params.TFMfreqbin_;
+  int Imin = params.Imin_, Imax = params.Imax_, Istep = params.Istep_; 
+  int prtlev = params.prtlev_;
+  bool FgTFMAC = true;
+  bool FgTFMCX = true;
+  string desctfmap;
+  bool FgTFM = params.fgTFM_;   // true -> create time-frequency maps
+
+  params.Print(cout);
+  cout <<"visiavg/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
+       <<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
+       <<"Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<" DeltaIAvg="<<deltaIavg<<"\n"
+       <<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
+
+  if (!FgTFM) {
+    cout<<" visiavg/parameter error : specify Time-Frequency map parameter with -tfm "<<endl;
+    return 5;
+  }
+
+  P4AVisiNumEncoder  visiencod;
+  vector<sa_size_t> KVAC = visiencod.getAllAutoCor();
+  vector<sa_size_t> KVCXHH = visiencod.getAllHCrossCor();
+  cout << " List of AutoCorrelation rows:"<<endl;
+  for(size_t k=0; k<KVAC.size(); k++) {
+    cout << "KVAC["<<k<<"]="<<KVAC[k]<<"  ->"<<visiencod.Convert2VisiName(KVAC[k])<<endl;
+  }
+  cout << " List of HH X-cor rows:"<<endl;
+  for(size_t k=0; k<KVCXHH.size(); k++) {
+    cout << "KVCXHH["<<k<<"]="<<KVCXHH[k]<<"  ->"<<visiencod.Convert2VisiName(KVCXHH[k])<<endl;
+  }
+  // --- Open file to store visibility matrices if requested
+
+  // ---
+  HiStatsInitiator _inia;
+  int rc = 0;
+  try {
+    ResourceUsage resu;
+
+    // Gain correction class
+    P4gnuGain p4g(params.gain_gnu_file_);
+
+    // setting up input visi reader
+    vector<string> paths; 
+    paths.push_back(params.inpath5_); 
+    paths.push_back(params.inpath6_);
+
+    
+    VisiP4ReaderBase * reader = VisiP4ReaderBase::getReader(paths);
+    VisiP4ReaderBase & vreader = (*reader);
+    vreader.setFreqReordering(params.fgreorderfreq_);
+    if (!params.fgserall_ && !params.fgtmsel_) {
+      cout << " vreader.SelectSerialNum(Imin="<<Imin<<" ,Imax="<<Imax<<" ,Istep="<<Istep<<")"<<endl;
+      vreader.SelectSerialNum(Imin,Imax,Istep);
+    }
+    else if (params.fgserall_) {
+      cout << " vreader.SelectAll() ... " << endl;
+      vreader.SelectAll();
+    }
+    else {
+      TimeStamp tustart = params.tmsel_tu_;  tustart.ShiftSeconds(-params.tmsel_duration_*30.);
+      TimeStamp tuend = params.tmsel_tu_;  tuend.ShiftSeconds(params.tmsel_duration_*30.);
+      cout << " vreader.SelectTimeFrame(TUStart="<<tustart.ToString()<<" ,TUEnd="<<tuend.ToString()
+	   <<" ,Istep="<<Istep<<")"<<endl;
+      vreader.SelectTimeFrame(tustart, tuend, Istep);
+    }
+    vreader.setPrintLevel(prtlev);
+
+    Imin = vreader.getSerialFirst(); Imax =  vreader.getSerialLast();   Istep =  vreader.getSerialStep();
+    cout << "visiavg/Info: processing visibility matrix serial/sequence number range "
+	 <<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
+
+    bool fgok=true;
+    // un vecteur avec les temps 
+    TVector< double > timevec((Imax-Imin)/Istep/deltaIavg); 
+    TMatrix< complex<r_4> > vismtx;
+    TMatrix< complex<r_4> > acsum;
+    TMatrix< complex<r_4> > cxsum;
+
+    TimeStamp dateobs, cfdate,datestart;
+    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 H-H cross-cor TimeFrequency maps 
+    vector< TArray< complex<r_4> > > vtfm;
+    //----- 8 auto-corr TimeFrequency maps 
+    vector< TArray< r_4 > > vtfmac;
+    
+    //---- 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)  break;
+      // Apply gain g(nu)
+      p4g.applyGain(vismtx);
+      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());
+	tfmSX=(Imax-Imin)/Istep/deltaIavg;
+	tfmSY=vismtx.NCols()/TFMfbin;
+	//allocating 8 Auto-Corr time-frequency maps 
+	cout<<"visiavg/Info: allocating 8 AutoCor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
+	for(int k=0; k<8; k++) vtfmac.push_back( TArray< r_4 >(tfmSX, tfmSY) ); 
+	//allocating 6 Cross-Corr H-H time-frequency maps 
+	cout<<"visiavg/Info: allocating H-H cross-cor Time-Frequency 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) );
+
+	// recupere le jour de depart @ 0h
+	datestart = TimeStamp(cfdate.DaysPart(),0.);
+ 
+      }
+      if (I==0) {   // start filling a new time bin 
+	dateobs=cfdate;
+	if (prtlev>0) 
+	  cout<<"visiavg/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
+	acsum = complex<r_4>(0.,0.);
+	cxsum = complex<r_4>(0.,0.);
+      }
+      
+      //   sum (integration) along the time axis 
+      for(size_t k=0; k<KVAC.size(); k++)      acsum.Row(k) += vismtx.Row(KVAC[k]);     // Les auto-correlations 
+      for(size_t k=0; k<KVCXHH.size(); 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 
+	for(size_t k=0; k<KVAC.size(); k++) {  // Loop over the 8 auto-correlations 
+	  TVector<r_4> vac = real(acsum.Row(k));
+	  if (TFMtmidx<tfmSX) {  // we check that our time index did not go beyond the allocated array size (might not be necessary)
+	    TArray< r_4 > & tfmap = vtfmac[k];
+	    for(sa_size_t jy=0; jy<tfmSY; jy++) {  // frequency binning 
+	      tfmap(TFMtmidx, jy) = vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
+	    } 
+	  }  
+	}  //----- end of loop over the 8 AutoCor
+	//---- On s'occupe des cross-correlations  1H-2H ... 3H-4H 
+ 	for(size_t k=0; k<KVCXHH.size(); k++)   {   // loop over the 6 Xcor 	  
+	  TVector< complex<r_4> > vcx = cxsum.Row(k);
+	  if (TFMtmidx<tfmSX) {  // we check that our time index did not go beyond the allocated array size (might not be necessary) 
+	    TArray< complex<r_4> > & tfmap = vtfm[k];
+	    for(sa_size_t jy=0; jy<tfmSY; jy++) {
+	      tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
+	    } 
+	  }  
+	}  //----- end of loop over the 6 Xcor 
+	timevec(TFMtmidx) = cfdate.TimeDifferenceSeconds(cfdate,datestart);
+	TFMtmidx++;
+	//  ... done 
+	I=0;  cntnt++;
+      }
+      cnt++;
+      if ((cnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
+	cout<<"visiavg/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt<<" /Max="<<Imax<<" DateObs="<<dateobs<<endl;
+	pcntnt=cntnt;
+      }
+    }
+    cout<<"visiavg/Info: count="<<cnt<<" visimtx read "<<endl;
+
+    // --- Sauvegarde cartes temps-frequence 
+    POutPersist potfm(outfile);
+    // --- renormalizing and saving AutoCorr time-frequency maps 
+    cout<<"  visiavg/Info: Saving 8 AutoCorr time-frequency maps to PPF file "<<outfile<<endl;
+    const char* tfm_names[8]={"TFM_1H", "TFM_2H", "TFM_3H", "TFM_4H", "TFM_1V", "TFM_2V", "TFM_3V", "TFM_4V"};
+    for(int k=0; k<8; k++)  {   // loop over the 8 AutoCorr 
+      TArray< r_4 > & tfmap = vtfmac[k];
+      tfmap *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
+      potfm << PPFNameTag(tfm_names[k]) << tfmap;
+    }
+    // --- renormalizing and saving H-H Cross-Corr time-frequency maps 
+    cout<<"  visiavg/Info: Saving 6 H-H cross-corr time-frequency maps to PPF file "<<outfile<<endl;
+    const char* tfmCC_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(tfmCC_names[k]) << tfmap;
+    }
+    potfm << PPFNameTag("TimeVec") << timevec ;
+    // --- FIN sauvegarde cartes temps-frequence 
+    //    resu.Update();
+    cout << resu;   // Update est fait lors du print
+  }
+  catch (PException& exc) {
+    cerr << " visiavg.cc catched PException " << exc.Msg() << endl;
+    rc = 77;
+  }  
+  catch (std::exception& sex) {
+    cerr << "\n visiavg.cc std::exception :" 
+         << (string)typeid(sex).name() << "\n msg= " 
+         << sex.what() << endl;
+    rc = 78;
+  }
+  catch (...) {
+    cerr << " visiavg.cc catched unknown (...) exception  " << endl; 
+    rc = 79; 
+  } 
+
+  cout << ">>>> visiavg.cc ------- END ----------- RC=" << rc << endl;
+  return rc;
+
+}
+