Line data    Source code

       1             : //# FiltrationTVI.tcc: Template class for data filtering TVI
       2             : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
       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             : #ifndef _MSVIS_FILTRATIONTVI_TCC_
      27             : #define _MSVIS_FILTRATIONTVI_TCC_
      28             : 
      29             : #include <synthesis/MeasurementComponents/FiltrationTVI.h>
      30             : 
      31             : #include <climits>
      32             : 
      33             : #include <casacore/casa/Exceptions/Error.h>
      34             : #include <casacore/casa/Arrays/Array.h>
      35             : 
      36             : #include <msvis/MSVis/VisibilityIterator2.h>
      37             : 
      38             : using namespace casacore;
      39             : 
      40             : namespace {
      41             : template<class T>
      42       27036 : inline void FiltrateVector(Vector<T> const &feed,
      43             :     Vector<bool> const &is_filtrate, Vector<T> &filtrate) {
      44       27036 :   AlwaysAssert(feed.nelements() == is_filtrate.nelements(), AipsError);
      45             :   // filter_flag: true --> filtrate, false --> residue
      46       27036 :   auto const num_filtrates = ntrue(is_filtrate);
      47       27036 :   if (num_filtrates == feed.nelements()) {
      48       27036 :     filtrate.resize();
      49       27036 :     filtrate.reference(feed);
      50             :   } else {
      51           0 :     filtrate.resize(ntrue(is_filtrate));
      52           0 :     Int k = 0;
      53           0 :     for (size_t i = 0; i < feed.nelements(); ++i) {
      54           0 :       if (is_filtrate[i]) {
      55           0 :         filtrate[k] = feed[i];
      56           0 :         ++k;
      57             :       }
      58             :     }
      59             :   }
      60       27036 : }
      61             : 
      62             : //template<class T>
      63             : //inline void FiltrateMatrix(Matrix<T> const &feed,
      64             : //    Vector<bool> const &is_filtrate, Matrix<T> &filtrate) {
      65             : //  AlwaysAssert(feed.conform(is_filtrate), AipsError);
      66             : //  // filter_flag: true --> filtrate, false --> residue
      67             : //  filtrate.resize(ntrue(is_filtrate));
      68             : //  Int k = 0;
      69             : //  for (size_t i = 0; i < feed.nelements(); ++i) {
      70             : //    if (is_filtrate[i]) {
      71             : //      filtrate.column(k) = feed.column(i);
      72             : //      ++k;
      73             : //    }
      74             : //  }
      75             : //}
      76             : //
      77             : //template<class T>
      78             : //inline void FiltrateCube(Cube<T> const &feed, Vector<bool> const &is_filtrate,
      79             : //    Cube<T> &filtrate) {
      80             : //  AlwaysAssert(feed.conform(is_filtrate), AipsError);
      81             : //  // filter_flag: true --> filtrate, false --> residue
      82             : //  filtrate.resize(ntrue(is_filtrate));
      83             : //  Int k = 0;
      84             : //  for (size_t i = 0; i < feed.nelements(); ++i) {
      85             : //    if (is_filtrate[i]) {
      86             : //      filtrate.xyPlane(k) = feed.xyPlane(i);
      87             : //      ++k;
      88             : //    }
      89             : //  }
      90             : //}
      91             : 
      92             : template<class T>
      93        6632 : inline void FiltrateArray(Array<T> const &feed, Vector<bool> const &is_filtrate,
      94             :     Array<T> &filtrate) {
      95             :   // filter_flag: true --> filtrate, false --> residue
      96        6632 :   ssize_t const num_filtrates = ntrue(is_filtrate);
      97       13264 :   IPosition shape = feed.shape();
      98        6632 :   uInt const ndim = feed.ndim();
      99        6632 :   if (num_filtrates == shape[ndim - 1]) {
     100        6632 :     filtrate.resize();
     101        6632 :     filtrate.reference(feed);
     102             :   } else {
     103           0 :     shape[ndim - 1] = num_filtrates;
     104           0 :     filtrate.resize(shape);
     105           0 :     AlwaysAssert((size_t )feed.shape()[ndim - 1] == is_filtrate.nelements(),
     106             :         AipsError);
     107           0 :     IPosition iter_axis(1, ndim - 1);
     108           0 :     ArrayIterator<T> from_iter(feed, iter_axis, False);
     109           0 :     ArrayIterator<T> to_iter(filtrate, iter_axis, False);
     110           0 :     for (size_t i = 0; i < is_filtrate.nelements(); ++i) {
     111           0 :       if (is_filtrate[i]) {
     112           0 :         to_iter.array() = from_iter.array();
     113           0 :         to_iter.next();
     114             :       }
     115           0 :       from_iter.next();
     116             :     }
     117             :   }
     118        6632 : }
     119             : 
     120             : }
     121             : 
     122             : namespace casa { //# NAMESPACE CASA - BEGIN
     123             : 
     124             : namespace vi { //# NAMESPACE vi - BEGIN
     125             : 
     126             : // forward declaration
     127             : template<class Filter>
     128             : class FiltrationTVI;
     129             : 
     130             : // FiltrationTVI implementation
     131             : template<class Filter>
     132           9 : FiltrationTVI<Filter>::FiltrationTVI(ViImplementation2 * inputVi,
     133             :     Record const &configuration) :
     134             :     TransformingVi2(inputVi), configuration_p(configuration), filter_p(0), num_filtrates_p(
     135           9 :         0), is_filtrate_p(), is_valid_subchunk_p() {
     136             :   // new filter
     137             : //  MeasurementSet const &ms = getVii()->ms();
     138           9 :   filter_p = new Filter(configuration_p);
     139             : 
     140             :   // Initialize attached VisBuffer
     141           9 :   setVisBuffer(createAttachedVisBuffer(VbRekeyable));
     142           9 : }
     143             : 
     144             : template<class Filter>
     145          18 : FiltrationTVI<Filter>::~FiltrationTVI() {
     146           9 :   if (filter_p) {
     147           9 :     delete filter_p;
     148             :   }
     149          27 : }
     150             : 
     151             : template<class Filter>
     152        1667 : void FiltrationTVI<Filter>::origin() {
     153        1667 :   getVii()->origin();
     154             : //  cout << __func__ << ": subchunkId = " << getSubchunkId().subchunk() << endl;
     155             : 
     156             :   // Synchronize own VisBuffer -- is it required?
     157             :   //configureNewSubchunk();
     158             : 
     159             :   // filtration
     160        1667 :   filter();
     161             : 
     162             : //  if (more()) {
     163             : //    cout << __func__ << ": there is a subchunk passed through the filter ";
     164             : //    Vector<uInt> rowIds;
     165             : //    getVii()->getRowIds(rowIds);
     166             : //    cout << " is_filtrate_p = " << is_filtrate_p << ", rowIds = " << rowIds << endl;
     167             : //  } else {
     168             : //    cout << __func__ << ": no subchunk remaining after filtration" << endl;
     169             : //    cout << __func__ << ": is_filtrate_p = " << is_filtrate_p << endl;
     170             : //  }
     171        1667 : }
     172             : 
     173             : template<class Filter>
     174        1658 : void FiltrationTVI<Filter>::next() {
     175             :   // next subchunk
     176             :   // must call next at least one time
     177        1658 :   getVii()->next();
     178             : 
     179             :   // filtration
     180        1658 :   filter();
     181             : 
     182             : //  cout << __func__ << ": subchunkId = " << getSubchunkId().subchunk() << endl;
     183             : 
     184             : //  if (more()) {
     185             : //    cout << __func__ << ": there is a subchunk passed through the filter ";
     186             : //    Vector<uInt> rowIds;
     187             : //    getVii()->getRowIds(rowIds);
     188             : //    cout << " is_filtrate_p = " << is_filtrate_p << ", rowIds = " << rowIds << endl;
     189             : //  } else {
     190             : //    cout << __func__ << ": no subchunk remaining after filtration" << endl;
     191             : //    cout << __func__ << ": is_filtrate_p = " << is_filtrate_p << endl;
     192             : //  }
     193        1658 : }
     194             : 
     195             : template<class Filter>
     196          27 : void FiltrationTVI<Filter>::originChunks(Bool forceRewind) {
     197          27 :   TransformingVi2::originChunks(forceRewind);
     198             : 
     199             :   // sync
     200          27 :   filter_p->syncWith(this);
     201             : 
     202          27 :   filterChunk();
     203             : 
     204          27 :   getVii()->origin();
     205             : 
     206             : //  cout << __func__ << ": chunkId = " << getSubchunkId().chunk() << endl;
     207          27 : }
     208             : 
     209             : template<class Filter>
     210        3316 : void FiltrationTVI<Filter>::nextChunk() {
     211        3316 :   TransformingVi2::nextChunk();
     212             : 
     213             :   // sync
     214        3316 :   filter_p->syncWith(this);
     215             : 
     216        3316 :   filterChunk();
     217             : 
     218        3316 :   getVii()->origin();
     219             : 
     220             : //  cout << __func__ << ": chunkId = " << getSubchunkId().chunk() << endl;
     221        3316 : }
     222             : 
     223             : template<class Filter>
     224           0 : rownr_t FiltrationTVI<Filter>::nRows() const {
     225           0 :   return num_filtrates_p;
     226             : }
     227             : 
     228             : template<class Filter>
     229           0 : void FiltrationTVI<Filter>::getRowIds(Vector<rownr_t> &rowids) const {
     230           0 :   Vector<rownr_t> org;
     231           0 :   getVii()->getRowIds(org);
     232           0 :   ::FiltrateVector(org, is_filtrate_p, rowids);
     233           0 : }
     234             : 
     235             : template<class Filter>
     236        2474 : void FiltrationTVI<Filter>::antenna1(Vector<Int> &ant1) const {
     237        4948 :   Vector<Int> org;
     238        2474 :   getVii()->antenna1(org);
     239        2474 :   ::FiltrateVector(org, is_filtrate_p, ant1);
     240        2474 : }
     241             : 
     242             : template<class Filter>
     243        2474 : void FiltrationTVI<Filter>::antenna2(Vector<Int> &ant2) const {
     244        4948 :   Vector<Int> org;
     245        2474 :   getVii()->antenna2(org);
     246        2474 :   ::FiltrateVector(org, is_filtrate_p, ant2);
     247        2474 : }
     248             : 
     249             : template<class Filter>
     250         612 : void FiltrationTVI<Filter>::exposure(Vector<double> &expo) const {
     251        1224 :   Vector<double> org;
     252         612 :   getVii()->exposure(org);
     253         612 :   ::FiltrateVector(org, is_filtrate_p, expo);
     254         612 : }
     255             : 
     256             : template<class Filter>
     257           0 : void FiltrationTVI<Filter>::feed1(Vector<Int> &fd1) const {
     258           0 :   Vector<Int> org;
     259           0 :   getVii()->feed1(org);
     260           0 :   ::FiltrateVector(org, is_filtrate_p, fd1);
     261           0 : }
     262             : 
     263             : template<class Filter>
     264           0 : void FiltrationTVI<Filter>::feed2(Vector<Int> &fd2) const {
     265           0 :   Vector<Int> org;
     266           0 :   getVii()->feed2(org);
     267           0 :   ::FiltrateVector(org, is_filtrate_p, fd2);
     268           0 : }
     269             : 
     270             : template<class Filter>
     271        2474 : void FiltrationTVI<Filter>::fieldIds(Vector<Int> &fld) const {
     272        4948 :   Vector<Int> org;
     273        2474 :   getVii()->fieldIds(org);
     274        2474 :   ::FiltrateVector(org, is_filtrate_p, fld);
     275        2474 : }
     276             : 
     277             : template<class Filter>
     278        1658 : void FiltrationTVI<Filter>::arrayIds(Vector<Int> &arr) const {
     279        3316 :   Vector<Int> org;
     280        1658 :   getVii()->arrayIds(org);
     281        1658 :   ::FiltrateVector(org, is_filtrate_p, arr);
     282        1658 : }
     283             : 
     284             : template<class Filter>
     285        1658 : void FiltrationTVI<Filter>::flag(Cube<Bool> &flags) const {
     286        3316 :   Cube<Bool> org;
     287        1658 :   getVii()->flag(org);
     288        1658 :   ::FiltrateArray(org, is_filtrate_p, flags);
     289        1658 : }
     290             : 
     291             : template<class Filter>
     292           0 : void FiltrationTVI<Filter>::flag(Matrix<Bool> &flags) const {
     293           0 :   Matrix<Bool> org;
     294           0 :   getVii()->flag(org);
     295           0 :   ::FiltrateArray(org, is_filtrate_p, flags);
     296           0 : }
     297             : 
     298             : template<class Filter>
     299           0 : void FiltrationTVI<Filter>::flagCategory(Array<Bool> &flagCategories) const {
     300           0 :   Array<Bool> org;
     301           0 :   getVii()->flagCategory(org);
     302           0 :   ::FiltrateArray(org, is_filtrate_p, flagCategories);
     303           0 : }
     304             : 
     305             : template<class Filter>
     306        3316 : void FiltrationTVI<Filter>::flagRow(Vector<Bool> &rowflags) const {
     307        6632 :   Vector<Bool> org;
     308        3316 :   getVii()->flagRow(org);
     309        3316 :   ::FiltrateVector(org, is_filtrate_p, rowflags);
     310        3316 : }
     311             : 
     312             : template<class Filter>
     313        2474 : void FiltrationTVI<Filter>::observationId(Vector<Int> &obsids) const {
     314        4948 :   Vector<Int> org;
     315        2474 :   getVii()->observationId(org);
     316        2474 :   ::FiltrateVector(org, is_filtrate_p, obsids);
     317        2474 : }
     318             : 
     319             : template<class Filter>
     320           0 : void FiltrationTVI<Filter>::processorId(Vector<Int> &procids) const {
     321           0 :   Vector<Int> org;
     322           0 :   getVii()->processorId(org);
     323           0 :   ::FiltrateVector(org, is_filtrate_p, procids);
     324           0 : }
     325             : 
     326             : template<class Filter>
     327        2474 : void FiltrationTVI<Filter>::scan(Vector<Int> &scans) const {
     328        4948 :   Vector<Int> org;
     329        2474 :   getVii()->scan(org);
     330        2474 :   ::FiltrateVector(org, is_filtrate_p, scans);
     331        2474 : }
     332             : 
     333             : template<class Filter>
     334         816 : void FiltrationTVI<Filter>::stateId(Vector<Int> &stateids) const {
     335        1632 :   Vector<Int> org;
     336         816 :   getVii()->stateId(org);
     337         816 :   ::FiltrateVector(org, is_filtrate_p, stateids);
     338         816 : }
     339             : 
     340             : template<class Filter>
     341           0 : void FiltrationTVI<Filter>::jonesC(
     342             :     Vector<SquareMatrix<Complex, 2> > &cjones) const {
     343           0 :   Vector<SquareMatrix<Complex, 2>> org;
     344           0 :   getVii()->jonesC(org);
     345           0 :   ::FiltrateVector(org, is_filtrate_p, cjones);
     346           0 : }
     347             : 
     348             : template<class Filter>
     349        1658 : void FiltrationTVI<Filter>::sigma(Matrix<Float> &sigmat) const {
     350        3316 :   Matrix<Float> org;
     351        1658 :   getVii()->sigma(org);
     352        1658 :   ::FiltrateArray(org, is_filtrate_p, sigmat);
     353        1658 : }
     354             : 
     355             : template<class Filter>
     356        2474 : void FiltrationTVI<Filter>::spectralWindows(Vector<Int> &spws) const {
     357        4948 :   Vector<Int> org;
     358        2474 :   getVii()->spectralWindows(org);
     359        2474 :   ::FiltrateVector(org, is_filtrate_p, spws);
     360        2474 : }
     361             : 
     362             : template<class Filter>
     363        2474 : void FiltrationTVI<Filter>::time(Vector<double> &t) const {
     364        4948 :   Vector<double> org;
     365        2474 :   getVii()->time(org);
     366        2474 :   ::FiltrateVector(org, is_filtrate_p, t);
     367        2474 : }
     368             : 
     369             : template<class Filter>
     370        1658 : void FiltrationTVI<Filter>::timeCentroid(Vector<double> &t) const {
     371        3316 :   Vector<double> org;
     372        1658 :   getVii()->timeCentroid(org);
     373        1658 :   ::FiltrateVector(org, is_filtrate_p, t);
     374        1658 : }
     375             : 
     376             : template<class Filter>
     377           0 : void FiltrationTVI<Filter>::timeInterval(Vector<double> &ti) const {
     378           0 :   Vector<double> org;
     379           0 :   getVii()->timeInterval(org);
     380           0 :   ::FiltrateVector(org, is_filtrate_p, ti);
     381           0 : }
     382             : 
     383             : template<class Filter>
     384           0 : void FiltrationTVI<Filter>::uvw(Matrix<double> &uvwmat) const {
     385           0 :   Matrix<double> org;
     386           0 :   getVii()->uvw(org);
     387           0 :   ::FiltrateArray(org, is_filtrate_p, uvwmat);
     388           0 : }
     389             : 
     390             : template<class Filter>
     391           0 : void FiltrationTVI<Filter>::visibilityCorrected(Cube<Complex> &vis) const {
     392           0 :   Cube<Complex> org;
     393           0 :   getVii()->visibilityCorrected(org);
     394           0 :   ::FiltrateArray(org, is_filtrate_p, vis);
     395           0 : }
     396             : 
     397             : template<class Filter>
     398        1658 : void FiltrationTVI<Filter>::visibilityModel(Cube<Complex> &vis) const {
     399        3316 :   Cube<Complex> org;
     400        1658 :   getVii()->visibilityModel(org);
     401        1658 :   ::FiltrateArray(org, is_filtrate_p, vis);
     402        1658 : }
     403             : 
     404             : template<class Filter>
     405           0 : void FiltrationTVI<Filter>::visibilityObserved(Cube<Complex> &vis) const {
     406           0 :   Cube<Complex> org;
     407           0 :   getVii()->visibilityObserved(org);
     408           0 :   ::FiltrateArray(org, is_filtrate_p, vis);
     409           0 : }
     410             : 
     411             : template<class Filter>
     412        1658 : void FiltrationTVI<Filter>::floatData(Cube<Float> &fcube) const {
     413        3316 :   Cube<Float> org;
     414        1658 :   getVii()->floatData(org);
     415        1658 :   ::FiltrateArray(org, is_filtrate_p, fcube);
     416        1658 : }
     417             : 
     418             : template<class Filter>
     419           0 : IPosition FiltrationTVI<Filter>::visibilityShape() const {
     420           0 :   IPosition shape = getVii()->visibilityShape();
     421           0 :   shape[shape.nelements() - 1] = num_filtrates_p;
     422           0 :   return shape;
     423             : }
     424             : 
     425             : template<class Filter>
     426           0 : void FiltrationTVI<Filter>::weight(Matrix<Float> &wtmat) const {
     427           0 :   Matrix<Float> org;
     428           0 :   getVii()->weight(org);
     429           0 :   ::FiltrateArray(org, is_filtrate_p, wtmat);
     430           0 : }
     431             : 
     432             : template<class Filter>
     433           0 : void FiltrationTVI<Filter>::weightSpectrum(Cube<Float> &wtsp) const {
     434           0 :   Cube<Float> org;
     435           0 :   getVii()->weightSpectrum(org);
     436           0 :   ::FiltrateArray(org, is_filtrate_p, wtsp);
     437           0 : }
     438             : 
     439             : template<class Filter>
     440           0 : void FiltrationTVI<Filter>::sigmaSpectrum(Cube<Float> &sigmasp) const {
     441           0 :   Cube<Float> org;
     442           0 :   getVii()->sigmaSpectrum(org);
     443           0 :   ::FiltrateArray(org, is_filtrate_p, sigmasp);
     444           0 : }
     445             : 
     446             : template<class Filter>
     447        1658 : void FiltrationTVI<Filter>::dataDescriptionIds(Vector<Int> &ddids) const {
     448        3316 :   Vector<Int> org;
     449        1658 :   getVii()->dataDescriptionIds(org);
     450        1658 :   ::FiltrateVector(org, is_filtrate_p, ddids);
     451        1658 : }
     452             : 
     453             : template<class Filter>
     454        3325 : void FiltrationTVI<Filter>::filter() {
     455        3325 :   auto const vii = getVii();
     456        3325 :   auto const vb = vii->getVisBuffer();
     457        3325 :   for (; vii->more() && !is_valid_subchunk_p(vii->getSubchunkId().subchunk()); vii->next()) {
     458             :     // Synchronize own VisBuffer -- is it required to do inside the loop?
     459             :     //configureNewSubchunk();
     460             :   }
     461             : 
     462             :   // update filter information
     463        3325 :   num_filtrates_p = filter_p->isFiltratePerRow(vb, is_filtrate_p);
     464             : 
     465             :   // Synchronize own VisBuffer
     466        3325 :   if(vii->more())
     467        1667 :     configureNewSubchunk();
     468        3325 : }
     469             : 
     470             : template<class Filter>
     471        3343 : void FiltrationTVI<Filter>::filterChunk() {
     472        3343 :   size_t const block_size = max(getVii()->nRowsInChunk(), (size_t)100);
     473        3343 :   if (is_valid_subchunk_p.nelements() < block_size) {
     474           9 :     is_valid_subchunk_p.resize(block_size);
     475             :   }
     476        3343 :   is_valid_subchunk_p = false;
     477             : 
     478             :   // increment iterator while the chunk doesn't contain valid subchunk
     479      162271 :   for (auto vii = getVii(); vii->moreChunks(); vii->nextChunk()) {
     480      162253 :     is_valid_subchunk_p = false;
     481      324506 :     for (vii->origin(); vii->more(); vii->next()) {
     482      162253 :       size_t const subchunk_id = vii->getSubchunkId().subchunk();
     483      162253 :       size_t const n = is_valid_subchunk_p.nelements();
     484      162253 :       size_t m = n;
     485      162253 :       while (m < subchunk_id) {
     486           0 :         is_valid_subchunk_p.resize(m + block_size);
     487           0 :         m = is_valid_subchunk_p.nelements();
     488             :       }
     489      162253 :       is_valid_subchunk_p(Slice(n, m - n)) = false;
     490      162253 :       is_valid_subchunk_p[subchunk_id] = filter_p->isFiltrate(vii->getVisBuffer());
     491             : //      cout << __func__ << ": chunk " << vii->getSubchunkId().chunk() << " subchunk " << subchunk_id
     492             : //          << " is_valid " << is_valid_subchunk_p[subchunk_id] << endl;
     493             :     }
     494             : 
     495      162253 :     if (anyTrue(is_valid_subchunk_p)) {
     496        3325 :       break;
     497             :     }
     498             :   }
     499             : 
     500        3343 : }
     501             : 
     502             : } //# NAMESPACE vi - END
     503             : 
     504             : } //# NAMESPACE CASA - END
     505             : 
     506             : #endif // _MSVIS_FILTRATIONTVI_TCC_