Line data Source code
1 : //# SimpleComponentFTMachine.cc: Implementation of SimpleComponentFTMachine class
2 : //# Copyright (C) 1997,1998,1999,2000,2001
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/TransformMachines2/SimpleComponentFTMachine.h>
29 : #include <casacore/scimath/Mathematics/RigidVector.h>
30 : #include <components/ComponentModels/ComponentShape.h>
31 : #include <components/ComponentModels/ComponentList.h>
32 : #include <components/ComponentModels/ComponentType.h>
33 : #include <components/ComponentModels/Flux.h>
34 : #include <components/ComponentModels/SkyComponent.h>
35 : #include <components/ComponentModels/SpectralModel.h>
36 : #include <msvis/MSVis/VisBuffer2.h>
37 : #include <casacore/ms/MeasurementSets/MSColumns.h>
38 : #include <casacore/ms/MeasurementSets/MSIter.h>
39 : #include <casacore/casa/Arrays/ArrayMath.h>
40 : #include <casacore/casa/Arrays/ArrayLogical.h>
41 : #include <casacore/casa/Arrays/Cube.h>
42 : #include <casacore/casa/Arrays/Matrix.h>
43 : #include <casacore/casa/Arrays/Vector.h>
44 : #include <casacore/casa/Arrays/IPosition.h>
45 : #include <casacore/casa/BasicSL/Complex.h>
46 : #include <casacore/casa/BasicSL/Constants.h>
47 : #include <casacore/measures/Measures/UVWMachine.h>
48 : #ifdef _OPENMP
49 : #include <omp.h>
50 : #endif
51 :
52 :
53 : using namespace casacore;
54 : namespace casa { //# NAMESPACE CASA - BEGIN
55 :
56 : namespace refim { // namespace for refactor
57 : using namespace casacore;
58 : using namespace casa;
59 : using namespace casacore;
60 : using namespace casa::refim;
61 : using namespace casacore;
62 : using namespace casa::vi;
63 :
64 0 : void SimpleComponentFTMachine::get(VisBuffer2& vb, SkyComponent& component,
65 : Int row)
66 : {
67 : // If row is -1 then we pass through all rows
68 : uInt startRow, endRow;
69 0 : if (row < 0) {
70 0 : startRow = 0;
71 0 : endRow = vb.nRows() - 1;
72 : } else {
73 0 : startRow = endRow = row;
74 : }
75 0 : const uInt nRow = endRow - startRow + 1;
76 :
77 0 : component.flux().convertUnit(Unit("Jy"));
78 : // Rotate the uvw
79 0 : Matrix<Double> uvw(3, nRow);
80 0 : Vector<Double> dphase(nRow);
81 : {
82 : //const Vector<RigidVector<Double,3> >& uvwBuff = vb.uvw();
83 0 : for (uInt i = startRow, n = 0; i <= endRow; i++, n++) {
84 : //const RigidVector<Double,3>& uvwValue = uvwBuff(i);
85 : //NEGATING THE UV to be consitent with the GridFT::get
86 0 : for (Int idim = 0; idim < 2; idim++) {
87 0 : uvw(idim, n) = -vb.uvw()(idim, n);
88 : }
89 0 : uvw(2, n) = vb.uvw()(2,n);
90 : }
91 0 : rotateUVW(uvw, dphase, vb, component.shape().refDirection());
92 0 : dphase *= -C::_2pi;
93 : }
94 0 : uInt npol=vb.nCorrelations();
95 0 : uInt nChan=vb.nChannels();
96 0 : Cube<Complex> modelData;
97 0 : modelData.reference(vb.visCubeModel());
98 0 : Vector<Complex> visibility(4);
99 : //UNUSED: Double phase;
100 : //UNUSED: Double phaseMult;
101 0 : Vector<Double> frequency = vb.getFrequencies(0);
102 0 : Vector<Double> invLambda = frequency/C::c;
103 :
104 : // Find the offsets in polarization.
105 0 : Vector<Int> corrType(vb.getCorrelationTypesSelected().nelements());
106 0 : convertArray(corrType, vb.getCorrelationTypesSelected());
107 : {
108 0 : Int startPol = corrType(0);
109 0 : if((startPol > 4) && (startPol < 9)) {
110 0 : component.flux().convertPol(ComponentType::CIRCULAR);
111 0 : corrType -= 5;
112 : }
113 0 : else if((startPol > 8) && (startPol < 13)) {
114 0 : component.flux().convertPol(ComponentType::LINEAR);
115 0 : corrType -= 9;
116 : }
117 : else {
118 0 : component.flux().convertPol(ComponentType::STOKES);
119 0 : corrType -= startPol;
120 : }
121 : }
122 :
123 :
124 0 : uInt npart=1;
125 : #ifdef _OPENMP
126 0 : npart= numthreads_p <0 ? omp_get_max_threads() : min(numthreads_p, omp_get_max_threads());
127 : #endif
128 0 : if((nRow/npart)==0) npart=1;
129 0 : Block<Cube<DComplex> > dVisp(npart);
130 0 : Vector<uInt> nRowp(npart);
131 0 : Block<Matrix<Double> > uvwp(npart);
132 0 : Block<SkyComponent> compp(npart);
133 0 : Block<Int> startrow(npart);
134 0 : startrow[0]=0;
135 :
136 0 : nRowp.set(nRow/npart);
137 0 : nRowp(npart-1) += nRow%npart;
138 0 : Int sumrow=0;
139 0 : Record lala;
140 0 : String err;
141 0 : component.toRecord(err, lala);
142 0 : for (uInt k=0; k < npart; ++k){
143 0 : dVisp[k].resize(4, nChan, nRowp(k));
144 0 : uvwp[k].resize(3, nRowp(k));
145 0 : uvwp[k]=uvw(IPosition(2,0, sumrow ), IPosition(2, 2, sumrow+nRowp(k)-1));
146 0 : compp[k].fromRecord(err, lala);
147 0 : sumrow+=nRowp(k);
148 0 : startrow[k]=0;
149 0 : for (uInt j=0; j < k; ++j) startrow[k]+=nRowp(j);
150 : }
151 :
152 :
153 : //#pragma omp parallel default(none) firstprivate(npart) shared(frequency,dVisp, uvwp, compp) num_threads(npart)
154 0 : #pragma omp parallel firstprivate(npart) shared(dVisp, uvwp, compp) num_threads(npart)
155 : {
156 :
157 : #pragma omp for
158 : for (Int k=0; k < Int(npart); ++k){
159 : //Cube<DComplex> dVis(4, nChan, nRow);
160 : compp[k].visibility(dVisp[k], uvwp[k], frequency);
161 : }
162 :
163 : }
164 :
165 : // Loop over all rows
166 0 : sumrow=0;
167 :
168 : Bool isCopy;
169 0 : Complex *modData=modelData.getStorage(isCopy);
170 :
171 0 : #pragma omp parallel default(none) firstprivate(npart, npol, nChan, modData, corrType, nRowp, dphase, invLambda) shared(startrow, dVisp) num_threads(npart)
172 : {
173 : #pragma omp for
174 : for (Int k = 0; k < Int(npart); ++k) {
175 :
176 :
177 :
178 : applyPhasor(k, startrow, nRowp, dphase, invLambda, npol, nChan, corrType, dVisp[k], modData);
179 :
180 : }
181 :
182 : }
183 :
184 0 : modelData.putStorage(modData, isCopy);
185 :
186 0 : }
187 :
188 0 : void SimpleComponentFTMachine::applyPhasor(Int part, const Block<Int>& startrow, const Vector<uInt>& nRowp, const Vector<Double>& dphase, const Vector<Double>& invLambda, const Int npol, const Int nchan, const Vector<Int>& corrType, const Cube<DComplex>& dVis, Complex*& modData){
189 :
190 : Int r;
191 : Int rowoff;
192 : Double phaseMult;
193 : Double phase;
194 0 : for (uInt j=0; j< nRowp[part]; ++j){
195 0 : r=startrow[part]+j;
196 0 : rowoff=r*nchan*npol;
197 0 : phaseMult = dphase(r);
198 0 : for (Int chn = 0; chn < nchan; chn++) {
199 0 : phase = phaseMult * invLambda(chn);
200 0 : Complex phasor(cos(phase), sin(phase));
201 0 : for (Int pol=0; pol < npol; pol++) {
202 0 : const DComplex& val = dVis(corrType(pol), chn, j);
203 0 : modData[rowoff+chn*npol+pol] = Complex(val.real(), val.imag()) * conj(phasor);
204 : }
205 : }
206 : }
207 :
208 :
209 :
210 0 : }
211 :
212 :
213 640 : void SimpleComponentFTMachine::get(VisBuffer2& vb, const ComponentList& compList,
214 : Int row)
215 : {
216 : // If row is -1 then we pass through all rows
217 : uInt startRow, endRow;
218 640 : if (row < 0) {
219 640 : startRow = 0;
220 640 : endRow = vb.nRows() - 1;
221 : } else {
222 0 : startRow = endRow = row;
223 : }
224 640 : const uInt nRow = endRow - startRow + 1;
225 :
226 : // Rotate the uvw
227 640 : Matrix<Double> uvw0(3, nRow);
228 640 : Matrix<Double> uvw(3, nRow);
229 640 : Vector<Double> dphase(nRow);
230 :
231 : //const Vector<RigidVector<Double,3> >& uvwBuff = vb.uvw();
232 226915 : for (uInt i = startRow, n = 0; i <= endRow; i++, n++) {
233 : //const RigidVector<Double,3>& uvwValue = uvwBuff(i);
234 : //NEGATING the UV to be consitent with GridFT::get
235 678825 : for (Int idim = 0; idim < 2; idim++) {
236 452550 : uvw(idim, n) = -vb.uvw()(idim,n);
237 : }
238 226275 : uvw(2, n) = vb.uvw()(2,n);
239 : }
240 640 : uvw0=uvw.copy();
241 :
242 640 : uInt ncomponents=compList.nelements();
243 640 : uInt npol=vb.nCorrelations();
244 640 : uInt nChan=vb.nChannels();
245 640 : Cube<Complex> modelData;
246 640 : modelData.reference(vb.visCubeModel());
247 640 : modelData=0.0;
248 :
249 640 : Vector<Complex> visibility(4);
250 : //UNUSED: Double phase;
251 : //UNUSED: Double phaseMult;
252 640 : Vector<Double> frequency;
253 640 : frequency= vb.getFrequencies(0);
254 640 : Vector<Double> invLambda = frequency/C::c;
255 : // Find the offsets in polarization.
256 640 : Vector<Int> corrTypeL(vb.getCorrelationTypesSelected().nelements());
257 640 : convertArray(corrTypeL, vb.getCorrelationTypesSelected());
258 640 : Vector<Int> corrTypeC = corrTypeL.copy();
259 640 : Vector<Int> corrType = corrTypeL.copy();
260 640 : corrTypeL -= 9;
261 640 : corrTypeC -= 5;
262 640 : Cube<DComplex> dVis(4, nChan, nRow);
263 640 : ComponentType::Polarisation poltype=ComponentType::CIRCULAR;
264 640 : if(anyGT(Int(Stokes::RR), corrType)){
265 0 : poltype=ComponentType::STOKES;
266 0 : corrType = corrType-1;
267 : }
268 640 : else if(vb.polarizationFrame()==MSIter::Linear) {
269 0 : poltype=ComponentType::LINEAR;
270 0 : corrType = corrTypeL;
271 : } else {
272 640 : poltype=ComponentType::CIRCULAR;
273 640 : corrType = corrTypeC;
274 : }
275 :
276 640 : uInt npart=1;
277 : #ifdef _OPENMP
278 640 : npart= numthreads_p <0 ? omp_get_max_threads() : min(numthreads_p, omp_get_max_threads());
279 : #endif
280 :
281 :
282 :
283 640 : if((nRow/npart)==0) npart=1;
284 640 : Block<Cube<DComplex> > dVisp(npart);
285 640 : Vector<uInt> nRowp(npart);
286 640 : Block<Matrix<Double> > uvwp(npart);
287 640 : Block<ComponentList> compp(npart);
288 640 : Block<MeasFrame> mFramep(npart);
289 640 : Block<MDirection> mDir(npart);
290 640 : Block<uInt> startrow(npart);
291 :
292 640 : nRowp.set(nRow/npart);
293 640 : nRowp(npart-1) += nRow%npart;
294 :
295 :
296 640 : Block<Vector<Double> > dphasecomp(ncomponents);
297 640 : Block<Matrix<Double> > uvwcomp(ncomponents);
298 640 : Block<Block<Vector<Double> > > dphasecomps(npart);
299 640 : Block<Block<Matrix<Double> > > uvwcomps(npart);
300 21120 : for (uInt jj=0; jj < npart; ++jj){
301 20480 : dphasecomps[jj].resize(ncomponents);
302 20480 : uvwcomps[jj].resize(ncomponents);
303 : }
304 : ///Have to do this in one thread as MeasFrame and thus uvwmachine is not thread safe
305 1280 : for (uInt k=0; k < ncomponents; ++k){
306 640 : uvwcomp[k]=uvw;
307 640 : dphasecomp[k].resize(nRow);
308 640 : rotateUVW(uvwcomp[k], dphasecomp[k], vb, compList.component(k).shape().refDirection());
309 640 : dphasecomp[k] *= -C::_2pi;
310 640 : Int sumrow=0;
311 21120 : for (uInt jj=0; jj < npart; ++jj){
312 20480 : dphasecomps[jj][k]=dphasecomp[k](IPosition(1,sumrow), IPosition(1, sumrow+nRowp(jj)-1));
313 20480 : uvwcomps[jj][k]=uvwcomp[k](IPosition(2,0, sumrow ), IPosition(2, 2, sumrow+nRowp(jj)-1));
314 20480 : sumrow+=nRowp(jj);
315 : }
316 :
317 : }
318 :
319 :
320 :
321 : //UNUSED: Int sumrow=0;
322 640 : Record lala;
323 640 : String err;
324 640 : compList.toRecord(err, lala);
325 21120 : for (uInt k=0; k < npart; ++k){
326 20480 : compp[k].fromRecord(err, lala);
327 20480 : startrow[k]=0;
328 337920 : for (uInt j=0; j < k; ++j) startrow[k]+=nRowp(j);
329 : }
330 : Bool isCopy;
331 640 : Complex *modData=modelData.getStorage(isCopy);
332 640 : MDirection dataMDir=vb.phaseCenter();
333 :
334 640 : #pragma omp parallel default(none) firstprivate(npart, npol, nChan, modData, corrType, nRowp, invLambda, frequency, poltype) shared(startrow, compp, uvwcomps, dphasecomps) num_threads(npart)
335 : {
336 : #pragma omp for
337 :
338 : for (Int k=0; k < Int(npart); ++k){
339 :
340 : predictVis(modData, invLambda, frequency, compp[k], poltype, corrType,
341 : startrow[k], nRowp[k], nChan, npol, uvwcomps[k], dphasecomps[k]);
342 :
343 : }
344 :
345 :
346 : }
347 640 : modelData.putStorage(modData, isCopy);
348 :
349 640 : }
350 :
351 :
352 20480 : void SimpleComponentFTMachine::predictVis(Complex*& modData, const Vector<Double>& invLambda,
353 : const Vector<Double>& frequency, const ComponentList& compList,
354 : ComponentType::Polarisation poltype, const Vector<Int>& corrType,
355 : const uInt startrow, const uInt nrows, const uInt nchan, const uInt npol, const Block<Matrix<Double> > & uvwcomp, const Block<Vector<Double> > & dphasecomp){
356 20480 : Cube<DComplex> dVis(4, nchan, nrows);
357 20480 : uInt ncomponents=compList.nelements();
358 20480 : Vector<Double> dphase(nrows);
359 :
360 40960 : for (uInt icomp=0;icomp<ncomponents; ++icomp) {
361 20480 : SkyComponent component=compList.component(icomp);
362 20480 : component.flux().convertUnit(Unit("Jy"));
363 20480 : component.flux().convertPol(poltype);
364 20480 : MDirection compdir=component.shape().refDirection();
365 20480 : Vector<Double> thisRow(3);
366 20480 : thisRow=0.0;
367 : //UNUSED: uInt i;
368 20480 : Matrix<Double> uvw;
369 20480 : uvw=uvwcomp[icomp];
370 20480 : Vector<Double> dphase;
371 20480 : dphase=dphasecomp[icomp];
372 20480 : component.visibility(dVis, uvw, frequency);
373 : Double phaseMult;
374 : Double phase;
375 : uInt realrow;
376 : //// Loop over all rows
377 246755 : for (uInt r = 0; r < nrows ; r++) {
378 226275 : phaseMult = dphase(r);
379 226275 : realrow=r+startrow;
380 22275948 : for (uInt chn = 0; chn < nchan; chn++) {
381 22049673 : phase = phaseMult * invLambda(chn);
382 22049673 : Complex phasor(cos(phase), sin(phase));
383 108543963 : for (uInt pol=0; pol < npol; pol++) {
384 86494290 : const DComplex& val = dVis(corrType(pol), chn, r);
385 86494290 : modData[realrow*npol*nchan+chn*npol + pol] += Complex(val.real(), val.imag()) * conj(phasor);
386 : }
387 : }
388 : }
389 :
390 :
391 20480 : }
392 :
393 20480 : }
394 :
395 :
396 : }// end namespace refim
397 : } //# NAMESPACE CASA - END
398 :
|
