Line data Source code
1 : //# VisChunkAverager.cc: Implementation of VisChunkAverager.h
2 : //# Copyright (C) 2010
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 : //----------------------------------------------------------------------------
27 :
28 : #include <msvis/MSVis/VisChunkAverager.h>
29 :
30 : #include <casacore/casa/Arrays/ArrayPartMath.h>
31 :
32 : using namespace casacore;
33 : namespace casa { //# NAMESPACE CASA - BEGIN
34 :
35 : //----------------------------------------------------------------------------
36 :
37 0 : VisChunkAverager::VisChunkAverager(const Vector<MS::PredefinedColumns>& dataCols,
38 : const Bool doSpWeight,
39 0 : const Vector<Matrix<Int> >& chBounds)
40 0 : : colEnums_p(dataCols),
41 0 : doSpWeight_p(doSpWeight),
42 0 : chanAveBounds_p(chBounds),
43 0 : readyToHash_p(false),
44 0 : haveHashMap_p(false)
45 : {
46 0 : reset();
47 0 : }
48 :
49 : //----------------------------------------------------------------------------
50 :
51 0 : VisChunkAverager::~VisChunkAverager()
52 : {
53 : // Null default destructor
54 0 : }
55 :
56 : //----------------------------------------------------------------------------
57 :
58 0 : void VisChunkAverager::reset()
59 : {
60 0 : haveHashMap_p = false;
61 0 : sphash_to_inprows_p.clear();
62 0 : readyToHash_p = false;
63 0 : }
64 :
65 : // read all the columns into vb.
66 0 : void VisChunkAverager::fill_vb(VisBuffer& vb)
67 : {
68 0 : vb.time();
69 : // timeExtraPrec?
70 0 : vb.antenna1();
71 0 : vb.antenna2();
72 : // antenna3?
73 0 : vb.feed1();
74 0 : vb.feed2();
75 : // feed3?
76 : // ddid implicitly read during the implicit sort.
77 : // So shouldn't fieldId have been read as well?
78 0 : vb.processorId();
79 : // phaseId?
80 0 : vb.fieldId();
81 0 : vb.timeInterval();
82 0 : vb.exposure();
83 0 : vb.timeCentroid();
84 : // pulsar_bin?
85 : // pulsar_gate_id?
86 0 : vb.scan();
87 0 : vb.arrayId();
88 0 : vb.observationId();
89 0 : vb.stateId();
90 : // baseline_ref?
91 0 : vb.uvwMat();
92 : // vb.uvw2();
93 0 : for(Int i = colEnums_p.nelements(); i--;){
94 0 : if(colEnums_p[i] == MS::CORRECTED_DATA)
95 0 : vb.correctedVisCube();
96 0 : else if(colEnums_p[i] == MS::MODEL_DATA)
97 0 : vb.modelVisCube();
98 : // else if(colEnums_p[i] == MS::LAG_DATA)
99 : // VisBuffer doesn't handle LAG_DATA.
100 0 : else if(colEnums_p[i] == MS::FLOAT_DATA)
101 0 : vb.floatDataCube();
102 0 : else if(colEnums_p[i] == MS::DATA)
103 0 : vb.visCube();
104 : }
105 : // video_point?
106 0 : vb.sigmaMat();
107 0 : vb.weightMat();
108 0 : if(doSpWeight_p)
109 0 : vb.weightSpectrum();
110 0 : vb.flagCube();
111 0 : vb.flagCategory();
112 0 : vb.flagRow();
113 0 : }
114 :
115 0 : Bool VisChunkAverager::setupHashFunction(ROVisibilityIterator& vi)
116 : {
117 0 : if(readyToHash_p) // No-op.
118 0 : return true;
119 :
120 0 : VisBuffer vb(vi);
121 :
122 0 : maxant1p1_p = 1;
123 0 : maxant2p1_p = 1;
124 0 : maxfeed1p1_p = 1;
125 0 : maxfeed2p1_p = 1;
126 : // maxprocp1_p = max(vb.processorId()) + 1; Not needed
127 :
128 0 : for(vi.origin(); vi.more(); ++vi){
129 0 : maxant1p1_p = uIntMax(maxant1p1_p, 1 + max(vb.antenna1()));
130 0 : maxant2p1_p = uIntMax(maxant2p1_p, max(vb.antenna2()) + 1);
131 0 : maxfeed1p1_p = uIntMax(maxfeed1p1_p, max(vb.feed1()) + 1);
132 0 : maxfeed2p1_p = uIntMax(maxfeed2p1_p, max(vb.feed2()) + 1);
133 : }
134 0 : vi.origin(); // Go back to the start of the chunk.
135 0 : return true;
136 0 : }
137 :
138 0 : VisBuffer& VisChunkAverager::average(ROVisibilityIterator& vi)
139 : {
140 : using casacore::operator*;
141 :
142 : // Just in case findCollision() returned before doing it.
143 : // makeHashMap will return right away if it can.
144 0 : haveHashMap_p = makeHashMap(vi);
145 :
146 0 : VisBuffer vb(vi);
147 0 : uInt chunkletNum = 0;
148 0 : Bool firstValidOutRowInChunk = true;
149 0 : Vector<Bool> firstrowinslot(sphash_to_inprows_p.size());
150 0 : firstrowinslot.set(true);
151 :
152 0 : initialize(vb);
153 :
154 : Double minTime, maxTime, firstinterval, lastinterval;
155 0 : Int nrows = vb.nRow();
156 : // Paranoid initialization.
157 0 : minTime = vb.time()[0];
158 0 : firstinterval = vb.timeInterval()[0];
159 0 : maxTime = vb.time()[nrows - 1];
160 0 : lastinterval = vb.timeInterval()[nrows - 1];
161 :
162 0 : for(vi.origin(); vi.more(); ++vi){
163 : // Preload the things that need to be channel averaged.
164 0 : for(uInt colind = 0; colind < colEnums_p.nelements(); ++colind){
165 0 : if(colEnums_p[colind] == MS::DATA)
166 0 : vb.visCube();
167 0 : else if(colEnums_p[colind] == MS::MODEL_DATA)
168 0 : vb.modelVisCube();
169 0 : else if(colEnums_p[colind] == MS::CORRECTED_DATA)
170 0 : vb.correctedVisCube();
171 0 : else if(colEnums_p[colind] == MS::DATA)
172 0 : vb.floatDataCube();
173 : }
174 :
175 : // The flags and weights are already loaded by this point, UNLESS the
176 : // row flag was true for all the rows. Make sure they're loaded, or
177 : // they could end up with the wrong shape.
178 0 : vb.flagCube();
179 0 : if(vi.existsWeightSpectrum())
180 0 : vb.weightSpectrum();
181 0 : vb.weightMat();
182 0 : vb.sigmaMat();
183 0 : if(nCat_p > 0)
184 0 : vb.flagCategory();
185 :
186 0 : vb.channelAve(chanAveBounds_p[vi.spectralWindow()]);
187 :
188 : // Handle zeros
189 0 : vb.sigmaMat()(vb.sigmaMat()<FLT_EPSILON)=FLT_MAX;
190 0 : vb.weightMat()(vb.weightMat()<FLT_EPSILON)=FLT_EPSILON;
191 :
192 : // First iterate through the *unflagged* rows of the current VisBuffer.
193 0 : Int outrow = 0;
194 0 : Bool firstValidOutRowInIntegration = true;
195 :
196 0 : mapuIvIType::iterator sphend = sphash_to_inprows_p.end();
197 0 : for(mapuIvIType::iterator sphit = sphash_to_inprows_p.begin();
198 0 : sphit != sphend; ++sphit){
199 0 : Int inrow = sphit->second[chunkletNum];
200 :
201 0 : if(inrow >= 0){
202 0 : if(firstValidOutRowInIntegration){
203 0 : Double time = vb.time()[inrow];
204 :
205 0 : firstValidOutRowInIntegration = false;
206 :
207 : // Some initialization for the chunk
208 0 : if(firstValidOutRowInChunk){
209 0 : firstValidOutRowInChunk = false;
210 0 : minTime = time;
211 0 : maxTime = minTime;
212 0 : firstinterval = vb.timeInterval()[inrow];
213 0 : lastinterval = firstinterval;
214 : }
215 :
216 0 : if(time < minTime){
217 0 : minTime = time;
218 0 : firstinterval = max(vb.timeInterval()[inrow], 0.0);
219 : }
220 0 : else if(time > maxTime){
221 0 : maxTime = time;
222 0 : lastinterval = max(vb.timeInterval()[inrow], 0.0);
223 : }
224 : }
225 :
226 : // Add the VisBuffer row to the current accumulation
227 :
228 0 : if(firstrowinslot[outrow]){
229 0 : firstrowinslot[outrow] = false;
230 :
231 0 : avBuf_p.antenna1()[outrow] = vb.antenna1()[inrow];
232 0 : avBuf_p.antenna2()[outrow] = vb.antenna2()[inrow];
233 0 : avBuf_p.feed1()[outrow] = vb.feed1()[inrow];
234 0 : avBuf_p.feed2()[outrow] = vb.feed2()[inrow];
235 0 : avBuf_p.observationId()[outrow] = vb.observationId()[inrow];
236 0 : avBuf_p.processorId()[outrow] = vb.processorId()[inrow];
237 0 : avBuf_p.scan()[outrow] = vb.scan()[inrow];
238 0 : avBuf_p.stateId()[outrow] = vb.stateId()[inrow];
239 : }
240 :
241 0 : Vector<Float> wtM;
242 0 : if (anyEQ(colEnums_p,MS::CORRECTED_DATA))
243 : // use weightMat for corrected data
244 0 : wtM.reference(vb.weightMat().column(inrow));
245 : else
246 : // otherwise use sigmaMat (data, model)
247 0 : wtM=Float(1.)/square(vb.sigmaMat().column(inrow));
248 :
249 : // Accumulate the visibilities and weights, and set the flags.
250 : // For better or worse, calibration may adjust WEIGHT without adjusting
251 : // WEIGHT_SPECTRUM. This means that at least for now weightSpectrum
252 : // should be treated as a unnormalized relative quantity between
253 : // channels, and the overall weight is given by WEIGHT.
254 : // Note that avBuf_p.weightSpectrum() is used to accumulate the weights
255 : // regardless of vb.existsWeightSpectrum(), so its absolute magnitude
256 : // does matter (i.e. it holds the relative weights between one
257 : // integration and the next).
258 0 : Double totwt = 0.0; // Total weight for inrow.
259 0 : for(Int cor = 0; cor < nCorr_p; ++cor){
260 0 : Double constwtperchan = wtM[cor]; // Already adjusted by flagging.
261 :
262 0 : Double totwtsp = wtM[cor];
263 0 : Bool useWtSp = doSpWeight_p;
264 0 : if(doSpWeight_p){
265 0 : totwtsp = 0.0;
266 0 : for(Int ochan = 0; ochan < nChan_p; ++ochan)
267 0 : totwtsp += vb.weightSpectrum()(cor, ochan, inrow);
268 0 : useWtSp = totwtsp > 0.0;
269 : }
270 :
271 0 : Double unflaggedWt = 0.0;
272 0 : for(Int ochan = 0; ochan < nChan_p; ++ochan){
273 0 : if(totwtsp > 0.0){
274 0 : for(Int cat = 0; cat < nCat_p; ++cat){
275 0 : avBuf_p.flagCategory()(IPosition(4, cor, ochan, cat, outrow)) &=
276 0 : vb.flagCategory()(IPosition(4, cor, ochan, cat, inrow));
277 : }
278 : }
279 0 : if(!vb.flagCube()(cor, ochan, inrow)){
280 0 : Double wt = useWtSp ?
281 0 : wtM[cor] * vb.weightSpectrum()(cor, ochan, inrow) / totwtsp
282 0 : : constwtperchan;
283 :
284 0 : if(wt > 0.0){
285 0 : avBuf_p.weightSpectrum()(cor, ochan, outrow) += wt;
286 0 : avBuf_p.flagCube()(cor, ochan, outrow) = false;
287 :
288 0 : for(Int i = colEnums_p.nelements(); i--;){
289 0 : if(colEnums_p[i] == MS::CORRECTED_DATA)
290 0 : avBuf_p.correctedVisCube()(cor, ochan, outrow) +=
291 0 : (wt * vb.correctedVisCube()(cor, ochan, inrow));
292 0 : else if(colEnums_p[i] == MS::MODEL_DATA)
293 0 : avBuf_p.modelVisCube()(cor, ochan, outrow) +=
294 0 : (wt * vb.modelVisCube()(cor, ochan, inrow));
295 : // else if(colEnums_p[i] == MS::LAG_DATA)
296 : // VisBuffer doesn't handle LAG_DATA
297 0 : else if(colEnums_p[i] == MS::FLOAT_DATA)
298 0 : avBuf_p.floatDataCube()(cor, ochan, outrow) +=
299 0 : (wt * vb.floatDataCube()(cor, ochan, inrow));
300 0 : else if(colEnums_p[i] == MS::DATA)
301 0 : avBuf_p.visCube()(cor, ochan, outrow) +=
302 0 : (wt * vb.visCube()(cor, ochan, inrow));
303 : }
304 0 : unflaggedWt += wt;
305 : }
306 : }
307 : }
308 0 : avBuf_p.weightMat()(cor, outrow) += unflaggedWt;
309 0 : Double sig = vb.sigmaMat()(cor, inrow);
310 0 : avBuf_p.sigmaMat()(cor, outrow) += unflaggedWt * unflaggedWt * sig * sig;
311 0 : totwt += unflaggedWt;
312 : }
313 :
314 0 : Double totslotwt = 0.0;
315 :
316 0 : for(Int cor = 0; cor < nCorr_p; ++cor)
317 0 : totslotwt += avBuf_p.weightMat()(cor, outrow);
318 :
319 : // totwt > 0.0 implies totslotwt > 0.0, as required by
320 : // the running averages (mandatory for timeCentroid!).
321 0 : if(totwt > 0.0){
322 0 : Double leverage = avBuf_p.timeCentroid()[outrow] == 0.0 ? 1.0 :
323 0 : totwt / totslotwt; // Technically == 1 for the 1st input
324 : // integration, but numerical error is a real problem.
325 :
326 : // UVW (weighted only by weight for now)
327 : // // Iterate through slots.
328 : // for(ui2vmap::iterator slotit = bin_slots_p[tbn].begin();
329 : // slotit != bin_slots_p[tbn].end(); ++slotit){
330 : // for(uivector::iterator toikit = slotv.begin();
331 : // toikit != slotv.end(); ++toikit){
332 : // Double wv = 0.0;
333 : // Array<Complex>::const_iterator dataEnd = data_toikit.end();
334 : // for(Array<Complex>::const_iterator dit = data_toikit.begin();
335 : // dit != dataEnd; ++dit)
336 : // wv += fabs(*dit);
337 : // if(wv > 0.0){
338 : // swv += wv;
339 : // outUVW += (wv / swv) * (inUVW(*toikit) - outUVW);
340 : // }
341 : // } // End of loop through the slot's rows.
342 0 : for(uInt ax = 0; ax < 3; ++ax)
343 0 : avBuf_p.uvwMat()(ax, outrow) += leverage * (vb.uvwMat()(ax, inrow) -
344 0 : avBuf_p.uvwMat()(ax, outrow));
345 :
346 0 : avBuf_p.timeCentroid()[outrow] += leverage *
347 0 : (vb.timeCentroid()[inrow] - avBuf_p.timeCentroid()[outrow]);
348 0 : avBuf_p.exposure()[outrow] += totwt * vb.exposure()[inrow];
349 : }
350 0 : else if(totslotwt == 0.0){ // Put in a representative UVW.
351 0 : for(uInt ax = 0; ax < 3; ++ax)
352 0 : avBuf_p.uvwMat()(ax, outrow) = vb.uvwMat()(ax, inrow);
353 : }
354 0 : }
355 0 : ++outrow;
356 : } // End of loop over sphit for chunkletNum.
357 0 : ++chunkletNum;
358 : } // End of loop over chunkletNums (integrations) in vi's current chunk.
359 :
360 :
361 : // Now go back and see if there were any flagged earlier minTimes or later
362 : // maxTimes. Remember that although TIME and INTERVAL should not care
363 : // about flagging, sphash_to_inprows_p does. BUT, don't let RIDICULOUS
364 : // flagged times and intervals contaminate the average.
365 : //
366 : // ridiculous means off from the unflagged values by > 40 years.
367 : // That's a generous but maybe too wide margin. The most common bad values
368 : // seem to be 0s (= 1858 A.D.) and twice the current epoch.
369 : // const Double ridiculous = 1262304000.0;
370 : // Double acc_minTime = minTime - ridiculous;
371 : // Double acc_maxTime = maxTime + ridiculous;
372 : // for(vi.origin(); vi.more(); ++vi){
373 : // Double time = vb.time()[0];
374 :
375 : // if(time < minTime && time > acc_minTime){
376 : // minTime = time;
377 : // firstinterval = max(vb.timeInterval()[0], 0.0);
378 : // }
379 : // else if(time > maxTime && time < acc_maxTime){
380 : // maxTime = time;
381 : // lastinterval = max(vb.timeInterval()[0], 0.0);
382 : // }
383 : // }
384 :
385 0 : normalize(minTime, maxTime, firstinterval, lastinterval);
386 :
387 0 : return avBuf_p;
388 0 : }
389 :
390 0 : void VisChunkAverager::initialize(VisBuffer& vb)
391 : {
392 : // Initialize the averaging buffer.
393 : // Output to private data:
394 : // avBuf_p VisBuffer Averaging buffer
395 : // readyToHash_p Bool sets to true
396 : // nChan_p Int # of channels
397 : // nCorr_p Int # of polarizations
398 :
399 0 : fill_vb(vb);
400 :
401 : // Assign to vb to establish a two-way connection to the underlying VisIter.
402 : // Assign main meta info only
403 : // avBuf_p.assign(vb,true);
404 0 : avBuf_p.assign(vb, false);
405 0 : avBuf_p.updateCoordInfo(); // This is the (simplified) CalVisBuffer version!
406 :
407 : // Immutables:
408 0 : nChan_p = vb.nChannel();
409 0 : nCorr_p = vb.corrType().nelements();
410 0 : nCat_p = vb.flagCategory().nelements() > 0 ? vb.flagCategory().shape()(2) : 0;
411 :
412 : // More or less VisBuffAccumulator::initialize() up to
413 : // "// Fill in the antenna numbers for all rows"
414 0 : Int nRow = sphash_to_inprows_p.size();
415 :
416 0 : avBuf_p.nRow() = nRow;
417 0 : avBuf_p.nChannel() = nChan_p;
418 :
419 : // Resize and initialize the VisBuffer columns used here
420 0 : avBuf_p.antenna1().resize(nRow);
421 0 : avBuf_p.antenna2().resize(nRow);
422 : // avBuf_p.arrayId() is an Int, not a Vector.
423 0 : for(Int i = colEnums_p.nelements(); i--;){
424 0 : if(colEnums_p[i] == MS::CORRECTED_DATA)
425 0 : avBuf_p.correctedVisCube().resize(nCorr_p, nChan_p, nRow);
426 0 : else if(colEnums_p[i] == MS::MODEL_DATA)
427 0 : avBuf_p.modelVisCube().resize(nCorr_p, nChan_p, nRow);
428 : // else if(colEnums_p[i] == MS::LAG_DATA)
429 : // VisBuffer doesn't handle LAG_DATA
430 0 : else if(colEnums_p[i] == MS::FLOAT_DATA)
431 0 : avBuf_p.floatDataCube().resize(nCorr_p, nChan_p, nRow);
432 0 : else if(colEnums_p[i] == MS::DATA)
433 0 : avBuf_p.visCube().resize(nCorr_p, nChan_p, nRow);
434 : }
435 : // avBuf_p.dataDescriptionId() is an Int, not a Vector.
436 0 : avBuf_p.exposure().resize(nRow);
437 0 : avBuf_p.feed1().resize(nRow);
438 0 : avBuf_p.feed2().resize(nRow);
439 : // avBuf_p.fieldId() is an Int, not a Vector.
440 0 : avBuf_p.flagCube().resize(nCorr_p, nChan_p, nRow);
441 0 : avBuf_p.flagRow().resize(nRow);
442 :
443 0 : if(nCat_p > 0){
444 0 : avBuf_p.flagCategory().resize(IPosition(4, nCorr_p, nChan_p, nCat_p, nRow));
445 0 : avBuf_p.flagCategory() = true;
446 : }
447 :
448 0 : avBuf_p.observationId().resize(nRow);
449 0 : avBuf_p.processorId().resize(nRow);
450 0 : avBuf_p.scan().resize(nRow);
451 :
452 : // Use sigmaMat() throughout, not the correlation-averaged sigma().
453 0 : avBuf_p.sigmaMat().resize(nCorr_p, nRow);
454 :
455 0 : avBuf_p.stateId().resize(nRow);
456 0 : avBuf_p.time().resize(nRow);
457 0 : avBuf_p.timeCentroid().resize(nRow);
458 0 : avBuf_p.timeInterval().resize(nRow);
459 0 : avBuf_p.uvwMat().resize(3, nRow);
460 :
461 : // Use weightMat() throughout, not the correlation-averaged weight().
462 0 : avBuf_p.weightMat().resize(nCorr_p, nRow);
463 :
464 : // Because of channel dependent flagging, avBuf_p's weight spectrum gets used
465 : // for normalization regardless of vb.existsWeightSpectrum().
466 0 : avBuf_p.weightSpectrum().resize(nCorr_p, nChan_p, nRow);
467 :
468 : // Fill in the antenna numbers for all rows.
469 : // Do not assume a constant number of antennas like VisBuffAccumulator.
470 :
471 : // Initialize everything else
472 0 : const Complex czero(0.0);
473 0 : for(Int i = colEnums_p.nelements(); i--;){
474 0 : if(colEnums_p[i] == MS::CORRECTED_DATA)
475 0 : avBuf_p.correctedVisCube() = czero;
476 0 : else if(colEnums_p[i] == MS::MODEL_DATA)
477 0 : avBuf_p.modelVisCube() = czero;
478 : // else if(colEnums_p[i] == MS::LAG_DATA)
479 : // VisBuffer doesn't handle LAG_DATA.
480 0 : else if(colEnums_p[i] == MS::FLOAT_DATA)
481 0 : avBuf_p.floatDataCube() = 0.0;
482 0 : else if(colEnums_p[i] == MS::DATA)
483 0 : avBuf_p.visCube() = czero;
484 : }
485 0 : avBuf_p.exposure() = 0.0;
486 0 : avBuf_p.flagCube() = true;
487 0 : avBuf_p.flagRow() = false;
488 0 : avBuf_p.scan() = 0; // This should be unnecessary.
489 0 : avBuf_p.time() = 0.0; // This should be unnecessary.
490 0 : avBuf_p.timeInterval() = 0.0; // This should be unnecessary.
491 0 : avBuf_p.timeCentroid() = 0.0;
492 0 : avBuf_p.uvwMat() = 0.0;
493 0 : avBuf_p.weightMat() = 0.0;
494 0 : avBuf_p.sigmaMat() = 0.0;
495 0 : avBuf_p.weightSpectrum() = 0.0;
496 0 : }
497 :
498 : //----------------------------------------------------------------------------
499 :
500 0 : void VisChunkAverager::normalize(const Double minTime, const Double maxTime,
501 : const Double firstinterval,
502 : const Double lastinterval)
503 : {
504 : // Normalize the current accumulation timeBin
505 : // Output to private data:
506 : // avBuf_p VisBuffer& Averaged buffer
507 : //
508 :
509 : // These columns are independently normalized:
510 : // DATA_DESC_ID, UVW, TIME_CENTROID, WEIGHT,
511 : // WEIGHT_SPECTRUM, FLAG
512 :
513 : // Things that are the same throughout avBuf_p.
514 : // A short running average (for numerical stability).
515 0 : avBuf_p.time() = minTime + 0.5 * (maxTime - minTime);
516 :
517 0 : avBuf_p.timeInterval() = maxTime - minTime +
518 0 : 0.5 * (firstinterval + lastinterval);
519 :
520 : // Already done by initialize()'s avBuf_p.assign(vb, false);
521 : // avBuf_p.dataDescriptionId() = vb.dataDescriptionId();
522 :
523 0 : uInt outrow = 0;
524 0 : mapuIvIType::iterator sphend = sphash_to_inprows_p.end();
525 :
526 0 : for(mapuIvIType::iterator sphit = sphash_to_inprows_p.begin();
527 0 : sphit != sphend; ++sphit){
528 : // Divide by the weights
529 0 : Vector<Float> wtM(avBuf_p.weightMat().column(outrow));
530 0 : Double wt = 0.0; // total weight for the output row.
531 :
532 0 : for(Int cor = 0; cor < nCorr_p; ++cor){
533 0 : Float orw = wtM[cor];
534 :
535 : // sigma is the channel-averaged sigma for a single chan.
536 0 : avBuf_p.sigmaMat()(cor, outrow) = orw > 0.0 ?
537 0 : sqrt(avBuf_p.sigmaMat()(cor, outrow)) / orw
538 : : -1.0; // Seems safer than 0.0.
539 0 : wt += orw;
540 : }
541 :
542 0 : if(wt == 0.0){
543 0 : avBuf_p.flagRow()[outrow] = true;
544 :
545 : // Looks better than 0 (-> 1858).
546 0 : avBuf_p.timeCentroid()[outrow] = avBuf_p.time()[outrow];
547 : }
548 : else{
549 0 : uInt nrowsinslot = 0;
550 0 : std::vector<Int>::iterator vend = sphit->second.end();
551 :
552 0 : for(std::vector<Int>::iterator vit = sphit->second.begin();
553 0 : vit != vend; ++vit){
554 0 : if(*vit >= 0)
555 0 : ++nrowsinslot;
556 : }
557 :
558 0 : avBuf_p.exposure()[outrow] *= nrowsinslot / wt;
559 : }
560 :
561 0 : if(!avBuf_p.flagRow()[outrow]){
562 :
563 0 : Matrix<Float> rowWeightSpectrum(avBuf_p.weightSpectrum().xyPlane(outrow));
564 0 : Vector<Float> rowWeightMat(avBuf_p.weightMat().column(outrow));
565 0 : Vector<Float> rowSigmaMat(avBuf_p.sigmaMat().column(outrow));
566 :
567 : // extract chan-indep weight/sigma
568 0 : rowWeightMat = partialMedians(rowWeightSpectrum,IPosition(1,1));
569 :
570 0 : for(Int cor = 0; cor < nCorr_p; ++cor){
571 :
572 0 : if (rowWeightMat[cor]<FLT_EPSILON)
573 0 : rowWeightMat[cor]=FLT_EPSILON;
574 0 : rowSigmaMat[cor]=1/sqrt(rowWeightMat[cor]);
575 :
576 0 : for(Int ochan = 0; ochan < avBuf_p.nChannel(); ++ochan){
577 0 : if(!avBuf_p.flagCube()(cor, ochan, outrow)){
578 0 : Double w = avBuf_p.weightSpectrum()(cor, ochan, outrow);
579 :
580 : // The 2nd choice begs the question of why it isn't flagged.
581 0 : Double norm = w > 0.0 ? 1.0 / w : 0.0;
582 :
583 0 : for(Int i = colEnums_p.nelements(); i--;){
584 0 : if(colEnums_p[i] == MS::CORRECTED_DATA)
585 0 : avBuf_p.correctedVisCube()(cor, ochan, outrow) *= norm;
586 0 : else if(colEnums_p[i] == MS::MODEL_DATA)
587 0 : avBuf_p.modelVisCube()(cor, ochan, outrow) *= norm;
588 : // else if(colEnums_p[i] == MS::LAG_DATA)
589 : // VisBuffer doesn't handle LAG_DATA.
590 0 : else if(colEnums_p[i] == MS::FLOAT_DATA)
591 0 : avBuf_p.floatDataCube()(cor, ochan, outrow) *= norm;
592 0 : else if(colEnums_p[i] == MS::DATA)
593 0 : avBuf_p.visCube()(cor, ochan, outrow) *= norm;
594 : }
595 : } // ends if(!avBuf_p.flagCube()(cor, ochan, outrow))
596 : } // ends chan loop
597 : } // ends cor loop
598 0 : } // ends if(!avBuf_p.flagRow()(outrow))
599 0 : ++outrow;
600 0 : } // ends sphit loop
601 0 : }
602 :
603 0 : uInt VisChunkAverager::hashFunction(const Int ant1, const Int ant2,
604 : const Int feed1, const Int feed2,
605 : const Int procid) const
606 : {
607 0 : return feed2 +
608 0 : maxfeed2p1_p * (feed1 +
609 0 : maxfeed1p1_p * (ant2 +
610 0 : maxant2p1_p * (ant1 +
611 0 : maxant1p1_p * procid)));
612 : };
613 :
614 0 : Bool VisChunkAverager::makeHashMap(ROVisibilityIterator& vi)
615 : {
616 0 : if(haveHashMap_p) // No-op
617 0 : return true;
618 :
619 0 : readyToHash_p = setupHashFunction(vi);
620 0 : if(!readyToHash_p)
621 0 : return false;
622 :
623 0 : sphash_to_inprows_p.clear();
624 :
625 0 : uInt nChunklets = vi.nSubInterval(); // The # of vb's in vi's current chunk.
626 0 : uInt chunkletNum = 0;
627 0 : VisBuffer vb(vi);
628 :
629 0 : for(vi.origin(); vi.more(); ++vi){
630 0 : uInt nrow = vb.nRow();
631 :
632 0 : for(uInt rn = 0; rn < nrow; ++rn){
633 0 : if(!vb.flagRow()[rn] && !allTrue(vb.flagCube().xyPlane(rn))){
634 : // sphash = sparse hash
635 0 : uInt sphash = hashFunction(vb.antenna1()[rn], vb.antenna2()[rn],
636 0 : vb.feed1()[rn], vb.feed2()[rn],
637 0 : vb.processorId()[rn]);
638 :
639 : /* jagonzal: These insertions cause a seg. fault in MAC 10.00 using llvm and C++11
640 : if(sphash_to_inprows_p.count(sphash) < 1){
641 : sphash_to_inprows_p.insert(poshint,
642 : std::pair<uInt, std::vector<Int> >(sphash,
643 : std::vector<Int>()));
644 : ++poshint;
645 :
646 : // Fill the new vector of input rows with "the skip value".
647 : sphash_to_inprows_p[sphash].assign(nChunklets, -1);
648 : }
649 : */
650 :
651 0 : if(sphash_to_inprows_p.count(sphash) < 1){
652 0 : sphash_to_inprows_p[sphash]=std::vector<Int>(nChunklets,-1);
653 : }
654 :
655 0 : sphash_to_inprows_p[sphash][chunkletNum] = rn;
656 :
657 : }
658 : }
659 0 : ++chunkletNum;
660 : }
661 0 : vi.origin(); // Go back to the start of the chunk.
662 0 : return true;
663 0 : }
664 :
665 0 : Bool VisChunkAverager::check_chunk(ROVisibilityIterator& vi, Double& time_to_break,
666 : const Bool watch_obs, const Bool watch_scan,
667 : const Bool watch_state)
668 : {
669 : Bool foundbreak;
670 0 : Vector<MS::PredefinedColumns> dataCols;
671 : //VisChunkAverager vba(time_to_break, dataCols, false);
672 0 : VisChunkAverager vba(dataCols, false);
673 :
674 0 : foundbreak = vba.findCollision(vi, time_to_break, watch_obs,
675 : watch_scan, watch_state);
676 0 : return foundbreak;
677 0 : }
678 :
679 0 : Bool VisChunkAverager::findCollision(ROVisibilityIterator& vi, Double& time_to_break,
680 : const Bool watchObs, const Bool watchScan,
681 : const Bool watchState)
682 : {
683 0 : if(!(watchScan || watchState || watchObs))
684 0 : return false; // Nothing would count as a collision.
685 :
686 0 : if(vi.nRowChunk() < 2) // It takes two to collide, so
687 0 : return false; // return before accessing nonexistent rows.
688 :
689 0 : haveHashMap_p = makeHashMap(vi);
690 :
691 0 : if(!haveHashMap_p)
692 0 : throw(AipsError("VisChunkAverager::findCollision(): could not make a hash map"));
693 :
694 0 : Bool foundCollision = false;
695 0 : VisBuffer vb(vi);
696 0 : uInt chunkletNum = 0;
697 0 : uInt nslots = sphash_to_inprows_p.size();
698 0 : Vector<Int> firstScan; // The first scan for each sphash
699 0 : Vector<Int> firstState; // The first state for each sphash
700 0 : Vector<Int> firstObs; // The first obsID for each sphash
701 :
702 0 : if(watchScan)
703 0 : firstScan.resize(nslots);
704 0 : if(watchState)
705 0 : firstState.resize(nslots);
706 0 : if(watchObs)
707 0 : firstObs.resize(nslots);
708 :
709 0 : for(vi.origin(); vi.more(); ++vi){
710 0 : if(foundCollision) // The last chunklet did it.
711 0 : break;
712 :
713 : // We want the soonest time_to_break in all the keys. Default
714 : // time_to_break to the last time in the latest chunklet, assuming vb is
715 : // time sorted.
716 0 : time_to_break = vb.time()[vb.nRow() - 1];
717 :
718 0 : uInt slotnum = 0;
719 0 : for(mapuIvIType::const_iterator keyit = sphash_to_inprows_p.begin();
720 0 : keyit != sphash_to_inprows_p.end(); ++keyit){
721 0 : Int row = keyit->second[chunkletNum];
722 :
723 0 : if(row >= 0){ // If the chunklet has this sphash...
724 : // Don't collapse this into an || chain - time_to_break needs to be evaluated
725 : // even if a collision has already been found.
726 0 : if(watchScan)
727 0 : foundCollision |= checkForBreak(firstScan, vb.scan()[row], slotnum,
728 0 : chunkletNum, keyit->second);
729 0 : if(watchState)
730 0 : foundCollision |= checkForBreak(firstState, vb.stateId()[row], slotnum,
731 0 : chunkletNum, keyit->second);
732 : //if(watchArray) arrayIDs are already separated by the chunking.
733 0 : if(watchObs)
734 0 : foundCollision |= checkForBreak(firstObs, vb.observationId()[row],
735 0 : slotnum, chunkletNum, keyit->second);
736 :
737 0 : if(foundCollision){
738 0 : time_to_break = vb.time()[row];
739 0 : break;
740 : }
741 : }
742 0 : ++slotnum;
743 : }
744 0 : ++chunkletNum;
745 : }
746 0 : vi.origin(); // Go back to the start of the chunk.
747 0 : return foundCollision;
748 0 : }
749 :
750 0 : Bool VisChunkAverager::checkForBreak(Vector<Int>& firstVals, const Int i,
751 : const uInt slotnum,
752 : const uInt chunkletNum,
753 : const std::vector<Int>& inrows_for_slot) const
754 : {
755 0 : if(chunkletNum == 0 || inrows_for_slot[chunkletNum - 1] < 0)
756 0 : firstVals[slotnum] = i;
757 0 : return i != firstVals[slotnum];
758 : }
759 :
760 : } //# NAMESPACE CASA - END
761 :
|
