source: Sophya/trunk/AddOn/TAcq/brviscalc.cc@ 4012

Last change on this file since 4012 was 4012, checked in by ansari, 14 years ago

Codage du mode de calcul de visibilites par intervalle de temps (BRVisibiliyCalculator) et propagation des modifs (ajouts de parametres) ds mfacq.cc et vismfib.cc, Reza 02/08/2011

File size: 26.0 KB
RevLine 
[3872]1//----------------------------------------------------------------
2// Projet BAORadio - (C) LAL/IRFU 2008-2010
3// Classes de threads de calcul de visibilites pour BAORadio
4//----------------------------------------------------------------
5
6#include <stdlib.h>
7#include <string.h>
8#include <unistd.h>
9#include <fstream>
10#include <signal.h>
11
12#include "pexceptions.h"
13#include "fioarr.h"
[3956]14#include "fitsarrhand.h"
[3872]15#include "matharr.h"
[3956]16#include "arrctcast.h"
[3872]17#include "timestamp.h"
18#include "ctimer.h"
19
20#include "brviscalc.h"
21
22
23//---------------------------------------------------------------------
24// Classe de traitement - calcul de visibilite pour n fibres
25//---------------------------------------------------------------------
26
27/* --Methode-- */
28BRVisibilityCalculator::BRVisibilityCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, size_t nthr)
29 : BRBaseProcessor(memgr), paralex_(*this, nthr), nparthr_(nthr),
[4012]30 outpath_(outpath), nbcalc_(1), calcid_(0),
[3956]31 vpdata_(2*memgr.NbFibres()), vpdatar_(2*memgr.NbFibres())
[3872]32 // , dtfos_(outpath+"visdt.fits", Fits_Create), visdt_(dtfos_, 1024, true);
33{
[3956]34 SetFFTData();
[4012]35 SetNPaqIntervalMode(nmean);
[3872]36 DefineRank(1,0);
[3956]37 SetPPFOutput();
[3872]38
[4012]39 npaqcumul_=0;
40 moyfc_=moytt_=0.;
41 first_fc_=first_tt_=0;
42
[3872]43 uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();
[3909]44 chanids_.SetSize(2*memgr_.NbFibres());
45 chanpairnumall_.SetSize(maxnpairs);
46 chanpairsall_.SetSize(maxnpairs,2);
[3956]47 for(size_t i=0; i<2*memgr_.NbFibres(); i++) {
48 vpdata_[i]=NULL; vpdatar_[i]=NULL;
49 }
[3957]50 // Pas necessaire - appele par DefineRank- SelectPairs();
[3872]51
52 // visdt_.AddFloatColumn("mfc");
53 visdt_.AddDoubleColumn("mfc");
54 visdt_.AddDoubleColumn("mtt");
55 visdt_.AddIntegerColumn("jfreq");
56 visdt_.AddIntegerColumn("numch");
57 visdt_.AddFloatColumn("vre");
58 visdt_.AddFloatColumn("vim");
[3920]59 ActivateVisDTable();
[3872]60
61 if (nmean_ < 1) nmean_=memgr_.NbPaquets();
62 if (nmean_ < 1) nmean_=1;
63
64 totnbpaq_=0;
65 numfile_=0;
66 nb_flop_=0.;
67
68 fgallfibok=NULL;
69 fgcktt_=false;
70 setNameId("viscalc", 0);
71}
72
73/* --Methode-- */
74BRVisibilityCalculator::~BRVisibilityCalculator()
75{
[3876]76 if (totnbpaq_<1) return;
[3920]77 if (fgvisdt_) {
78 cout << " ~BRVisibilityCalculator - Visibility Datatable : " << endl;
79 cout << visdt_;
[3876]80
[3920]81 string filename;
82 filename = outpath_+"visdt.ppf";
83 if (nbcalc_>1) {
84 char sbuff[32];
85 sprintf(sbuff,"visdt_%d.ppf",(int)calcid_);
86 filename = outpath_+sbuff;
87 }
88 POutPersist po(filename);
89 po << visdt_;
[3872]90 }
[3920]91 if (calcid_ == 0) {
[3872]92 if (fgcktt_) {
93 cout << " BRVisibilityCalculator - Check TimeTag done: TotNPaqProc= " << totnbpaq_ << endl;
94 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
95 cout << " BRTTCheck-Fiber[" << fib << "] NBadTT=" << vbadtt_[fib] << " NDiffTT>5="
96 << vndiff5tt_[fib] << " NotSameTT=" << vnsamett_[fib] << endl;
97 }
98 POutPersist pott(outpath_+"ttfcmtx.ppf");
99 pott << PPFNameTag("FC") << fcmtx_;
100 pott << PPFNameTag("TT") << ttmtx_;
101 }
102 }
103}
104
105/* --Methode-- */
[3915]106void BRVisibilityCalculator::DefineRank(uint_4 nbc, uint_4 cid, uint_4 pair1, uint_4 nbpairs, bool fgpimp)
[3872]107{
108 if ((nbc>6)||(cid>=nbc))
109 throw ParmError("BRVisibilityCalculator::DefineRank() NbCalc > 6 !");
[3915]110 uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();
111 if ((pair1>=maxnpairs)||(nbpairs<1))
112 throw ParmError("BRVisibilityCalculator::DefineRank() (pair1>=maxnpairs)||(nbpairs<1) !");
113 maxnpairs-=pair1;
114 if (maxnpairs>nbpairs) maxnpairs=nbpairs;
115
[3872]116 nbcalc_=nbc;
117 calcid_=cid;
118 if (nbcalc_>1) {
[3915]119 uint_4 npairspth=maxnpairs/nbcalc_; // nb de paires par thread
[3872]120 if (calcid_==(nbcalc_-1))
[3915]121 SelectPairs(calcid_*npairspth+pair1, maxnpairs-calcid_*npairspth, fgpimp);
[3872]122 else
[3915]123 SelectPairs(calcid_*npairspth+pair1, npairspth, fgpimp);
[3956]124 MemZaction mmzas[10]={MemZA_ProcA,MemZA_ProcB,MemZA_ProcC,MemZA_ProcD,MemZA_ProcE,MemZA_ProcF,
125 MemZA_ProcG,MemZA_ProcH,MemZA_ProcI,MemZA_ProcJ};
[3872]126 SetMemZAction(mmzas[calcid_]);
127 setNameId("viscalc_grp", calcid_);
128 }
[3915]129 else SelectPairs(pair1,maxnpairs, fgpimp);
[3872]130 return ;
131}
132
133/* --Methode-- */
[3915]134uint_4 BRVisibilityCalculator::SelectPairs(uint_4 pair1, uint_4 nbpairs, bool fgpimp)
[3872]135{
136 BRPaquet paq(memgr_.PaqSize());
137 uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();
138
139 if (pair1 >= maxnpairs) pair1=maxnpairs-1;
140 if (nbpairs > maxnpairs-pair1) nbpairs=maxnpairs-pair1;
[3915]141
[3872]142 pairst_=pair1;
143 nbpairs_=nbpairs;
[3915]144 fgpimp_=fgpimp;
[3909]145
[3915]146 uint_4 nbvisicomp=0;
147 sa_size_t kpair=0;
148 for(size_t i=0; i<vpdata_.size(); i++) {
149 for(size_t j=i; j<vpdata_.size(); j++) {
150 kpair++;
151 if (kpair<(pairst_+1)) continue;
152 if (kpair>=(pairst_+nbpairs_+1)) break;
153 if (fgpimp_&&(i!=j)&&((i+j)%2==0)) continue; // calcul des visib avec numero pair-impair + autocorrel
154 nbvisicomp++;
155 }
156 }
[3909]157
[3915]158 vismtx_.SetSize(nbvisicomp, paq.DataSize()/4);
159
160 chanpairnum_.SetSize(nbvisicomp);
161 chanpairs_.SetSize(nbvisicomp,2);
162
163 return nbvisicomp;
[3872]164}
165
166/* --Methode-- */
167int BRVisibilityCalculator::SelectFreqBinning(uint_4 freq1, uint_4 freq2, uint_4 nbfreq)
168{
[3993]169 if (freq2<freq1) {
170 freq1=0; freq2=vismtx_.NCols()-1; nbfreq=1;
171 }
[3872]172 jf1_=freq1; jf2_=freq2;
173 if ((jf1_<1)||(jf1_>=vismtx_.NCols())) jf1_=1;
174 if ((jf2_<1)||(jf2_>=vismtx_.NCols())||(jf2_<jf1_)) jf2_=vismtx_.NCols()-1;
[3876]175 if (nbfreq<1) nbfreq=(jf2_-jf1_);
[3872]176 djf_=(jf2_-jf1_)/nbfreq;
[3876]177 if (djf_<1) djf_=1;
[3872]178 cout << " BRVisibilityCalculator::SelectFreqBinning/Info JF1=" << jf1_
179 << " JF2=" << jf2_ << " DJF=" << djf_ << endl;
180
181}
182
[3909]183/* --Methode-- */
184void BRVisibilityCalculator::UpdateChanIds()
185{
[3872]186
[3909]187 for(size_t i=0; i<memgr_.NbFibres(); i++) {
188 chanids_(2*i)=memgr_.FiberId(i)*2-1;
189 chanids_(2*i+1)=memgr_.FiberId(i)*2;
190 }
191 sa_size_t k=0; // numero de ligne dans la matrice des visibilites
192 for(size_t i=0; i<vpdata_.size(); i++) {
193 for(size_t j=i; j<vpdata_.size(); j++) {
194 chanpairnumall_(k)=chanids_(i)*CHANPAIRCONVFAC+chanids_(j);
195 chanpairsall_(k,0)=chanids_(i); chanpairsall_(k,1)=chanids_(j); k++;
196 }
197 }
198 sa_size_t kpair=0;
199 k=0; // numero de ligne dans la matrice des visibilites
200 for(size_t i=0; i<vpdata_.size(); i++) {
201 for(size_t j=i; j<vpdata_.size(); j++) {
202 kpair++;
203 if (kpair<(pairst_+1)) continue;
204 if (kpair>=(pairst_+nbpairs_+1)) break;
[3915]205 if (fgpimp_&&(i!=j)&&((i+j)%2==0)) continue; // calcul des visib avec numero pair-impair + autocorrel
[3909]206 chanpairnum_(k)=chanids_(i)*CHANPAIRCONVFAC+chanids_(j);
207 chanpairs_(k,0)=chanids_(i); chanpairs_(k,1)=chanids_(j); k++;
208 }
209 }
210
211 string filename;
[3956]212 filename = (outpath_+"chanum.")+OutFileExtension();
[3909]213 if (nbcalc_>1) {
214 char sbuff[32];
[3956]215 sprintf(sbuff,"chanum_%d.%s",(int)calcid_,OutFileExtension());
[3909]216 filename = outpath_+sbuff;
217 }
[3956]218 if (fgfitsout_) { // Ecriture au format FITS
219 FitsInOutFile foc(filename, FitsInOutFile::Fits_Create);
220 foc.SetNextExtensionName("chanids");
221 foc << chanids_;
222 foc.SetNextExtensionName("chanpairs");
223 foc << chanpairs_;
224 foc.SetNextExtensionName("chanpairnum");
225 foc << chanpairnum_;
226 foc.SetNextExtensionName("chanpairsall");
227 foc << chanpairsall_;
228 foc.SetNextExtensionName("chanpairnumall");
229 foc << chanpairnumall_;
230 }
231 else { // Format PPF
232 POutPersist poc(filename);
233 poc << PPFNameTag("chanids") << chanids_;
234 poc << PPFNameTag("chanpairs") << chanpairs_;
235 poc << PPFNameTag("chanpairnum") << chanpairnum_;
236 poc << PPFNameTag("chanpairsall") << chanpairsall_;
237 poc << PPFNameTag("chanpairnumall") << chanpairnumall_;
238 }
239 cout << "BRVisibilityCalculator[" << calcid_ << "]::UpdateChanIds() Channel Ids/Pairs saved to file "
[3910]240 << filename << endl;
[3915]241 cout << " ... NbVisib=NbChanPairs=" << chanpairs_.NRows() << " ChannelPairs= " ;
[3910]242 for(sa_size_t ir=0; ir<chanpairs_.NRows(); ir++) {
[3911]243 if (ir%10==0) cout << endl;
[3910]244 cout << "(" << chanpairs_(ir,0) << "," << chanpairs_(ir,1) << ") ";
245 }
246 cout << endl;
[3909]247 return;
248}
249
[3872]250/* --Methode-- */
251int BRVisibilityCalculator::ActivateTimeTagCheck(uint_8 maxnpaq)
252{
253 mindeltatt_=memgr_.PaqSize()/2;
254 if (mindeltatt_<1) mindeltatt_=1;
255 fcmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
256 ttmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
257 vlasttt_.resize(memgr_.NbFibres(), 0);
258 vbadtt_.resize(memgr_.NbFibres(), 0);
259 vnsamett_.resize(memgr_.NbFibres(), 0);
260 vndiff5tt_.resize(memgr_.NbFibres(), 0);
261
262 fgcktt_=true;
263 cout << " BRVisibilityCalculator::ActivateTimeTagCheck() - TT/Fc matrix NCols=" << maxnpaq
264 << " MinDeltaTT=" << mindeltatt_ << endl;
265
266 return 0;
267}
268
269
270/* --Methode-- */
271void BRVisibilityCalculator::run()
272{
273 cout << " BRVisibilityCalculator[" << calcid_ << "/" << nbcalc_
274 << "]::run() - Starting " << " NFibers=" << memgr_.NbFibres()
275 << " NChan=" << 2*memgr_.NbFibres() << " NPairs=" << nbpairs_ << " First:" << pairst_ << endl;
[4012]276
277 if (fgnpaq_time)
278 cout << " BRVisibilityCalculator::run()/Info NPaq=" << nmean_ << " Interval mode for visibility mtrix save" << endl;
279 else
280 cout << " BRVisibilityCalculator::run()/Info Time Tnterval (dt=" << dtimemean_ << ") for visibility mtrix save" << endl;
281
[3872]282 if (nparthr_ < 2) return BRBaseProcessor::run();
283 // Execution multithread parallele
284 setRC(1);
285 int rc=0;
286 try {
287 if ((nmean_%memgr_.NbPaquets())!=0) {
288 uint_4 mnmean = (nmean_/memgr_.NbPaquets()+1)*memgr_.NbPaquets();
289 cout << " BRVisibilityCalculator::run()/Info changing nmean=" << nmean_ << " to multiple of"
290 << " memgr_.NbPaquets() -> " << mnmean << endl;
291 nmean_=mnmean;
292 }
293 paralex_.SetParallelTask(*this);
294 cout << " BRVisibilityCalculator::run()/Info : starting ParallelExecutor with nThreads="
295 << paralex_.nThreads() << " ... " << endl;
296 paralex_.start();
297
298 fgallfibok = new bool[memgr_.NbPaquets()];
299
300 size_t paqsz=memgr_.PaqSize();
[3956]301 size_t procpaqsz=memgr_.ProcPaqSize();
[3872]302 bool fgrun=true;
303 while (fgrun) {
304 if (stop_) break;
305 if (memgr_.GetRunState() == MemZR_Stopped) break;
306 int mid = memgr_.FindMemZoneId(mmact_); // (MemZA_ProcA);
307 // Byte* buffg = memgr_.GetMemZone(mid);
308 // if (buffg == NULL) {
309 if (mid < 0) {
310 cout << "BRVisibilityCalculator[" << calcid_ << "]::run()/ERROR FindMemZoneId("
311 << (int)mmact_ << ") ->" << mid << ") -> NULL" << endl;
312 setRC(7); fgrun=false;
313 break;
314 }
[3909]315 cts_=memgr_.GetAuxData(mid)->FillTime(); // get associated date/time (DATEOBS)
[3910]316 if (totnbpaq_==0) UpdateChanIds(); // Appele ici pour etre sur que le thread de remplissage a mis l'info a jour.
[3909]317
[3872]318 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
319 fbuff_[fib] = memgr_.GetMemZone(mid,fib);
320 if (fbuff_[fib] == NULL) { // cela ne devrait pas arriver
321 cout << "BRVisibilityCalculator[" << calcid_ << "]::run()/ERROR memgr.GetMemZone(" << mid << "," << fib << ") -> NULL" << endl;
322 setRC(9); fgrun=false;
323 break;
324 }
[3956]325 if ((procpaqsz>0)&&((fprocbuff_[fib]=memgr_.GetProcMemZone(mid,fib))==NULL)) { // cela ne devrait pas arriver non plus
326 cout << "BRVisibilityCalculator[" << calcid_ << "]::run()/ERROR memgr.GetProcMemZone("
327 << mid << "," << fib << ") -> NULL" << endl;
328 setRC(9); fgrun=false;
329 break;
330 }
[3872]331 }
332
[4012]333 if ( CheckInterval4VisMtxSave() ) {
334 if ((prtlev_>0)&&(numfile_%prtmodulo_==0)) {
335 cout << numfile_ << "-BRVisCalc[" << calcid_ << "/" << nbcalc_ << "]::run() NPaqProc="
336 << totnbpaq_ << " TotMegaFLOP=" << (uint_8)TotNbMegaFLOP() << " -> VisibMtx FileNum=" << numfile_ << endl;
[3872]337 }
[4012]338 SaveVisibilityMatrix();
[3872]339 }
340
341 for(size_t jp=0; jp<memgr_.NbPaquets(); jp++) { // boucle sur les paquets d'une zone
342 fgallfibok[jp]=fgokallfibers_=true;
343 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
344 vpaq_[fib].Set(fbuff_[fib]+jp*paqsz);
345 vfgok_[fib] = vpchk_[fib].Check(vpaq_[fib],curfc_[fib]);
346 if (!vfgok_[fib]) fgallfibok[jp]=fgokallfibers_=false;
[3967]347 //Pas utile if (procpaqsz>0) vprocpaq_[fib] = fprocbuff_[fib]+jp*procpaqsz;
[3872]348 }
349 if (fgokallfibers_) {
350 if (totprocnpaq_==0) {
351 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
352 fcfirst_[fib]=curfc_[fib];
353 ttfirst_[fib]=vpaq_[fib].TimeTag();
354 }
355 }
356 totprocnpaq_++;
357 moyfc_ += curfc_[0];
358 moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
359 if ((fgcktt_)&&(calcid_==0)) CheckTimeTag();
[4012]360 if (npaqcumul_==0) {
[3895]361 first_fc_=curfc_[0];
[4012]362 first_tt_= vpaq_[0].TimeTag();
[3895]363 }
[4012]364 npaqcumul_++;
[3872]365 totnbpaq_++;
366 }
367 } // Fin de boucle sur les paquets
[3895]368
[3872]369 // Execution parallele pour calcul des visibilites par bandes de frequence
370 int rcpex=paralex_.execute();
371 if (rcpex!=0) cout << " BRVisibilityCalculator[" << calcid_ << "]::run() / Error Rc[paralex_.execute()]=" << rcpex << endl;
372
373 memgr_.FreeMemZone(mid, mmsta_); // (MemZS_ProcA);
374 } // Fin de boucle sur les zones a traiter
375 //------------------------------------
376 cout << " --------- END BRVisibilityCalculator[" << calcid_ << "]::run() , TotNbProcPaq=" << totprocnpaq_ << endl;
377 /*
378 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) vpchk_[fib].Print();
379 cout << " ------------------------------------ " << endl;
380 */
381 delete[] fgallfibok;
382 }
383 catch (std::exception& exc) {
384 cout << " BRVisibilityCalculator[" << calcid_ << "]::run()/catched std::exception " << exc.what() << endl;
385 setRC(98);
386 return;
387 }
388 catch(...) {
389 cout << " BRVisibilityCalculator[" << calcid_ << "]::run()/catched unknown ... exception " << endl;
390 setRC(99);
391 return;
392 }
393
394}
395
[3967]396
[3872]397/* --Methode-- */
398int BRVisibilityCalculator::Process()
399{
[3910]400 if (totnbpaq_==0) UpdateChanIds(); // Appele ici pour etre sur que le thread de remplissage a mis l'info a jour.
[3956]401 if (fgdataraw_) { // Donnees firmware RAW apres TF soft
402 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
403 vpdatar_[2*fib] = reinterpret_cast< complex<r_4>* > (vprocpaq_[fib]);
404 vpdatar_[2*fib+1] = reinterpret_cast< complex<r_4>* >(vprocpaq_[fib]+memgr_.ProcPaqSize()/2) ;
405 }
406 }
407 else { // donnees firmware FFT
408 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
409 vpdata_[2*fib] = vpaq_[fib].Data1C();
410 vpdata_[2*fib+1] = vpaq_[fib].Data2C();
411 }
[3872]412 }
413
[4012]414 if ( CheckInterval4VisMtxSave() ) {
415 if ((prtlev_>0)&&(numfile_%prtmodulo_==0)) {
416 cout << numfile_ << "-BRVisCalc[" << calcid_ << "/" << nbcalc_ << "]::Process() NPaqProc="
417 << totnbpaq_ << " TotMegaFLOP=" << (uint_8)TotNbMegaFLOP() << " -> VisibMtx FileNum= " << numfile_ << endl;
[3872]418 }
[4012]419 SaveVisibilityMatrix();
[3872]420 }
421
[3893]422// kpair=numero sequentiel de la paire: 0->(0,0), 1->(0,1), 2->(0,2), 3->(0,3), 4->(1,1), 5->(1,2) ...
[3895]423 sa_size_t kpair=0;
424 sa_size_t k=0; // numero de ligne dans la matrice des visibilites
[3872]425 for(size_t i=0; i<vpdata_.size(); i++) {
426 for(size_t j=i; j<vpdata_.size(); j++) {
[3895]427 kpair++;
428 if (kpair<(pairst_+1)) continue;
429 if (kpair>=(pairst_+nbpairs_+1)) break;
[3915]430 if (fgpimp_&&(i!=j)&&((i+j)%2==0)) continue; // calcul des visib avec numero pair-impair + autocorrel
[3922]431 TVector< complex<r_4> > vis = vismtx_.Row(k); k++;
[3956]432
433 if (fgdataraw_) { // Donnees firmware RAW apres TF soft
434 for(sa_size_t f=1; f<vis.Size(); f++) {
[3963]435 vis(f) += vpdatar_[i][f] * conj(vpdatar_[j][f]);
[3956]436 }
[3872]437 }
[3956]438 else { // donnees firmware FFT
439 for(sa_size_t f=1; f<vis.Size(); f++) {
440 vis(f) += complex<r_4>((r_4)vpdata_[i][f].realB(), (r_4)vpdata_[i][f].imagB()) *
441 complex<r_4>((r_4)vpdata_[j][f].realB(), -(r_4)vpdata_[j][f].imagB());
442 }
443 }
[3895]444 nb_flop_ += (8.*(r_8)(vis.Size()-1));
[3872]445 }
446 }
447
[4012]448 npaqcumul_++;
[3872]449 moyfc_ += curfc_[0];
450 moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
451 if ((fgcktt_)&&(calcid_==0)) CheckTimeTag();
452 totnbpaq_++;
453 return 0;
454}
[4012]455/* --Methode-- */
456void BRVisibilityCalculator::SaveVisibilityMatrix()
457{
458 if (totnbpaq_ > 0) {
459 moyfc_/=npaqcumul_;
460 moytt_/=npaqcumul_;
461 UpdateVisMtxInfo(); // add/update keywords in the Info DVList
462 // ATTENTION : Matrice visibilites non moyennee
463 char nfile[48];
464 if (nbcalc_==1)
465 sprintf(nfile,"vismtx%d.%s",numfile_,OutFileExtension());
466 else
467 sprintf(nfile,"vismtx_%d_%d.%s",(int)calcid_,numfile_,OutFileExtension());
468 string flnm=outpath_+nfile;
469 if (fgfitsout_) { // Ecriture au format FITS
470 FitsInOutFile fo(flnm, FitsInOutFile::Fits_Create);
471 TArray<r_4> arvismtx = ArrCastC2R(vismtx_);
472 arvismtx.Info()=vismtx_.Info();
473 fo << arvismtx;
474 }
475 else { // Format PPF
476 POutPersist po(flnm);
477 po << vismtx_;
478 }
479 if (fgvisdt_) FillVisibTable(moyfc_, moytt_);
480 numfile_++;
481 }
482 vismtx_ = complex<r_4>((r_4)0.,(r_4)0.);
483 npaqcumul_=0;
484 moyfc_=moytt_=0.;
485 first_fc_=curfc_[0];
486 first_tt_=vpaq_[0].TimeTag();
487 // first_tmstamp_.SetNow(); // Current date and time
488 first_tmstamp_=cts_; // Current date and time
489
490 return;
491}
[3872]492
493/* --Methode-- */
[3878]494void BRVisibilityCalculator::UpdateVisMtxInfo()
495{
496 string ikey,ikdesc;
497 ikey="DATEOBS"; ikdesc=" Date/Time corresponding to TimeTagFirst";
498 vismtx_.Info().SetS(ikey,first_tmstamp_.ToString());
499 vismtx_.Info().SetComment(ikey,ikdesc);
500 ikey="FirstFC"; ikdesc="First FrameCounter";
501 vismtx_.Info().SetI(ikey,first_fc_);
502 vismtx_.Info().SetComment(ikey,ikdesc);
503 ikey="FirstTT"; ikdesc="First TimeTag";
[4012]504 vismtx_.Info().SetI(ikey,(first_tt_-ttfirst_[0]));
[3878]505 vismtx_.Info().SetComment(ikey,ikdesc);
506 ikey="LastFC"; ikdesc="Last FrameCounter";
507 vismtx_.Info().SetI(ikey,curfc_[0]);
508 vismtx_.Info().SetComment(ikey,ikdesc);
509 ikey="LastTT"; ikdesc="Last TimeTag";
510 vismtx_.Info().SetI(ikey,vpaq_[0].TimeTag()-ttfirst_[0]);
511 vismtx_.Info().SetComment(ikey,ikdesc);
512 ikey="MeanFC"; ikdesc="Mean FrameCounter";
513 vismtx_.Info().SetD(ikey,moyfc_);
514 vismtx_.Info().SetComment(ikey,ikdesc);
515 ikey="MeanTT"; ikdesc="Mean TimeTag";
516 vismtx_.Info().SetD(ikey,moytt_);
517 vismtx_.Info().SetComment(ikey,ikdesc);
518 ikey="NPAQSUM"; ikdesc="Number of paquets summed";
[4012]519 vismtx_.Info().SetI(ikey,npaqcumul_);
[3878]520}
521
522
523/* --Methode-- */
[3872]524int BRVisibilityCalculator::execute(int tid)
525{
526 vector<TwoByteComplex*> pvpdata(2*memgr_.NbFibres());
[3956]527 vector< complex<r_4>* > pvpdatar(2*memgr_.NbFibres());
[3872]528 size_t paqsz=memgr_.PaqSize();
[3967]529 size_t procpaqsz=memgr_.ProcPaqSize();
[3872]530 BRPaquet ppaq(paqsz);
531
[3895]532
[3872]533 sa_size_t fdelt = vismtx_.NCols()/nparthr_;
534 sa_size_t fdeb = tid*fdelt;
535 sa_size_t ffin = (tid+1)*fdelt;
536
537 if (fdeb<1) fdeb=1;
538 if ((ffin>vismtx_.NCols())||(tid==(nparthr_-1))) ffin=vismtx_.NCols();
539
540 for(size_t jp=0; jp<memgr_.NbPaquets(); jp++) { // boucle sur les paquets d'une zone
541 if (!fgallfibok[jp]) continue;
[3956]542 if (fgdataraw_) { // Donnees firmware RAW apres TF soft
543 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
[3967]544 Byte* procdatap=fprocbuff_[fib]+jp*procpaqsz;
545 pvpdatar[2*fib] = reinterpret_cast< complex<r_4>* > (procdatap);
546 pvpdatar[2*fib+1] = reinterpret_cast< complex<r_4>* >(procdatap+procpaqsz/2) ;
[3956]547 }
548 }
549 else { // donnees firmware FFT
550 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
551 ppaq.Set(fbuff_[fib]+jp*paqsz);
552 pvpdata[2*fib] = ppaq.Data1C();
553 pvpdata[2*fib+1] = ppaq.Data2C();
554 }
[3872]555 }
[3895]556
557
[3893]558// kpair=numero sequentiel de la paire: 0->(0,0), 1->(0,1), 2->(0,2), 3->(0,3), 4->(1,1), 5->(1,2) ...
[3895]559 sa_size_t kpair=0;
560 sa_size_t k=0; // numero de ligne dans la matrice des visibilites
[3872]561 for(size_t i=0; i<vpdata_.size(); i++) {
562 for(size_t j=i; j<vpdata_.size(); j++) {
[3895]563 kpair++;
564 if (kpair<(pairst_+1)) continue;
565 if (kpair>=(pairst_+nbpairs_+1)) break;
[3915]566 if (fgpimp_&&(i!=j)&&((i+j)%2==0)) continue; // calcul des visib avec numero pair-impair + autocorrel
[3895]567 TVector< complex<r_4> > vis = vismtx_.Row(k); k++;
[3956]568 if (fgdataraw_) { // Donnees firmware RAW apres TF soft
569 for(sa_size_t f=fdeb; f<ffin; f++) {
[3963]570 vis(f) += pvpdatar[i][f] * conj(pvpdatar[j][f]);
[3956]571 }
[3872]572 }
[3956]573 else { // donnees firmware FFT
574 for(sa_size_t f=fdeb; f<ffin; f++) {
575 vis(f) += complex<r_4>((r_4)pvpdata[i][f].realB(), (r_4)pvpdata[i][f].imagB()) *
576 complex<r_4>((r_4)pvpdata[j][f].realB(), -(r_4)pvpdata[j][f].imagB());
577 }
578 }
[3895]579 nb_flop_ += (8.*(r_8)(ffin-fdeb));
[3872]580 }
581 }
582
583 } // Fin de boucle sur les paquets
584
585 return 0;
586}
587
588/* --Methode-- */
589int BRVisibilityCalculator::FillVisibTable(double fcm, double ttm)
590{
591 double xnt[10];
592 xnt[0]=fcm; xnt[1]=ttm/1.25e8;
593
594 if (djf_<2) {
595 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
596 for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
597 xnt[2]=jf;
[3909]598 xnt[3]=chanpairnumall_(rv+pairst_);
[3872]599 xnt[4]=vismtx_(rv,jf).real()/(r_4)(nmean_);
600 xnt[5]=vismtx_(rv,jf).imag()/(r_4)(nmean_);
601 visdt_.AddRow(xnt);
602 }
603 }
604 }
605 else {
606 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
607 for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
608 r_4 moyreal=0.;
609 r_4 moyimag=0.;
610 sa_size_t jjfmx=jf+djf_;
611 if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
612 for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
613 moyreal+=vismtx_(rv,jjf).real();
614 moyimag+=vismtx_(rv,jjf).imag();
615 }
616 xnt[2]=jf+djf_/2;
[3909]617 xnt[3]=chanpairnumall_(rv+pairst_);
[3872]618 xnt[4]=moyreal/(r_4)(nmean_*djf_);
619 xnt[5]=moyimag/(r_4)(nmean_*djf_);
620 visdt_.AddRow(xnt);
621 }
622 }
623 }
624 return 0;
625}
626
627/* --Methode-- */
628int BRVisibilityCalculator::CheckTimeTag()
629{
630 if (totnbpaq_==0) {
631 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
632 vlasttt_[fib]=ttfirst_[fib];
633 if (ttmtx_.NCols()>0) {
634 fcmtx_(fib,totnbpaq_) = curfc_[fib];
635 ttmtx_(fib,totnbpaq_) = vlasttt_[fib];
636 }
637 }
638 return 0;
639 }
640 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
641 int_8 ld = (int_8)vpaq_[fib].TimeTag()-(int_8)vlasttt_[fib];
642 int_8 fd = (int_8)vpaq_[fib].TimeTag()-(int_8)ttfirst_[fib]-(int_8)vpaq_[0].TimeTag()+(int_8)ttfirst_[0];
643 /* if ( (ld < mindeltatt_) || (fd<-5) || (fd>5)) { vbadtt_[fib]++; vnsamett_[fib]++; }
644 else {
645 if (fd!=0) vnsamett_[fib]++;
646 }
647 */
648 if (ld < mindeltatt_) vbadtt_[fib]++;
649 else {
650 if (fd != 0) vnsamett_[fib]++;
651 if ((fd<-5)||(fd>5)) vndiff5tt_[fib]++;
652 }
653 vlasttt_[fib]=vpaq_[fib].TimeTag();
654 if (totnbpaq_<ttmtx_.NCols()) {
655 fcmtx_(fib,totnbpaq_) = curfc_[fib];
656 ttmtx_(fib,totnbpaq_) = vlasttt_[fib];
657 }
658 }
659 return 0;
660}
661
662//-------------------------------------------------------------------------------
663// Classe Groupe (ensemble) de Calculateur de Visibilites, tournant en parallele
664//-------------------------------------------------------------------------------
665
666/* --Methode-- */
[3915]667BRVisCalcGroup::BRVisCalcGroup(size_t nbcalc, RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean,
668 uint_4 pair1, uint_4 nbpairs, bool fgpimp, size_t nthr)
[3872]669 : tm_(false)
670{
671 if ((nbcalc<1)||(nbcalc>6))
672 throw ParmError("BRVisCalcGroup::BRVisCalcGroup NbCalc > 6 !");
673 for(size_t i=0; i<nbcalc; i++) {
674 BRVisibilityCalculator * viscp=new BRVisibilityCalculator(memgr, outpath, nmean, nthr);
[3915]675 viscp->DefineRank(nbcalc, i, pair1, nbpairs, fgpimp);
[3872]676 viscalcp_.push_back(viscp);
677 }
678}
679/* --Methode-- */
680BRVisCalcGroup::~BRVisCalcGroup()
681{
682 for(size_t i=0; i<viscalcp_.size(); i++)
683 delete viscalcp_[i];
684}
685/* --Methode-- */
[3956]686MemZStatus BRVisCalcGroup::SetMemZAction(MemZaction mmact)
687{
688 MemZaction mmzas[10]={MemZA_ProcA,MemZA_ProcB,MemZA_ProcC,MemZA_ProcD,MemZA_ProcE,MemZA_ProcF,
689 MemZA_ProcG,MemZA_ProcH,MemZA_ProcI,MemZA_ProcJ};
690 MemZStatus mmzss[10]={MemZS_ProcA,MemZS_ProcB,MemZS_ProcC,MemZS_ProcD,MemZS_ProcE,MemZS_ProcF,
691 MemZS_ProcG,MemZS_ProcH,MemZS_ProcI,MemZS_ProcJ};
692 size_t ia=9999;
693 for(int i=0; i<10; i++) {
694 if (mmact==mmzas[i]) { ia=i; break; }
695 }
696 if (ia>=10)
697 throw ParmError("BRVisCalcGroup::SetMemZAction() Bad MemZaction mmact !");
698 if ((ia+viscalcp_.size())>10)
699 throw ParmError("BRVisCalcGroup::SetMemZAction() MemZaction mmact too high !");
700 for(size_t i=0; i<viscalcp_.size(); i++) {
701 viscalcp_[i]->SetMemZAction(mmzas[ia]); ia++;
702 }
703 return mmzss[ia-1];
704}
705/* --Methode-- */
[3872]706int BRVisCalcGroup::SelectFreqBinning(uint_4 freq1, uint_4 freq2, uint_4 nbfreq)
707{
708 int rc=0;
709 for(size_t i=0; i<viscalcp_.size(); i++)
710 rc=viscalcp_[i]->SelectFreqBinning(freq1, freq2, nbfreq);
711 return rc;
712}
713/* --Methode-- */
[3920]714void BRVisCalcGroup::ActivateVisDTable(bool fgfdt)
715{
716 for(size_t i=0; i<viscalcp_.size(); i++)
717 viscalcp_[i]->ActivateVisDTable(fgfdt);
718}
719/* --Methode-- */
[3956]720void BRVisCalcGroup::SetFitsOutput()
721{
722 for(size_t i=0; i<viscalcp_.size(); i++)
723 viscalcp_[i]->SetFitsOutput();
724}
725/* --Methode-- */
726void BRVisCalcGroup::SetPPFOutput()
727{
728 for(size_t i=0; i<viscalcp_.size(); i++)
729 viscalcp_[i]->SetPPFOutput();
[3967]730}
[3956]731/* --Methode-- */
[3967]732void BRVisCalcGroup::SetFFTData()
733{
734 for(size_t i=0; i<viscalcp_.size(); i++)
735 viscalcp_[i]->SetFFTData();
736}
737/* --Methode-- */
738void BRVisCalcGroup::SetRawData()
739{
740 for(size_t i=0; i<viscalcp_.size(); i++)
741 viscalcp_[i]->SetRawData();
742}
743/* --Methode-- */
[4012]744void BRVisCalcGroup::SetNPaqIntervalMode(uint_4 nmean)
745{
746 for(size_t i=0; i<viscalcp_.size(); i++)
747 viscalcp_[i]->SetNPaqIntervalMode(nmean);
748}
749/* --Methode-- */
750void BRVisCalcGroup::SetTimeIntervalMode(double dtime)
751{
752 for(size_t i=0; i<viscalcp_.size(); i++)
753 viscalcp_[i]->SetTimeIntervalMode(dtime);
754}
755/* --Methode-- */
[3923]756void BRVisCalcGroup::SetPrintLevel(int lev, uint_8 prtmodulo)
757{
758 for(size_t i=0; i<viscalcp_.size(); i++)
759 viscalcp_[i]->SetPrintLevel(lev,prtmodulo);
760}
761/* --Methode-- */
[3872]762void BRVisCalcGroup::start()
763{
764 for(size_t i=0; i<viscalcp_.size(); i++)
765 viscalcp_[i]->start();
766 tm_.SplitQ();
767}
768/* --Methode-- */
769void BRVisCalcGroup::join()
770{
771 r_8 totflop=0.;
772 for(size_t i=0; i<viscalcp_.size(); i++) {
773 viscalcp_[i]->join();
[3895]774 // cout << " BRVisCalcGroup::join()/ VisibCalc[" << i << "]->TotNbMegaFLOP()="
775 // << viscalcp_[i]->TotNbFLOP()/1024e3 << endl;
[3872]776 totflop += viscalcp_[i]->TotNbFLOP();
777 }
778 tm_.SplitQ();
779 cout << " ----------------------------------------------------------" << endl;
[3895]780 cout << " BRVisCalcGroup::join() : Finished " << viscalcp_.size() << " VisibilityCalculator(s)" << endl;
781 cout << " ... Elaspsed time: " << tm_.PartialElapsedTimems()
782 << " ms (total:" << tm_.TotalElapsedTimems() << ")" << endl;
783 double mflopsrate=totflop/(r_8)tm_.PartialElapsedTimems()/(1024.);
[3872]784 cout << " ... TotalMegaFLOP= " << totflop/(1024.e3) << " @ "
[3895]785 << mflopsrate << " MFLOP/s" << " (=" << mflopsrate/(r_8)viscalcp_.size() << "/VisibCalcObject)" << endl;
[3872]786 cout << " ----------------------------------------------------------" << endl;
787 return;
788}
Note: See TracBrowser for help on using the repository browser.