Line data Source code
1 : //# CalVisBuffer.cc: Extends VisBuffer for calibration purposes
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 <msvis/MSVis/CalVisBuffer.h>
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <casacore/casa/Exceptions/Error.h>
31 : #include <casacore/casa/Utilities/Assert.h>
32 :
33 : using namespace casacore;
34 : namespace casa { //# NAMESPACE CASA - BEGIN
35 :
36 250 : CalVisBuffer::CalVisBuffer() :
37 : VisBuffer(),
38 250 : focusChan_p(-1),
39 : // infocusFlagCube_p(),
40 250 : infocusFlag_p(),
41 250 : infocusVisCube_p(),
42 250 : infocusModelVisCube_p(),
43 250 : residuals_p(),
44 250 : residFlag_p(),
45 500 : diffResiduals_p()
46 250 : {}
47 :
48 0 : CalVisBuffer::CalVisBuffer(ROVisibilityIterator& iter) :
49 : VisBuffer(iter),
50 0 : focusChan_p(-1),
51 : // infocusFlagCube_p(),
52 0 : infocusFlag_p(),
53 0 : infocusVisCube_p(),
54 0 : infocusModelVisCube_p(),
55 0 : residuals_p(),
56 0 : residFlag_p(),
57 0 : diffResiduals_p()
58 :
59 0 : {}
60 :
61 0 : CalVisBuffer::CalVisBuffer(const CalVisBuffer& vb) :
62 : VisBuffer(vb),
63 0 : focusChan_p(-1),
64 : // infocusFlagCube_p(),
65 0 : infocusFlag_p(),
66 0 : infocusVisCube_p(),
67 0 : infocusModelVisCube_p(),
68 0 : residuals_p(),
69 0 : residFlag_p(),
70 0 : diffResiduals_p()
71 0 : {}
72 :
73 250 : CalVisBuffer::~CalVisBuffer()
74 250 : {}
75 :
76 0 : CalVisBuffer& CalVisBuffer::operator=(const VisBuffer& other)
77 : {
78 :
79 0 : if (this!=&other) {
80 : // delete workspaces
81 0 : this->cleanUp();
82 : // call parent
83 0 : VisBuffer::operator=(other);
84 : }
85 :
86 0 : return *this;
87 : }
88 :
89 250 : CalVisBuffer& CalVisBuffer::assign(const VisBuffer& other, Bool copy)
90 : {
91 : // delete workspaces
92 250 : this->cleanUp();
93 :
94 : // Call parent:
95 250 : VisBuffer::assign(other,copy);
96 :
97 250 : return *this;
98 :
99 : }
100 :
101 0 : void CalVisBuffer::updateCoordInfo(const VisBuffer * /*vb=NULL*/, const Bool /*dirDependent=true*/)
102 : {
103 : // Just do the nominally non-row-dep values
104 0 : arrayId();
105 0 : fieldId();
106 0 : spectralWindow();
107 0 : nCorr();
108 0 : nChannel();
109 0 : frequency();
110 :
111 0 : }
112 :
113 124 : void CalVisBuffer::enforceAPonData(const String& apmode)
114 : {
115 :
116 : // ONLY if something to do
117 124 : if (apmode=="A" || apmode=="P") {
118 0 : Int nCor(nCorr());
119 0 : Float amp(1.0);
120 0 : Complex cor(1.0);
121 0 : Bool *flR=flagRow().data();
122 0 : Bool *fl =flag().data();
123 0 : Vector<Float> ampCorr(nCor);
124 0 : Vector<Int> n(nCor,0);
125 0 : for (Int irow=0;irow<nRow();++irow,++flR) {
126 0 : if (!flagRow()(irow)) {
127 0 : ampCorr=0.0f;
128 0 : n=0;
129 0 : for (Int ich=0;ich<nChannel();++ich,++fl) {
130 0 : if (!flag()(ich,irow)) {
131 0 : for (Int icorr=0;icorr<nCor;icorr++) {
132 :
133 0 : amp=abs(visCube()(icorr,ich,irow));
134 0 : if (amp>0.0f) {
135 :
136 0 : if (apmode=="P") {
137 : // we will scale by amp to make data phase-only
138 0 : cor=Complex(amp,0.0);
139 : // keep track for weight adjustment
140 0 : ampCorr(icorr)+=abs(cor);
141 0 : n(icorr)++;
142 : }
143 0 : else if (apmode=="A")
144 : // we will scale by "phase" to make data amp-only
145 0 : cor=visCube()(icorr,ich,irow)/amp;
146 :
147 : // Apply the complex scaling
148 0 : visCube()(icorr,ich,irow)/=cor;
149 : }
150 : } // icorr
151 : } // !*fl
152 : } // ich
153 : // Make appropriate weight adjustment
154 : // (only required for phase-only, since only it rescales data)
155 0 : if (apmode=="P") {
156 0 : for (Int icorr=0;icorr<nCor;icorr++)
157 0 : if (n(icorr)>0)
158 : // weights adjusted by square of the mean(amp)
159 0 : weightMat()(icorr,irow)*=square(ampCorr(icorr)/Float(n(icorr)));
160 : else
161 : // weights now zero
162 0 : weightMat()(icorr,irow)=0.0f;
163 : }
164 : } // !*flR
165 : } // irow
166 :
167 0 : } // phase- or amp-only
168 :
169 : // cout << "amp(visCube())=" << amplitude(visCube().reform(IPosition(1,visCube().nelements()))) << endl;
170 :
171 124 : }
172 :
173 :
174 :
175 0 : void CalVisBuffer::setFocusChan(const Int focusChan)
176 : {
177 : // Nominally focus on the whole data array
178 0 : IPosition focusblc(3,0,0,0);
179 0 : IPosition focustrc(visCube().shape());
180 0 : focustrc-=1;
181 :
182 : // if focusChan non-negative, select the single channel
183 0 : if (focusChan>-1)
184 0 : focusblc(1)=focustrc(1)=focusChan;
185 :
186 : // infocusFlagCube_p.reference(flagCube()(focusblc,focustrc));
187 0 : infocusFlag_p.reference(flag()(focusblc.getLast(2),focustrc.getLast(2)));
188 0 : infocusVisCube_p.reference(visCube()(focusblc,focustrc));
189 0 : infocusModelVisCube_p.reference(modelVisCube()(focusblc,focustrc));
190 :
191 : // Remember current in-focus channel
192 0 : focusChan_p=focusChan;
193 :
194 0 : }
195 :
196 0 : void CalVisBuffer::sizeResiduals(const Int& nPar,
197 : const Int& nDiff)
198 : {
199 :
200 0 : IPosition ip1(visCube().shape());
201 0 : if (focusChan_p>-1)
202 0 : ip1(1)=1;
203 0 : residuals_p.resize(ip1);
204 0 : residuals_p.set(0.0);
205 :
206 0 : if (nPar>0 && nDiff>0) {
207 0 : IPosition ip2(5,ip1(0),nPar,ip1(1),ip1(2),nDiff);
208 0 : diffResiduals_p.resize(ip2);
209 0 : diffResiduals_p.set(0.0);
210 0 : }
211 :
212 0 : }
213 :
214 0 : void CalVisBuffer::initResidWithModel()
215 : {
216 : // Copy (literally) the in-focus model to the residual workspace
217 : // TBD: obey flags here!?
218 : // TBD: weights and flagCube...
219 0 : residuals_p = infocusModelVisCube_p;
220 :
221 : // TBD: should probably move this to setFocusChan, so that
222 : // we can optimize handling of flags better (e.g., prior to repeated
223 : // calls to SVC.differentiate, etc.) (initResidWithModel should
224 : // be viewed as just _refreshing_ the residuals with the model
225 : // data for a new trial corrupt)
226 0 : if (focusChan_p>-1) {
227 : // copy flags so they are contiguous
228 0 : residFlag_p.resize(infocusFlag_p.shape());
229 0 : residFlag_p=infocusFlag_p;
230 : }
231 : else
232 : // just reference the whole (contiguous) infocusFlag array
233 0 : residFlag_p.reference(infocusFlag_p);
234 :
235 : // Ensure contiguity, because CalSolver will depend on this
236 0 : AlwaysAssert(residFlag_p.contiguousStorage(),AipsError);
237 0 : AlwaysAssert(residuals_p.contiguousStorage(),AipsError);
238 :
239 0 : }
240 :
241 0 : void CalVisBuffer::finalizeResiduals()
242 : {
243 : // Subtract in-focus obs data from residuals workspace
244 0 : residuals_p -= infocusVisCube_p;
245 :
246 : // TBD: zero flagged samples here?
247 :
248 0 : }
249 :
250 250 : void CalVisBuffer::cleanUp()
251 : {
252 :
253 : // Zero-size all workspaces
254 : // infocusFlagCube_p.resize();
255 250 : infocusFlag_p.resize();
256 250 : infocusVisCube_p.resize();
257 250 : infocusModelVisCube_p.resize();
258 :
259 250 : residuals_p.resize();
260 250 : residFlag_p.resize();
261 250 : diffResiduals_p.resize();
262 :
263 250 : }
264 :
265 :
266 :
267 : } //# NAMESPACE CASA - END
268 :
|
