visiavg.cc 13.8 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
// 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 
40
#include "visp4winreader.h"
OP's avatar
OP committed
41

OP's avatar
OP committed
42 43 44
#include "fitsioserver.h"
#include "fiosinit.h"

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

OP's avatar
OP committed
64 65
  string outfile = params.outfile_;
  if (outfile.length()<1)  outfile = "visavg.ppf";
OP's avatar
OP committed
66

OP's avatar
OP committed
67 68 69 70
  string fitsoutfile = params.fitsoutfile_;
  if (fitsoutfile.length()>=1) {
    fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
  }
OP's avatar
modifs  
OP committed
71

OP's avatar
OP committed
72 73 74 75 76 77 78 79 80 81 82
  int deltaIavg = params.TFMtimebin_;
  sa_size_t TFMfbin = params.TFMfreqbin_;
  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"
83
       <<" DeltaIAvg="<<deltaIavg<<"\n"
OP's avatar
OP committed
84 85 86 87 88 89 90 91 92 93
       <<"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();
94 95

  
OP's avatar
OP committed
96 97 98 99 100 101 102 103
  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
modifs  
OP committed
104 105


OP's avatar
OP committed
106 107 108
  // ---
  HiStatsInitiator _inia;
  int rc = 0;
OP's avatar
modifs  
OP committed
109

OP's avatar
OP committed
110 111 112 113
  try {
    ResourceUsage resu;


114
    VisiP4WindowReader wreader(params);
115

116 117 118
    long Imin = wreader.getReader().getSerialFirst();
    long Imax =  wreader.getReader().getSerialLast();
    long Istep =  wreader.getReader().getSerialStep();
OP's avatar
OP committed
119 120
    cout << "visiavg/Info: processing visibility matrix serial/sequence number range "
	 <<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
121
    cout << " WindowSize="<<wreader.getWindowSize()<<"  -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;
OP's avatar
OP committed
122 123

    bool fgok=true;
OP's avatar
OP committed
124 125
    // vecteur de noms 
    vector <string> ext_names;
OP's avatar
OP committed
126
    // un vecteur avec les temps 
127 128
    TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg); 
    TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg); 
OP's avatar
OP committed
129 130
    TMatrix< complex<r_4> > vismtx;
    TMatrix< complex<r_4> > acsum;
OP's avatar
OP committed
131
    TMatrix< r_4 > acsum_sq; // wil sum only the real part ^2
OP's avatar
OP committed
132
    TMatrix< complex<r_4> > cxsum;
OP's avatar
OP committed
133 134 135
    // for sums of real and imag parts 
    TMatrix< r_4 > cxsum_sq_rp;
    TMatrix< r_4 > cxsum_sq_ip;
OP's avatar
OP committed
136 137 138 139 140 141

    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
142
    // for
OP's avatar
OP committed
143 144
    //----- 6 H-H cross-cor TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm;
OP's avatar
OP committed
145 146 147 148
    //----- 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
149 150
    //----- 8 auto-corr TimeFrequency maps 
    vector< TArray< r_4 > > vtfmac;
OP's avatar
OP committed
151 152
    //----- 8 auto-corr TimeFrequency variance maps 
    vector< TArray< r_4 > > vtfmac_sq;
OP's avatar
OP committed
153 154 155 156 157 158
    
    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

    while (fgok) {
159
      //reads next visibility matrix window 
160 161
      fgok = wreader.Shift();
      
162 163
      if (!fgok)  break;
     
164
      vismtx = wreader.getAverageVisMtx(cfdate);
165
      
166
      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
OP's avatar
OP committed
167
	acsum.SetSize(8, vismtx.NCols());
OP's avatar
OP committed
168
	acsum_sq.SetSize(8, vismtx.NCols());
OP's avatar
OP committed
169
	cxsum.SetSize(6, vismtx.NCols());
OP's avatar
OP committed
170 171
	cxsum_sq_rp.SetSize(6, vismtx.NCols());
	cxsum_sq_ip.SetSize(6, vismtx.NCols());
172
	
173
	tfmSX=wreader.getTotalNbWindows()/deltaIavg;
OP's avatar
OP committed
174 175 176 177
	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
178 179 180
	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
181 182 183
	//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) );
184
	cout << "and the 2x6 for squares of real & imaginary parts " << endl;
OP's avatar
OP committed
185 186
	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
187 188 189 190 191 192 193 194 195 196

	// 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
197
	acsum_sq = 0.;
OP's avatar
OP committed
198
	cxsum = complex<r_4>(0.,0.);
OP's avatar
OP committed
199 200
	cxsum_sq_rp = 0.;
	cxsum_sq_ip = 0.;
OP's avatar
OP committed
201 202 203 204
      }
      
      //   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
205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
      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
222 223 224
      I++;    // incrementing DeltaTime counter 
      
      if (I==deltaIavg) {
OP's avatar
OP committed
225
 
OP's avatar
OP committed
226 227 228
	//---- 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
229
	  TVector<r_4> vacsq = acsum_sq.Row(k);
OP's avatar
OP committed
230 231
	  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
232
	    TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
OP's avatar
OP committed
233 234
	    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
235 236
	      tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();

OP's avatar
OP committed
237
	    } 
OP's avatar
OP committed
238
	  }
OP's avatar
OP committed
239
	}  //----- end of loop over the 8 AutoCor
OP's avatar
OP committed
240

OP's avatar
OP committed
241
	//---- On s'occupe des 6 cross-correlations  1H-2H ... 3H-4H 
OP's avatar
OP committed
242 243
 	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
244 245
	  TVector<r_4>  vcxprsq = cxsum_sq_rp.Row(k);
	  TVector<r_4>  vcxpisq = cxsum_sq_ip.Row(k);
OP's avatar
OP committed
246 247
	  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
248 249
	    TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
	    TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
250 251
	    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
252 253
	      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
254 255 256
	    } 
	  }  
	}  //----- end of loop over the 6 Xcor 
OP's avatar
OP committed
257 258 259
	
	double tdif =  cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
	timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart);	// centre du bin 
260
	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs.ShiftSeconds (tdif ));
OP's avatar
OP committed
261 262 263 264 265 266 267 268 269 270 271 272
	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
273
    // --- Sauvegarde cartes temps-frequence en fits 
OP's avatar
OP committed
274 275 276
    //FitsABTWriter * fbtw = NULL;
    FitsInOutFile  * fos = NULL ;

OP's avatar
OP committed
277
    if (fitsoutfile.length()>=1){
278
      cout << " fitsoutfile :" <<fitsoutfile<<":"<< endl;
OP's avatar
OP committed
279 280 281
      fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
    }
    
OP's avatar
OP committed
282

OP's avatar
OP committed
283
    POutPersist potfm(outfile);
OP's avatar
OP committed
284 285
    char bufnam[10];
    int numkey =0;
OP's avatar
OP committed
286 287 288
    // --- 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
289
    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
290 291 292 293
    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;
OP's avatar
OP committed
294

OP's avatar
OP committed
295 296 297
      if (fos != NULL) {
	ext_names.push_back(tfm_names[k]);
	(*fos)<<  tfmap;
OP's avatar
OP committed
298 299 300 301 302 303
      }
      TArray< r_4 > & tfmapsq = vtfmac_sq[k];
      tfmapsq *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
      tfmapsq = tfmapsq - tfmap.MulElt(tfmap,tfmap) ; // tfmap ->tfmap*tfmap !! 
      potfm << PPFNameTag(tfmsq_names[k]) << tfmapsq;
      if (fos != NULL) {
OP's avatar
OP committed
304
	ext_names.push_back(tfmsq_names[k]);
OP's avatar
OP committed
305
	(*fos)<<  tfmapsq;
OP's avatar
OP committed
306
      }
OP's avatar
OP committed
307
      
OP's avatar
OP committed
308
    }
OP's avatar
OP committed
309

OP's avatar
OP committed
310 311
    // --- 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
312

OP's avatar
OP committed
313
    const char* tfmCC_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
OP's avatar
OP committed
314 315 316
    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
317 318
    for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
      TArray< complex<r_4> > & tfmap = vtfm[k];
OP's avatar
OP committed
319 320 321

      TArray< r_4 > & tfmap_sqpr = vtfm_rp_sq[k];
      TArray< r_4 > & tfmap_sqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
322
      tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
OP's avatar
OP committed
323 324 325 326 327 328 329
      potfm << PPFNameTag(tfmCC_names[k]) << tfmap;
      if (fos != NULL) {
	ext_names.push_back(string(tfmCC_names[k])+"_real");
	(*fos)<<  real(tfmap);
	ext_names.push_back(string(tfmCC_names[k])+"_imag");
	(*fos)<<  imag(tfmap);
      }
OP's avatar
OP committed
330 331 332 333 334 335 336 337 338 339 340
      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;
      potfm << PPFNameTag(vrtfmCC_names[k]) << tfmap_sqpr;
      potfm << PPFNameTag(vitfmCC_names[k]) << tfmap_sqpi;
OP's avatar
OP committed
341 342 343 344 345 346
      if (fos != NULL) {
	ext_names.push_back(vrtfmCC_names[k]);
	(*fos)<<  tfmap_sqpr;
	ext_names.push_back(vitfmCC_names[k]);
	(*fos)<<  tfmap_sqpi;
      }
OP's avatar
OP committed
347

OP's avatar
OP committed
348
    }
OP's avatar
OP committed
349

OP's avatar
OP committed
350
    potfm << PPFNameTag("TimeVec") << timevec ;
OP's avatar
OP committed
351
    potfm << PPFNameTag("RAVec") << ravec ;
OP's avatar
OP committed
352 353
    // --- FIN sauvegarde cartes temps-frequence 
    //    resu.Update();
OP's avatar
OP committed
354 355
    P4FreqBand myp4fre;
    TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
356
    
OP's avatar
OP committed
357
    double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
OP's avatar
OP committed
358
    for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
OP's avatar
OP committed
359
      avg_freqs(kf) = frbase ;
360
    
OP's avatar
OP committed
361
    potfm << PPFNameTag("FreqVec") << avg_freqs;
OP's avatar
OP committed
362 363 364 365 366 367 368 369 370 371 372
    if (fos != NULL) {
      ext_names.push_back("Frequences");
      (*fos)<< avg_freqs ;
      ext_names.push_back("RAs");
      (*fos)<<  ravec;
      ext_names.push_back("Times");
      (*fos)<<  timevec;
      cout << " number of objs in fits "<< ext_names.size() << endl;
      cout << ext_names << endl;
      delete(fos);
    }
373
    cout << " return code "<<rc<<endl;
OP's avatar
OP committed
374
    cout << resu;   // Update est fait lors du print
375

OP's avatar
OP committed
376 377 378 379 380 381 382 383 384 385 386 387 388 389
  }
  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; 
390
  }
OP's avatar
OP committed
391 392 393 394 395 396

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

}