Changeset 2852

Timestamp:

09/13/13 19:09:11 (11 years ago)

Author:

Kana Sugimoto

Message:

New Development: Yes

JIRA Issue: Yes (CAS-5139/TRAC-290)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: A couple of private functions are added to STSideBandSep class.

Test Programs:

Put in Release Notes: Yes

Module(s): asap.sbseparator

Description: Implemented a way to transfer channel based flag of imput tables to output tables. Channel flag is accumulated by logical SUM.

Location:

trunk/src

Files:

• STSideBandSep.cpp (modified) (10 diffs)
• STSideBandSep.h (modified) (1 diff)

trunk/src/STSideBandSep.cpp

@@ -261 +261 @@
   remRowIds.resize(0);
   Matrix<float> specMat(nchan_, nshift_);
+  Matrix<bool> flagMat(nchan_, nshift_);
   vector<float> sigSpec(nchan_), imgSpec(nchan_);
+  Vector<bool> flagVec(nchan_);
   vector<uInt> tabIdvec;
 
@@ -276 +278 @@
     // Get a set of spectra to solve
     if (!getSpectraToSolve(polId, beamId, dir[0], dir[1],
-                           specMat, tabIdvec)){
+                           specMat, flagMat, tabIdvec)){
       remRowIds.push_back(irow);
 #ifdef KS_DEBUG
@@ -289 +291 @@
       // unflag the spectrum since there should be some valid data
       sigTab_p->flagRow(vector<uInt>(irow), true);
+      // need to unflag whole channels anyway
       sigTab_p->flag(irow, vector<bool>(), true);
     }
+    // apply channel flag
+    flagVec = collapseFlag(flagMat, tabIdvec, true);
+    //boolVec = !boolVec; // flag
+    sigTab_p->flag(irow, flagVec.tovector(), false);
 
     // Solve image sideband
@@ -299 +306 @@
         // unflag the spectrum since there should be some valid data
         imgTab_p->flagRow(vector<uInt>(irow), true);
+        // need to unflag whole channels anyway
         imgTab_p->flag(irow, vector<bool>(), true);
       }
+      // apply channel flag
+      flagVec = collapseFlag(flagMat, tabIdvec, false);
+      //boolVec = !boolVec; // flag
+      imgTab_p->flag(irow, flagVec.tovector(), true);
     }
   } // end of row loop
@@ -594 +606 @@
 };
 
+// bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
+//                                    const double dirX, const double dirY,
+//                                    Matrix<float> &specMat, vector<uInt> &tabIdvec)
 bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
                                       const double dirX, const double dirY,
-                                      Matrix<float> &specMat, vector<uInt> &tabIdvec)
+                                      Matrix<float> &specMat,
+                                      Matrix<bool> &flagMat,
+                                      vector<uInt> &tabIdvec)
 {
   LogIO os(LogOrigin("STSideBandSep","getSpectraToSolve()", WHERE));
@@ -603 +620 @@
   specMat.resize(nchan_, nshift_);
   Vector<float> spec;
+  Vector<bool> boolVec;
   uInt nspec = 0;
   STMath stm(false); // insitu has no effect for average.
@@ -675 +693 @@
     // Interpolate flagged channels of the spectrum.
     Vector<Float> tmpSpec = avetab_p->getSpectrum(0);
+    // Mask is true if the data is valid (!flag)
     vector<bool> mask = avetab_p->getMask(0);
     mathutil::doZeroOrderInterpolation(tmpSpec, mask);
     spec.reference(specMat.column(nspec));
     spec = tmpSpec;
+    boolVec.reference(flagMat.column(nspec));
+    boolVec = mask; // cast std::vector to casa::Vector
+    boolVec = !boolVec;
     tabIdvec.push_back((uInt) itab);
     nspec++;
+    //Liberate from reference
+    spec.unique();
+    boolVec.unique();
   } // end of table loop
   // Check the number of selected spectra and resize matrix.
@@ -686 +711 @@
     //shiftSpecmat.resize(nchan_, nspec, true);
     specMat.resize(nchan_, nspec, true);
+    flagMat.resize(nchan_, nspec, true);
 #ifdef KS_DEBUG
       cout << "Could not find corresponding rows in some tables."
@@ -701 +727 @@
   return true;
 };
+
+
+Vector<bool> STSideBandSep::collapseFlag(const Matrix<bool> &flagMat,
+                                         const vector<uInt> &tabIdvec,
+                                         const bool signal)
+{
+  LogIO os(LogOrigin("STSideBandSep","collapseFlag()", WHERE));
+  if (tabIdvec.size() == 0)
+    throw(AipsError("Internal error. Table index is not defined."));
+  if (flagMat.ncolumn() != tabIdvec.size())
+    throw(AipsError("Internal error. The row number of input matrix is not conformant."));
+  if (flagMat.nrow() != nchan_)
+    throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
+  
+  const size_t nspec = tabIdvec.size();
+  vector<double> *thisShift;
+  if (signal == otherside_) {
+    // (solve signal && solveother = T) OR (solve image && solveother = F)
+    thisShift = &imgShift_;
+  } else {
+    // (solve signal && solveother = F) OR (solve image && solveother = T)
+    thisShift =  &sigShift_;
+ }
+
+  Vector<bool> outflag(nchan_, false);
+  double tempshift;
+  Vector<bool> shiftvec(nchan_, false);
+  Vector<bool> accflag(nchan_, false);
+  uInt shiftId;
+  for (uInt i = 0 ; i < nspec; ++i) {
+    shiftId = tabIdvec[i];
+    tempshift = - thisShift->at(shiftId);
+    shiftFlag(flagMat.column(i), tempshift, shiftvec);
+    // Now accumulate Flag
+    for (uInt j = 0 ; j < nchan_ ; ++j)
+      accflag[j] |= shiftvec[j];
+  }
+  // Shift back Flag
+  shiftFlag(accflag, thisShift->at(0), outflag);
+
+  return outflag;
+}
+
 
 vector<float> STSideBandSep::solve(const Matrix<float> &specmat,
@@ -815 +884 @@
   }
 };
+
+
+void STSideBandSep::shiftFlag(const Vector<bool> &invec,
+                                  double shift,
+                                  Vector<bool> &outvec)
+{
+  LogIO os(LogOrigin("STSideBandSep","shiftFlag()", WHERE));
+  if (invec.size() != nchan_)
+    throw(AipsError("Internal error. The length of input vector differs from nchan_"));
+  if (outvec.size() != nchan_)
+    throw(AipsError("Internal error. The length of output vector differs from nchan_"));
+
+#ifdef KS_DEBUG
+  cout << "Start shifting flag for " << shift << "channels" << endl;
+#endif
+
+  // shift is almost integer think it as int.
+  // tolerance should be in 0 - 1
+  double tolerance = 0.01;
+  // tweak shift to be in 0 ~ nchan_-1
+  if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
+  if (shift < 0.) shift += nchan_;
+  double rweight = fmod(shift, 1.);
+  bool ruse(true), luse(true);
+  if (rweight < 0.) rweight += 1.;
+  if (rweight < tolerance){
+    ruse = true;
+    luse = false;
+  }
+  if (rweight > 1-tolerance){
+    ruse = false;
+    luse = true;
+  }
+  uInt lchan, rchan;
+
+  outvec = false;
+  for (uInt i = 0 ; i < nchan_ ; i++) {
+    lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
+    if (lchan < 0.) lchan += nchan_;
+    rchan = ( (lchan + 1) % nchan_ );
+    outvec(lchan) |= (invec(i) && luse);
+    outvec(rchan) |= (invec(i) && ruse);
+  }
+};
+
 
 void STSideBandSep::deconvolve(Matrix<float> &specmat,

trunk/src/STSideBandSep.h

@@ -130 +130 @@
                          const double refval, const double refpix,
                          const double increment, const int nChan);
+  //  bool getSpectraToSolve(const int polId, const int beamId,
+  //                     const double dirX, const double dirY,
+  //                     Matrix<float> &specmat, vector<uInt> &tabIdvec);
   bool getSpectraToSolve(const int polId, const int beamId,
                          const double dirX, const double dirY,
-                         Matrix<float> &specmat, vector<uInt> &tabIdvec);
+                         Matrix<float> &specMat, Matrix<bool> &flagMat,
+                         vector<uInt> &tabIdvec);
 
-  vector<float> solve(const Matrix<float> &specmat,
+  vector<float> solve(const Matrix<float> &specMat,
                       const vector<uInt> &tabIdvec,
                       const bool signal = true);
 
+  Vector<bool> collapseFlag(const Matrix<bool> &flagMat,
+                            const vector<uInt> &tabIdvec,
+                            const bool signal = true);
+
   void shiftSpectrum(const Vector<float> &invec, double shift,
                      Vector<float> &outvec);
+
+  void shiftFlag(const Vector<bool> &invec, double shift,
+                     Vector<bool> &outvec);
 
   void deconvolve(Matrix<float> &specmat, const vector<double> shiftvec,