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: 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$
      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           0 : 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           0 :   if (row < 0) {
     219           0 :     startRow = 0;
     220           0 :     endRow = vb.nRows() - 1;
     221             :   } else {
     222           0 :     startRow = endRow = row;
     223             :   }
     224           0 :   const uInt nRow = endRow - startRow + 1;
     225             :    
     226             :    // Rotate the uvw
     227           0 :   Matrix<Double> uvw0(3, nRow);
     228           0 :   Matrix<Double> uvw(3, nRow); 
     229           0 :   Vector<Double> dphase(nRow);
     230             : 
     231             :   //const Vector<RigidVector<Double,3> >& uvwBuff = vb.uvw();
     232           0 :   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           0 :     for (Int idim = 0; idim < 2; idim++) {
     236           0 :       uvw(idim, n) = -vb.uvw()(idim,n);
     237             :     }
     238           0 :     uvw(2, n) = vb.uvw()(2,n);
     239             :   }
     240           0 :   uvw0=uvw.copy();
     241             :   
     242           0 :   uInt ncomponents=compList.nelements();
     243           0 :   uInt npol=vb.nCorrelations();
     244           0 :   uInt nChan=vb.nChannels();
     245           0 :   Cube<Complex> modelData;
     246           0 :   modelData.reference(vb.visCubeModel());
     247           0 :   modelData=0.0;
     248             : 
     249           0 :   Vector<Complex> visibility(4);
     250             :   //UNUSED: Double phase;
     251             :   //UNUSED: Double phaseMult;
     252           0 :   Vector<Double> frequency;
     253           0 :   frequency= vb.getFrequencies(0);
     254           0 :   Vector<Double> invLambda = frequency/C::c;
     255             :   // Find the offsets in polarization. 
     256           0 :   Vector<Int> corrTypeL(vb.getCorrelationTypesSelected().nelements());
     257           0 :   convertArray(corrTypeL, vb.getCorrelationTypesSelected());
     258           0 :   Vector<Int> corrTypeC = corrTypeL.copy();
     259           0 :   Vector<Int> corrType = corrTypeL.copy();
     260           0 :   corrTypeL -= 9;
     261           0 :   corrTypeC -= 5;
     262           0 :   Cube<DComplex> dVis(4, nChan, nRow);
     263           0 :   ComponentType::Polarisation poltype=ComponentType::CIRCULAR;
     264           0 :   if(anyGT(Int(Stokes::RR), corrType)){
     265           0 :     poltype=ComponentType::STOKES;
     266           0 :     corrType = corrType-1;
     267             :   }
     268           0 :   else if(vb.polarizationFrame()==MSIter::Linear) {
     269           0 :     poltype=ComponentType::LINEAR;
     270           0 :     corrType = corrTypeL;
     271             :     } else {
     272           0 :     poltype=ComponentType::CIRCULAR;
     273           0 :     corrType = corrTypeC;
     274             :   }
     275             : 
     276           0 :  uInt npart=1;
     277             : #ifdef _OPENMP
     278           0 :  npart= numthreads_p <0 ? omp_get_max_threads() : min(numthreads_p, omp_get_max_threads());
     279             : #endif
     280             : 
     281             : 
     282             : 
     283           0 :   if((nRow/npart)==0) npart=1;
     284           0 :   Block<Cube<DComplex> > dVisp(npart);
     285           0 :   Vector<uInt> nRowp(npart);
     286           0 :   Block<Matrix<Double> > uvwp(npart);
     287           0 :   Block<ComponentList> compp(npart);
     288           0 :   Block<MeasFrame> mFramep(npart);
     289           0 :   Block<MDirection> mDir(npart);
     290           0 :   Block<uInt> startrow(npart);
     291             : 
     292           0 :   nRowp.set(nRow/npart);
     293           0 :   nRowp(npart-1) += nRow%npart;
     294             : 
     295             :   
     296           0 :   Block<Vector<Double> > dphasecomp(ncomponents);
     297           0 :   Block<Matrix<Double> > uvwcomp(ncomponents);
     298           0 :   Block<Block<Vector<Double> > > dphasecomps(npart);
     299           0 :   Block<Block<Matrix<Double> > > uvwcomps(npart);
     300           0 :   for (uInt jj=0; jj < npart; ++jj){
     301           0 :     dphasecomps[jj].resize(ncomponents);
     302           0 :     uvwcomps[jj].resize(ncomponents);
     303             :   }
     304             :   ///Have to do this in one thread as MeasFrame and thus uvwmachine is not thread safe
     305           0 :   for (uInt k=0; k < ncomponents; ++k){
     306           0 :     uvwcomp[k]=uvw;
     307           0 :     dphasecomp[k].resize(nRow);
     308           0 :     rotateUVW(uvwcomp[k], dphasecomp[k],  vb, compList.component(k).shape().refDirection());
     309           0 :     dphasecomp[k] *= -C::_2pi;
     310           0 :     Int sumrow=0;
     311           0 :     for (uInt jj=0; jj < npart; ++jj){
     312           0 :       dphasecomps[jj][k]=dphasecomp[k](IPosition(1,sumrow), IPosition(1, sumrow+nRowp(jj)-1));
     313           0 :       uvwcomps[jj][k]=uvwcomp[k](IPosition(2,0, sumrow  ), IPosition(2, 2, sumrow+nRowp(jj)-1));
     314           0 :       sumrow+=nRowp(jj);
     315             :     }
     316             :    
     317             :   }
     318             : 
     319             : 
     320             : 
     321             :   //UNUSED: Int sumrow=0;
     322           0 :   Record lala;
     323           0 :   String err;
     324           0 :   compList.toRecord(err, lala);
     325           0 :   for (uInt k=0; k < npart; ++k){
     326           0 :     compp[k].fromRecord(err, lala);
     327           0 :     startrow[k]=0;
     328           0 :     for (uInt j=0; j < k; ++j) startrow[k]+=nRowp(j);
     329             :   }
     330             :   Bool isCopy;
     331           0 :   Complex *modData=modelData.getStorage(isCopy);
     332           0 :   MDirection dataMDir=vb.phaseCenter();
     333             : 
     334           0 : #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           0 :   modelData.putStorage(modData, isCopy);
     348             :  
     349           0 : }
     350             : 
     351             : 
     352           0 :   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           0 :     Cube<DComplex> dVis(4, nchan, nrows);
     357           0 :     uInt ncomponents=compList.nelements();
     358           0 :     Vector<Double> dphase(nrows);
     359             :     
     360           0 :     for (uInt icomp=0;icomp<ncomponents; ++icomp) {
     361           0 :       SkyComponent component=compList.component(icomp);
     362           0 :       component.flux().convertUnit(Unit("Jy"));
     363           0 :       component.flux().convertPol(poltype);
     364           0 :       MDirection compdir=component.shape().refDirection();
     365           0 :       Vector<Double> thisRow(3);
     366           0 :       thisRow=0.0;
     367             :       //UNUSED: uInt i;
     368           0 :       Matrix<Double> uvw;
     369           0 :       uvw=uvwcomp[icomp];
     370           0 :       Vector<Double> dphase;
     371           0 :       dphase=dphasecomp[icomp];
     372           0 :       component.visibility(dVis, uvw, frequency);
     373             :       Double phaseMult;
     374             :       Double phase;
     375             :       uInt realrow;
     376             :       //// Loop over all rows
     377           0 :       for (uInt r = 0; r < nrows ; r++) {
     378           0 :         phaseMult = dphase(r);
     379           0 :         realrow=r+startrow;
     380           0 :         for (uInt chn = 0; chn < nchan; chn++) {
     381           0 :           phase = phaseMult * invLambda(chn);   
     382           0 :           Complex phasor(cos(phase), sin(phase));
     383           0 :           for (uInt pol=0; pol < npol; pol++) {
     384           0 :             const DComplex& val = dVis(corrType(pol), chn, r);
     385           0 :             modData[realrow*npol*nchan+chn*npol + pol] += Complex(val.real(), val.imag()) * conj(phasor);
     386             :           }
     387             :         }
     388             :       }
     389             : 
     390             : 
     391             :     }
     392             :     
     393           0 :   }  
     394             : 
     395             : 
     396             : }// end namespace refim
     397             : } //# NAMESPACE CASA - END
     398             :