visiavg.cc 13.6 KB
Newer Older
OP's avatar
OP committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
// 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"

OP's avatar
OP committed
44 45 46
#include "fitsioserver.h"
#include "fiosinit.h"

OP's avatar
OP committed
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
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();
OP's avatar
OP committed
62
  FitsIOServerInit();
OP's avatar
OP committed
63 64
  P4AnaParams params;
  params.DecodeArgs(narg, arg);
OP's avatar
OP committed
65

OP's avatar
OP committed
66 67
  string outfile = params.outfile_;
  if (outfile.length()<1)  outfile = "visavg.ppf";
OP's avatar
OP committed
68 69 70 71 72

  string fitsoutfile = params.fitsoutfile_;
  if (fitsoutfile.length()>=1) {
    fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
  }
OP's avatar
OP committed
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
  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;
  }
OP's avatar
OP committed
104
  
OP's avatar
OP committed
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
  // ---
  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); 
OP's avatar
OP committed
147
    TVector< double > ravec((Imax-Imin)/Istep/deltaIavg); 
OP's avatar
OP committed
148 149
    TMatrix< complex<r_4> > vismtx;
    TMatrix< complex<r_4> > acsum;
OP's avatar
OP committed
150
    TMatrix< r_4 > acsum_sq; // wil sum only the real part ^2
OP's avatar
OP committed
151
    TMatrix< complex<r_4> > cxsum;
OP's avatar
OP committed
152 153 154
    // for sums of real and imag parts 
    TMatrix< r_4 > cxsum_sq_rp;
    TMatrix< r_4 > cxsum_sq_ip;
OP's avatar
OP committed
155 156 157 158 159 160

    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; 
OP's avatar
OP committed
161
    // for
OP's avatar
OP committed
162 163
    //----- 6 H-H cross-cor TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm;
OP's avatar
OP committed
164 165 166 167
    //----- 6 H-H cross-cor TimeFrequency maps for the variances of real and imag parts 
    vector< TArray< r_4 > > vtfm_rp_sq;
    vector< TArray< r_4 > > vtfm_ip_sq;

OP's avatar
OP committed
168 169
    //----- 8 auto-corr TimeFrequency maps 
    vector< TArray< r_4 > > vtfmac;
OP's avatar
OP committed
170 171
    //----- 8 auto-corr TimeFrequency variance maps 
    vector< TArray< r_4 > > vtfmac_sq;
OP's avatar
OP committed
172 173 174 175 176 177
    
    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

    while (fgok) {
OP's avatar
OP committed
178
      //reads next visimtx
OP's avatar
OP committed
179
      fgok=vreader.ReadNext(vismtx, cfdate, mttag);
OP's avatar
OP committed
180

OP's avatar
OP committed
181 182 183
      if (!fgok)  break;
      // Apply gain g(nu)
      p4g.applyGain(vismtx);
OP's avatar
OP committed
184

OP's avatar
OP committed
185 186
      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
	acsum.SetSize(8, vismtx.NCols());
OP's avatar
OP committed
187
	acsum_sq.SetSize(8, vismtx.NCols());
OP's avatar
OP committed
188
	cxsum.SetSize(6, vismtx.NCols());
OP's avatar
OP committed
189 190 191
	cxsum_sq_rp.SetSize(6, vismtx.NCols());
	cxsum_sq_ip.SetSize(6, vismtx.NCols());

OP's avatar
OP committed
192 193 194 195 196
	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) ); 
OP's avatar
OP committed
197 198 199
	cout <<" and 8 for the the variance maps "<<endl;
	for(int k=0; k<8; k++) vtfmac_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) ); 

OP's avatar
OP committed
200 201 202
	//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) );
OP's avatar
OP committed
203 204 205
	cout << "and the 2x6 for squares of rela & imag parts " << endl;
	for(int k=0; k<6; k++) vtfm_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
	for(int k=0; k<6; k++) vtfm_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
OP's avatar
OP committed
206 207 208 209 210 211 212 213 214 215

	// 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.);
OP's avatar
OP committed
216
	acsum_sq = 0.;
OP's avatar
OP committed
217
	cxsum = complex<r_4>(0.,0.);
OP's avatar
OP committed
218 219
	cxsum_sq_rp = 0.;
	cxsum_sq_ip = 0.;
OP's avatar
OP committed
220 221 222 223
      }
      
      //   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 
OP's avatar
OP committed
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
      for(size_t k=0; k<KVAC.size(); k++)      {
	TVector<r_4> tmp = real(vismtx.Row(KVAC[k]));
	acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ;
      }

      for(size_t k=0; k<KVCXHH.size(); k++){
	cxsum.Row(k) += vismtx.Row(KVCXHH[k]);   // les cross-correlations 

	TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]));
	cxsum_sq_rp.Row(k) +=  tmp.MulElt(tmp,tmp) ;

	tmp = imag(vismtx.Row(KVCXHH[k]));
	cxsum_sq_ip.Row(k) +=  tmp.MulElt(tmp,tmp) ;

      }


OP's avatar
OP committed
241 242 243 244 245 246
      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));
OP's avatar
OP committed
247
	  TVector<r_4> vacsq = acsum_sq.Row(k);
OP's avatar
OP committed
248 249
	  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];
OP's avatar
OP committed
250
	    TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
OP's avatar
OP committed
251 252
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {  // frequency binning 
	      tfmap(TFMtmidx, jy) = vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
OP's avatar
OP committed
253 254
	      tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();

OP's avatar
OP committed
255
	    } 
OP's avatar
OP committed
256 257
	  }
  
OP's avatar
OP committed
258
	}  //----- end of loop over the 8 AutoCor
OP's avatar
OP committed
259
	//---- On s'occupe des 6 cross-correlations  1H-2H ... 3H-4H 
OP's avatar
OP committed
260 261
 	for(size_t k=0; k<KVCXHH.size(); k++)   {   // loop over the 6 Xcor 	  
	  TVector< complex<r_4> > vcx = cxsum.Row(k);
OP's avatar
OP committed
262 263
	  TVector<r_4>  vcxprsq = cxsum_sq_rp.Row(k);
	  TVector<r_4>  vcxpisq = cxsum_sq_ip.Row(k);
OP's avatar
OP committed
264 265
	  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];
OP's avatar
OP committed
266 267
	    TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
	    TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
268 269
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {
	      tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
OP's avatar
OP committed
270 271
	      tfmapsqpr(TFMtmidx, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	      tfmapsqpi(TFMtmidx, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
OP's avatar
OP committed
272 273 274
	    } 
	  }  
	}  //----- end of loop over the 6 Xcor 
OP's avatar
OP committed
275 276 277 278
	
	double tdif =  cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
	timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart);	// centre du bin 
	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs);
OP's avatar
OP committed
279 280 281 282 283 284 285 286 287 288 289 290
	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;

OP's avatar
OP committed
291 292 293 294 295 296
    // --- Sauvegarde cartes temps-frequence en fits 
    if (fitsoutfile.length()>=1)
      FitsInOutFile fos(fitsoutfile, FitsInOutFile::Fits_Create);

    //    DataTable dt;

OP's avatar
OP committed
297 298 299 300
    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"};
OP's avatar
OP committed
301
    const char* tfmsq_names[8]={"VARTFM_1H", "VARTFM_2H", "VARTFM_3H", "VARTFM_4H", "VARTFM_1V", "VARTFM_2V", "VARTFM_3V", "VARTFM_4V"};
OP's avatar
OP committed
302 303 304
    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));
OP's avatar
OP committed
305 306 307
      TArray< r_4 > & tfmapsq = vtfmac_sq[k];
      tfmapsq *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
      tfmapsq = tfmapsq - tfmap.MulElt(tfmap,tfmap) ; 
OP's avatar
OP committed
308
      potfm << PPFNameTag(tfm_names[k]) << tfmap;
OP's avatar
OP committed
309
      potfm << PPFNameTag(tfmsq_names[k]) << tfmapsq;
OP's avatar
OP committed
310
      
OP's avatar
OP committed
311 312 313
    }
    // --- 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;
OP's avatar
OP committed
314

OP's avatar
OP committed
315
    const char* tfmCC_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
OP's avatar
OP committed
316 317 318
    const char* vrtfmCC_names[6]={"RVARTFM_1H2H", "RVARTFM_1H3H", "RVARTFM_1H4H", "RVARTFM_2H3H", "RVARTFM_2H4H", "RVARTFM_3H4H"};
    const char* vitfmCC_names[6]={"IVARTFM_1H2H", "IVARTFM_1H3H", "IVARTFM_1H4H", "IVARTFM_2H3H", "IVARTFM_2H4H", "IVARTFM_3H4H"};

OP's avatar
OP committed
319 320
    for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
      TArray< complex<r_4> > & tfmap = vtfm[k];
OP's avatar
OP committed
321 322 323

      TArray< r_4 > & tfmap_sqpr = vtfm_rp_sq[k];
      TArray< r_4 > & tfmap_sqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
324
      tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
OP's avatar
OP committed
325 326 327 328 329 330 331 332 333
      TArray< r_4 >  tfmapr = real(tfmap);
      TArray< r_4 >  tfmapi = imag(tfmap);

      tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
      tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
      tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
      tfmap_sqpr -= tfmapr ;
      tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
      tfmap_sqpi -= tfmapi;
OP's avatar
OP committed
334
      potfm << PPFNameTag(tfmCC_names[k]) << tfmap;
OP's avatar
OP committed
335 336
      potfm << PPFNameTag(vrtfmCC_names[k]) << tfmap_sqpr;
      potfm << PPFNameTag(vitfmCC_names[k]) << tfmap_sqpi;
OP's avatar
OP committed
337 338
    }
    potfm << PPFNameTag("TimeVec") << timevec ;
OP's avatar
OP committed
339
    potfm << PPFNameTag("RAVec") << ravec ;
OP's avatar
OP committed
340 341
    // --- FIN sauvegarde cartes temps-frequence 
    //    resu.Update();
OP's avatar
OP committed
342 343 344 345 346 347 348 349
    P4FreqBand myp4fre;
    TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
    double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;

    for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase +=myp4fre.getP4FreqResolution() )
      avg_freqs(kf) = frbase ;
    potfm << PPFNameTag("FreqVec") << avg_freqs;

OP's avatar
OP committed
350
    cout << resu;   // Update est fait lors du print
OP's avatar
OP committed
351 352
    

OP's avatar
OP committed
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373
  }
  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;

}