Changeset 2727

Timestamp:

01/11/13 18:48:37 (11 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: Yes CAS-4770, CAS-4774

Ready for Test: Yes

Interface Changes: 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...

Updated STApplyCal to be able to specify interpolation method.
The method can be specified in time and frequency axes independently.
Possible options are nearest, linear (default), (natural) cubic spline,
and polynomial with arbitrary order.

Location:

trunk/src

Files:

• BisectionLocator.cpp (modified) (2 diffs)
• BisectionLocator.h (modified) (1 diff)
• BufferedBisectionLocator.cpp (added)
• BufferedBisectionLocator.h (added)
• BufferedLinearInterpolator1D.cpp (added)
• BufferedLinearInterpolator1D.h (added)
• CMakeLists.txt (modified) (1 diff)
• CubicSplineInterpolator1D.cpp (added)
• CubicSplineInterpolator1D.h (added)
• HuntLocator.cpp (added)
• HuntLocator.h (added)
• Interpolator1D.cpp (modified) (6 diffs)
• Interpolator1D.h (modified) (3 diffs)
• LinearInterpolator1D.cpp (added)
• LinearInterpolator1D.h (added)
• Locator.cpp (modified) (1 diff)
• Locator.h (modified) (1 diff)
• PolynomialInterpolator1D.cpp (added)
• PolynomialInterpolator1D.h (added)
• STApplyCal.cpp (modified) (4 diffs)
• STApplyCal.h (modified) (1 diff)
• STCalEnum.h (modified) (1 diff)
• STCalSkyTable.cpp (modified) (1 diff)

trunk/src/BisectionLocator.cpp

@@ -12 +12 @@
 #include <assert.h>
 
-#include <casa/Arrays/Vector.h>
-#include <casa/Arrays/ArrayIO.h>
-#include <casa/Exceptions/Error.h>
-
 #include "BisectionLocator.h"
-
-using namespace casa;
 
 namespace asap {
@@ -39 +33 @@
     else if (x > x_[n_-1])
       return n_;
+
+    unsigned int jl = 0;
+    unsigned int ju = n_;
+    unsigned int jm;
+    while (ju - jl > 1) {
+      jm = (ju + jl) >> 1;
+      if (x > x_[jm])
+        jl = jm;
+      else
+        ju = jm;
+    }
+    return ju;
   }
   else {
-    if (x > x_[0])
+    if (x >= x_[0])
       return 0;
-    else if (x <= x_[n_-1])
+    else if (x < x_[n_-1])
       return n_;
+
+    unsigned int jl = 0;
+    unsigned int ju = n_;
+    unsigned int jm;
+    while (ju - jl > 1) {
+      jm = (ju + jl) >> 1;
+      if (x < x_[jm])
+        jl = jm;
+      else
+        ju = jm;
+    }
+    return ju;
   }
-  unsigned int jl = 0;
-  unsigned int ju = n_;
-  unsigned int jm;
-
-  while (ju - jl > 1) {
-    jm = (ju + jl) >> 1;
-    if ((x >= x_[jm]) == ascending)
-      jl = jm;
-    else
-      ju = jm;
-  }
-  return ju;
 }
 

trunk/src/BisectionLocator.h

@@ -12 +12 @@
 #ifndef ASAP_BISECTION_LOCATOR_H
 #define ASAP_BISECTION_LOCATOR_H
-
-#include <memory>
-#include <vector>
-
-#include <casa/aips.h>
-#include <casa/Arrays/Vector.h>
-#include <casa/Arrays/Matrix.h>
-#include <casa/BasicSL/String.h>
-#include <casa/Utilities/CountedPtr.h>
 
 #include "Locator.h"

trunk/src/CMakeLists.txt

@@ -82 +82 @@
      ${SRCDIR}/Locator.cpp
      ${SRCDIR}/BisectionLocator.cpp
+     ${SRCDIR}/HuntLocator.cpp
+     ${SRCDIR}/BufferedBisectionLocator.cpp
      ${SRCDIR}/Interpolator1D.cpp
-     ${SRCDIR}/NearestInterpolator1D.cpp )
+     ${SRCDIR}/NearestInterpolator1D.cpp
+     ${SRCDIR}/LinearInterpolator1D.cpp
+     ${SRCDIR}/BufferedLinearInterpolator1D.cpp
+     ${SRCDIR}/CubicSplineInterpolator1D.cpp
+     ${SRCDIR}/PolynomialInterpolator1D.cpp )
 
 set( ASAP_PYSRCS

trunk/src/Interpolator1D.cpp

@@ -17 +17 @@
 #include "Interpolator1D.h"
 #include "BisectionLocator.h"
+#include "HuntLocator.h"
 
 using namespace casa;
@@ -23 +24 @@
 
 Interpolator1D::Interpolator1D()
-  : order_(0),
+  : order_(1),
     n_(0),
     x_(0),
-    y_(0)
+    y_(0),
+    locator_(0)
 {
-  locator_ = new BisectionLocator();
 }
 
@@ -42 +43 @@
   y_ = y;
   n_ = n;
+  createLocator();
   locator_->set(x, n);
 }
@@ -50 +52 @@
   x_ = x;
   n_ = n;
+  createLocator();
   locator_->set(x, n);
 }
@@ -65 +68 @@
   x_ = 0;
   y_ = 0;
+  if (locator_) {
+    delete locator_;
+    locator_ = 0;
+  }
 }
 
@@ -77 +84 @@
 }
 
+void Interpolator1D::createLocator()
+{
+  if (!locator_) {
+    if (n_ > 1000)
+      locator_ = new HuntLocator();
+    else
+      locator_ = new BisectionLocator();
+  }
 }
+
+}

trunk/src/Interpolator1D.h

@@ -12 +12 @@
 #ifndef ASAP_INTERPOLATOR_1D_H
 #define ASAP_INTERPOLATOR_1D_H
-
-#include <memory>
-#include <vector>
-
-#include <casa/aips.h>
-#include <casa/Containers/Block.h>
-#include <casa/Arrays/Vector.h>
-#include <casa/Arrays/Matrix.h>
-#include <casa/BasicSL/String.h>
-#include <casa/Utilities/CountedPtr.h>
 
 #include "Locator.h"
@@ -41 +31 @@
   void setY(float *y, unsigned int n);
   void reset();
+
+  // currently only effective for polynomial interpolation
   void setOrder(unsigned int order) {order_ = order;}
 
@@ -48 +40 @@
   int locate(double x);
   bool isready();
+  void createLocator();
 
   unsigned int order_;

trunk/src/Locator.cpp

@@ -12 +12 @@
 #include <assert.h>
 
-#include <casa/Arrays/Vector.h>
-#include <casa/Exceptions/Error.h>
-
 #include "Locator.h"
-
-using namespace casa;
 
 namespace asap {

trunk/src/Locator.h

@@ -12 +12 @@
 #ifndef ASAP_LOCATOR_H
 #define ASAP_LOCATOR_H
-
-#include <memory>
-#include <vector>
-
-#include <casa/aips.h>
-#include <casa/Arrays/Vector.h>
-#include <casa/Arrays/Matrix.h>
-#include <casa/BasicSL/String.h>
-#include <casa/Utilities/CountedPtr.h>
 
 namespace asap {

trunk/src/STApplyCal.cpp

@@ -34 +34 @@
 #include "PSAlmaCalibrator.h"
 #include "NearestInterpolator1D.h"
+#include "BufferedLinearInterpolator1D.h"
+#include "PolynomialInterpolator1D.h"
+#include "CubicSplineInterpolator1D.h"
 #include <atnf/PKSIO/SrcType.h>
 
@@ -62 +65 @@
   doTsys_ = False;
   interp_.resize((int)STCalEnum::NumAxis);
+  // default is linear interpolation
+  for (unsigned int i = 0; i < interp_.size(); i++) {
+    interp_[i] = STCalEnum::LinearInterpolation;
+  }
 }
 
@@ -125 +132 @@
 
   // interpolator
-  interpolatorS_ = new NearestInterpolator1D();
-  interpolatorT_ = new NearestInterpolator1D();
-  interpolatorF_ = new NearestInterpolator1D();
+  initInterpolator();
 
   // select data
@@ -450 +455 @@
 }
 
-}
+void STApplyCal::initInterpolator()
+{
+  int ta = (int)STCalEnum::TimeAxis;
+  int fa = (int)STCalEnum::FrequencyAxis;
+  int order = (order_ > 0) ? order_ : 1;
+  switch (interp_[ta]) {
+  case STCalEnum::NearestInterpolation:
+    {
+      os_ << "use NearestInterpolator in time axis" << LogIO::POST;
+      interpolatorS_ = new NearestInterpolator1D();
+      interpolatorT_ = new NearestInterpolator1D();
+      break;
+    }
+  case STCalEnum::LinearInterpolation:
+    {
+      os_ << "use BufferedLinearInterpolator in time axis" << LogIO::POST;
+      interpolatorS_ = new BufferedLinearInterpolator1D();
+      interpolatorT_ = new BufferedLinearInterpolator1D();
+      break;      
+    }
+  case STCalEnum::CubicSplineInterpolation:
+    {
+      os_ << "use CubicSplineInterpolator in time axis" << LogIO::POST;
+      interpolatorS_ = new CubicSplineInterpolator1D();
+      interpolatorT_ = new CubicSplineInterpolator1D();
+      break;
+    }
+  case STCalEnum::PolynomialInterpolation:
+    {
+      os_ << "use PolynomialInterpolator in time axis" << LogIO::POST;
+      if (order == 0) {
+        interpolatorS_ = new NearestInterpolator1D();
+        interpolatorT_ = new NearestInterpolator1D();
+      }
+      else {
+        interpolatorS_ = new PolynomialInterpolator1D();
+        interpolatorT_ = new PolynomialInterpolator1D();
+        interpolatorS_->setOrder(order);
+        interpolatorT_->setOrder(order);
+      }
+      break;
+    }
+  default:
+    {
+      os_ << "use BufferedLinearInterpolator in time axis" << LogIO::POST;
+      interpolatorS_ = new BufferedLinearInterpolator1D();
+      interpolatorT_ = new BufferedLinearInterpolator1D();
+      break;      
+    }
+  }
+   
+  switch (interp_[fa]) {
+  case STCalEnum::NearestInterpolation:
+    {
+      os_ << "use NearestInterpolator in frequency axis" << LogIO::POST;
+      interpolatorF_ = new NearestInterpolator1D();
+      break;
+    }
+  case STCalEnum::LinearInterpolation:
+    {
+      os_ << "use BufferedLinearInterpolator in frequency axis" << LogIO::POST;
+      interpolatorF_ = new BufferedLinearInterpolator1D();
+      break;      
+    }
+  case STCalEnum::CubicSplineInterpolation:
+    {
+      os_ << "use CubicSplineInterpolator in frequency axis" << LogIO::POST;
+      interpolatorF_ = new CubicSplineInterpolator1D();
+      break;
+    }
+  case STCalEnum::PolynomialInterpolation:
+    {
+      os_ << "use PolynomialInterpolator in frequency axis" << LogIO::POST;
+      if (order == 0) {
+        interpolatorF_ = new NearestInterpolator1D();
+      }
+      else {
+        interpolatorF_ = new PolynomialInterpolator1D();
+        interpolatorF_->setOrder(order);
+      }
+      break;
+    }
+  default:
+    {
+      os_ << "use LinearInterpolator in frequency axis" << LogIO::POST;
+      interpolatorF_ = new BufferedLinearInterpolator1D();
+      break;      
+    }
+  }
+}
+}

trunk/src/STApplyCal.h

@@ -73 +73 @@
   void init();
 
+  // setup interpolator
+  void initInterpolator();
+
   // single loop element in apply()
   void doapply(casa::uInt beamno, casa::uInt ifno, casa::uInt polno, 

trunk/src/STCalEnum.h

@@ -33 +33 @@
                           LinearInterpolation,
                           PolynomialInterpolation,
-                          CSplineInterpolation};
+                          CubicSplineInterpolation};
   enum InterpolationAxis {TimeAxis = 0,
                           FrequencyAxis,

trunk/src/STCalSkyTable.cpp

@@ -10 +10 @@
 //
 //
+#include <assert.h>
+
 #include <casa/Exceptions/Error.h>
 #include <casa/Logging/LogIO.h>