Changeset 2928

Timestamp:

04/16/14 12:19:06 (11 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: test_sdcal2, test_tsdcal2

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Improved performance of apply mode by revisiting array structure.

File:

• trunk/src/STApplyCal.cpp (modified) (7 diffs)

trunk/src/STApplyCal.cpp

@@ -206 +206 @@
 
   // list of apply tables for sky calibration
-  Vector<uInt> skycalList;
+  Vector<uInt> skycalList(skytable_.size());
   uInt numSkyCal = 0;
-  uInt nrowSky = 0;
 
   // list of apply tables for Tsys calibration
@@ -214 +213 @@
     STCalEnum::CalType caltype = STApplyTable::getCalType(skytable_[i]);
     if (caltype == caltype_) {
-      skycalList.resize(numSkyCal+1, True);
       skycalList[numSkyCal] = i;
       numSkyCal++;
-      nrowSky += skytable_[i]->nrow();
-    }
-  }
+    }
+  }
+  skycalList.resize(numSkyCal, True);
 
 
@@ -303 +301 @@
 
   uInt nchanSp = skytable_[skylist[0]]->nchan(ifno);
-  Vector<Double> timeSky(nrowSky);
-  Matrix<Float> spoff(nrowSky, nchanSp);
-  Vector<Float> iOff(nchanSp);
-  nrowSky = 0;
-  for (uInt i = 0 ; i < skylist.nelements(); i++) {
-    STCalSkyTable *p = skytable_[skylist[i]];
-    Vector<Double> t = p->getTime();
-    Matrix<Float> sp = p->getSpectra();
-    for (uInt j = 0; j < t.nelements(); j++) {
-      timeSky[nrowSky] = t[j];
-      spoff.row(nrowSky) = sp.column(j);
-      nrowSky++;
-    }
-  }
-
-  Vector<uInt> skyIdx = timeSort(timeSky);
-
-  Double *xa = new Double[skyIdx.nelements()];
-  Float *ya = new Float[skyIdx.nelements()];
-  IPosition ipos(1, skyIdx.nelements());
-  Vector<Double> timeSkySorted(ipos, xa, TAKE_OVER);
-  Vector<Float> tmpOff(ipos, ya, TAKE_OVER);
-  for (uInt i = 0 ; i < skyIdx.nelements(); i++) {
-    timeSkySorted[i] = timeSky[skyIdx[i]];
-  }
-
-  interpolatorS_->setX(xa, skyIdx.nelements());
-
-  Vector<uInt> tsysIdx;
-  Vector<Double> timeTsys(nrowTsys);
-  Matrix<Float> tsys;
+  uInt nrowSkySorted = nrowSky;
+  Vector<Double> timeSkySorted;
+  Matrix<Float> spoffSorted;
+  {
+    Vector<Double> timeSky(nrowSky);
+    Matrix<Float> spoff(nrowSky, nchanSp);
+    nrowSky = 0;
+    for (uInt i = 0 ; i < skylist.nelements(); i++) {
+      STCalSkyTable *p = skytable_[skylist[i]];
+      Vector<Double> t = p->getTime();
+      Matrix<Float> sp = p->getSpectra();
+      for (uInt j = 0; j < t.nelements(); j++) {
+        timeSky[nrowSky] = t[j];
+        spoff.row(nrowSky) = sp.column(j);
+        nrowSky++;
+      }
+    }
+    
+    Vector<uInt> skyIdx = timeSort(timeSky);
+    nrowSkySorted = skyIdx.nelements();
+    
+    timeSkySorted.takeStorage(IPosition(1, nrowSkySorted),
+                              new Double[nrowSkySorted],
+                              TAKE_OVER);
+    for (uInt i = 0 ; i < nrowSkySorted; i++) {
+      timeSkySorted[i] = timeSky[skyIdx[i]];
+    }
+    interpolatorS_->setX(timeSkySorted.data(), nrowSkySorted);
+    
+    spoffSorted.takeStorage(IPosition(2, nrowSky, nchanSp),
+                            new Float[nrowSky * nchanSp],
+                            TAKE_OVER);
+    for (uInt i = 0 ; i < nrowSky; i++) {
+      spoffSorted.row(i) = spoff.row(skyIdx[i]);
+    }
+  }
+
+  uInt nrowTsysSorted = nrowTsys;
+  Matrix<Float> tsysSorted;
   Vector<Double> timeTsysSorted;
-  Vector<Float> tmpTsys;
   if (doTsys) {
     //os_ << "doTsys" << LogIO::POST;
-    timeTsys.resize(nrowTsys);
-    tsys.resize(nrowTsys, nchanTsys);
+    Vector<Double> timeTsys(nrowTsys);
+    Matrix<Float> tsys(nrowTsys, nchanTsys);
+    tsysSorted.takeStorage(IPosition(2, nrowTsys, nchanTsys),
+                           new Float[nrowTsys * nchanTsys],
+                           TAKE_OVER);
     nrowTsys = 0;
     for (uInt i = 0 ; i < tsystable_.size(); i++) {
@@ -351 +359 @@
       }
     }
-    tsysIdx = timeSort(timeTsys);
-
-    Double *xb = new Double[tsysIdx.nelements()];
-    Float *yb = new Float[tsysIdx.nelements()];
-    IPosition ipos(1, tsysIdx.nelements());
-    timeTsysSorted.takeStorage(ipos, xb, TAKE_OVER);
-    tmpTsys.takeStorage(ipos, yb, TAKE_OVER);
-    for (uInt i = 0 ; i < tsysIdx.nelements(); i++) {
+    Vector<uInt> tsysIdx = timeSort(timeTsys);
+    nrowTsysSorted = tsysIdx.nelements();
+
+    timeTsysSorted.takeStorage(IPosition(1, nrowTsysSorted),
+                               new Double[nrowTsysSorted],
+                               TAKE_OVER);
+    for (uInt i = 0 ; i < nrowTsysSorted; i++) {
       timeTsysSorted[i] = timeTsys[tsysIdx[i]];
     }
-    interpolatorT_->setX(xb, tsysIdx.nelements());
+    interpolatorT_->setX(timeTsysSorted.data(), nrowTsysSorted);
+
+    for (uInt i = 0; i < nrowTsys; ++i) {
+      tsysSorted.row(i) = tsys.row(tsysIdx[i]);
+    }
   }
 
@@ -371 +382 @@
 
   // Array for scaling factor (aka Tsys)
-  Vector<Float> iTsys(IPosition(1,nchanSp), new Float[nchanSp], TAKE_OVER);
+  Vector<Float> iTsys(IPosition(1, nchanSp), new Float[nchanSp], TAKE_OVER);
+  // Array for Tsys interpolation
+  // This is empty array and is never referenced if doTsys == false
+  // (i.e. nchanTsys == 0)
+  Vector<Float> iTsysT(IPosition(1, nchanTsys), new Float[nchanTsys], TAKE_OVER);
+
+  // Array for interpolated off spectrum
+  Vector<Float> iOff(IPosition(1, nchanSp), new Float[nchanSp], TAKE_OVER);
   
   for (uInt i = 0; i < rows.nelements(); i++) {
@@ -385 +403 @@
     Double t0 = timeCol(irow);
     for (uInt ichan = 0; ichan < nchanSp; ichan++) {
-      for (uInt j = 0; j < skyIdx.nelements(); j++) {
-        tmpOff[j] = spoff(skyIdx[j], ichan);
-      }
-      interpolatorS_->setY(ya, skyIdx.nelements());
+      Float *tmpY = &(spoffSorted.data()[ichan * nrowSkySorted]); 
+      interpolatorS_->setY(tmpY, nrowSkySorted);
       iOff[ichan] = interpolatorS_->interpolate(t0);
     }
@@ -396 +412 @@
     if (doTsys) {
       // Tsys correction
-      Float *yt = new Float[nchanTsys];
-      Vector<Float> iTsysT(IPosition(1,nchanTsys), yt, TAKE_OVER);
-      Float *yb = tmpTsys.data();
+      Float *yt = iTsysT.data();
       for (uInt ichan = 0; ichan < nchanTsys; ichan++) {
-        for (uInt j = 0; j < tsysIdx.nelements(); j++) {
-          tmpTsys[j] = tsys(tsysIdx[j], ichan);
-        }
-        interpolatorT_->setY(yb, tsysIdx.nelements());
+        Float *tmpY = &(tsysSorted.data()[ichan * nrowTsysSorted]);
+        interpolatorT_->setY(tmpY, nrowTsysSorted);
         iTsysT[ichan] = interpolatorT_->interpolate(t0);
       }