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

Last change on this file since 3876 was 3876, checked in by ansari, 15 years ago

Ajout de la fonctionalite de calcul des visibilites (sur donnees firmware FFT uniquement) dans le programme d'acquisition mfacq.cc , Reza 09/09/2010
File size: 16.0 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 "matharr.h"
15#include "timestamp.h"
16#include "ctimer.h"
17
18#include "brviscalc.h"
19
20
21//---------------------------------------------------------------------
22// Classe de traitement - calcul de visibilite pour n fibres
23//---------------------------------------------------------------------
24
25/* --Methode-- */
26BRVisibilityCalculator::BRVisibilityCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, size_t nthr)
27 : BRBaseProcessor(memgr), paralex_(*this, nthr), nparthr_(nthr),
28 outpath_(outpath), nmean_(nmean), nbcalc_(1), calcid_(0), vpdata_(2*memgr.NbFibres())
29 // , dtfos_(outpath+"visdt.fits", Fits_Create), visdt_(dtfos_, 1024, true);
30{
31 DefineRank(1,0);
32
33 uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();
34 chanum_.SetSize(maxnpairs);
35 sa_size_t k=0;
36 for(size_t i=0; i<2*memgr_.NbFibres(); i++) vpdata_[i]=NULL;
37 for(size_t i=0; i<2*memgr_.NbFibres(); i++) {
38 for(size_t j=i; j<2*memgr_.NbFibres(); j++) {
39 chanum_(k) = (i+1)*100+(j+1); k++;
40 }
41 }
42 SelectPairs();
43
44 // visdt_.AddFloatColumn("mfc");
45 visdt_.AddDoubleColumn("mfc");
46 visdt_.AddDoubleColumn("mtt");
47 visdt_.AddIntegerColumn("jfreq");
48 visdt_.AddIntegerColumn("numch");
49 visdt_.AddFloatColumn("vre");
50 visdt_.AddFloatColumn("vim");
51
52 if (nmean_ < 1) nmean_=memgr_.NbPaquets();
53 if (nmean_ < 1) nmean_=1;
54
55 cout << " BRVisibilityCalculator::/Info nmean=" << nmean_ << endl;
56
57 totnbpaq_=0;
58 numfile_=0;
59 nb_flop_=0.;
60 moyfc_=moytt_=0.;
61
62 fgallfibok=NULL;
63 fgcktt_=false;
64 setNameId("viscalc", 0);
65}
66
67/* --Methode-- */
68BRVisibilityCalculator::~BRVisibilityCalculator()
69{
70 if (totnbpaq_<1) return;
71 cout << " BRVisibilityCalculator - Visibility Datatable : " << endl;
72 cout << visdt_;
73
74 string filename;
75 filename = outpath_+"visdt.ppf";
76 if (nbcalc_>1) {
77 char sbuff[32];
78 sprintf(sbuff,"visdt_%d",(int)calcid_);
79 filename = outpath_+sbuff;
80 }
81 POutPersist po(filename);
82 po << visdt_;
83 if (calcid_ == 0) {
84 POutPersist poc(outpath_+"chanum.ppf");
85 poc << chanum_;
86
87 if (fgcktt_) {
88 cout << " BRVisibilityCalculator - Check TimeTag done: TotNPaqProc= " << totnbpaq_ << endl;
89 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
90 cout << " BRTTCheck-Fiber[" << fib << "] NBadTT=" << vbadtt_[fib] << " NDiffTT>5="
91 << vndiff5tt_[fib] << " NotSameTT=" << vnsamett_[fib] << endl;
92 }
93 POutPersist pott(outpath_+"ttfcmtx.ppf");
94 pott << PPFNameTag("FC") << fcmtx_;
95 pott << PPFNameTag("TT") << ttmtx_;
96 }
97 }
98}
99
100/* --Methode-- */
101void BRVisibilityCalculator::DefineRank(uint_4 nbc, uint_4 cid)
102{
103 if ((nbc>6)||(cid>=nbc))
104 throw ParmError("BRVisibilityCalculator::DefineRank() NbCalc > 6 !");
105 nbcalc_=nbc;
106 calcid_=cid;
107 if (nbcalc_>1) {
108 uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();
109 uint_4 npairs=maxnpairs/nbcalc_;
110 if (calcid_==(nbcalc_-1))
111 SelectPairs(calcid_*npairs, maxnpairs-calcid_*npairs);
112 else
113 SelectPairs(calcid_*npairs, npairs);
114 MemZaction mmzas[6]={MemZA_ProcA,MemZA_ProcB,MemZA_ProcC,MemZA_ProcD,MemZA_ProcE,MemZA_ProcF};
115 SetMemZAction(mmzas[calcid_]);
116 setNameId("viscalc_grp", calcid_);
117 }
118 return ;
119}
120
121/* --Methode-- */
122uint_4 BRVisibilityCalculator::SelectPairs(uint_4 pair1, uint_4 nbpairs)
123{
124 BRPaquet paq(memgr_.PaqSize());
125 uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();
126
127 if (pair1 >= maxnpairs) pair1=maxnpairs-1;
128 if (nbpairs > maxnpairs-pair1) nbpairs=maxnpairs-pair1;
129 pairst_=pair1;
130 nbpairs_=nbpairs;
131 vismtx_.SetSize(nbpairs_, paq.DataSize()/4);
132 return nbpairs_;
133}
134
135/* --Methode-- */
136int BRVisibilityCalculator::SelectFreqBinning(uint_4 freq1, uint_4 freq2, uint_4 nbfreq)
137{
138 jf1_=freq1; jf2_=freq2;
139 if ((jf1_<1)||(jf1_>=vismtx_.NCols())) jf1_=1;
140 if ((jf2_<1)||(jf2_>=vismtx_.NCols())||(jf2_<jf1_)) jf2_=vismtx_.NCols()-1;
141 if (nbfreq<1) nbfreq=(jf2_-jf1_);
142 djf_=(jf2_-jf1_)/nbfreq;
143 if (djf_<1) djf_=1;
144 cout << " BRVisibilityCalculator::SelectFreqBinning/Info JF1=" << jf1_
145 << " JF2=" << jf2_ << " DJF=" << djf_ << endl;
146
147}
148
149
150/* --Methode-- */
151int BRVisibilityCalculator::ActivateTimeTagCheck(uint_8 maxnpaq)
152{
153 mindeltatt_=memgr_.PaqSize()/2;
154 if (mindeltatt_<1) mindeltatt_=1;
155 fcmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
156 ttmtx_.SetSize(memgr_.NbFibres(), maxnpaq);
157 vlasttt_.resize(memgr_.NbFibres(), 0);
158 vbadtt_.resize(memgr_.NbFibres(), 0);
159 vnsamett_.resize(memgr_.NbFibres(), 0);
160 vndiff5tt_.resize(memgr_.NbFibres(), 0);
161
162 fgcktt_=true;
163 cout << " BRVisibilityCalculator::ActivateTimeTagCheck() - TT/Fc matrix NCols=" << maxnpaq
164 << " MinDeltaTT=" << mindeltatt_ << endl;
165
166 return 0;
167}
168
169
170/* --Methode-- */
171void BRVisibilityCalculator::run()
172{
173 cout << " BRVisibilityCalculator[" << calcid_ << "/" << nbcalc_
174 << "]::run() - Starting " << " NFibers=" << memgr_.NbFibres()
175 << " NChan=" << 2*memgr_.NbFibres() << " NPairs=" << nbpairs_ << " First:" << pairst_ << endl;
176
177 if (nparthr_ < 2) return BRBaseProcessor::run();
178 // Execution multithread parallele
179 setRC(1);
180 int rc=0;
181 try {
182 if ((nmean_%memgr_.NbPaquets())!=0) {
183 uint_4 mnmean = (nmean_/memgr_.NbPaquets()+1)*memgr_.NbPaquets();
184 cout << " BRVisibilityCalculator::run()/Info changing nmean=" << nmean_ << " to multiple of"
185 << " memgr_.NbPaquets() -> " << mnmean << endl;
186 nmean_=mnmean;
187 }
188 paralex_.SetParallelTask(*this);
189 cout << " BRVisibilityCalculator::run()/Info : starting ParallelExecutor with nThreads="
190 << paralex_.nThreads() << " ... " << endl;
191 paralex_.start();
192
193 fgallfibok = new bool[memgr_.NbPaquets()];
194
195 size_t paqsz=memgr_.PaqSize();
196 bool fgrun=true;
197 while (fgrun) {
198 if (stop_) break;
199 if (memgr_.GetRunState() == MemZR_Stopped) break;
200 int mid = memgr_.FindMemZoneId(mmact_); // (MemZA_ProcA);
201 // Byte* buffg = memgr_.GetMemZone(mid);
202 // if (buffg == NULL) {
203 if (mid < 0) {
204 cout << "BRVisibilityCalculator[" << calcid_ << "]::run()/ERROR FindMemZoneId("
205 << (int)mmact_ << ") ->" << mid << ") -> NULL" << endl;
206 setRC(7); fgrun=false;
207 break;
208 }
209 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
210 fbuff_[fib] = memgr_.GetMemZone(mid,fib);
211 if (fbuff_[fib] == NULL) { // cela ne devrait pas arriver
212 cout << "BRVisibilityCalculator[" << calcid_ << "]::run()/ERROR memgr.GetMemZone(" << mid << "," << fib << ") -> NULL" << endl;
213 setRC(9); fgrun=false;
214 break;
215 }
216 }
217
218 if (totnbpaq_%nmean_ == 0) {
219 if (totnbpaq_ > 0) {
220 moyfc_/=nmean_;
221 moytt_/=nmean_;
222 vismtx_.Info()["MeanFC"] = moyfc_;
223 vismtx_.Info()["MeanTT"] = moytt_;
224 vismtx_.Info()["NPAQSUM"] = nmean_;
225
226 // ATTENTION : Matrice visibilites non moyennee
227 char nfile[48];
228 if (nbcalc_==1)
229 sprintf(nfile,"vismtx%d.ppf",numfile_);
230 else
231 sprintf(nfile,"vismtx_%d_%d.ppf",(int)calcid_,numfile_);
232 string flnm=outpath_+nfile;
233 POutPersist po(flnm);
234 po << vismtx_;
235 cout << numfile_ << "-BRVisibilityCalculator[" << calcid_ << "]::run() NPaqProc="
236 << totnbpaq_ << " -> Visibility Matrix in " << flnm << endl;
237 FillVisibTable(moyfc_, moytt_);
238 numfile_++;
239 }
240 vismtx_ = complex<r_4>((r_4)0.,(r_4)0.);
241 moyfc_=moytt_=0.;
242 }
243
244 for(size_t jp=0; jp<memgr_.NbPaquets(); jp++) { // boucle sur les paquets d'une zone
245 fgallfibok[jp]=fgokallfibers_=true;
246 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
247 vpaq_[fib].Set(fbuff_[fib]+jp*paqsz);
248 vfgok_[fib] = vpchk_[fib].Check(vpaq_[fib],curfc_[fib]);
249 if (!vfgok_[fib]) fgallfibok[jp]=fgokallfibers_=false;
250 }
251 if (fgokallfibers_) {
252 if (totprocnpaq_==0) {
253 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
254 fcfirst_[fib]=curfc_[fib];
255 ttfirst_[fib]=vpaq_[fib].TimeTag();
256 }
257 }
258 totprocnpaq_++;
259 moyfc_ += curfc_[0];
260 moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
261 if ((fgcktt_)&&(calcid_==0)) CheckTimeTag();
262 totnbpaq_++;
263 }
264 } // Fin de boucle sur les paquets
265
266 // Execution parallele pour calcul des visibilites par bandes de frequence
267 int rcpex=paralex_.execute();
268 if (rcpex!=0) cout << " BRVisibilityCalculator[" << calcid_ << "]::run() / Error Rc[paralex_.execute()]=" << rcpex << endl;
269
270 memgr_.FreeMemZone(mid, mmsta_); // (MemZS_ProcA);
271 } // Fin de boucle sur les zones a traiter
272 //------------------------------------
273 cout << " --------- END BRVisibilityCalculator[" << calcid_ << "]::run() , TotNbProcPaq=" << totprocnpaq_ << endl;
274 /*
275 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) vpchk_[fib].Print();
276 cout << " ------------------------------------ " << endl;
277 */
278 delete[] fgallfibok;
279 }
280 catch (std::exception& exc) {
281 cout << " BRVisibilityCalculator[" << calcid_ << "]::run()/catched std::exception " << exc.what() << endl;
282 setRC(98);
283 return;
284 }
285 catch(...) {
286 cout << " BRVisibilityCalculator[" << calcid_ << "]::run()/catched unknown ... exception " << endl;
287 setRC(99);
288 return;
289 }
290
291}
292
293/* --Methode-- */
294int BRVisibilityCalculator::Process()
295{
296
297 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
298 vpdata_[2*fib] = vpaq_[fib].Data1C();
299 vpdata_[2*fib+1] = vpaq_[fib].Data2C();
300 }
301
302 if (totnbpaq_%nmean_ == 0) {
303 if (totnbpaq_ > 0) {
304 moyfc_/=nmean_;
305 moytt_/=nmean_;
306 vismtx_.Info()["MeanFC"] = moyfc_;
307 vismtx_.Info()["MeanTT"] = moytt_;
308 vismtx_.Info()["NPAQSUM"] = nmean_;
309
310 // ATTENTION : Matrice visibilites non moyennee
311 char nfile[48];
312 if (nbcalc_==1)
313 sprintf(nfile,"vismtx%d.ppf",numfile_);
314 else
315 sprintf(nfile,"vismtx_%d_%d.ppf",(int)calcid_,numfile_);
316 string flnm=outpath_+nfile;
317 POutPersist po(flnm);
318 po << vismtx_;
319 cout << numfile_ << "-BRVisibilityCalculator::Process() NPaqProc="
320 << totnbpaq_ << " -> Visibility Matrix in " << flnm << endl;
321 FillVisibTable(moyfc_, moytt_);
322 numfile_++;
323 }
324 vismtx_ = complex<r_4>((r_4)0.,(r_4)0.);
325 moyfc_=moytt_=0.;
326 }
327
328 sa_size_t k=0;
329 for(size_t i=0; i<vpdata_.size(); i++) {
330 for(size_t j=i; j<vpdata_.size(); j++) {
331 size_t kpair=i*vpdata_.size()+j;
332 if (kpair<pairst_) continue;
333 if (kpair>=(pairst_+nbpairs_)) break;
334 TVector< complex<r_4> > vis = vismtx_.Row(k); k++;
335 for(sa_size_t f=1; f<vis.Size(); f++) {
336 vis(f) += complex<r_4>((r_4)vpdata_[i][f].realB(), (r_4)vpdata_[i][f].imagB()) *
337 complex<r_4>((r_4)vpdata_[j][f].realB(), -(r_4)vpdata_[j][f].imagB());
338 }
339 nb_flop_ += (7.*(r_8)vis.Size());
340 }
341 }
342
343 moyfc_ += curfc_[0];
344 moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
345 if ((fgcktt_)&&(calcid_==0)) CheckTimeTag();
346 totnbpaq_++;
347 return 0;
348}
349
350/* --Methode-- */
351int BRVisibilityCalculator::execute(int tid)
352{
353 vector<TwoByteComplex*> pvpdata(2*memgr_.NbFibres());
354 size_t paqsz=memgr_.PaqSize();
355 BRPaquet ppaq(paqsz);
356
357 sa_size_t fdelt = vismtx_.NCols()/nparthr_;
358 sa_size_t fdeb = tid*fdelt;
359 sa_size_t ffin = (tid+1)*fdelt;
360
361 if (fdeb<1) fdeb=1;
362 if ((ffin>vismtx_.NCols())||(tid==(nparthr_-1))) ffin=vismtx_.NCols();
363
364 for(size_t jp=0; jp<memgr_.NbPaquets(); jp++) { // boucle sur les paquets d'une zone
365 if (!fgallfibok[jp]) continue;
366 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
367 ppaq.Set(fbuff_[fib]+jp*paqsz);
368 pvpdata[2*fib] = ppaq.Data1C();
369 pvpdata[2*fib+1] = ppaq.Data2C();
370 }
371 sa_size_t k=0;
372 for(size_t i=0; i<vpdata_.size(); i++) {
373 for(size_t j=i; j<vpdata_.size(); j++) {
374 size_t kpair=i*vpdata_.size()+j;
375 if (kpair<pairst_) continue;
376 if (kpair>=(pairst_+nbpairs_)) break;
377 TVector< complex<r_4> > vis = vismtx_.Row(k); k++;
378 for(sa_size_t f=fdeb; f<ffin; f++) {
379 vis(f) += complex<r_4>((r_4)pvpdata[i][f].realB(), (r_4)pvpdata[i][f].imagB()) *
380 complex<r_4>((r_4)pvpdata[j][f].realB(), -(r_4)pvpdata[j][f].imagB());
381 }
382 nb_flop_ += (7.*(r_8)(ffin-fdeb));
383 }
384 }
385
386 } // Fin de boucle sur les paquets
387
388 return 0;
389}
390
391/* --Methode-- */
392int BRVisibilityCalculator::FillVisibTable(double fcm, double ttm)
393{
394 double xnt[10];
395 xnt[0]=fcm; xnt[1]=ttm/1.25e8;
396
397 if (djf_<2) {
398 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
399 for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
400 xnt[2]=jf;
401 xnt[3]=chanum_(rv+pairst_);
402 xnt[4]=vismtx_(rv,jf).real()/(r_4)(nmean_);
403 xnt[5]=vismtx_(rv,jf).imag()/(r_4)(nmean_);
404 visdt_.AddRow(xnt);
405 }
406 }
407 }
408 else {
409 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
410 for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
411 r_4 moyreal=0.;
412 r_4 moyimag=0.;
413 sa_size_t jjfmx=jf+djf_;
414 if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
415 for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
416 moyreal+=vismtx_(rv,jjf).real();
417 moyimag+=vismtx_(rv,jjf).imag();
418 }
419 xnt[2]=jf+djf_/2;
420 xnt[3]=chanum_(rv+pairst_);
421 xnt[4]=moyreal/(r_4)(nmean_*djf_);
422 xnt[5]=moyimag/(r_4)(nmean_*djf_);
423 visdt_.AddRow(xnt);
424 }
425 }
426 }
427 return 0;
428}
429
430/* --Methode-- */
431int BRVisibilityCalculator::CheckTimeTag()
432{
433 if (totnbpaq_==0) {
434 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
435 vlasttt_[fib]=ttfirst_[fib];
436 if (ttmtx_.NCols()>0) {
437 fcmtx_(fib,totnbpaq_) = curfc_[fib];
438 ttmtx_(fib,totnbpaq_) = vlasttt_[fib];
439 }
440 }
441 return 0;
442 }
443 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
444 int_8 ld = (int_8)vpaq_[fib].TimeTag()-(int_8)vlasttt_[fib];
445 int_8 fd = (int_8)vpaq_[fib].TimeTag()-(int_8)ttfirst_[fib]-(int_8)vpaq_[0].TimeTag()+(int_8)ttfirst_[0];
446 /* if ( (ld < mindeltatt_) || (fd<-5) || (fd>5)) { vbadtt_[fib]++; vnsamett_[fib]++; }
447 else {
448 if (fd!=0) vnsamett_[fib]++;
449 }
450 */
451 if (ld < mindeltatt_) vbadtt_[fib]++;
452 else {
453 if (fd != 0) vnsamett_[fib]++;
454 if ((fd<-5)||(fd>5)) vndiff5tt_[fib]++;
455 }
456 vlasttt_[fib]=vpaq_[fib].TimeTag();
457 if (totnbpaq_<ttmtx_.NCols()) {
458 fcmtx_(fib,totnbpaq_) = curfc_[fib];
459 ttmtx_(fib,totnbpaq_) = vlasttt_[fib];
460 }
461 }
462 return 0;
463}
464
465//-------------------------------------------------------------------------------
466// Classe Groupe (ensemble) de Calculateur de Visibilites, tournant en parallele
467//-------------------------------------------------------------------------------
468
469/* --Methode-- */
470BRVisCalcGroup::BRVisCalcGroup(size_t nbcalc, RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, size_t nthr)
471 : tm_(false)
472{
473 if ((nbcalc<1)||(nbcalc>6))
474 throw ParmError("BRVisCalcGroup::BRVisCalcGroup NbCalc > 6 !");
475 for(size_t i=0; i<nbcalc; i++) {
476 BRVisibilityCalculator * viscp=new BRVisibilityCalculator(memgr, outpath, nmean, nthr);
477 viscp->DefineRank(nbcalc, i);
478 viscalcp_.push_back(viscp);
479 }
480}
481/* --Methode-- */
482BRVisCalcGroup::~BRVisCalcGroup()
483{
484 for(size_t i=0; i<viscalcp_.size(); i++)
485 delete viscalcp_[i];
486}
487/* --Methode-- */
488int BRVisCalcGroup::SelectFreqBinning(uint_4 freq1, uint_4 freq2, uint_4 nbfreq)
489{
490 int rc=0;
491 for(size_t i=0; i<viscalcp_.size(); i++)
492 rc=viscalcp_[i]->SelectFreqBinning(freq1, freq2, nbfreq);
493 return rc;
494}
495/* --Methode-- */
496void BRVisCalcGroup::start()
497{
498 for(size_t i=0; i<viscalcp_.size(); i++)
499 viscalcp_[i]->start();
500 tm_.SplitQ();
501}
502/* --Methode-- */
503void BRVisCalcGroup::join()
504{
505 r_8 totflop=0.;
506 for(size_t i=0; i<viscalcp_.size(); i++) {
507 viscalcp_[i]->join();
508 totflop += viscalcp_[i]->TotNbFLOP();
509 }
510 tm_.SplitQ();
511 cout << " ----------------------------------------------------------" << endl;
512 cout << " BRVisCalcGroup::join() : Finished , Elaspsed time: " << tm_.PartialElapsedTimems()
513 << " ms (total:" << tm_.TotalElapsedTimems() << ")" << endl;
514 cout << " ... TotalMegaFLOP= " << totflop/(1024.e3) << " @ "
515 << totflop/(r_8)tm_.PartialElapsedTimems()/(1024) << " MFLOP/s" << endl;
516 cout << " ----------------------------------------------------------" << endl;
517 return;
518}
Note: See TracBrowser for help on using the repository browser.