Line data Source code
1 : //# SolveDataBuffer.cc: Implementation of SolveDataBuffer.h
2 : //# Copyright (C) 2008
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
29 : #include <msvis/MSVis/VisBuffer2.h>
30 : #include <msvis/MSVis/VisBufferComponents2.h>
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/casa/Arrays/ArrayPartMath.h>
33 : #include <casacore/casa/Arrays/MaskArrMath.h>
34 : #include <casacore/casa/IO/ArrayIO.h>
35 : #include <casacore/casa/iomanip.h>
36 : #include <casacore/casa/Exceptions/Error.h>
37 : #include <casacore/casa/Containers/Block.h>
38 : #include <casacore/casa/Utilities/Assert.h>
39 : #include <casacore/casa/Quanta/MVTime.h>
40 : #include <casacore/casa/BasicSL/Constants.h>
41 : #include <casacore/casa/BasicMath/Math.h>
42 :
43 : using namespace casacore;
44 : using namespace casa::vi;
45 :
46 :
47 : using namespace casacore;
48 : namespace casa { //# NAMESPACE CASA - BEGIN
49 :
50 0 : SolveDataBuffer::SolveDataBuffer() :
51 0 : vb_(0),
52 0 : nAnt_(0),
53 0 : freqs_(0),
54 0 : centroidFreq_(0.0),
55 0 : corrs_(0),
56 0 : feedPa_(0),
57 0 : focusChan_p(-1),
58 0 : infocusFlagCube_p(),
59 0 : infocusWtSpec_p(),
60 0 : infocusVisCube_p(),
61 0 : infocusModelVisCube_p(),
62 0 : workingFlagCube_p(),
63 0 : workingWtSpec_p(),
64 0 : residuals_p(),
65 0 : residFlagCube_p(),
66 0 : diffResiduals_p()
67 0 : {}
68 :
69 47196 : SolveDataBuffer::SolveDataBuffer(const vi::VisBuffer2& vb) :
70 47196 : vb_(0),
71 47196 : nAnt_(0),
72 47196 : freqs_(0),
73 47196 : centroidFreq_(0.0),
74 47196 : corrs_(0),
75 47196 : feedPa_(0),
76 47196 : focusChan_p(-1),
77 47196 : infocusFlagCube_p(),
78 47196 : infocusWtSpec_p(),
79 47196 : infocusVisCube_p(),
80 47196 : infocusModelVisCube_p(),
81 47196 : workingFlagCube_p(),
82 47196 : workingWtSpec_p(),
83 47196 : residuals_p(),
84 47196 : residFlagCube_p(),
85 94392 : diffResiduals_p()
86 :
87 : {
88 :
89 47196 : vb_= vi::VisBuffer2::factory(VbRekeyable);
90 :
91 47196 : initFromVB(vb);
92 47196 : }
93 :
94 0 : SolveDataBuffer::SolveDataBuffer(const SolveDataBuffer& sdb) :
95 0 : vb_(),
96 0 : nAnt_(0),
97 0 : freqs_(0),
98 0 : centroidFreq_(0.0),
99 0 : corrs_(0),
100 0 : feedPa_(0),
101 0 : focusChan_p(-1),
102 0 : infocusFlagCube_p(),
103 0 : infocusWtSpec_p(),
104 0 : infocusVisCube_p(),
105 0 : infocusModelVisCube_p(),
106 0 : workingFlagCube_p(),
107 0 : workingWtSpec_p(),
108 0 : residuals_p(),
109 0 : residFlagCube_p(),
110 0 : diffResiduals_p()
111 : {
112 : // Copy from the other's VB2
113 : // NB: that VB2 won't be attached!!
114 0 : initFromVB(*(sdb.vb_));
115 :
116 : // copy over freqs_, corrs_,feedPa, nAnt_
117 : // (things that normally require being attached)
118 0 : freqs_.assign(sdb.freqs_);
119 0 : centroidFreq_=sdb.centroidFreq_;
120 0 : corrs_.assign(sdb.corrs_);
121 0 : feedPa_.assign(sdb.feedPa_);
122 0 : nAnt_=sdb.nAnt_;
123 0 : }
124 :
125 47196 : SolveDataBuffer::~SolveDataBuffer()
126 : {
127 47196 : if (vb_) delete vb_;
128 47196 : vb_=0;
129 :
130 47196 : }
131 :
132 : //SolveDataBuffer& SolveDataBuffer::operator=(const VisBuffer& other)
133 :
134 :
135 748805 : Bool SolveDataBuffer::Ok() {
136 : // Ok if net unflagged weight is positive
137 748805 : if (nfalse(flagCube())>0) {
138 748002 : Float wtsum=sum(weightSpectrum()(!flagCube()));
139 748002 : return (wtsum>0.0f);
140 : }
141 : else
142 803 : return false;
143 : }
144 :
145 : /*
146 : // Divide corrected by model
147 : void SolveDataBuffer::divideCorrByModel() {
148 :
149 : Int nCor(nCorrelations());
150 : Int nChan(nChannels());
151 : Int nRow(nRows());
152 : Float amp(1.0);
153 : Complex cor(1.0);
154 :
155 : Cube<Complex> vC; vC.reference(visCubeCorrected());
156 : Cube<Complex> vM; vM.reference(visCubeModel());
157 : Cube<Float> wS; wS.reference(weightSpectrum());
158 :
159 : for (Int irow=0;irow<nRow;++irow) {
160 : if (!flagRow()(irow)) {
161 : for (Int ich=0;ich<nChan;++ich) {
162 : for (Int icorr=0;icorr<nCor;icorr++) {
163 : if (!flagCube()(icorr,ich,irow)) {
164 : amp=abs(vM(icorr,ich,irow));
165 : if (amp>0.0f) {
166 : // Divide corr by model
167 : vC(icorr,ich,irow)/=vM(icorr,ich,irow);
168 : // Adjust weight by square of model amp
169 : wS(icorr,ich,irow)*=square(amp);
170 : }
171 : } // !*fl
172 : else {
173 : // zero data and weight
174 : vC(icorr,ich,irow)=Complex(0.0);
175 : wS(icorr,ich,irow)=0.0;
176 : }
177 : // model always unity after division
178 : vM(icorr,ich,irow)=Complex(1.0);
179 : } // icorr
180 : } // ich
181 : } // !*flR
182 : } // irow
183 : }
184 : */
185 :
186 47161 : void SolveDataBuffer::enforceAPonData(const String& apmode)
187 : {
188 :
189 : // TBD: apply to model, too? (never matters, since we prob /modelVis already?
190 :
191 : // ONLY if something to do
192 47161 : if (apmode=="A" || apmode=="P") {
193 8 : Int nCor(nCorrelations());
194 8 : Int nChan(nChannels());
195 8 : Int nRow(nRows());
196 8 : Float amp(1.0);
197 8 : Complex cor(1.0);
198 :
199 8 : Cube<Complex> vC; vC.reference(visCubeCorrected());
200 8 : Cube<Float> wS; wS.reference(weightSpectrum());
201 :
202 368 : for (Int irow=0;irow<nRow;++irow) {
203 360 : if (!flagRow()(irow)) {
204 720 : for (Int ich=0;ich<nChan;++ich) {
205 1800 : for (Int icorr=0;icorr<nCor;icorr++) {
206 1440 : if (!flagCube()(icorr,ich,irow)) {
207 :
208 1440 : amp=abs(vC(icorr,ich,irow));
209 1440 : if (amp>0.0f) {
210 :
211 1440 : if (apmode=="P") {
212 : // we will scale by amp to make data phase-only
213 720 : cor=Complex(amp,0.0);
214 : // Adjust weight by square of amp
215 720 : wS(icorr,ich,irow)*=square(amp);
216 : }
217 720 : else if (apmode=="A")
218 : // we will scale by "phase" to make data amp-only
219 : // no weight adjustment
220 720 : cor=vC(icorr,ich,irow)/amp;
221 :
222 : // Apply the complex scaling
223 1440 : vC(icorr,ich,irow)/=cor;
224 : }
225 : } // !*fl
226 : else {
227 0 : vC(icorr,ich,irow)=Complex(0.0);
228 0 : wS(icorr,ich,irow)=0.0;
229 : }
230 : } // icorr
231 : } // ich
232 : } // !*flR
233 : } // irow
234 :
235 8 : } // phase- or amp-only
236 :
237 47161 : }
238 :
239 47161 : void SolveDataBuffer::enforceSolveWeights(const Bool phandonly)
240 : {
241 : // If requested set cross-hand weights to zero (also---or only---flags?)
242 47161 : if(phandonly && this->nCorrelations()>2)
243 24272 : this->weightSpectrum()(Slice(1, 2, 1), Slice(), Slice()).set(0.0);
244 :
245 : // Set flagged weights to zero, ensuring they don't get used in accumulations
246 47161 : this->weightSpectrum()(this->flagCube())=0.0f;
247 :
248 : // Set flagged data to zero (some solve types don't look at flags)
249 47161 : Cube<Complex> vCC(this->visCubeCorrected());
250 47161 : vCC(this->flagCube())=Complex(0.0f);
251 :
252 47161 : }
253 :
254 :
255 :
256 348816 : void SolveDataBuffer::setFocusChan(const Int focusChan)
257 : {
258 :
259 : // Nominally focus on the whole data array
260 348816 : IPosition focusblc(3,0,0,0);
261 348816 : IPosition focustrc(vb_->getShape());
262 348816 : focustrc-=1;
263 :
264 : // if focusChan non-negative, select the single channel
265 348816 : if (focusChan>-1)
266 348816 : focusblc(1)=focustrc(1)=focusChan;
267 :
268 348816 : infocusFlagCube_p.reference(flagCube()(focusblc,focustrc));
269 348816 : infocusWtSpec_p.reference(weightSpectrum()(focusblc,focustrc));
270 348816 : infocusVisCube_p.reference(visCubeCorrected()(focusblc,focustrc));
271 348816 : infocusModelVisCube_p.reference(visCubeModel()(focusblc,focustrc));
272 :
273 : // Remember current in-focus channel
274 348816 : focusChan_p=focusChan;
275 :
276 348816 : }
277 :
278 329767 : void SolveDataBuffer::sizeResiduals(const Int& nPar,
279 : const Int& nDiff)
280 : {
281 :
282 329767 : IPosition ip1(vb_->getShape());
283 329767 : if (focusChan_p>-1)
284 329767 : ip1(1)=1;
285 329767 : residuals_p.resize(ip1);
286 329767 : residuals_p.set(0.0);
287 329767 : residFlagCube_p.resize(ip1);
288 329767 : residFlagCube_p.set(false);
289 :
290 329767 : if (nPar>0 && nDiff>0) {
291 329767 : IPosition ip2(5,ip1(0),nPar,ip1(1),ip1(2),nDiff);
292 329767 : diffResiduals_p.resize(ip2);
293 329767 : diffResiduals_p.set(0.0);
294 329767 : }
295 :
296 329767 : }
297 :
298 329767 : void SolveDataBuffer::initResidWithModel()
299 : {
300 :
301 : // Copy (literally) the in-focus model to the residual workspace
302 : // TBD: weights?
303 329767 : residuals_p = infocusModelVisCube_p;
304 329767 : residFlagCube_p = infocusFlagCube_p;
305 :
306 : // Ensure contiguity, because CalSolver will depend on this
307 329767 : AlwaysAssert(residFlagCube_p.contiguousStorage(),AipsError);
308 329767 : AlwaysAssert(residuals_p.contiguousStorage(),AipsError);
309 :
310 :
311 329767 : }
312 :
313 329767 : void SolveDataBuffer::finalizeResiduals()
314 : {
315 :
316 : // Subtract in-focus obs data from residuals workspace
317 329767 : residuals_p -= infocusVisCube_p;
318 :
319 : // TBD: zero flagged samples here?
320 :
321 329767 : }
322 :
323 : // Manage working weights
324 0 : void SolveDataBuffer::updateWorkingFlags()
325 : {
326 :
327 0 : workingFlagCube_p|=residFlagCube_p;
328 :
329 0 : }
330 :
331 :
332 : // Manage working weights
333 0 : void SolveDataBuffer::updateWorkingWeights(Bool doL1, Float L1clamp)
334 : {
335 :
336 0 : workingWtSpec_p.resize(0,0,0); // nominally forces implicit use of infocusWtSpec
337 0 : if (doL1) {
338 :
339 : //cout << "****Using L1 weights!!! (" << L1clamp << ") ******" << endl;
340 :
341 0 : workingWtSpec_p.assign(infocusWtSpec_p); // init from nominal weights
342 0 : const Cube<Bool>& wFC(this->const_workingFlagCube()); // adaptive access to
343 :
344 0 : Cube<Float> Ramp(amplitude(this->residuals()));
345 :
346 0 : IPosition shRamp(Ramp.shape());
347 0 : Int nCorr=shRamp(0);
348 0 : Int nChan=shRamp(1);
349 0 : Int nRow=shRamp(2);
350 :
351 : // TBD: assert common shapes!
352 :
353 0 : for (Int irow=0;irow<nRow;++irow) {
354 0 : if (!flagRow()(irow)) {
355 0 : for (Int ich=0;ich<nChan;++ich) {
356 0 : for (Int icorr=0;icorr<nCorr;++icorr) {
357 0 : Float& nomWt=workingWtSpec_p(icorr,ich,irow);
358 0 : Float& Ra=Ramp(icorr,ich,irow);
359 0 : if (!wFC(icorr,ich,irow) &&
360 0 : nomWt>0.0 &&
361 0 : Ra>0.0 ) {
362 0 : nomWt/=sqrt(Ra*Ra+L1clamp);
363 : }
364 : else {
365 0 : nomWt=0.0;
366 : }
367 : }
368 : }
369 : }
370 : }
371 :
372 0 : }
373 : // else {
374 : // cout << "****Using nominal weights!!!******" << endl;
375 : //infocusWorkingWtSpec_p.reference(infocusWtSpec_p);
376 : //}
377 :
378 0 : }
379 :
380 0 : void SolveDataBuffer::reportData()
381 : {
382 :
383 0 : Slice corrs(0,2,3), sl; // p-hands only
384 :
385 0 : Vector<Int> a1(this->antenna1()), a2(this->antenna2());
386 0 : Cube<Complex> vCC(this->visCubeCorrected()(corrs,sl,sl));
387 0 : Cube<Complex> vCM(this->visCubeModel()(corrs,sl,sl));
388 0 : Cube<Float> wtS(this->weightSpectrum()(corrs,sl,sl));
389 0 : Cube<Bool> flg(this->flagCube()(corrs,sl,sl));
390 :
391 : /*
392 : cout << "******Making weights globally uniform!!!!!!!!!!!!!!!!" << endl;
393 : // this makes VI2 match VI1
394 : wtS.set(1.0f);
395 : //*/
396 :
397 : /*
398 : cout << "*******Zeroing zero-wt data (for K solutions)!!!!!!!!!!!!!" << endl;
399 : // this makes VI2 match VI1 solint='int' when some channels flagged
400 : vCC(wtS==0.0f)=Complex(0.0);
401 : //*/
402 :
403 :
404 : /*
405 : cout << "*******making weights uniform (for K solution testing)!!!!!!!!!!!!!" << endl;
406 : wtS.set(max(wtS));
407 : //*/
408 :
409 : /*
410 : cout << "*******zero off-wt data (for K solution testing)!!!!!!!!!!!!!" << endl;
411 : vCC(wtS!=max(wtS))=0.0f;
412 : wtS(wtS!=max(wtS))=0.0f;
413 : //*/
414 :
415 : /*
416 : // this makes VI2 match VI1 for solint='inf' (and > 'int')
417 : cout << "*******renormalizing off-wt data (for K solution testing)!!!!!!!!!!!!!" << endl;
418 : vCC*=wtS;
419 : vCC/=max(wtS);
420 : //*/
421 :
422 :
423 0 : cout << "Time=" << MVTime(time()(0)/C::day).string(MVTime::YMD,7) << endl;
424 0 : cout.precision(8);
425 0 : for (Int irow=0;irow<nRows();++irow) {
426 0 : for (Int ich=0;ich<nChannels();++ich) {
427 0 : cout << std::setw(2) << a1(irow) << "-" << std::setw(2) << a2(irow) << " ";
428 0 : if (nChannels()>1) cout << "ich=" << ich << " ";
429 0 : cout << "fl=" << flg(Slice(),Slice(ich),Slice(irow)).nonDegenerate() << " ";
430 0 : cout << "wt=" << wtS(Slice(),Slice(ich),Slice(irow)).nonDegenerate() << " ";
431 0 : cout << "cM=" << vCM(Slice(),Slice(ich),Slice(irow)).nonDegenerate() << " ";
432 0 : cout << "cC=" << vCC(Slice(),Slice(ich),Slice(irow)).nonDegenerate() << " ";
433 0 : cout << endl;
434 : }
435 : }
436 0 : cout << "*****************************************************************************" << endl;
437 0 : }
438 :
439 47196 : void SolveDataBuffer::initFromVB(const vi::VisBuffer2& vb)
440 : {
441 :
442 : // The required VB2 components
443 : vi::VisBufferComponents2 comps =
444 : vi::VisBufferComponents2::these({VisBufferComponent2::ObservationId,
445 : VisBufferComponent2::ArrayId,
446 : VisBufferComponent2::Scan,
447 : VisBufferComponent2::FieldId,
448 : VisBufferComponent2::DataDescriptionIds,
449 : VisBufferComponent2::SpectralWindows,
450 : VisBufferComponent2::Antenna1,
451 : VisBufferComponent2::Antenna2,
452 : VisBufferComponent2::Time,
453 : VisBufferComponent2::TimeCentroid,
454 : VisBufferComponent2::NCorrelations,
455 : VisBufferComponent2::NChannels,
456 : VisBufferComponent2::NRows,
457 : VisBufferComponent2::FlagRow,
458 : VisBufferComponent2::FlagCube,
459 : VisBufferComponent2::WeightSpectrum,
460 : VisBufferComponent2::VisibilityCubeCorrected,
461 47196 : VisBufferComponent2::VisibilityCubeModel});
462 :
463 : // Copy required components from the supplied VB2:
464 47196 : vb_->copyComponents(vb,comps,True,True); // will fetch things, if needed
465 :
466 : // Set weights for flagged data to zero
467 47196 : Cube<Float> wtsp(vb_->weightSpectrum());
468 47196 : wtsp(vb_->flagCube())=0.0;
469 :
470 : // Store the frequeny info
471 : // TBD: also need bandwidth info....
472 47196 : if (vb.isAttached()) {
473 47196 : freqs_.assign(vb.getFrequencies(0));
474 47196 : corrs_.assign(vb.correlationTypes());
475 47196 : nAnt_=vb.nAntennas();
476 : }
477 : else {
478 : // Probably only needed in testing.... (gmoellen, 2016Aug04, 2019Jan07)
479 0 : cout << "The supplied VisBuffer2 is not attached to a ViImplementation2," << endl
480 0 : << " which is necessary to generate accurate frequency info." << endl
481 0 : << " This is probably just a test with a naked VisBuffer2." << endl
482 0 : << " Spoofing freq axis with 1 MHz channels at 100+10*ispw GHz." << endl
483 0 : << " Spoofing corr axis with [5,6,7,8] (circulars)" << endl;
484 :
485 0 : freqs_.resize(vb.nChannels());
486 0 : indgen(freqs_);
487 0 : freqs_*=1e6;
488 0 : freqs_+=100.0005e9; // _edge_ of first channel at 100 GHz.
489 0 : freqs_+=(10.0e9*vb.spectralWindows()(0)); // 10 GHz spacing of spws
490 :
491 0 : Int nC=vb.nCorrelations();
492 0 : corrs_.resize(nC);
493 0 : corrs_[0]=5;
494 0 : if (nC>1) corrs_[nC-1]=8;
495 0 : if (nC==4) {
496 0 : corrs_[1]=6;
497 0 : corrs_[2]=7;
498 : }
499 :
500 : // nAnt is last a2 index +1
501 : // Assumes simple sorting of these (which is how test data works)
502 0 : Int nR=vb.nRows();
503 0 : nAnt_=vb.antenna2()(nR-1)+1;
504 : }
505 :
506 : // Store the centroid freq (use mean, for now)
507 47196 : centroidFreq_ = mean(freqs_);
508 :
509 : // Store the feedPa info
510 47196 : if (vb.isAttached())
511 : // Assumes vb.time() is constant!
512 47196 : feedPa_.assign(vb.feedPa(vb.time()(0)));
513 :
514 47196 : }
515 0 : void SolveDataBuffer::cleanUp()
516 : {
517 :
518 : // Zero-size all workspaces
519 0 : infocusFlagCube_p.resize();
520 0 : infocusWtSpec_p.resize();
521 0 : infocusVisCube_p.resize();
522 0 : infocusModelVisCube_p.resize();
523 :
524 0 : residuals_p.resize();
525 0 : residFlagCube_p.resize();
526 0 : diffResiduals_p.resize();
527 :
528 0 : }
529 :
530 552 : String SolveDataBuffer::polBasis() const
531 : {
532 :
533 : // UNKNOWN, nominally
534 552 : String polBas("UNKNOWN");
535 :
536 552 : Int nC(corrs_.nelements());
537 552 : if (nC<1)
538 : // Can't tell, so return UNKNOWN
539 0 : return polBas;
540 :
541 552 : if (corrs_(0)==5)
542 480 : polBas=String("CIRC");
543 552 : if (corrs_(0)==9)
544 72 : polBas=String("LIN");
545 :
546 552 : return polBas;
547 :
548 0 : }
549 :
550 :
551 19605 : SDBList::SDBList() :
552 19605 : nSDB_(0),
553 19605 : SDB_(),
554 19605 : freqs_(),
555 19605 : aggCentroidFreq_(0),
556 19605 : aggCentroidFreqOK_(false)
557 19605 : {}
558 :
559 19605 : SDBList::~SDBList()
560 : {
561 : // Delete the SDBs
562 66801 : for (Int i=0;i<nSDB_;++i)
563 47196 : delete SDB_[i];
564 19605 : SDB_.resize(0);
565 19605 : nSDB_=0;
566 19605 : }
567 :
568 47196 : void SDBList::add(const vi::VisBuffer2& vb)
569 : {
570 :
571 : // Enlarge the list, copying existing SDB pointers
572 47196 : SDB_.resize(nSDB_+1,true);
573 :
574 : // Generate the new SolveDataBuffer
575 47196 : SDB_[nSDB_] = new SolveDataBuffer(vb);
576 :
577 : // increment the count
578 47196 : nSDB_++;
579 :
580 : // Clear the freqs_ info (forces recalculation)
581 47196 : freqs_.resize(0);
582 47196 : aggCentroidFreq_=0.0;
583 47196 : aggCentroidFreqOK_=false;
584 :
585 47196 : }
586 :
587 1811555 : SolveDataBuffer& SDBList::operator()(Int i)
588 : {
589 1811555 : if (i<nSDB_)
590 1811555 : return *SDB_[i];
591 : else
592 0 : throw(AipsError("SDBList::operator(): requests non-existent SolveDataBuffer."));
593 :
594 : }
595 :
596 19605 : Int SDBList::aggregateObsId() const {
597 19605 : if (nSDB_>0)
598 : // Obs Id from first SDB
599 19605 : return SDB_[0]->observationId()(0);
600 0 : throw(AipsError("SDBList::aggregateObsId(): No SDBs in this SDBList yet."));
601 : }
602 :
603 19605 : Int SDBList::aggregateScan() const {
604 19605 : if (nSDB_>0)
605 : // Scan number from first SDB
606 19605 : return SDB_[0]->scan()(0);
607 0 : throw(AipsError("SDBList::aggregateScan(): No SDBs in this SDBList yet."));
608 : }
609 :
610 39231 : Int SDBList::aggregateSpw() const {
611 39231 : if (nSDB_>0)
612 : // from first SDB
613 39231 : return SDB_[0]->spectralWindow()(0);
614 0 : throw(AipsError("SDBList::aggregateSpw(): No SDBs in this SDBList yet."));
615 : }
616 :
617 19605 : Int SDBList::aggregateFld() const {
618 19605 : if (nSDB_>0)
619 : // from first SDB
620 19605 : return SDB_[0]->fieldId()(0);
621 0 : throw(AipsError("SDBList::aggregateFld(): No SDBs in this SDBList yet."));
622 : }
623 :
624 35 : Double SDBList::aggregateTime() const {
625 :
626 : // Simple average of the mean times in each SDB
627 : // (TBD: Improve with attention to flags/weights?)
628 35 : if (nSDB_>0) {
629 35 : Double aTime(0.0);
630 70 : for (Int isdb=0;isdb<nSDB_;++isdb)
631 35 : aTime+=mean(SDB_[isdb]->time());
632 35 : aTime/=nSDB_;
633 35 : return aTime;
634 : }
635 : else
636 0 : throw(AipsError("SDBList::aggregateTime(): No SDBs in this SDBList yet."));
637 : }
638 :
639 19605 : Double SDBList::aggregateTimeCentroid() const {
640 :
641 : // Weighted average of SDBs' timeCentroids
642 19605 : if (nSDB_>0) {
643 19605 : Double aTime(0.0);
644 19605 : Double aWt(0.0);
645 19605 : Vector<Double> wtvD;
646 66801 : for (Int isdb=0;isdb<nSDB_;++isdb) {
647 47196 : SolveDataBuffer& sdb(*SDB_[isdb]);
648 94392 : Vector<Float> wtv(partialSums(sdb.weightSpectrum(),IPosition(2,0,1)));
649 47196 : wtvD.resize(wtv.nelements());
650 47196 : convertArray(wtvD,wtv);
651 47196 : aTime+=sum(wtvD*sdb.timeCentroid());
652 47196 : aWt+=sum(wtvD);
653 47196 : }
654 19605 : if (aWt>0.0)
655 19570 : aTime/=aWt;
656 : else
657 : // Use aggregateTime if no unflagged data
658 35 : aTime=aggregateTime();
659 :
660 19605 : return aTime;
661 19605 : }
662 : else
663 0 : throw(AipsError("SDBList::aggregateFld(): No SDBs in this SDBList yet."));
664 : }
665 :
666 :
667 : // How many data chans?
668 : // Currently, this insists on uniformity over all SDBs
669 : // In future, we may _sum_ the SDBs nChans, and
670 : // enable forming aggregate spectra (e.g., for common normalization)
671 : // This will require focusChan loop over SDBs...
672 3816 : Int SDBList::nChannels() const {
673 :
674 3816 : Int nChan=SDB_[0]->nChannels(); // from first
675 :
676 : // Trap non-uniformity, for now
677 49491 : for (Int isdb=1;isdb<nSDB_;++isdb)
678 45675 : AlwaysAssert((SDB_[isdb]->nChannels()==nChan),AipsError);
679 :
680 : // Reach here, then ok
681 3816 : return nChan;
682 :
683 : }
684 :
685 1848 : const Vector<Double>& SDBList::freqs() const {
686 :
687 1848 : if (nSDB_==0)
688 0 : throw(AipsError("SDBList::freqs(): No SDBs in this SDBList yet."));
689 :
690 :
691 1848 : if (nSDB_==1) {
692 : // Only one SDB, just return that one's freqs
693 1768 : const Vector<Double>& f(SDB_[0]->freqs()); // from first SDB
694 1768 : return f;
695 : }
696 :
697 : // Reach here, more than one SDB, need to gather info
698 80 : if (freqs_.nelements()==0) {
699 : // Haven't accumumlated yet
700 :
701 : // How many channels in aggregate?
702 : // (This will insist on uniformity over SDBs)
703 80 : Int nchan(this->nChannels());
704 :
705 : // Will accumulate mean freqs here
706 80 : freqs_.resize(nchan);
707 80 : freqs_.set(0.0);
708 :
709 : // Average over SDBs, counting each spw exactly _once_
710 : // Map to keep track of unique spws
711 80 : std::set<Int> spws;
712 80 : Int nSpw(0);
713 1240 : for (Int isdb=0;isdb<nSDB_;++isdb) {
714 1160 : Int ispw=SDB_[isdb]->spectralWindow()(0);
715 1160 : if (spws.count(ispw)<1) {
716 83 : freqs_+=SDB_[isdb]->freqs();
717 83 : spws.insert(ispw); // Record that we got this spw
718 83 : ++nSpw;
719 : }
720 : }
721 : // Divide by nSpw
722 80 : freqs_/=Double(nSpw);
723 :
724 80 : }
725 : // else: freqs_ already filled previously...
726 :
727 80 : return freqs_;
728 :
729 : }
730 :
731 :
732 : // ~Centroid frequency over all SDBs
733 438 : casacore::Double SDBList::centroidFreq() const {
734 :
735 : // Calculate a simple mean frequency, for now
736 438 : Double fsum(0.0f);
737 438 : Int nf(0);
738 53508 : for (Int isdb=0;isdb<nSDB_;++isdb) {
739 53070 : const Vector<Double> f(SDB_[isdb]->freqs());
740 1366590 : for (uInt ich=0;ich<f.nelements();++ich) {
741 1313520 : fsum+=f(ich);
742 1313520 : ++nf;
743 : }
744 53070 : }
745 438 : return fsum/Double(nf);
746 : }
747 :
748 : // ~Centroid frequency over all SDBs
749 17757 : casacore::Double SDBList::aggregateCentroidFreq() const {
750 :
751 : // Trap no data case
752 17757 : if (nSDB_==0)
753 0 : throw(AipsError("SDBList::aggregateCentroidFreq(): No SDBs in this SDBList yet."));
754 :
755 : // Need to calculate?
756 17757 : if (!aggCentroidFreqOK_) {
757 :
758 17757 : if (nSDB_==1) {
759 : // from first and only SDB
760 17488 : aggCentroidFreq_=SDB_[0]->centroidFreq();
761 : }
762 : else {
763 : // More than one SDB, need to gather simple mean
764 : // TBD: weight this by per-SDB bandwidth
765 269 : aggCentroidFreq_=0.0;
766 27049 : for (Int isdb=0;isdb<nSDB_;++isdb)
767 26780 : aggCentroidFreq_+=SDB_[isdb]->centroidFreq();
768 269 : aggCentroidFreq_/=Double(nSDB_);
769 : }
770 : // We've calculated it
771 17757 : aggCentroidFreqOK_=true;
772 : }
773 :
774 : // Reach here, one way or another we have a good value, so return it
775 17757 : return aggCentroidFreq_;
776 :
777 : }
778 :
779 63 : String SDBList::polBasis() const
780 : {
781 63 : String polBas(SDB_[0]->polBasis());
782 :
783 : // Trap non-uniformity, for now
784 552 : for (Int isdb=1;isdb<nSDB_;++isdb)
785 489 : AlwaysAssert((SDB_[isdb]->polBasis()==polBas),AipsError);
786 :
787 63 : return polBas;
788 0 : }
789 :
790 :
791 : // How many antennas?
792 : // Currently, this insists on uniformity over all SDBs
793 21 : Int SDBList::nAntennas() const {
794 :
795 21 : Int nAnt=SDB_[0]->nAntennas();
796 :
797 : // Trap non-uniformity, for now
798 184 : for (Int isdb=1;isdb<nSDB_;++isdb)
799 163 : AlwaysAssert((SDB_[isdb]->nAntennas()==nAnt),AipsError);
800 :
801 : // Reach here, then ok
802 21 : return nAnt;
803 :
804 : }
805 :
806 : // How many correlations?
807 : // Currently, this insists on uniformity over all SDBs
808 1832 : Int SDBList::nCorrelations() const {
809 :
810 1832 : Int nCorr=SDB_[0]->nCorrelations();
811 :
812 : // Trap non-uniformity, for now
813 54527 : for (Int isdb=1;isdb<nSDB_;++isdb)
814 52695 : AlwaysAssert((SDB_[isdb]->nCorrelations()==nCorr),AipsError);
815 :
816 : // Reach here, then ok
817 1832 : return nCorr;
818 :
819 : }
820 :
821 :
822 :
823 :
824 :
825 19605 : Bool SDBList::Ok() {
826 :
827 19672 : for (Int i=0;i<nSDB_;++i)
828 19637 : if (SDB_[i]->Ok()) return True;
829 :
830 : // If we get here, either no SDBs, or none have non-zero weight.
831 35 : return False;
832 :
833 : }
834 :
835 : /*
836 : void SDBList::divideCorrByModel()
837 : {
838 : for (Int i=0;i<nSDB_;++i)
839 : SDB_[i]->divideCorrByModel();
840 : }
841 : */
842 :
843 19570 : void SDBList::enforceAPonData(const String& apmode)
844 : {
845 66731 : for (Int i=0;i<nSDB_;++i)
846 47161 : SDB_[i]->enforceAPonData(apmode);
847 19570 : }
848 19570 : void SDBList::enforceSolveWeights(const Bool pHandOnly)
849 : {
850 66731 : for (Int i=0;i<nSDB_;++i)
851 47161 : SDB_[i]->enforceSolveWeights(pHandOnly);
852 19570 : }
853 :
854 0 : void SDBList::sizeResiduals(const Int& nPar, const Int& nDiff)
855 : {
856 0 : for (Int i=0;i<nSDB_;++i)
857 0 : SDB_[i]->sizeResiduals(nPar,nDiff);
858 0 : }
859 0 : void SDBList::initResidWithModel()
860 : {
861 0 : for (Int i=0;i<nSDB_;++i)
862 0 : SDB_[i]->initResidWithModel();
863 0 : }
864 0 : void SDBList::finalizeResiduals()
865 : {
866 0 : for (Int i=0;i<nSDB_;++i)
867 0 : SDB_[i]->finalizeResiduals();
868 0 : }
869 :
870 : // Manage working flags
871 0 : void SDBList::updateWorkingFlags()
872 : {
873 0 : for (Int i=0;i<nSDB_;++i)
874 0 : SDB_[i]->updateWorkingFlags();
875 0 : }
876 :
877 : // Manage working weights
878 0 : void SDBList::updateWorkingWeights(casacore::Bool doL1,casacore::Float clamp)
879 : {
880 0 : for (Int i=0;i<nSDB_;++i)
881 0 : SDB_[i]->updateWorkingWeights(doL1,clamp);
882 0 : }
883 :
884 :
885 0 : void SDBList::reportData()
886 : {
887 0 : cout << "nSDB=" << nSDB_ << endl;
888 0 : for (Int i=0;i<nSDB_;++i) {
889 0 : cout << "isdb=" << i << endl;
890 0 : SDB_[i]->reportData();
891 : }
892 0 : }
893 :
894 :
895 21 : Int SDBList::extendCrossHandBaselineFlags(String& message)
896 : {
897 :
898 :
899 : // This enforces uniform nAnt in all SDBs
900 21 : Int nAnt(this->nAntennas());
901 21 : Int nBln=nAnt*(nAnt-1)/2;
902 :
903 : // channelized flags per baseline, cross-corr pair
904 21 : Cube<Bool> aggFlag(nChannels(),nAnt,nAnt,true);
905 :
906 :
907 : // Initialize aggregate flags
908 : // If an antenna pair is not avail here, it won't be used (default flagged)
909 : // Subsequent SDBs might flag more baselines
910 :
911 21 : Int isdb0(0);
912 : {
913 :
914 : // Step forward to first SDB that has some unflagged data
915 : // TBD: ...for now; in future, may wish to use least-flagged one, if we can also apply a
916 : // non-zero threshold in accumulation below...
917 21 : while (isdb0<nSDB_ && nfalse(SDB_[isdb0]->flagCube())<1)
918 0 : ++isdb0;
919 :
920 : // Trap all bad data case
921 21 : if (isdb0==nSDB_)
922 0 : throw(AipsError("Accumulated data entirely flagged."));
923 :
924 21 : Cube<Bool>& flc(SDB_[isdb0]->flagCube());
925 :
926 21 : const Vector<Int>& a1(SDB_[isdb0]->antenna1());
927 21 : const Vector<Int>& a2(SDB_[isdb0]->antenna2());
928 :
929 21 : Int nRows(SDB_[isdb0]->nRows());
930 21 : Int nChan(SDB_[isdb0]->nChannels());
931 :
932 966 : for (Int irow=0;irow<nRows;++irow) {
933 8505 : for (Int ichan=0;ichan<nChan;++ichan) {
934 7560 : aggFlag(ichan,a1(irow),a2(irow))=(flc(1,ichan,irow)||flc(2,ichan,irow)); // Either cross-hand flagged
935 : }
936 : }
937 : }
938 :
939 :
940 : // Accumulate from other SDBs
941 21 : Int nBadSDB(0);
942 205 : for (Int isdb=0;isdb<nSDB_;++isdb) {
943 :
944 : // Skip the one we initialized with
945 184 : if (isdb==isdb0)
946 21 : continue;
947 :
948 163 : Cube<Bool>& flc(SDB_[isdb]->flagCube());
949 163 : const Vector<Int>& a1(SDB_[isdb]->antenna1());
950 163 : const Vector<Int>& a2(SDB_[isdb]->antenna2());
951 :
952 163 : Int ngood=nfalse(flc);
953 163 : if (ngood==0) {
954 0 : nBadSDB+=1;
955 0 : continue;
956 : }
957 :
958 163 : Int nRows(SDB_[isdb]->nRows());
959 163 : Int nChan(SDB_[isdb]->nChannels());
960 :
961 7498 : for (Int irow=0;irow<nRows;++irow) {
962 66015 : for (Int ichan=0;ichan<nChan;++ichan) {
963 58680 : aggFlag(ichan,a1(irow),a2(irow))|=(flc(1,ichan,irow)||flc(2,ichan,irow)); // Either cross-hand flagged
964 : }
965 : }
966 : }
967 :
968 : // cout << "aggFlag: " << ntrue(aggFlag) << "/" << nfalse(aggFlag) << " sh=" << aggFlag.shape() << endl;
969 :
970 : // Count good baselines, by channel
971 42 : Vector<size_t> goodBL(partialNFalse(aggFlag,IPosition(2,1,2)));
972 :
973 : // Median number of good baselines, over channel
974 21 : size_t medGoodBL=median(goodBL);
975 :
976 : // Count good channels (those with at least median number of good baselines)
977 21 : Int goodChan=ntrue(goodBL>=medGoodBL);
978 :
979 : // Apply aggregate flags to each SDB
980 : //cout << "Flag sum = ";
981 21 : Int nGoodSDB(0);
982 205 : for (Int isdb=0;isdb<nSDB_;++isdb) {
983 184 : Cube<Bool>& flc(SDB_[isdb]->flagCube());
984 184 : const Vector<Int>& a1(SDB_[isdb]->antenna1());
985 184 : const Vector<Int>& a2(SDB_[isdb]->antenna2());
986 :
987 184 : Int nRows(SDB_[isdb]->nRows());
988 184 : Int nChan(SDB_[isdb]->nChannels());
989 :
990 8464 : for (Int irow=0;irow<nRows;++irow) {
991 74520 : for (Int ichan=0;ichan<nChan;++ichan) {
992 66240 : if (aggFlag(ichan,a1(irow),a2(irow))) {
993 0 : flc(Slice(1,2,1),Slice(ichan),Slice(irow))=true;
994 : }
995 : }
996 : }
997 184 : if (nfalse(flc)>0) ++nGoodSDB;
998 :
999 : //cout << ntrue(flc) << " ";
1000 : }
1001 : //cout << endl;
1002 :
1003 : // Assemble message
1004 21 : ostringstream ostr;
1005 21 : ostr << "There is usable (unflagged) data in " << nGoodSDB << " (of " << nSDB_ << " total) data segments, with "
1006 21 : << goodChan << " good (of " << nChannels() << " total) channels (" << floor(10000.0*Double(goodChan)/Double(nChannels()))/100.0 << "%) having at least "
1007 21 : << medGoodBL << " good (of " << nBln << " total) baselines (" << floor(10000.0*Double(medGoodBL)/Double(nBln))/100.0 <<"%).";
1008 : // cout << ostr << endl;
1009 21 : message=ostr;
1010 :
1011 :
1012 : // Return number of surviving SDBs...
1013 21 : return nGoodSDB;
1014 :
1015 21 : }
1016 :
1017 : } //# NAMESPACE CASA - END
1018 :
|
