LCOV - code coverage report
Current view: top level - flagging/Flagging - RFASpectralRej.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 0 186 0.0 %
Date: 2024-11-06 17:42:47 Functions: 0 11 0.0 %

          Line data    Source code
       1             : //# RFASpectralRej.cc: this defines RFASpectralRej
       2             : //# Copyright (C) 2000,2001,2002
       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 <flagging/Flagging/RFASpectralRej.h> 
      29             : #include <casacore/scimath/Functionals/Polynomial.h>
      30             : #include <msvis/MSVis/VisibilityIterator.h>
      31             : #include <msvis/MSVis/VisBuffer.h>
      32             : #include <casacore/casa/Arrays/ArrayMath.h>
      33             : #include <casacore/casa/Arrays/ArrayLogical.h>
      34             : #include <casacore/casa/Arrays/Slice.h>
      35             : #include <casacore/casa/System/PGPlotterInterface.h>
      36             :     
      37             : using namespace casacore;
      38             : namespace casa { //# NAMESPACE CASA - BEGIN
      39             : 
      40           0 : void RFASpectralRej::addSegment ( Int spwid,Double fq0,Double fq1,Int ch0,Int ch1 )
      41             : {
      42           0 :   Segment seg = { spwid,fq0,fq1,ch0,ch1 };
      43           0 :   uInt n = segments.nelements();
      44           0 :   segments.resize(n+1);
      45           0 :   segments[n] = seg;
      46           0 : }
      47             : 
      48             : // -----------------------------------------------------------------------
      49             : // parseRegion
      50             : // Helper function to parse one region specification. Returns
      51             : // the number of ranges parsed.
      52             : // -----------------------------------------------------------------------
      53           0 : void RFASpectralRej::parseRegion ( const RecordInterface &parm )
      54             : {
      55           0 :   Bool parsed = false;
      56           0 :   Int spwid = -1;
      57           0 :   if( isFieldSet(parm,RF_SPWID) && fieldType(parm,RF_SPWID,TpInt) )
      58           0 :     spwid = parm.asInt(RF_SPWID) - indexingBase();
      59             : // figure out how channel ranges are specified - frequencies or channel #'s
      60             : // First version is frequencies
      61           0 :   if( fieldSize(parm,RF_FREQS)>1 )
      62             :   {
      63           0 :     Array<Double> freqarr;
      64             :     try {
      65           0 :       freqarr = parm.toArrayDouble(RF_FREQS);
      66             :       // make sure array is of the right shape (can be reformed into 2xN)
      67           0 :       if( freqarr.ndim()>2 || (freqarr.nelements()%2) !=0 )
      68           0 :         throw( AipsError("") );
      69           0 :     } catch( AipsError x ) {
      70           0 :       os<<"Illegal \""<<RF_FREQS<<"\" array\n"<<LogIO::EXCEPTION;
      71           0 :     }
      72           0 :     Matrix<Double> fq( freqarr.reform(IPosition(2,2,freqarr.nelements()/2)) );
      73             :     // enqueue the region specs
      74           0 :     for( uInt i=0; i<fq.ncolumn(); i++ )
      75             :     {
      76           0 :       if( fq(0,i) >= fq(1,i) )
      77             :       {
      78             :         char s[128];
      79           0 :         sprintf(s,"Illegal spectral region %0.2f-%0.2f\n",fq(0,i),fq(1,i));
      80           0 :         os<<s<<LogIO::EXCEPTION;
      81             :       }
      82           0 :       addSegment(spwid,fq(0,i),fq(1,i),-1,-1);
      83             :     }
      84           0 :     parsed=true;
      85           0 :   }
      86             : // second option is channel numbers
      87           0 :   if( fieldSize(parm,RF_CHANS)>1 )
      88             :   {
      89           0 :     Array<Int> arr;
      90             :     try {
      91           0 :       arr = parm.toArrayInt(RF_CHANS);
      92             :       // make sure array is of the right shape (can be reformed into 2xN)
      93           0 :       if( arr.ndim()>2 || (arr.nelements()%2) !=0 )
      94           0 :         throw( AipsError("") );
      95           0 :     } catch( AipsError x ) {
      96           0 :       os<<"Illegal \""<<RF_CHANS<<"\" array\n"<<LogIO::EXCEPTION;
      97           0 :     }
      98           0 :     Matrix<Int> ch( arr.reform(IPosition(2,2,arr.nelements()/2)) );
      99             :     // enqueue the region specs
     100           0 :     for( uInt i=0; i<ch.ncolumn(); i++ )
     101             :     {
     102           0 :       if( ch(0,i) >= ch(1,i) )
     103             :       {
     104             :         char s[128];
     105           0 :         sprintf(s,"Illegal spectral region #%d-%d\n",ch(0,i),ch(1,i));
     106           0 :         os<<s<<LogIO::EXCEPTION;
     107             :       }
     108           0 :       ch -= (Int)indexingBase();
     109           0 :       addSegment(spwid,0,0,ch(0,i),ch(1,i));
     110             :     }
     111           0 :     parsed=true;
     112           0 :   }
     113           0 :   if( !parsed )
     114           0 :     os<<"\""<<RF_FREQS<<"\" or \""<<RF_CHANS<<"\" must be specified\n"<<LogIO::EXCEPTION;
     115           0 : }
     116             :     
     117           0 : RFASpectralRej::RFASpectralRej  ( RFChunkStats &ch,const RecordInterface &parm ) : 
     118             :     RFAFlagCubeBase(ch,parm),
     119             :     RFDataMapper(parm.asArrayString(RF_EXPR),parm.asString(RF_COLUMN)),
     120           0 :     ndeg( parm.asInt(RF_NDEG) ),
     121           0 :     halfwin( parm.asInt(RF_ROW_HW) ),
     122           0 :     threshold( parm.asDouble(RF_ROW_THR) ),
     123           0 :     rowclipper(chunk,flag,threshold,halfwin)
     124             : {
     125             : // figure out how channel ranges are specified
     126             : // if a full region record is specified, parse each element
     127             : // otherwise just parse the parameter record itself
     128           0 :   if( isValidRecord(parm,RF_REGION) ) // full region record
     129             :   {
     130           0 :     const RecordInterface &reg( parm.asRecord(RF_REGION) );
     131           0 :     for( uInt i=0; i<reg.nfields(); i++ )
     132             :     {
     133           0 :       if( reg.type(i) != TpRecord )
     134           0 :         os<<"\""<<RF_REGION<<"\" must be a record of records\n"<<LogIO::EXCEPTION;
     135           0 :       parseRegion(reg.asRecord(i));
     136             :     }
     137             :   }
     138             :   else // else assume only one region specified
     139           0 :     parseRegion(parm);
     140             :   
     141           0 :   if( !segments.nelements() )
     142           0 :     os<<"No spectral region has been specified\n"<<LogIO::EXCEPTION;
     143             :     
     144             : // set up fitter
     145           0 :   Polynomial<AutoDiff<Float> > poly(ndeg);
     146           0 :   fitter.setFunction(poly);
     147             :   
     148             : // set up debugging info
     149           0 :   if( fieldType(parm,RF_DEBUG,TpArrayInt) )
     150             :   {
     151           0 :     Vector<Int> dbg;
     152           0 :     parm.get(RF_DEBUG,dbg);
     153           0 :     if(dbg.nelements() != 2 && dbg.nelements() != 3)
     154             :     {
     155           0 :       os<<"\""<<RF_DEBUG<<"\" parameter must be [NIFR,NTIME] or [ANT1,ANT2,NTIME]"<<LogIO::EXCEPTION;
     156             :     }
     157           0 :   }
     158           0 : }
     159             : 
     160             : // -----------------------------------------------------------------------
     161             : // newChunk
     162             : // At each new chunk, figure out which channels fit into the
     163             : // specified fitting regions.
     164             : // -----------------------------------------------------------------------
     165           0 :   Bool RFASpectralRej::newChunk (Int &maxmem)
     166             : {
     167             : // compute correlations mask, return false if fails
     168           0 :   corrmask = RFDataMapper::corrMask(chunk.visIter());
     169           0 :   if( !corrmask )
     170             :   {
     171           0 :     os<<LogIO::WARN<<"missing selected correlations, ignoring this chunk\n"<<LogIO::POST;
     172           0 :     return active=false;
     173             :   }
     174             : // figure out active channels (i.e. within specified segments)
     175           0 :   Int spwid = chunk.visBuf().spectralWindow();
     176           0 :   fitchan.resize(num(CHAN));
     177           0 :   fitchan.set(false);
     178           0 :   const Vector<Double> & fq( chunk.frequency()*1e-6 );
     179           0 :   for( uInt i=0; i<segments.nelements(); i++)
     180             :   {
     181           0 :     const Segment &seg ( segments[i] );
     182             :     // compare spectral windows, if specified
     183           0 :     if( seg.spwid >= 0 && seg.spwid != spwid )
     184           0 :       continue;
     185           0 :     if( seg.ch0 >= 0 )  // use channel numbers
     186             :     {
     187           0 :       if( (uInt)seg.ch0 < num(CHAN) )
     188             :       {
     189           0 :         Int ch1 = num(CHAN)-1;
     190           0 :         if( seg.ch1 < ch1 )
     191           0 :           ch1 = seg.ch1;
     192           0 :         fitchan(Slice(seg.ch0,ch1-seg.ch0+1)) = true;
     193             :       }
     194             :     }
     195             :     else // use frequencies
     196             :     {
     197           0 :       fitchan = fitchan || ( fq >= seg.fq0 && fq <= seg.fq1 );
     198             :     }
     199             :   }
     200             : // count number of fitted channels  
     201           0 :   num_fitchan = 0;
     202           0 :   for( uInt i=0; i<num(CHAN); i++ )
     203             :   {
     204           0 :     if( fitchan(i) )
     205             :     {
     206           0 :       xnorm = i;
     207           0 :       num_fitchan++;
     208             :     }
     209             :   }
     210             : // return if none 
     211           0 :   os<<num_fitchan<<" channels will be fitted in this chunk\n"<<LogIO::POST;
     212           0 :   if( num_fitchan<ndeg+2 )
     213             :   {
     214           0 :     os<<LogIO::WARN<<"not enough channels, ignoring chunk\n"<<LogIO::POST;
     215           0 :     return active=false;
     216             :   }
     217             : // finish with init  
     218           0 :   RFAFlagCubeBase::newChunk(maxmem-=1);
     219           0 :   rowclipper.init(num(IFR),num(TIME));
     220           0 :   return active=true;
     221           0 : }
     222             : 
     223           0 : void RFASpectralRej::endChunk ()
     224             : {
     225           0 :   RFAFlagCubeBase::endChunk();
     226           0 :   flag.cleanup();
     227           0 :   rowclipper.cleanup();
     228           0 : }
     229             : 
     230           0 : void RFASpectralRej::startData (bool verbose)
     231             : {
     232           0 :   RFAFlagCubeBase::startData(verbose);
     233           0 :   rowclipper.reset();
     234           0 : }
     235             : 
     236           0 : RFA::IterMode RFASpectralRej::iterTime (uInt it)
     237             : {
     238           0 :   RFAFlagCubeBase::iterTime(it);
     239           0 :   RFDataMapper::setVisBuffer(chunk.visBuf());
     240           0 :   return RFA::CONT;
     241             : }
     242             : 
     243           0 : RFA::IterMode RFASpectralRej::iterRow ( uInt irow )
     244             : {
     245             : // during first pass, compute diff-median. Also keep track of the AAD.
     246           0 :   uInt iifr = chunk.ifrNum(irow);
     247           0 :   uInt it = chunk.iTime();
     248           0 :   Bool rowfl = chunk.npass() ? flag.rowFlagged(iifr,it) 
     249           0 :                             : flag.rowPreFlagged(iifr,it);
     250           0 :   if( !rowfl )   
     251             :   {
     252           0 :     Vector<Float> x(num_fitchan),y(num_fitchan);
     253           0 :     Vector<uInt>  chan(num_fitchan);
     254           0 :     uInt np=0;
     255             : // loop over channels, collect valid (non-flagged) pixels into the x and y
     256             : // vectors
     257           0 :     for( uInt ich=0; ich<num(CHAN); ich++ )
     258             :     {
     259           0 :       if( fitchan(ich) && !(chunk.npass() ? flag.anyFlagged(ich,iifr) : flag.preFlagged(ich,iifr)) )
     260             :       {
     261           0 :         x(np) = ich/xnorm;
     262           0 :         y(np) = mapValue(ich,irow);
     263           0 :         chan(np++) = ich;
     264             :       }
     265             :     }
     266             : // check that we have enough points to constrain the fit
     267           0 :     if( np > ndeg+2 )
     268             :     {
     269             :   // resize x/y vectors, and do the fit
     270           0 :       Slice S(0,np);
     271           0 :       Vector<Float> sigma(np,1.0f),x1(x(S)),y1(y(S));
     272           0 :       Vector<Float> c = fitter.fit(x1,y1,sigma);
     273           0 :       Float chisq = fitter.chiSquare();
     274           0 :       rowclipper.setSigma(iifr,it,chisq);
     275           0 :     } // endif( np>ndeg+2)
     276           0 :   } // endif( !rowfl )
     277             :   
     278           0 :   return RFA::CONT;
     279             : }
     280             : 
     281             : // -----------------------------------------------------------------------
     282             : // endData
     283             : // Ends data pass
     284             : // -----------------------------------------------------------------------
     285           0 : RFA::IterMode RFASpectralRej::endData ()
     286             : {
     287           0 :   RFAFlagCubeBase::endData();
     288             :   uInt dum;
     289           0 :   rowclipper.updateSigma(dum,dum);
     290           0 :   return RFA::STOP;
     291             : }
     292             : 
     293             : // -----------------------------------------------------------------------
     294             : // RFASpectralRej::getDesc
     295             : // Returns description of parameters
     296             : // -----------------------------------------------------------------------
     297           0 : String RFASpectralRej::getDesc ()
     298             : {
     299           0 :   String desc( RFDataMapper::description()+";" );
     300             :   char s[256];
     301             : // build up description of spectral segments
     302           0 :   for( uInt i=0; i<segments.nelements(); i++)
     303             :   {
     304           0 :     const Segment &seg ( segments[i] );
     305             :     // is spwid specified?
     306             :     char s1[32];
     307           0 :     if( seg.spwid >= 0 )
     308           0 :       sprintf(s1,"%d:",seg.spwid);
     309             :     else
     310           0 :       s1[0]=0;
     311           0 :     if( seg.ch0 >= 0 )  // use channel numbers
     312           0 :       sprintf(s, " %s#%d-%d",s1,seg.ch0,seg.ch1);
     313             :     else
     314           0 :       sprintf(s, " %s%.2f-%.2fMHz",s1,seg.fq0,seg.fq1);
     315           0 :     desc+=s;
     316             :   }
     317           0 :   desc += "; ";
     318           0 :   sprintf(s,
     319             :           "%s=%d %s=%.1f %s=%d",
     320             :           RF_NDEG,
     321             :           ndeg,
     322             :           RF_ROW_THR,
     323             :           threshold,
     324             :           RF_ROW_HW,
     325             :           halfwin);
     326           0 :   desc += s;
     327           0 :   desc += RFAFlagCubeBase::getDesc();
     328           0 :   return desc;
     329           0 : }
     330             : 
     331             : // -----------------------------------------------------------------------
     332             : // RFASpectralRej::getDefaults
     333             : // Returns record of default parameters
     334             : // -----------------------------------------------------------------------
     335           0 : const RecordInterface & RFASpectralRej::getDefaults ()
     336             : {
     337           0 :   static Record rec;
     338             : // create record description on first entry
     339           0 :   if( !rec.nfields() )
     340             :   {
     341           0 :     rec = RFAFlagCubeBase::getDefaults();
     342           0 :     rec.define(RF_NAME,"SpectralRej");
     343           0 :     rec.define(RF_COLUMN,"DATA");
     344           0 :     rec.define(RF_EXPR,"+ ABS XX YY");
     345           0 :     rec.define(RF_NDEG,(Int)2);
     346           0 :     rec.define(RF_ROW_THR,(Double)5);
     347           0 :     rec.define(RF_ROW_HW,(Int)6);
     348           0 :     rec.define(RF_REGION,false);
     349           0 :     rec.define(RF_SPWID,false);
     350           0 :     rec.define(RF_FREQS,false);
     351           0 :     rec.define(RF_CHANS,false);
     352           0 :     rec.define(RF_DEBUG,false);
     353             :     
     354           0 :     rec.setComment(RF_COLUMN,"Use column: [DATA|MODEL|CORRected]");
     355           0 :     rec.setComment(RF_EXPR,"Expression for deriving value (e.g. \"ABS XX\", \"+ ABS XX YY\")");
     356           0 :     rec.setComment(RF_NDEG,"Number of degrees for polynomial fit");
     357           0 :     rec.setComment(RF_ROW_THR,"Row flagging threshold, in AADs");
     358           0 :     rec.setComment(RF_ROW_HW,"Row flagging, half-window of sliding median (0 to use overall median)");
     359           0 :     rec.setComment(RF_SPWID,"Spectral window ID (F or -1 for all)");
     360           0 :     rec.setComment(RF_FREQS,"Range(s) of frequencies (in MHz)");
     361           0 :     rec.setComment(RF_CHANS,"Range(s) of channel numbers");
     362           0 :     rec.setComment(RF_REGION,"For several spectral regions, record of records");
     363           0 :     rec.setComment(RF_DEBUG,"Set to [NIFR,NTIME] or [ANT1,ANT2,NTIME] to produce debugging plots");
     364             :   }
     365           0 :   return rec;
     366             : }
     367             : 
     368             :   
     369             : 
     370             : 
     371             : } //# NAMESPACE CASA - END
     372             : 

Generated by: LCOV version 1.16