Changeset 2737

Timestamp:

01/18/13 19:09:41 (11 years ago)

Author:

WataruKawasaki

Message:

New Development: No

JIRA Issue: Yes CAS-4794

Ready for Test: No

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sd

Description:

Location:

trunk/src

Files:

• STBaselineEnum.h (modified) (1 diff)
• STBaselineParamTable.cpp (deleted)
• STBaselineParamTable.h (deleted)
• Scantable.cpp (modified) (15 diffs)
• Scantable.h (modified) (2 diffs)

trunk/src/STBaselineEnum.h

@@ -26 +26 @@
 @version $Revision:$
 */
-class STBaselineEnum  {
+class STBaselineFunc  {
 public:
-  enum BaselineType {Polynomial = 0,
-                     CSpline,
-                     Sinusoid,
-                     Chebyshev};
+  enum FuncName {Polynomial = 0,
+                 CSpline,
+                 Sinusoid,
+                 Chebyshev};
 };
 

trunk/src/Scantable.cpp

@@ -66 +66 @@
 #include "MathUtils.h"
 #include "STAttr.h"
+#include "STBaselineTable.h"
 #include "STLineFinder.h"
 #include "STPolCircular.h"
@@ -2637 +2638 @@
     std::vector<bool> chanMask;
 
+    std::vector<float> sect(0);
+    std::vector<bool> finalChanMask;
+    float rms;
+
+    //---
+    bool outBaselineParamTable = true;
+    STBaselineTable* bt = new STBaselineTable(*this);
+    //---
+
     for (int whichrow = 0; whichrow < nRow; ++whichrow) {
       std::vector<float> sp = getSpectrum(whichrow);
@@ -2642 +2652 @@
       std::vector<float> params(order+1);
       int nClipped = 0;
-      std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, nClipped, thresClip, nIterClip, getResidual);
-
-      setSpectrum(res, whichrow);
+      std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, finalChanMask, rms, nClipped, thresClip, nIterClip, getResidual);
+
+      //---
+      if (outBaselineParamTable) {
+        bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)), 
+                       uInt(0), 
+                       Double(0.0), // <-- Double(getTime(whichrow, false)), 
+                       uInt(nchan(whichrow)), STBaselineFunc::Chebyshev, 
+                       uInt(order), uInt(nIterClip), Float(thresClip), Vector<Float>(sect), Vector<Float>(params), 
+                       Vector<Float>(params), // <-- Vector<Float>(finalChanMask), 
+                       Float(rms));
+      } else {
+        setSpectrum(res, whichrow);
+      }
+      //---
+
       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "chebyshevBaseline()", params, nClipped);
       showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     }
     
+    //---
+    if (outBaselineParamTable) {
+      bt->save("chebyparamtable");
+    }
+    //---
+    delete bt;
+
     if (outTextFile) ofs.close();
 
@@ -2697 +2727 @@
   int nparam;
   std::vector<float> params;
+  std::vector<bool> finalChanMask;
+  float rms;
   int nClipped = 0;
   float thresClip = 0.0;
@@ -2703 +2735 @@
   if (blfunc == "chebyshev") {
     nparam = order + 1;
-    res = doChebyshevFitting(spec, mask, order, params, nClipped, thresClip, nIterClip, true);
+    res = doChebyshevFitting(spec, mask, order, params, finalChanMask, rms, nClipped, thresClip, nIterClip, true);
   } else if (blfunc == "sinusoid") {
     std::vector<int> nWaves;
@@ -2798 +2830 @@
     int minNRow;
     parseProgressInfo(progressInfo, showProgress, minNRow);
+
+    std::vector<bool> finalChanMask;
+    float rms;
 
     for (int whichrow = 0; whichrow < nRow; ++whichrow) {
@@ -2831 +2866 @@
       std::vector<float> params(order+1);
       int nClipped = 0;
-      std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, nClipped, thresClip, nIterClip, getResidual);
+      std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, finalChanMask, rms, nClipped, thresClip, nIterClip, getResidual);
       setSpectrum(res, whichrow);
 
@@ -2901 +2936 @@
 }
 
-std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
+std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, std::vector<bool>& finalMask, float& rms, int& nClipped, float thresClip, int nIterClip, bool getResidual)
 {
   if (data.size() != mask.size()) {
@@ -2911 +2946 @@
 
   int nChan = data.size();
+
+  finalMask.clear();
+  finalMask.resize(nChan);
+
   std::vector<int> maskArray;
   std::vector<int> x;
@@ -2918 +2957 @@
       x.push_back(i);
     }
+    finalMask[i] = mask[i];
   }
 
@@ -3041 +3081 @@
     }
 
+    double stdDev = 0.0;
+    for (int i = 0; i < nChan; ++i) {
+      stdDev += residual[i]*residual[i]*(double)maskArray[i];
+    }
+    stdDev = sqrt(stdDev/(double)nData);
+    rms = (float)stdDev;
+
     if ((nClip == nIterClip) || (thresClip <= 0.0)) {
       break;
     } else {
-      double stdDev = 0.0;
-      for (int i = 0; i < nChan; ++i) {
-        stdDev += residual[i]*residual[i]*(double)maskArray[i];
-      }
-      stdDev = sqrt(stdDev/(double)nData);
-      
+
       double thres = stdDev * thresClip;
       int newNData = 0;
@@ -3055 +3097 @@
         if (abs(residual[i]) >= thres) {
           maskArray[i] = 0;
+          finalMask[i] = false;
         }
         if (maskArray[i] > 0) {
@@ -3065 +3108 @@
         nData = newNData;
       }
+
     }
   }
@@ -4315 +4359 @@
 {
   /****
-  double ms1TimeStart, ms1TimeEnd, ms2TimeStart, ms2TimeEnd;
+  double ms1TimeStart, ms1TimeEnd;
   double elapse1 = 0.0;
-  double elapse2 = 0.0;
-
   ms1TimeStart = mathutil::gettimeofday_sec();
   ****/
 
   Vector<Float> spec;
-
   specCol_.get(whichrow, spec);
 
@@ -4330 +4371 @@
   elapse1 = ms1TimeEnd - ms1TimeStart;
   std::cout << "rm1   : " << elapse1 << " (sec.)" << endl;
-  ms2TimeStart = mathutil::gettimeofday_sec();
   ****/
 
-  float mean = 0.0;
-  float smean = 0.0;
+  return (float)doGetRms(mask, spec);
+}
+
+double Scantable::doGetRms(const std::vector<bool>& mask, const Vector<Float>& spec)
+{
+  double mean = 0.0;
+  double smean = 0.0;
   int n = 0;
   for (uInt i = 0; i < spec.nelements(); ++i) {
     if (mask[i]) {
-      mean += spec[i];
-      smean += spec[i]*spec[i];
+      double val = (double)spec[i];
+      mean += val;
+      smean += val*val;
       n++;
     }
   }
 
-  mean /= (float)n;
-  smean /= (float)n;
-
-  /****
-  ms2TimeEnd = mathutil::gettimeofday_sec();
-  elapse2 = ms2TimeEnd - ms2TimeStart;
-  std::cout << "rm2   : " << elapse2 << " (sec.)" << endl;
-  ****/
+  mean /= (double)n;
+  smean /= (double)n;
 
   return sqrt(smean - mean*mean);
 }
-
 
 std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose, bool csvformat) const

trunk/src/Scantable.h

@@ -748 +748 @@
                                           int order,
                                           std::vector<float>& params,
+                                          std::vector<bool>& finalMask, 
+                                          float& rms, 
                                           int& nClipped,
                                           float thresClip=3.0, 
@@ -818 +820 @@
                                            const std::string& blfunc, 
                                            int order);
+  double doGetRms(const std::vector<bool>& mask, const casa::Vector<casa::Float>& spec);
 
 };