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 0 : SolveDataBuffer::SolveDataBuffer(const vi::VisBuffer2& vb) :
70 0 : vb_(0),
71 0 : nAnt_(0),
72 0 : freqs_(0),
73 0 : centroidFreq_(0.0),
74 0 : corrs_(0),
75 0 : feedPa_(0),
76 0 : focusChan_p(-1),
77 0 : infocusFlagCube_p(),
78 0 : infocusWtSpec_p(),
79 0 : infocusVisCube_p(),
80 0 : infocusModelVisCube_p(),
81 0 : workingFlagCube_p(),
82 0 : workingWtSpec_p(),
83 0 : residuals_p(),
84 0 : residFlagCube_p(),
85 0 : diffResiduals_p()
86 :
87 : {
88 :
89 0 : vb_= vi::VisBuffer2::factory(VbRekeyable);
90 :
91 0 : initFromVB(vb);
92 0 : }
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 0 : SolveDataBuffer::~SolveDataBuffer()
126 : {
127 0 : if (vb_) delete vb_;
128 0 : vb_=0;
129 :
130 0 : }
131 :
132 : //SolveDataBuffer& SolveDataBuffer::operator=(const VisBuffer& other)
133 :
134 :
135 0 : Bool SolveDataBuffer::Ok() {
136 : // Ok if net unflagged weight is positive
137 0 : if (nfalse(flagCube())>0) {
138 0 : Float wtsum=sum(weightSpectrum()(!flagCube()));
139 0 : return (wtsum>0.0f);
140 : }
141 : else
142 0 : 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 0 : 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 0 : if (apmode=="A" || apmode=="P") {
193 0 : Int nCor(nCorrelations());
194 0 : Int nChan(nChannels());
195 0 : Int nRow(nRows());
196 0 : Float amp(1.0);
197 0 : Complex cor(1.0);
198 :
199 0 : Cube<Complex> vC; vC.reference(visCubeCorrected());
200 0 : Cube<Float> wS; wS.reference(weightSpectrum());
201 :
202 0 : for (Int irow=0;irow<nRow;++irow) {
203 0 : if (!flagRow()(irow)) {
204 0 : for (Int ich=0;ich<nChan;++ich) {
205 0 : for (Int icorr=0;icorr<nCor;icorr++) {
206 0 : if (!flagCube()(icorr,ich,irow)) {
207 :
208 0 : amp=abs(vC(icorr,ich,irow));
209 0 : if (amp>0.0f) {
210 :
211 0 : if (apmode=="P") {
212 : // we will scale by amp to make data phase-only
213 0 : cor=Complex(amp,0.0);
214 : // Adjust weight by square of amp
215 0 : wS(icorr,ich,irow)*=square(amp);
216 : }
217 0 : else if (apmode=="A")
218 : // we will scale by "phase" to make data amp-only
219 : // no weight adjustment
220 0 : cor=vC(icorr,ich,irow)/amp;
221 :
222 : // Apply the complex scaling
223 0 : 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 0 : } // phase- or amp-only
236 :
237 0 : }
238 :
239 0 : void SolveDataBuffer::enforceSolveWeights(const Bool phandonly)
240 : {
241 : // If requested set cross-hand weights to zero (also---or only---flags?)
242 0 : if(phandonly && this->nCorrelations()>2)
243 0 : 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 0 : this->weightSpectrum()(this->flagCube())=0.0f;
247 :
248 : // Set flagged data to zero (some solve types don't look at flags)
249 0 : Cube<Complex> vCC(this->visCubeCorrected());
250 0 : vCC(this->flagCube())=Complex(0.0f);
251 :
252 0 : }
253 :
254 :
255 :
256 0 : void SolveDataBuffer::setFocusChan(const Int focusChan)
257 : {
258 :
259 : // Nominally focus on the whole data array
260 0 : IPosition focusblc(3,0,0,0);
261 0 : IPosition focustrc(vb_->getShape());
262 0 : focustrc-=1;
263 :
264 : // if focusChan non-negative, select the single channel
265 0 : if (focusChan>-1)
266 0 : focusblc(1)=focustrc(1)=focusChan;
267 :
268 0 : infocusFlagCube_p.reference(flagCube()(focusblc,focustrc));
269 0 : infocusWtSpec_p.reference(weightSpectrum()(focusblc,focustrc));
270 0 : infocusVisCube_p.reference(visCubeCorrected()(focusblc,focustrc));
271 0 : infocusModelVisCube_p.reference(visCubeModel()(focusblc,focustrc));
272 :
273 : // Remember current in-focus channel
274 0 : focusChan_p=focusChan;
275 :
276 0 : }
277 :
278 0 : void SolveDataBuffer::sizeResiduals(const Int& nPar,
279 : const Int& nDiff)
280 : {
281 :
282 0 : IPosition ip1(vb_->getShape());
283 0 : if (focusChan_p>-1)
284 0 : ip1(1)=1;
285 0 : residuals_p.resize(ip1);
286 0 : residuals_p.set(0.0);
287 0 : residFlagCube_p.resize(ip1);
288 0 : residFlagCube_p.set(false);
289 :
290 0 : if (nPar>0 && nDiff>0) {
291 0 : IPosition ip2(5,ip1(0),nPar,ip1(1),ip1(2),nDiff);
292 0 : diffResiduals_p.resize(ip2);
293 0 : diffResiduals_p.set(0.0);
294 0 : }
295 :
296 0 : }
297 :
298 0 : void SolveDataBuffer::initResidWithModel()
299 : {
300 :
301 : // Copy (literally) the in-focus model to the residual workspace
302 : // TBD: weights?
303 0 : residuals_p = infocusModelVisCube_p;
304 0 : residFlagCube_p = infocusFlagCube_p;
305 :
306 : // Ensure contiguity, because CalSolver will depend on this
307 0 : AlwaysAssert(residFlagCube_p.contiguousStorage(),AipsError);
308 0 : AlwaysAssert(residuals_p.contiguousStorage(),AipsError);
309 :
310 :
311 0 : }
312 :
313 0 : void SolveDataBuffer::finalizeResiduals()
314 : {
315 :
316 : // Subtract in-focus obs data from residuals workspace
317 0 : residuals_p -= infocusVisCube_p;
318 :
319 : // TBD: zero flagged samples here?
320 :
321 0 : }
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 0 : 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 0 : VisBufferComponent2::VisibilityCubeModel});
462 :
463 : // Copy required components from the supplied VB2:
464 0 : vb_->copyComponents(vb,comps,True,True); // will fetch things, if needed
465 :
466 : // Set weights for flagged data to zero
467 0 : Cube<Float> wtsp(vb_->weightSpectrum());
468 0 : wtsp(vb_->flagCube())=0.0;
469 :
470 : // Store the frequeny info
471 : // TBD: also need bandwidth info....
472 0 : if (vb.isAttached()) {
473 0 : freqs_.assign(vb.getFrequencies(0));
474 0 : corrs_.assign(vb.correlationTypes());
475 0 : 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 0 : centroidFreq_ = mean(freqs_);
508 :
509 : // Store the feedPa info
510 0 : if (vb.isAttached())
511 : // Assumes vb.time() is constant!
512 0 : feedPa_.assign(vb.feedPa(vb.time()(0)));
513 :
514 0 : }
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 0 : String SolveDataBuffer::polBasis() const
531 : {
532 :
533 : // UNKNOWN, nominally
534 0 : String polBas("UNKNOWN");
535 :
536 0 : Int nC(corrs_.nelements());
537 0 : if (nC<1)
538 : // Can't tell, so return UNKNOWN
539 0 : return polBas;
540 :
541 0 : if (corrs_(0)==5)
542 0 : polBas=String("CIRC");
543 0 : if (corrs_(0)==9)
544 0 : polBas=String("LIN");
545 :
546 0 : return polBas;
547 :
548 0 : }
549 :
550 :
551 0 : SDBList::SDBList() :
552 0 : nSDB_(0),
553 0 : SDB_(),
554 0 : freqs_(),
555 0 : aggCentroidFreq_(0),
556 0 : aggCentroidFreqOK_(false)
557 0 : {}
558 :
559 0 : SDBList::~SDBList()
560 : {
561 : // Delete the SDBs
562 0 : for (Int i=0;i<nSDB_;++i)
563 0 : delete SDB_[i];
564 0 : SDB_.resize(0);
565 0 : nSDB_=0;
566 0 : }
567 :
568 0 : void SDBList::add(const vi::VisBuffer2& vb)
569 : {
570 :
571 : // Enlarge the list, copying existing SDB pointers
572 0 : SDB_.resize(nSDB_+1,true);
573 :
574 : // Generate the new SolveDataBuffer
575 0 : SDB_[nSDB_] = new SolveDataBuffer(vb);
576 :
577 : // increment the count
578 0 : nSDB_++;
579 :
580 : // Clear the freqs_ info (forces recalculation)
581 0 : freqs_.resize(0);
582 0 : aggCentroidFreq_=0.0;
583 0 : aggCentroidFreqOK_=false;
584 :
585 0 : }
586 :
587 0 : SolveDataBuffer& SDBList::operator()(Int i)
588 : {
589 0 : if (i<nSDB_)
590 0 : return *SDB_[i];
591 : else
592 0 : throw(AipsError("SDBList::operator(): requests non-existent SolveDataBuffer."));
593 :
594 : }
595 :
596 0 : Int SDBList::aggregateObsId() const {
597 0 : if (nSDB_>0)
598 : // Obs Id from first SDB
599 0 : return SDB_[0]->observationId()(0);
600 0 : throw(AipsError("SDBList::aggregateObsId(): No SDBs in this SDBList yet."));
601 : }
602 :
603 0 : Int SDBList::aggregateScan() const {
604 0 : if (nSDB_>0)
605 : // Scan number from first SDB
606 0 : return SDB_[0]->scan()(0);
607 0 : throw(AipsError("SDBList::aggregateScan(): No SDBs in this SDBList yet."));
608 : }
609 :
610 0 : Int SDBList::aggregateSpw() const {
611 0 : if (nSDB_>0)
612 : // from first SDB
613 0 : return SDB_[0]->spectralWindow()(0);
614 0 : throw(AipsError("SDBList::aggregateSpw(): No SDBs in this SDBList yet."));
615 : }
616 :
617 0 : Int SDBList::aggregateFld() const {
618 0 : if (nSDB_>0)
619 : // from first SDB
620 0 : return SDB_[0]->fieldId()(0);
621 0 : throw(AipsError("SDBList::aggregateFld(): No SDBs in this SDBList yet."));
622 : }
623 :
624 0 : Double SDBList::aggregateTime() const {
625 :
626 : // Simple average of the mean times in each SDB
627 : // (TBD: Improve with attention to flags/weights?)
628 0 : if (nSDB_>0) {
629 0 : Double aTime(0.0);
630 0 : for (Int isdb=0;isdb<nSDB_;++isdb)
631 0 : aTime+=mean(SDB_[isdb]->time());
632 0 : aTime/=nSDB_;
633 0 : return aTime;
634 : }
635 : else
636 0 : throw(AipsError("SDBList::aggregateTime(): No SDBs in this SDBList yet."));
637 : }
638 :
639 0 : Double SDBList::aggregateTimeCentroid() const {
640 :
641 : // Weighted average of SDBs' timeCentroids
642 0 : if (nSDB_>0) {
643 0 : Double aTime(0.0);
644 0 : Double aWt(0.0);
645 0 : Vector<Double> wtvD;
646 0 : for (Int isdb=0;isdb<nSDB_;++isdb) {
647 0 : SolveDataBuffer& sdb(*SDB_[isdb]);
648 0 : Vector<Float> wtv(partialSums(sdb.weightSpectrum(),IPosition(2,0,1)));
649 0 : wtvD.resize(wtv.nelements());
650 0 : convertArray(wtvD,wtv);
651 0 : aTime+=sum(wtvD*sdb.timeCentroid());
652 0 : aWt+=sum(wtvD);
653 0 : }
654 0 : if (aWt>0.0)
655 0 : aTime/=aWt;
656 : else
657 : // Use aggregateTime if no unflagged data
658 0 : aTime=aggregateTime();
659 :
660 0 : return aTime;
661 0 : }
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 0 : Int SDBList::nChannels() const {
673 :
674 0 : Int nChan=SDB_[0]->nChannels(); // from first
675 :
676 : // Trap non-uniformity, for now
677 0 : for (Int isdb=1;isdb<nSDB_;++isdb)
678 0 : AlwaysAssert((SDB_[isdb]->nChannels()==nChan),AipsError);
679 :
680 : // Reach here, then ok
681 0 : return nChan;
682 :
683 : }
684 :
685 0 : const Vector<Double>& SDBList::freqs() const {
686 :
687 0 : if (nSDB_==0)
688 0 : throw(AipsError("SDBList::freqs(): No SDBs in this SDBList yet."));
689 :
690 :
691 0 : if (nSDB_==1) {
692 : // Only one SDB, just return that one's freqs
693 0 : const Vector<Double>& f(SDB_[0]->freqs()); // from first SDB
694 0 : return f;
695 : }
696 :
697 : // Reach here, more than one SDB, need to gather info
698 0 : 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 0 : Int nchan(this->nChannels());
704 :
705 : // Will accumulate mean freqs here
706 0 : freqs_.resize(nchan);
707 0 : freqs_.set(0.0);
708 :
709 : // Average over SDBs, counting each spw exactly _once_
710 : // Map to keep track of unique spws
711 0 : std::set<Int> spws;
712 0 : Int nSpw(0);
713 0 : for (Int isdb=0;isdb<nSDB_;++isdb) {
714 0 : Int ispw=SDB_[isdb]->spectralWindow()(0);
715 0 : if (spws.count(ispw)<1) {
716 0 : freqs_+=SDB_[isdb]->freqs();
717 0 : spws.insert(ispw); // Record that we got this spw
718 0 : ++nSpw;
719 : }
720 : }
721 : // Divide by nSpw
722 0 : freqs_/=Double(nSpw);
723 :
724 0 : }
725 : // else: freqs_ already filled previously...
726 :
727 0 : return freqs_;
728 :
729 : }
730 :
731 :
732 : // ~Centroid frequency over all SDBs
733 0 : casacore::Double SDBList::centroidFreq() const {
734 :
735 : // Calculate a simple mean frequency, for now
736 0 : Double fsum(0.0f);
737 0 : Int nf(0);
738 0 : for (Int isdb=0;isdb<nSDB_;++isdb) {
739 0 : const Vector<Double> f(SDB_[isdb]->freqs());
740 0 : for (uInt ich=0;ich<f.nelements();++ich) {
741 0 : fsum+=f(ich);
742 0 : ++nf;
743 : }
744 0 : }
745 0 : return fsum/Double(nf);
746 : }
747 :
748 : // ~Centroid frequency over all SDBs
749 0 : casacore::Double SDBList::aggregateCentroidFreq() const {
750 :
751 : // Trap no data case
752 0 : if (nSDB_==0)
753 0 : throw(AipsError("SDBList::aggregateCentroidFreq(): No SDBs in this SDBList yet."));
754 :
755 : // Need to calculate?
756 0 : if (!aggCentroidFreqOK_) {
757 :
758 0 : if (nSDB_==1) {
759 : // from first and only SDB
760 0 : 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 0 : aggCentroidFreq_=0.0;
766 0 : for (Int isdb=0;isdb<nSDB_;++isdb)
767 0 : aggCentroidFreq_+=SDB_[isdb]->centroidFreq();
768 0 : aggCentroidFreq_/=Double(nSDB_);
769 : }
770 : // We've calculated it
771 0 : aggCentroidFreqOK_=true;
772 : }
773 :
774 : // Reach here, one way or another we have a good value, so return it
775 0 : return aggCentroidFreq_;
776 :
777 : }
778 :
779 0 : String SDBList::polBasis() const
780 : {
781 0 : String polBas(SDB_[0]->polBasis());
782 :
783 : // Trap non-uniformity, for now
784 0 : for (Int isdb=1;isdb<nSDB_;++isdb)
785 0 : AlwaysAssert((SDB_[isdb]->polBasis()==polBas),AipsError);
786 :
787 0 : return polBas;
788 0 : }
789 :
790 :
791 : // How many antennas?
792 : // Currently, this insists on uniformity over all SDBs
793 0 : Int SDBList::nAntennas() const {
794 :
795 0 : Int nAnt=SDB_[0]->nAntennas();
796 :
797 : // Trap non-uniformity, for now
798 0 : for (Int isdb=1;isdb<nSDB_;++isdb)
799 0 : AlwaysAssert((SDB_[isdb]->nAntennas()==nAnt),AipsError);
800 :
801 : // Reach here, then ok
802 0 : return nAnt;
803 :
804 : }
805 :
806 : // How many correlations?
807 : // Currently, this insists on uniformity over all SDBs
808 0 : Int SDBList::nCorrelations() const {
809 :
810 0 : Int nCorr=SDB_[0]->nCorrelations();
811 :
812 : // Trap non-uniformity, for now
813 0 : for (Int isdb=1;isdb<nSDB_;++isdb)
814 0 : AlwaysAssert((SDB_[isdb]->nCorrelations()==nCorr),AipsError);
815 :
816 : // Reach here, then ok
817 0 : return nCorr;
818 :
819 : }
820 :
821 :
822 :
823 :
824 :
825 0 : Bool SDBList::Ok() {
826 :
827 0 : for (Int i=0;i<nSDB_;++i)
828 0 : if (SDB_[i]->Ok()) return True;
829 :
830 : // If we get here, either no SDBs, or none have non-zero weight.
831 0 : 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 0 : void SDBList::enforceAPonData(const String& apmode)
844 : {
845 0 : for (Int i=0;i<nSDB_;++i)
846 0 : SDB_[i]->enforceAPonData(apmode);
847 0 : }
848 0 : void SDBList::enforceSolveWeights(const Bool pHandOnly)
849 : {
850 0 : for (Int i=0;i<nSDB_;++i)
851 0 : SDB_[i]->enforceSolveWeights(pHandOnly);
852 0 : }
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 0 : Int SDBList::extendCrossHandBaselineFlags(String& message)
896 : {
897 :
898 :
899 : // This enforces uniform nAnt in all SDBs
900 0 : Int nAnt(this->nAntennas());
901 0 : Int nBln=nAnt*(nAnt-1)/2;
902 :
903 : // channelized flags per baseline, cross-corr pair
904 0 : 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 0 : 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 0 : while (isdb0<nSDB_ && nfalse(SDB_[isdb0]->flagCube())<1)
918 0 : ++isdb0;
919 :
920 : // Trap all bad data case
921 0 : if (isdb0==nSDB_)
922 0 : throw(AipsError("Accumulated data entirely flagged."));
923 :
924 0 : Cube<Bool>& flc(SDB_[isdb0]->flagCube());
925 :
926 0 : const Vector<Int>& a1(SDB_[isdb0]->antenna1());
927 0 : const Vector<Int>& a2(SDB_[isdb0]->antenna2());
928 :
929 0 : Int nRows(SDB_[isdb0]->nRows());
930 0 : Int nChan(SDB_[isdb0]->nChannels());
931 :
932 0 : for (Int irow=0;irow<nRows;++irow) {
933 0 : for (Int ichan=0;ichan<nChan;++ichan) {
934 0 : 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 0 : Int nBadSDB(0);
942 0 : for (Int isdb=0;isdb<nSDB_;++isdb) {
943 :
944 : // Skip the one we initialized with
945 0 : if (isdb==isdb0)
946 0 : continue;
947 :
948 0 : Cube<Bool>& flc(SDB_[isdb]->flagCube());
949 0 : const Vector<Int>& a1(SDB_[isdb]->antenna1());
950 0 : const Vector<Int>& a2(SDB_[isdb]->antenna2());
951 :
952 0 : Int ngood=nfalse(flc);
953 0 : if (ngood==0) {
954 0 : nBadSDB+=1;
955 0 : continue;
956 : }
957 :
958 0 : Int nRows(SDB_[isdb]->nRows());
959 0 : Int nChan(SDB_[isdb]->nChannels());
960 :
961 0 : for (Int irow=0;irow<nRows;++irow) {
962 0 : for (Int ichan=0;ichan<nChan;++ichan) {
963 0 : 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 0 : Vector<size_t> goodBL(partialNFalse(aggFlag,IPosition(2,1,2)));
972 :
973 : // Median number of good baselines, over channel
974 0 : size_t medGoodBL=median(goodBL);
975 :
976 : // Count good channels (those with at least median number of good baselines)
977 0 : Int goodChan=ntrue(goodBL>=medGoodBL);
978 :
979 : // Apply aggregate flags to each SDB
980 : //cout << "Flag sum = ";
981 0 : Int nGoodSDB(0);
982 0 : for (Int isdb=0;isdb<nSDB_;++isdb) {
983 0 : Cube<Bool>& flc(SDB_[isdb]->flagCube());
984 0 : const Vector<Int>& a1(SDB_[isdb]->antenna1());
985 0 : const Vector<Int>& a2(SDB_[isdb]->antenna2());
986 :
987 0 : Int nRows(SDB_[isdb]->nRows());
988 0 : Int nChan(SDB_[isdb]->nChannels());
989 :
990 0 : for (Int irow=0;irow<nRows;++irow) {
991 0 : for (Int ichan=0;ichan<nChan;++ichan) {
992 0 : if (aggFlag(ichan,a1(irow),a2(irow))) {
993 0 : flc(Slice(1,2,1),Slice(ichan),Slice(irow))=true;
994 : }
995 : }
996 : }
997 0 : if (nfalse(flc)>0) ++nGoodSDB;
998 :
999 : //cout << ntrue(flc) << " ";
1000 : }
1001 : //cout << endl;
1002 :
1003 : // Assemble message
1004 0 : ostringstream ostr;
1005 0 : ostr << "There is usable (unflagged) data in " << nGoodSDB << " (of " << nSDB_ << " total) data segments, with "
1006 0 : << goodChan << " good (of " << nChannels() << " total) channels (" << floor(10000.0*Double(goodChan)/Double(nChannels()))/100.0 << "%) having at least "
1007 0 : << medGoodBL << " good (of " << nBln << " total) baselines (" << floor(10000.0*Double(medGoodBL)/Double(nBln))/100.0 <<"%).";
1008 : // cout << ostr << endl;
1009 0 : message=ostr;
1010 :
1011 :
1012 : // Return number of surviving SDBs...
1013 0 : return nGoodSDB;
1014 :
1015 0 : }
1016 :
1017 : } //# NAMESPACE CASA - END
1018 :
|
