Changeset 2900

Timestamp:

03/12/14 12:23:31 (11 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: Yes CAS-5875

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdcoadd

Put in Release Notes: No

Module(s): sd

Description: Describe your changes here...

Added freq_tol parameter to asapmath.merge.
The freq_tol allows to specify frequency tolerance as
numeric value (1.0e6) or string ('1MHz'). The value will
be used to merge FREQUENCIES rows, FREQ_ID, and IFNO.

Location:

trunk

Files:

• python/asapmath.py (modified) (5 diffs)
• src/STFrequencies.cpp (modified) (2 diffs)
• src/STFrequencies.h (modified) (1 diff)
• src/STMath.cpp (modified) (4 diffs)
• src/STMath.h (modified) (1 diff)
• src/STMathWrapper.h (modified) (1 diff)

trunk/python/asapmath.py

@@ -823 +823 @@
 
 @asaplog_post_dec
-def merge(*args):
+def merge(*args, **kwargs):
     """
     Merge a list of scanatables, or comma-sperated scantables into one
@@ -829 +829 @@
     Parameters:
         A list [scan1, scan2] or scan1, scan2.
+        freq_tol: frequency tolerance for merging IFs. numeric values
+                  in units of Hz (1.0e6 -> 1MHz) and string ('1MHz')
+                  is allowed.
     Example:
         myscans = [scan1, scan2]
@@ -834 +837 @@
         # or equivalent
         sameallscans = merge(scan1, scan2)
+        # with freqtol
+        allscans = merge(scan1, scan2, freq_tol=1.0e6)
+        # or equivalently
+        allscans = merge(scan1, scan2, freq_tol='1MHz')
     """
     varlist = vars()
@@ -842 +849 @@
     else:
         lst = tuple(args)
+    if kwargs.has_key('freq_tol'):
+        freq_tol = str(kwargs['freq_tol'])
+    else:
+        freq_tol = ''
     varlist["args"] = "%d scantables" % len(lst)
     # need special formatting her for history...
@@ -850 +861 @@
             msg = "Please give a list of scantables"
             raise TypeError(msg)
-    s = scantable(stm._merge(lst))
+    s = scantable(stm._merge(lst, freq_tol))
     s._add_history("merge", varlist)
     return s

trunk/src/STFrequencies.cpp

@@ -16 +16 @@
 #include <casa/Arrays/IPosition.h>
 #include <casa/Logging/LogIO.h>
+#include <casa/BasicMath/Math.h>
 
 #include <tables/Tables/TableDesc.h>
@@ -450 +451 @@
 }
 
+bool STFrequencies::match( Double refpix, Double refval,
+                           Double inc, Double freqTolInHz,
+                           uInt &id)
+{
+  ROScalarColumn<uInt> idCol(table_, "ID");
+  ROScalarColumn<Double> refPixCol(table_, "REFPIX");
+  ROScalarColumn<Double> refValCol(table_, "REFVAL");
+  ROScalarColumn<Double> incCol(table_, "INCREMENT");
+  for (uInt irow = 0; irow < table_.nrow(); ++irow) {
+    Double refInc = incCol(irow);
+    Double refFreq = refValCol(irow) - refPixCol(irow) * refInc;
+    Double freq0 = refval - refpix * inc;
+    if (nearAbs(inc, refInc, freqTolInHz) &&
+        nearAbs(refFreq, freq0, freqTolInHz)) {
+      id = irow;
+      return true;
+    }
+  }
+  return false;
+}
+  
 } // namespace

trunk/src/STFrequencies.h

@@ -172 +172 @@
   const casa::String& name() const { return name_; }
 
+  /**
+   * Examine given set of refpix, refval, and increment matches
+   * any of the rows within a tolerance of freqTolInHz. If match,
+   * return true and id is filled properly. Otherwise, return false
+   * and id may have invalid value.
+   *
+   * @param[in] refpix
+   * @param[in] refval
+   * @param[in] inc
+   * @param[in] freqTolInHz
+   * @param[out] id
+   * @return boolean indicating match with any rows or not
+   */
+  bool match( casa::Double refpix, casa::Double refval, casa::Double inc,
+              casa::Double freqTolInHz, casa::uInt &id);
+
 private:
 

trunk/src/STMath.cpp

@@ -2541 +2541 @@
 
 CountedPtr< Scantable >
-STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
-//STMath::merge( const std::vector< CountedPtr < Scantable > >& in,
-//             const std::string &freqTol )
-{
-  /**
-  LogIO os;
-  Double freqTolInHz = 0.0;
+STMath::merge( const std::vector< CountedPtr < Scantable > >& in,
+               const std::string &freqTol )
+{
+  Double freqTolInHz = 1.0; // default is 1.0Hz according to concat task
   if (freqTol.size() > 0) {
     Quantum<Double> freqTolInQuantity;
@@ -2558 +2555 @@
     freqTolInHz = freqTolInQuantity.getValue("Hz");
   }
-  **/
   
   if ( in.size() < 2 ) {
@@ -2579 +2575 @@
   uInt newscanno = max(scannos)+1;
   ++it;
+
+  // new IFNO
+  uInt ifnoCounter = max(ifnocol.getColumn()) + 1;
+  
   while ( it != in.end() ){
-    /**
     // Check FREQUENCIES/BASEFRAME
     if ( out->frequencies().getFrame(true) != (*it)->frequencies().getFrame(true) ) {
       throw(AipsError("BASEFRAME is not identical"));
     }
-    **/
     
     if ( ! (*it)->conformant(*out) ) {
@@ -2621 +2619 @@
         (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
         uInt id;
-        id = out->frequencies().addEntry(rp, rv, inc);
-        //if ( !out->frequencies().match(rp, rv, inc, freqTolInHz, id) ) {
-        //  id = out->frequencies().addEntry(rp, rv, inc);
-        //}
+        
+        // default value is new unique IFNO
+        uInt newifno = ifnoCounter;
+        uInt nchan = rec.asArrayFloat("SPECTRA").shape()[0];
+        //id = out->frequencies().addEntry(rp, rv, inc);
+        if ( !out->frequencies().match(rp, rv, inc, freqTolInHz, id) ) {
+          // add new entry to FREQUENCIES table
+          id = out->frequencies().addEntry(rp, rv, inc);
+
+          // increment counter for IFNO
+          ifnoCounter++;
+        }
+        else {
+          // should renumber IFNO to be same as existing rows that have same FREQ_ID
+          LogIO os(LogOrigin("STMath", "merge", WHERE));
+          Table outFreqIdSelected = tout(tout.col("FREQ_ID") == id);
+          TableIterator _iter(outFreqIdSelected, "IFNO");
+          map<uInt, uInt> nchanMap;
+          while (!_iter.pastEnd()) {
+            const Table _table = _iter.table();
+            ROTableRow _row(_table);
+            const TableRecord &_rec = _row.get(0); 
+            uInt nchan = _rec.asArrayFloat("SPECTRA").shape()[0];
+            if (nchanMap.find(nchan) != nchanMap.end()) {
+              throw(AipsError("There are non-unique IFNOs assigned to spectra that have same FREQ_ID and same nchan. Something wrong."));
+            }
+            nchanMap[nchan] = _rec.asuInt("IFNO");
+            _iter.next();
+          }
+
+          os << LogIO::DEBUGGING << "nchanMap for " << id << ":" << LogIO::POST;
+          for (map<uInt, uInt>::iterator i = nchanMap.begin(); i != nchanMap.end(); ++i) {
+            os << LogIO::DEBUGGING << "nchanMap[" << i->first << "] = " << i->second << LogIO::POST;
+          }
+
+          if (nchanMap.find(nchan) == nchanMap.end()) {
+            // increment counter for IFNO
+            ifnoCounter++;
+          }
+          else {
+            // renumber IFNO to be same as existing value that corresponds to nchan
+            newifno = nchanMap[nchan];
+          }
+          os << LogIO::DEBUGGING << "newifno = " << newifno << LogIO::POST;
+        }
         thecolvals = id;
         freqidcol.putColumnRange(slice, thecolvals);
+
+        thecolvals = newifno;
+        ifnocol.putColumnRange(slice, thecolvals);
         
         // Set the proper MOLECULE_ID

trunk/src/STMath.h

@@ -310 +310 @@
 
   casa::CountedPtr<Scantable>
-    merge(const std::vector<casa::CountedPtr<Scantable> >& in);
+    merge(const std::vector<casa::CountedPtr<Scantable> >& in,
+          const std::string &freqTol = "");
 
   casa::CountedPtr<Scantable>

trunk/src/STMathWrapper.h

@@ -175 +175 @@
 
   ScantableWrapper
-    merge(const std::vector<ScantableWrapper >& in)
-
+    merge(const std::vector<ScantableWrapper >& in,
+          const std::string &freqTol)
   {
     std::vector<casa::CountedPtr<Scantable> > sts;
     for (unsigned int i=0; i<in.size(); ++i) sts.push_back(in[i].getCP());
-    return ScantableWrapper(STMath::merge(sts)); }
+    return ScantableWrapper(STMath::merge(sts, freqTol)); }
 
   ScantableWrapper rotateXYPhase( const ScantableWrapper& in, float angle)