Changeset 2726 for trunk/src

Timestamp:

01/10/13 20:37:28 (12 years ago)

Author:

Kana Sugimoto

Message:

New Development: No

JIRA Issue: Yes (CAS-4141)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: moved various set methods in python to c++

Test Programs:

Put in Release Notes: No

Module(s): sbseparator

Description:

Moved set methods in sbseparator module to c++ (STSideBandSep class).
Defined boost python interface methods, set_freq, set_dirtol, set_shift,
set_limit, solve_both, subtract_other to access them from python codes.

Location:

trunk/src

Files:

• STSideBandSep.cpp (modified) (3 diffs)
• STSideBandSep.h (modified) (4 diffs)
• python_STSideBandSep.cpp (modified) (2 diffs)

trunk/src/STSideBandSep.cpp

@@ -32 +32 @@
 using namespace asap ;
 
-//#ifndef KS_DEBUG
-//#define KS_DEBUG
-//#endif
+#ifndef KS_DEBUG
+#define KS_DEBUG
+#endif
 
 namespace asap {
 
 // constructors
-STSideBandSep::STSideBandSep()
-  : sigIfno_(0), ftol_(-1), 
-    loTime_(-1), lo1Freq_(-1), loDir_("")
-{
-#ifdef KS_DEBUG
-  cout << "Default constructor STSideBandSep()" << endl;
-#endif
-  // Default LO frame is TOPO
-  loFrame_ = MFrequency::TOPO;
+STSideBandSep::STSideBandSep(const vector<string> &names)
+{
+  LogIO os(LogOrigin("STSideBandSep","STSideBandSep()", WHERE));
+  os << "Setting scantable names to process." << LogIO::POST ;
+  // Set file names
+  ntable_ = names.size();
+  infileList_.resize(ntable_);
+  for (unsigned int i = 0; i < ntable_; i++){
+    if (!checkFile(names[i], "d"))
+      throw( AipsError("File does not exist") );
+    infileList_[i] = names[i];
+  }
+  intabList_.resize(0);
+
+  init();
+
+  {// Summary
+    os << ntable_ << " files are set: [";
+    for (unsigned int i = 0; i < ntable_; i++) {
+      os << " '" << infileList_[i] << "' ";
+      if (i != ntable_-1) os << ",";
+    }
+    os << "] " << LogIO::POST;
+  }
+};
+
+STSideBandSep::STSideBandSep(const vector<ScantableWrapper> &tables)
+{
+  LogIO os(LogOrigin("STSideBandSep","STSideBandSep()", WHERE));
+  os << "Setting list of scantables to process." << LogIO::POST ;
+  // Set file names
+  ntable_ = tables.size();
+  intabList_.resize(ntable_);
+  for (unsigned int i = 0; i < ntable_; i++){
+    intabList_[i] = tables[i].getCP();
+  }
+  infileList_.resize(0);
+
+  init();
+
+  os << ntable_ << " tables are set." << LogIO::POST;
 };
 
@@ -57 +89 @@
 };
 
-
-void STSideBandSep::setFrequency(const unsigned int ifno, const double freqtol, string frame)
-{
-#ifdef KS_DEBUG
-  cout << "STSideBandSep::setFrequency" << endl;
-  cout << "IFNO = " << ifno << endl;
-  cout << "freq tol = " << freqtol << " [Hz]" << endl;
-  cout << "frame = " << frame << endl;
-#endif
+void STSideBandSep::init()
+{
+  // frequency setup
+  sigIfno_= 0;
+  ftol_ = -1;
+  solFrame_ = MFrequency::N_Types;
+  // shifts
+  initshift();
+  // direction tolerance
+  xtol_ = ytol_ = 9.69627e-6; // 2arcsec
+  // solution parameters
+  otherside_ = false;
+  doboth_ = false;
+  rejlimit_ = 0.2;
+  // LO1 values
+  lo1Freq_ = -1;
+  loTime_ = -1;
+  loDir_ = "";
+  // Default LO frame is TOPO
+  loFrame_ = MFrequency::TOPO;
+};
+
+void STSideBandSep::initshift()
+{
+  // shifts
+  nshift_ = 0;
+  nchan_ = 0;
+  sigShift_.resize(0);
+  imgShift_.resize(0);
+  tableList_.resize(0);
+};
+
+void STSideBandSep::setFrequency(const unsigned int ifno,
+                                 const string freqtol,
+                                 const string frame)
+{
+  LogIO os(LogOrigin("STSideBandSep","setFrequency()", WHERE));
+
+  initshift();
+
+  // IFNO
+  sigIfno_ = ifno;
+
+  // Frequency tolerance
+  Quantum<Double> qftol;
+  readQuantity(qftol, String(freqtol));
+  if (!qftol.getUnit().empty()){
+    // make sure the quantity is frequency
+    if (qftol.getFullUnit().getValue() != Unit("Hz").getValue())
+      throw( AipsError("Invalid quantity for frequency tolerance.") );
+    qftol.convert("Hz");
+  }
+  ftol_ = qftol;
+
+  // Frequency Frame
+  if (!frame.empty()){
+    MFrequency::Types mft;
+    if (!MFrequency::getType(mft, frame))
+      throw( AipsError("Invalid frame type.") );
+    solFrame_ = mft;
+  } else {
+    solFrame_ = MFrequency::N_Types;
+  }
+
+  {// Summary
+    const String sframe = ( (solFrame_ == MFrequency::N_Types) ?
+                            "table frame" :
+                            MFrequency::showType(solFrame_) );
+    os << "Frequency setup to search IF group: "
+       << "IFNO of table[0] = " << sigIfno_ 
+       << " , Freq tolerance = " << ftol_.getValue() << " [ "
+       << (ftol_.getUnit().empty() ? "channel" : ftol_.getUnit() )
+       << " ] (in " << sframe <<")" << LogIO::POST;
+  }
+};
+
+
+void STSideBandSep::setDirTolerance(const vector<string> dirtol)
+{
+  LogIO os(LogOrigin("STSideBandSep","setDirTolerance()", WHERE));
+  Quantum<Double> qcell;
+  if ( (dirtol.size() == 1) && !dirtol[0].empty() ) {
+    readQuantity(qcell, String(dirtol[0]));
+    if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
+      xtol_ = ytol_ = qcell.getValue("rad");
+    else
+      throw( AipsError("Invalid unit for direction tolerance.") );
+  }
+  else if (dirtol.size() > 1) {
+    if ( dirtol[0].empty() && dirtol[1].empty() )
+      throw( AipsError("Direction tolerance is empty.") );
+    if ( !dirtol[0].empty() ) {
+      readQuantity(qcell, String(dirtol[0]));
+      if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
+        xtol_ = qcell.getValue("rad");
+      else
+        throw( AipsError("Invalid unit for direction tolerance.") );
+    }
+    if ( !dirtol[1].empty() ) {
+      readQuantity(qcell, String(dirtol[1]));
+      if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
+        ytol_ = qcell.getValue("rad");
+      else
+        throw( AipsError("Invalid unit for direction tolerance.") );
+    }
+    else {
+      ytol_ = xtol_;
+    }
+  }
+  else throw( AipsError("Invalid direction tolerance.") );
+
+  os << "Direction tolerance: ( "
+     << xtol_ << " , " << ytol_ << " ) [rad]" << LogIO::POST;
+};
+
+void STSideBandSep::setShift(const vector<double> &shift)
+{
+  LogIO os(LogOrigin("STSideBandSep","setShift()", WHERE));
+  imgShift_.resize(shift.size());
+  for (unsigned int i = 0; i < shift.size(); i++)
+    imgShift_[i] = shift[i];
+
+  if (imgShift_.size() == 0) {
+    os << "Channel shifts are cleared." << LogIO::POST;
+  } else {
+    os << "Channel shifts of image sideband are set: ( ";
+    for (unsigned int i = 0; i < imgShift_.size(); i++) {
+      os << imgShift_[i];
+      if (i != imgShift_.size()-1) os << " , ";
+    }
+    os << " ) [channel]" << LogIO::POST;
+  }
+};
+
+void STSideBandSep::setThreshold(const double limit)
+{
+  LogIO os(LogOrigin("STSideBandSep","setThreshold()", WHERE));
+  if (limit < 0)
+    throw( AipsError("Rejection limit should be positive number.") );
+
+  rejlimit_ = limit;
+  os << "Rejection limit is set to " << rejlimit_;
 };
 
@@ -81 +246 @@
 
 //
-void STSideBandSep::setLO1(const double lo1, string frame, double reftime, string refdir)
+void STSideBandSep::setLO1(const double lo1, const string frame,
+                           const double reftime, const string refdir)
 {
   lo1Freq_ = lo1;

trunk/src/STSideBandSep.h

@@ -36 +36 @@
    * constructors and a destructor
    **/
-  STSideBandSep();
-  //explicit STSideBandSep(const vector<string> infile);
-  //explicit STSideBandSep(const vector<ScantableWrapper> &tables);
+  STSideBandSep() { throw( AipsError("No data set to process") ); };
+  explicit STSideBandSep(const vector<string> &names);
+  explicit STSideBandSep(const vector<ScantableWrapper> &tables);
   virtual ~STSideBandSep();
 
   /**
-   * Set parameters for sideband separation
+   * Set IFNO and frequency tolerance to select data to process
    **/
-  void setFrequency(const unsigned int ifno, const double freqtol, string frame="");
+  void setFrequency(const unsigned int ifno, const string freqtol,
+                    const string frame="");
+
+  /**
+   * Set direction tolerance to group spectra.
+   * The spectra within this range will be averaged before procesing.
+   **/
+  void setDirTolerance(const vector<string> dirtol);
+
+  /**
+   * Set the number of channels shifted in image side band 
+   * of each of scantable.
+   **/
+  void setShift(const vector<double> &shift);
+
+  /**
+   * Set rejection limit of solution.
+   **/
+  void setThreshold(const double limit);
+
+  /**
+   * Resolve both image and signal sideband when true is set.
+   **/
+  void solveBoth(const bool flag) { doboth_ = flag; };
+
+  /**
+   * Obtain spectra by subtracting the solution of the other sideband.
+   **/
+  void solvefromOther(const bool flag) { otherside_ = flag; };
 
   /**
@@ -51 +79 @@
   void setImageTable(const ScantableWrapper &s);
   void setScanTb0(const ScantableWrapper &s);
+
   /**
    * Set additional information to fill frequencies of image sideband
    **/
-  void setLO1(const double lo1, string frame="TOPO", double reftime=-1, string refdir="");
+  void setLO1(const double lo1, const string frame="TOPO",
+              const double reftime=-1, string refdir="");
   void setLO1Root(const string name);
+
   /**
    * Actual calculation of frequencies of image sideband
@@ -62 +93 @@
 
 private:
+  /** Initialize member variables **/
+  void init();
+  void initshift();
+
+  /** Return if the path exists (optionally, check file type) **/
   Bool checkFile(const string name, string type="");
+
+  /** 
+   * Get LO1 frequency to solve the frequencies of image side band
+   **/
   bool getLo1FromAsdm(const string asdmname,
                       const double refval, const double refpix,
@@ -73 +113 @@
                          const double increment, const int nChan);
 
+  /** Member variables **/
+  // input tables
+  vector<string> infileList_;
+  vector< CountedPtr<Scantable> > intabList_;
+  unsigned int ntable_;
+  // frequency and direction setup to select data.
   unsigned int sigIfno_;
-  double ftol_;
+  Quantum<Double> ftol_;
+  MFrequency::Types solFrame_;
+  vector<double> sigShift_, imgShift_;
+  unsigned int nshift_, nchan_;
+  vector< CountedPtr<Scantable> > tableList_;
+  Double xtol_, ytol_;
+  // solution parameters
+  bool otherside_, doboth_;
+  double rejlimit_;
+  // LO1
   double lo1Freq_;
   MFrequency::Types loFrame_;

trunk/src/python_STSideBandSep.cpp

@@ -10 +10 @@
 #include <boost/python/args.hpp>
 
-
 #include "STSideBandSep.h"
 
@@ -21 +20 @@
   class_<STSideBandSep>("SBSeparator")
     .def( init <> () )
+    .def( init < const std::vector<std::string> > () )
+    .def( init < const std::vector<ScantableWrapper> > () )
     .def( "set_freq", &STSideBandSep::setFrequency,
           (boost::python::arg("frame")="") )
+    .def( "set_dirtol", &STSideBandSep::setDirTolerance )
+    .def( "set_shift", &STSideBandSep::setShift )
+    .def( "set_limit", &STSideBandSep::setThreshold )
+    .def( "solve_both", &STSideBandSep::solveBoth )
+    .def( "subtract_other", &STSideBandSep::solvefromOther )
     .def( "set_lo1", &STSideBandSep::setLO1,
           (boost::python::arg("frame")="TOPO",