Commit 9946addd authored by Reza Ansari's avatar Reza Ansari
Browse files

1) Modification de la classe VisiP4WindowReader pour qu'elle prenne en charge

les corrections de gain (g(nu), G(t)) et de phase
2) Constructeur specifique VisiP4WindowReader(P4AnaParams &)
3) Adaptation/simplification du programme visavg.cc pour utiliser ces conctionalites
de VisiP4WindowReader au lieu de gerer ces corrections ds visiavg.cc
                                                        Reza  17/01/2018
parent 4e5e744a
......@@ -197,8 +197,8 @@ std::vector<sa_size_t> P4AVisiNumEncoder::getAllVCrossCor()
//---------------------------------------------------------------------
P4AnaParams::P4AnaParams()
: inpath5_("bao5/"), inpath6_("bao6/"), Imin_(0), Imax_(10), Istep_(1),
fgserall_(false), fgtmsel_(false), inwsz_(8), fgreorderfreq_(false),
gain_gnu_file_("gain.ppf"), gain_var_file_(""), gain_blind_filt_file_(""), phase_corr_param_file_(""),
fgserall_(false), fgtmsel_(false), inwsz_(1), fgreorderfreq_(false),
gain_gnu_file_(""), gain_var_file_(""), gain_blind_filt_file_(""), phase_corr_param_file_(""),
outfile_("p4a.ppf"), fitsoutfile_(""), prtlev_(0), prtmodulo_(1),
Nmean_(10), fgTFM_(false), TFMtimebin_(30), TFMfreqbin_(4), fgdtable_(false)
{
......
......@@ -37,11 +37,7 @@
#include "timing.h"
// include lecteur de fichiers visibilites
#include "p4autils.h"
#include "visp4winreader.h"
#include "p4gnugain.h"
#include "p4gvcor.h"
#include "p4phcor.h"
#include "fitsioserver.h"
#include "fiosinit.h"
......@@ -75,7 +71,6 @@ int main(int narg, const char* arg[])
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;
......@@ -85,7 +80,7 @@ int main(int narg, const char* arg[])
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"
<<" DeltaIAvg="<<deltaIavg<<"\n"
<<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
if (!FgTFM) {
......@@ -105,72 +100,20 @@ int main(int narg, const char* arg[])
cout << "KVCXHH["<<k<<"]="<<KVCXHH[k]<<" ->"<<visiencod.Convert2VisiName(KVCXHH[k])<<endl;
}
string gvfile = params.gain_var_file_;
string fltfile = params.gain_blind_filt_file_;
P4gvCor * p4gvf = NULL ;
cout <<" gvfile "<< gvfile << endl;
cout <<"fltfile "<< fltfile << endl;
if ( gvfile.length()>1 && fltfile.length() > 1 ) p4gvf = new P4gvCor( gvfile, fltfile ) ;
// ---
HiStatsInitiator _inia;
int rc = 0;
// Gain correction class
P4gnuGain * p4g = NULL ;
try {
p4g = new P4gnuGain(params.gain_gnu_file_);
}
catch (PException& exc) {
cerr << " visiavg.cc catched PException " << exc.Msg() << endl;
cerr << " .... Continuing execution ... " << endl;
p4g = NULL ;
}
P4PhaseCor * p4phc = NULL;
if (params.phase_corr_param_file_.length() > 0) {
try {
p4phc = new P4PhaseCor(params.phase_corr_param_file_);
}
catch (PException& exc) {
cerr << " visiavg.cc catched PException when reading Phase Correction parameter file: " << params.phase_corr_param_file_ << endl;
cerr << " Exception message: " << exc.Msg() << endl;
cerr << " .... Continuing execution ... " << endl;
p4phc = NULL;
}
}
try {
ResourceUsage resu;
// setting up input visi reader
vector<string> paths;
paths.push_back(params.inpath5_);
paths.push_back(params.inpath6_);
VisiP4WindowReader wreader(params);
VisiP4WindowReader wreader(paths, params.inwsz_);
wreader.getReader().setFreqReordering(params.fgreorderfreq_);
if (!params.fgserall_ && !params.fgtmsel_) {
cout << " vreader.SelectSerialNum(Imin="<<Imin<<" ,Imax="<<Imax<<" ,Istep="<<Istep<<")"<<endl;
wreader.getReader().SelectSerialNum(Imin,Imax,Istep);
}
else if (params.fgserall_) {
cout << " vreader.SelectAll() ... " << endl;
wreader.getReader().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;
wreader.getReader().SelectTimeFrame(tustart, tuend, Istep);
}
wreader.getReader().setPrintLevel(prtlev);
Imin = wreader.getReader().getSerialFirst();
Imax = wreader.getReader().getSerialLast();
Istep = wreader.getReader().getSerialStep();
long Imin = wreader.getReader().getSerialFirst();
long Imax = wreader.getReader().getSerialLast();
long Istep = wreader.getReader().getSerialStep();
cout << "visiavg/Info: processing visibility matrix serial/sequence number range "
<<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
cout << " WindowSize="<<wreader.getWindowSize()<<" -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;
......@@ -215,26 +158,15 @@ int main(int narg, const char* arg[])
fgok = wreader.Shift();
if (!fgok) break;
vismtx = wreader.getAverageVisMtx();
cfdate = wreader.getAverageTimeTU();
// Apply gain g(nu)
if (p4g != NULL )
p4g->applyGain(vismtx);
// Apply time dependent (/ temperature) gain G(t)
if (p4gvf != NULL)
p4gvf->applyGain(vismtx,cfdate);
// Apply phase correction
if (p4phc != NULL)
p4phc->applyPhase(vismtx);
vismtx = wreader.getAverageVisMtx(cfdate);
if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
acsum.SetSize(8, vismtx.NCols());
acsum_sq.SetSize(8, vismtx.NCols());
cxsum.SetSize(6, vismtx.NCols());
cxsum_sq_rp.SetSize(6, vismtx.NCols());
cxsum_sq_ip.SetSize(6, vismtx.NCols());
tfmSX=wreader.getTotalNbWindows()/deltaIavg;
tfmSY=vismtx.NCols()/TFMfbin;
//allocating 8 Auto-Corr time-frequency maps
......@@ -436,9 +368,7 @@ int main(int narg, const char* arg[])
}
cout << resu; // Update est fait lors du print
if (p4g != NULL) delete(p4g);
}
catch (PException& exc) {
cerr << " visiavg.cc catched PException " << exc.Msg() << endl;
......
......@@ -9,17 +9,69 @@
/* --Methode-- */
VisiP4WindowReader::VisiP4WindowReader(vector<string> const& vpath, size_t wsz)
: wsz_(wsz) , vtu_(wsz)
: wsz_(wsz) , p4vreaderp_(NULL), vtu_(wsz),
p4gnu_(NULL), p4gvf_(NULL), p4phc_(NULL),
fgown_p4gnu_(false), fgown_p4gvf_(false), fgown_p4phc_(false)
{
cout << " VisiP4WindowReader - Creating Reader with window size= " << wsz_ << endl;
BaseArray::SetDefaultMemoryMapping(BaseArray::CMemoryMapping);
p4vreaderp_ = VisiP4ReaderBase::getReader(vpath);
}
/* --Methode-- */
VisiP4WindowReader::VisiP4WindowReader(P4AnaParams & params)
: wsz_(params.inwsz_) , p4vreaderp_(NULL), vtu_(params.inwsz_),
p4gnu_(NULL), p4gvf_(NULL), p4phc_(NULL),
fgown_p4gnu_(false), fgown_p4gvf_(false), fgown_p4phc_(false)
{
cout << " VisiP4WindowReader(P4AnaParams & params) - Creating Reader with window size= " << wsz_ << endl;
BaseArray::SetDefaultMemoryMapping(BaseArray::CMemoryMapping);
vector<string> paths;
paths.push_back(params.inpath5_);
paths.push_back(params.inpath6_);
p4vreaderp_ = VisiP4ReaderBase::getReader(paths);
getReader().setFreqReordering(params.fgreorderfreq_);
if (!params.fgserall_ && !params.fgtmsel_) {
cout << "VisiP4WindowReader: SelectSerialNum(Imin="<<params.Imin_<<" ,Imax="<<params.Imax_
<<" ,Istep="<<params.Istep_<<")"<<endl;
getReader().SelectSerialNum(params.Imin_,params.Imax_,params.Istep_);
}
else if (params.fgserall_) {
cout << "VisiP4WindowReader: SelectAll() ... " << endl;
getReader().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 << "VisiP4WindowReader: SelectTimeFrame(TUStart="<<tustart.ToString()<<" ,TUEnd="<<tuend.ToString()
<<" ,Istep="<<params.Istep_<<")"<<endl;
getReader().SelectTimeFrame(tustart, tuend, params.Istep_);
}
getReader().setPrintLevel(params.prtlev_);
if (params.gain_gnu_file_.length()>0) {
p4gnu_= new P4gnuGain(params.gain_gnu_file_);
fgown_p4gnu_=true;
}
if ( params.gain_var_file_.length()>1 && params.gain_blind_filt_file_.length() > 1 ) {
p4gvf_ = new P4gvCor(params.gain_var_file_, params.gain_blind_filt_file_);
fgown_p4gvf_=true;
}
if (params.phase_corr_param_file_.length() > 0) {
p4phc_ = new P4PhaseCor(params.phase_corr_param_file_);
fgown_p4phc_=true;
}
}
/* --Methode-- */
VisiP4WindowReader::~VisiP4WindowReader()
{
delete p4vreaderp_;
delete p4vreaderp_;
if (fgown_p4gnu_ && p4gnu_) delete p4gnu_;
if (fgown_p4gvf_ && p4gvf_) delete p4gvf_;
if (fgown_p4phc_ && p4phc_) delete p4phc_;
}
/* --Methode-- */
......@@ -28,6 +80,28 @@ size_t VisiP4WindowReader::getTotalNbWindows()
return ((size_t)((p4vreaderp_->getSerialLast()-p4vreaderp_->getSerialFirst())/(p4vreaderp_->getSerialStep()*wsz_)));
}
/* --Methode-- */
void VisiP4WindowReader::setGainCorrector(P4gnuGain& gnuc)
{
if (fgown_p4gnu_ && p4gnu_) delete p4gnu_;
p4gnu_=&gnuc; fgown_p4gnu_=false;
}
/* --Methode-- */
void VisiP4WindowReader::setGainCorrector(P4gvCor& gvc)
{
if (fgown_p4gvf_ && p4gvf_) delete p4gvf_;
p4gvf_=&gvc; fgown_p4gvf_=false;
}
/* --Methode-- */
void VisiP4WindowReader::setPhaseCorrector(P4PhaseCor& phc)
{
if (fgown_p4phc_ && p4phc_) delete p4phc_;
p4phc_=&phc; fgown_p4phc_=false;
}
/* --Methode-- */
bool VisiP4WindowReader::Shift()
{
......@@ -55,20 +129,31 @@ bool VisiP4WindowReader::Shift()
}
/* --Methode-- */
TMatrix< complex<r_4> > VisiP4WindowReader::getVisMtx(size_t k)
TMatrix< complex<r_4> > VisiP4WindowReader::getVisMtx(size_t k, bool fgcor)
{
if (k >= wsz_) throw RangeCheckError("VisiP4WindowReader::getVisMtx(k) k >= WindowSize");
TMatrix< complex<r_4> > vismtx = visarr_(Range::all(), Range::all(), Range((sa_size_t)k));
return vismtx;
if (!fgcor) return vismtx;
TMatrix< complex<r_4> > rvmx(vismtx, false); // copy the visibility matrix to avoid changing the array
if (p4gnu_) p4gnu_->applyGain(rvmx);
if (p4gvf_) p4gvf_->applyGain(rvmx,vtu_[k]);
if (p4phc_) p4phc_->applyPhase(rvmx);
return rvmx;
}
/* --Methode-- */
TMatrix< complex<r_4> > VisiP4WindowReader::getAverageVisMtx()
TMatrix< complex<r_4> > VisiP4WindowReader::getAverageVisMtx(TimeStamp& tu, bool fgcor)
{
TArray< complex<r_4> > avar(visarr_.SizeX(), visarr_.SizeY());
for(size_t k=0; k<wsz_; k++) avar += visarr_(Range::all(), Range::all(), Range((sa_size_t)k));
avar /= complex<r_4>((float)wsz_, 0.);
return TMatrix< complex<r_4> >(avar);
TMatrix< complex<r_4> > rvmx(avar);
tu = getAverageTimeTU();
if (!fgcor) return rvmx;
if (p4gnu_) p4gnu_->applyGain(rvmx);
if (p4gvf_) p4gvf_->applyGain(rvmx, tu);
if (p4phc_) p4phc_->applyPhase(rvmx);
return rvmx;
}
/* --Methode-- */
......
......@@ -17,6 +17,12 @@
#include "pexceptions.h"
#include "visip4reader.h"
#include "p4autils.h"
//-- correction de gains et de phase
#include "p4gnugain.h"
#include "p4gvcor.h"
#include "p4phcor.h"
using namespace std;
using namespace SOPHYA;
......@@ -36,7 +42,9 @@ public:
time window size \b wsz (specified as the number of visibility matrices)
*/
VisiP4WindowReader(vector<string> const& vpath, size_t wsz=8);
//! \brief Constructor all parameters defined through the provided P4AnaParams object
VisiP4WindowReader(P4AnaParams & params);
~VisiP4WindowReader();
//! Return the VisiP4ReaderBase object
......@@ -49,13 +57,24 @@ public:
//! return total number of windows to be read
size_t getTotalNbWindows();
//! \brief Set/change the frequency dependent gain g(nu) correction P4gnuGain object.
void setGainCorrector(P4gnuGain& gnuc);
//! \brief Set/change the time (or temperature) dependent gain G(t) correction P4gvCor object.
void setGainCorrector(P4gvCor& gvc);
//! \brief Set/change the phase correction P4PhaseCor object.
void setPhaseCorrector(P4PhaseCor& phc);
//! Reads the visibility matrices corresponding to the next time window - return true if OK, false otherwise
bool Shift();
//! Returns the visibiliy matrix with the index k in the window
TMatrix< complex<r_4> > getVisMtx(size_t k);
//! Returns the average visibiliy matrix of the window
TMatrix< complex<r_4> > getAverageVisMtx();
//! Returns the visibiliy matrix with the index k in the window, gain and phase corrected if \b fgcor is true
TMatrix< complex<r_4> > getVisMtx(size_t k, bool fgcor=true);
/*! \brief Returns the average visibiliy matrix of the window, gain and phase corrected
if \b fgcor is true and the corresponding TU time (getAverageTimeTU()) */
TMatrix< complex<r_4> > getAverageVisMtx(TimeStamp& tu, bool fgcor=true);
//! Returns the average visibiliy matrix of the window, gain and phase corrected if \b fgcor is true
inline TMatrix< complex<r_4> > getAverageVisMtx(bool fgcor=true)
{ TimeStamp ts; return getAverageVisMtx(ts, fgcor); }
//! Return the timestamp (TU time) of the visibility matrix with index k in the window
inline TimeStamp getTimeTU(size_t k) { return vtu_[k]; }
//! Return the average TU time of the window
......@@ -66,6 +85,12 @@ protected:
VisiP4ReaderBase * p4vreaderp_;
TArray< complex<r_4> > visarr_;
vector<TimeStamp> vtu_;
//---- corrections de gains et de phase
P4gnuGain * p4gnu_;
P4gvCor * p4gvf_;
P4PhaseCor * p4phc_;
bool fgown_p4gnu_, fgown_p4gvf_, fgown_p4phc_; // true->the corresponding pointer is owned by the VisiP4WindowReader object
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment