rdvisip4.cc 9.29 KB
Newer Older
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  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 
40
#include "visip4reader.h"
41

42 43 44 45 46 47 48
// ---- 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);
 
49 50
//--------------------------- Fonctions de ce fichier   ------------------- 
int Usage(bool fgshort=true);
51
// int DecodeArgs(int narg, char* arg[]);
52 53

/* --Fonction-- */
54
int Usage(bool fgea)
55 56
{
  cout << " --- rdvisip4.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
57
  if (fgea) cout << " rdvisip4: Missing or wrong argument ! " << endl;
58
  cout << " Usage: rdvisip4 [-mr] InPathBAO5 InPathBAO6 Imin,Imax,step OutPPFFile GainPPFFile [PrintLev=1] \n"
59
       << "  -m: read also the mean Fourier coeff matrix \n"
60 61
       << "  -r: reorder frequence (for FFT firmeware)  \n"
       << "  -mr: do both  \n"
62 63
       << "  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"
64 65
       << "  OutPPFFile: Output PPF file name (with the average visibility matrix) \n"
       << "  GainPPFFile: PPF file name containing the computed gain g(nu) \n"<<endl;
66
  /*
67 68 69 70
  if (fgshort) {
    cout << " rdvisip4 -h for detailed instructions" << endl;
    return 1;
  }
71
  */
72 73 74 75 76 77 78 79
  return 1;
}

//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
80
  if (narg<6) return Usage(true);
81
  if ((narg>1)&&((strcmp(arg[1],"-m")==0)||(strcmp(arg[1],"-r")==0)||(strcmp(arg[1],"-mr")==0)) && (narg<7))  return Usage(false);
82 83 84

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

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

99
  cout << " rdvisip4/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
100
       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"  OutFile:"<<outfile
101 102 103
       << " PrtLev="<<prtlev
       <<endl;
  if(fgrdmean) cout << "Trying to read also the FFT mean coef matrix" <<endl;
104
  
105 106 107 108 109 110
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    ResourceUsage resu;
111 112 113 114
    vector<string> paths; 
    paths.push_back(path5); 
    paths.push_back(path6);
 
115 116 117
    int nfiles=(Imax-Imin+1)/Istep;
    int nmod=nfiles/20;
    if (nmod<1) nmod=1;
118
    VisiP4Reader vreader(paths, Imin,Imax,Istep,fgreorderfreq);
119
    vreader.setPrintLevel(prtlev);
120
    bool fgok=true;
121 122 123
    //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;
124 125 126 127
    TimeStamp dateobs;
    double mttag;
    int cnt=0;
    while (fgok) {
128 129
      if(fgrdmean) {
	fgok=vreader.ReadNext(vismtx, meanmtx, dateobs, mttag);
130
	vismtx_center = vreader.SubtractOffset(vismtx, meanmtx);
131 132 133
      } else {
	fgok=vreader.ReadNext(vismtx, dateobs, mttag);
      }
134 135 136
      if (fgok)  {
	if (cnt==0)  vismtx_mean=vismtx;
	else vismtx_mean+=vismtx;
137 138 139
	if(fgrdmean) {
	  if (cnt==0)  meanmtx_mean = meanmtx;
	  else meanmtx_mean += meanmtx;
140 141
	  if (cnt==0)  vismtx_center_mean = vismtx_center;
	  else  vismtx_center_mean += vismtx_center;
142
	}
143
	cnt++;
144
	if (cnt%nmod==0)  cout<<"rdvisip4/Info file read count="<<cnt<<" / nfiles="<<nfiles<<endl;
145 146 147 148 149
      } 
    }
    cout << " rdvisip4/Info: count="<<cnt<<" visimtx read "<<endl;
    if (cnt>0)  { 
      vismtx_mean /= complex<r_4>((r_4)cnt, (r_4)0.);
150 151
      cout << " rdvisip4/Info: computing gains ..."<<endl;
      computeSaveGain(vismtx_mean, gainfile);
152 153 154
      cout<<"  rdvisip4/Info: Saving vismtx_mean to PPF file "<<outfile<<endl;
      POutPersist po(outfile);
      po<<vismtx_mean;
155 156 157 158
      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;
159 160 161 162

	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;
163 164
      }

165 166 167
    } 
    //    resu.Update();
    cout << resu;   // Update est fait lors du print
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188
  }
  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;

}

189 190 191 192 193 194 195 196 197 198 199 200 201 202 203
//---------------------  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);
204
  
205 206 207 208 209 210 211
  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;
}
212

213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
// ----- 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;
255
 return gain;
256 257 258 259 260 261 262 263 264
}

/*--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();
265 266
    double fcg=1000./gn;    // we decide arbitray to renormalize the gain to have the sum[gains(i)]=1000 
    gain.Row(r) *= fcg;
267 268
    sprintf(buff,"GNORM%d",(int)r);
    gain.Info()[buff]=gn;
269
    gnorm(r)=gn/1000.;
270 271 272 273 274 275 276 277 278 279
    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;
}