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

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

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