Changeset 2963 for trunk

Timestamp:

06/30/14 12:46:50 (10 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: Yes CAS-6585, CAS-6571

Ready for Test: Yes

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

STApplyCal is updated so that it takes into account channel flag information
when data is interpolated in time and in frequency.

Location:

trunk/src

Files:

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

trunk/src/STApplyCal.cpp

@@ -45 +45 @@
 using namespace std;
 
+namespace {
+// utility functions
+Vector<uInt> indexSort(Vector<Double> &t)
+{
+  Sort sort;
+  sort.sortKey(&t[0], TpDouble, 0, Sort::Ascending);
+  Vector<uInt> idx;
+  sort.sort(idx, t.nelements(), Sort::QuickSort|Sort::NoDuplicates);
+  return idx;
+}
+
+void getSortedData(const Vector<Double> key, const Vector<Float> data,
+                   const Vector<uChar> flag,
+                   Vector<Double> sortedKey, Vector<Float> sortedData,
+                   Vector<uChar> sortedFlag)
+{
+}
+}
+
 namespace asap {
-
 STApplyCal::STApplyCal()
 {
@@ -292 +310 @@
       uInt nrowThisTsys = tsystable_[i]->nrow();
       nrowTsys += nrowThisTsys;
-      if (nrowThisTsys > 0 and nchanTsys == 0) {
+      if (nrowThisTsys > 0 && nchanTsys == 0) {
         nchanTsys = tsystable_[i]->nchan(tsysifno);
         ftsys = tsystable_[i]->getBaseFrequency(0);
@@ -323 +341 @@
     }
     
-    Vector<uInt> skyIdx = timeSort(timeSky);
+    Vector<uInt> skyIdx = indexSort(timeSky);
     nrowSkySorted = skyIdx.nelements();
     
@@ -371 +389 @@
       }
     }
-    Vector<uInt> tsysIdx = timeSort(timeTsys);
+    Vector<uInt> tsysIdx = indexSort(timeTsys);
     nrowTsysSorted = tsysIdx.nelements();
 
@@ -407 +425 @@
   // Array for interpolated off spectrum
   Vector<Float> iOff(IPosition(1, nchanSp), new Float[nchanSp], TAKE_OVER);
+
+  // working array for interpolation with flags
+  uInt arraySize = max(max(nrowTsysSorted, nchanTsys), nrowSkySorted);
+  Vector<Double> xwork(IPosition(1, arraySize), new Double[arraySize], TAKE_OVER);
+  Vector<Float> ywork(IPosition(1, arraySize), new Float[arraySize], TAKE_OVER);
+  Vector<uChar> fwork(IPosition(1, nchanTsys), new uChar[nchanTsys], TAKE_OVER);
   
   for (uInt i = 0; i < rows.nelements(); i++) {
@@ -420 +444 @@
     // interpolation
     Double t0 = timeCol(irow);
+    Double *xwork_p = xwork.data();
+    Float *ywork_p = ywork.data();
     for (uInt ichan = 0; ichan < nchanSp; ichan++) {
       Float *tmpY = &(spoffSorted.data()[ichan * nrowSkySorted]);
-      if (allNE(flagoffSorted.column(ichan), (uChar)0)) {
+      uChar *tmpF = &(flagoffSorted.data()[ichan * nrowSkySorted]);
+      uInt wnrow = 0;
+      for (uInt ir = 0; ir < nrowSkySorted; ++ir) {
+        if (tmpF[ir] == 0) {
+          xwork_p[wnrow] = timeSkySorted.data()[ir];
+          ywork_p[wnrow] = tmpY[ir];
+          wnrow++;
+        }
+      }
+      //if (allNE(flagoffSorted.column(ichan), (uChar)0)) {
+      if (wnrow > 0) {
+        // any valid reference data
+        interpolatorS_->setData(xwork_p, ywork_p, wnrow);
+      }
+      else {
+        // no valid reference data
+        // interpolate data regardless of flag
+        interpolatorS_->setData(timeSkySorted.data(), tmpY, nrowSkySorted);
+        // flag this channel for calibrated data
         flag[ichan] = 1 << 7; // user flag
       }
-      interpolatorS_->setY(tmpY, nrowSkySorted);
+      //interpolatorS_->setY(tmpY, nrowSkySorted);
       iOff[ichan] = interpolatorS_->interpolate(t0);
     }
@@ -433 +477 @@
     if (doTsys) {
       // Tsys correction
+      // interpolation on time axis
       Float *yt = iTsysT.data();
+      uChar *fwork_p = fwork.data();
       for (uInt ichan = 0; ichan < nchanTsys; ichan++) {
         Float *tmpY = &(tsysSorted.data()[ichan * nrowTsysSorted]);
-        interpolatorT_->setY(tmpY, nrowTsysSorted);
-        iTsysT[ichan] = interpolatorT_->interpolate(t0);
+        uChar *tmpF = &(flagtsysSorted.data()[ichan * nrowTsysSorted]);
+        uInt wnrow = 0;
+        for (uInt ir = 0; ir < nrowTsysSorted; ++ir) {
+          if (tmpF[ir] == 0) {
+            xwork_p[wnrow] = timeTsysSorted.data()[ir];
+            ywork_p[wnrow] = tmpY[ir];
+            wnrow++;
+          }
+        }
+        if (wnrow > 0) {
+          // any valid value exists
+          //interpolatorT_->setY(tmpY, nrowTsysSorted);
+          interpolatorT_->setData(xwork_p, ywork_p, wnrow);
+          iTsysT[ichan] = interpolatorT_->interpolate(t0);
+          fwork_p[ichan] = 0;
+        }
+        else {
+          // no valid data
+          fwork_p[ichan] = 1 << 7; // user flag
+        }
       }
       if (nchanSp == 1) {
@@ -446 +510 @@
         // interpolation on frequency axis
         Vector<Double> fsp = getBaseFrequency(rows[i]);
-        interpolatorF_->setY(yt, nchanTsys);
+        uInt wnchan = 0;
+        for (uInt ichan = 0; ichan < nchanTsys; ++ichan) {
+          if (fwork_p[ichan] == 0) {
+            xwork_p[wnchan] = ftsys.data()[ichan];
+            ywork_p[wnchan] = yt[ichan];
+            ++wnchan;
+          }
+        }
+        //interpolatorF_->setY(yt, nchanTsys);
+        interpolatorF_->setData(xwork_p, ywork_p, wnchan);
         for (uInt ichan = 0; ichan < nchanSp; ichan++) {
           iTsys[ichan] = interpolatorF_->interpolate(fsp[ichan]);
@@ -488 +561 @@
   interpolatorF_->reset();
   interpolatorT_->reset();
-}
-
-Vector<uInt> STApplyCal::timeSort(Vector<Double> &t)
-{
-  Sort sort;
-  sort.sortKey(&t[0], TpDouble, 0, Sort::Ascending);
-  Vector<uInt> idx;
-  sort.sort(idx, t.nelements(), Sort::QuickSort|Sort::NoDuplicates);
-  return idx;
 }
 
@@ -629 +693 @@
   }
 }
-}
+
+}

trunk/src/STApplyCal.h

@@ -96 +96 @@
   casa::Vector<casa::Double> getBaseFrequency(casa::uInt whichrow);
 
-  // time sort
-  casa::Vector<casa::uInt> timeSort(casa::Vector<casa::Double> &t);
-
   // search spwmap_ to get IFNO for Tsys
   casa::uInt getIFForTsys(casa::uInt to);