Line data Source code
1 : //# VisBuffGroupAcc.cc: Implementation of VisBuffGroupAcc.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: aips2-request@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: VisBuffAccumulator.cc,v 19.7 2006/01/17 08:22:27 gvandiep Exp $
27 : //----------------------------------------------------------------------------
28 :
29 : #include <msvis/MSVis/VisBuffGroupAcc.h>
30 : #include <casacore/ms/MeasurementSets/MSColumns.h>
31 : #include <casacore/ms/MSSel/MSSelection.h>
32 : #include <casacore/casa/Exceptions/Error.h>
33 :
34 : using namespace casacore;
35 : namespace casa { //# NAMESPACE CASA - BEGIN
36 :
37 : //----------------------------------------------------------------------------
38 :
39 39 : VisBuffGroupAcc::VisBuffGroupAcc(const Int& nAnt,
40 : const Int& nSpw,
41 : const Int& nFld,
42 : const Double& subinterval,
43 39 : const Bool fillModel)
44 39 : : nAnt_p(nAnt),
45 39 : nSpw_p(nSpw),
46 39 : nFld_p(nFld),
47 : nBuf_p(0),
48 39 : subinterval_p(subinterval),
49 : fillModel_p(fillModel),
50 : prenorm_p(false),
51 : globalTimeStamp_p(0.0),
52 : VBA_p(),
53 : spwfldids_p(nSpw,nFld,-1),
54 39 : tvi_debug(False)
55 : {
56 : // Interval MUST be strictly greater than zero
57 39 : if (subinterval_p < DBL_EPSILON)
58 0 : subinterval_p=0.1; // TBD: is this reasonable?
59 :
60 : // NB: Enforcing no prenormalization here, for now
61 :
62 39 : }
63 :
64 : //----------------------------------------------------------------------------
65 :
66 39 : VisBuffGroupAcc::~VisBuffGroupAcc()
67 : {
68 : // Null default destructor
69 : //
70 :
71 : // Delete all VBAs.
72 : // TBD: encapsulate this for more general use?
73 163 : for (Int i=0;i<nBuf_p;++i)
74 124 : if (VBA_p[i]) delete VBA_p[i];
75 :
76 39 : VBA_p.resize(0);
77 39 : }
78 :
79 0 : void VisBuffGroupAcc::clearChanMask(std::map<Int, Vector<Bool>*>& chanmask)
80 : {
81 0 : for(std::map<Int, Vector<Bool>*>::iterator it = chanmask.begin();
82 0 : it != chanmask.end(); ++it)
83 0 : if(it->second)
84 0 : delete it->second;
85 0 : chanmask.clear();
86 0 : }
87 :
88 0 : Bool VisBuffGroupAcc::fillChanMask(std::map<Int, Vector<Bool>*>& chanmask,
89 : const String& spwstr,
90 : const MeasurementSet& ms)
91 : {
92 0 : clearChanMask(chanmask);
93 :
94 0 : MSSelection mssel;
95 0 : mssel.setSpwExpr(spwstr != "" ? spwstr : "*");
96 0 : Matrix<Int> chansel = mssel.getChanList(&ms, 1);
97 0 : uInt nranges = chansel.nrow();
98 0 : Bool didSel = nranges > 0;
99 :
100 0 : if(didSel){
101 0 : MSSpWindowColumns spwCols(ms.spectralWindow());
102 0 : Vector<Int> nChan0 = spwCols.numChan().getColumn();
103 :
104 0 : Vector<Int> uspw(chansel.column(0));
105 0 : Vector<Int> ustart(chansel.column(1));
106 0 : Vector<Int> uend(chansel.column(2));
107 0 : Vector<Int> ustep(chansel.column(3));
108 :
109 0 : for(uInt rangenum = 0; rangenum < nranges; ++rangenum){
110 0 : Int spw = uspw[rangenum];
111 :
112 : // Initialize this spw mask, if necessary (def = masked)
113 0 : if(chanmask.count(spw) < 1)
114 0 : chanmask[spw] = new Vector<Bool>(nChan0[spw], true);
115 :
116 0 : Int unchan = uend[rangenum] - ustart[rangenum] + 1;
117 :
118 : // Update the mask (false = selected)
119 0 : (*chanmask[spw])(Slice(ustart[rangenum], unchan, ustep[rangenum])) = false;
120 : } // rangenum
121 : } // non-triv spw selection
122 0 : return didSel;
123 : }
124 :
125 0 : uInt VisBuffGroupAcc::applyChanMask(std::map<Int, Vector<Bool>*>& chanmask)
126 : {
127 0 : uInt naffected = 0;
128 :
129 0 : for(Int i = 0; i < nBuf_p; ++i){
130 0 : Int spw = VBA_p[i]->aveCalVisBuff().spectralWindow();
131 0 : if(chanmask.count(spw) > 0){
132 0 : Int chan0 = VBA_p[i]->aveCalVisBuff().channel()(0);
133 0 : Int nchan = VBA_p[i]->aveCalVisBuff().nChannel();
134 0 : if(sum((*(chanmask[spw]))(Slice(chan0, nchan))) > 0){
135 : // There are some channels to mask...
136 0 : Vector<Bool> fr(VBA_p[i]->aveCalVisBuff().flagRow());
137 0 : Matrix<Bool> f(VBA_p[i]->aveCalVisBuff().flag());
138 0 : Vector<Bool> fc;
139 0 : Vector<Bool> chm((*(chanmask[spw]))(Slice(chan0, nchan)));
140 0 : uInt nr = VBA_p[i]->aveCalVisBuff().nRow();
141 :
142 0 : for(uInt irow = 0; irow < nr; ++irow){
143 0 : if(!fr(irow)){
144 0 : fc.reference(f.column(irow));
145 0 : fc = fc || chm;
146 : }
147 : }
148 0 : ++naffected;
149 : }
150 : }
151 : }
152 0 : return naffected;
153 : }
154 :
155 : //----------------------------------------------------------------------------
156 :
157 124 : void VisBuffGroupAcc::accumulate (const VisBuffer& vb)
158 : {
159 124 : Int spw = vb.spectralWindow();
160 124 : Int fld = vb.fieldId();
161 124 : Int ibuf = spwfldids_p(spw, fld);
162 :
163 : // Create the new VisBuffAccumulator, if needed
164 124 : if (ibuf<0) {
165 124 : ibuf=nBuf_p;
166 124 : ++nBuf_p;
167 124 : VBA_p.resize(nBuf_p, false, true);
168 124 : spwfldids_p(spw, fld) = ibuf;
169 124 : VBA_p[ibuf] = new VisBuffAccumulator(nAnt_p, subinterval_p, prenorm_p, fillModel_p);
170 124 : VBA_p[ibuf]->setTVIDebug(tvi_debug);
171 : }
172 :
173 : // ibuf should be non-negative now
174 124 : if(ibuf < 0)
175 0 : throw(AipsError("VisBuffGroupAcc: VisBuffAccumulator index failure."));
176 :
177 : // Accumulate the vb into the correct accumulator
178 124 : VBA_p[ibuf]->accumulate(vb);
179 124 : }
180 :
181 : //----------------------------------------------------------------------------
182 :
183 39 : void VisBuffGroupAcc::finalizeAverage()
184 : {
185 39 : globalTimeStamp_p=0.0;
186 39 : Double globalTimeStampWt(0.0);
187 39 : Double t0(-1.0); // avoid round-off problems in time average
188 :
189 163 : for (Int ibuf=0;ibuf<nBuf_p;++ibuf) {
190 : // tell a VBA to finalize itself
191 124 : VBA_p[ibuf]->finalizeAverage();
192 :
193 : // Accumulate weighted timestamps
194 124 : Double& thistimewt(VBA_p[ibuf]->timeStampWt());
195 124 : if (thistimewt>0.0) {
196 124 : globalTimeStampWt+=(thistimewt);
197 124 : Double& thistime(VBA_p[ibuf]->timeStamp());
198 124 : if (t0<0.0)
199 39 : t0=thistime;
200 : else
201 85 : globalTimeStamp_p+=(thistimewt*(thistime-t0));
202 : }
203 : }
204 :
205 : // Form global timestamp
206 39 : if (globalTimeStampWt>0.0)
207 39 : globalTimeStamp_p/=globalTimeStampWt;
208 :
209 : // Add offset back in!
210 39 : globalTimeStamp_p+=t0;
211 :
212 : // NB: the per-VBA timestamp weights are approximately (exactly?)
213 : // the relative weights of the _data_ going into the solution.
214 : // This could be a useful log message?....
215 :
216 39 : };
217 :
218 39 : void VisBuffGroupAcc::enforceAPonData(const String& apmode)
219 : {
220 :
221 : // Delegate to each CalVisBuffer in turn
222 163 : for (Int ibuf=0;ibuf<nBuf_p;++ibuf)
223 124 : if (VBA_p[ibuf])
224 124 : this->operator()(ibuf).enforceAPonData(apmode);
225 :
226 39 : }
227 :
228 39 : void VisBuffGroupAcc::enforceSolveCorrWeights(const Bool phandonly)
229 : {
230 :
231 : // If requested set cross-hand weights to zero (if they exist):
232 39 : if(phandonly)
233 0 : for(Int ibuf = 0; ibuf < nBuf_p; ++ibuf)
234 0 : if (VBA_p[ibuf]) {
235 0 : CalVisBuffer& cvb(this->operator()(ibuf));
236 0 : if (cvb.nCorr() > 2)
237 0 : cvb.weightMat()(Slice(1, 2, 1), Slice()).set(0.0);
238 : }
239 :
240 39 : }
241 :
242 4646 : CalVisBuffer& VisBuffGroupAcc::operator()(const Int& buf)
243 : {
244 4646 : if (buf > -1 && buf < nBuf_p)
245 4646 : return VBA_p[buf]->aveCalVisBuff();
246 : else
247 0 : throw(AipsError("VisBuffGroupAcc: operator(buf) index out-of-range."));
248 : }
249 :
250 0 : CalVisBuffer& VisBuffGroupAcc::operator()(const Int& spw,const Int& fld)
251 : {
252 0 : if (spw>-1 && fld > -1 && spwfldids_p(spw,fld) > 0)
253 0 : return this->operator()(spwfldids_p(spw,fld));
254 : else
255 0 : throw(AipsError("VisBuffGroupAcc: operator(spw,fld) index out-of-range."));
256 : }
257 :
258 0 : const Vector<Int>& VisBuffGroupAcc::outToInRow(const Int buf,
259 : const Bool hurl) const
260 : {
261 0 : if(buf > -1 && buf < nBuf_p)
262 0 : return VBA_p[buf]->outToInRow(hurl);
263 : else
264 0 : throw(AipsError("VisBuffGroupAcc outToInRow: buf index out of range."));
265 : }
266 :
267 0 : const Vector<Int>& VisBuffGroupAcc::outToInRow(const Int spw,
268 : const Int fld,
269 : const Bool hurl) const
270 : {
271 0 : if(spw > -1 && fld > -1 && spwfldids_p(spw, fld) > 0)
272 0 : return outToInRow(spwfldids_p(spw, fld), hurl);
273 : else
274 : throw(AipsError(
275 0 : "VisBuffGroupAcc outToInRow: (spw, fld) index out of range."));
276 : }
277 :
278 0 : void VisBuffGroupAcc::reportData()
279 : {
280 0 : cout << "nBuf=" << nBuf_p << endl;
281 0 : for(Int ibuf = 0; ibuf < nBuf_p; ++ibuf) {
282 0 : cout << "iBuf=" << ibuf << endl;
283 0 : if (VBA_p[ibuf]) VBA_p[ibuf]->reportData();
284 : }
285 :
286 0 : }
287 :
288 : } //# NAMESPACE CASA - END
289 :
|
