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

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

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