Changeset 2058 for trunk/src

Timestamp:

03/23/11 16:23:45 (14 years ago)

Author:

WataruKawasaki

Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sdbaseline

Description: bugfix for Scantable::doCubicSplineBaseline()

Location:

trunk/src

Files:

• Scantable.cpp (modified) (12 diffs)
• Scantable.h (modified) (2 diffs)

trunk/src/Scantable.cpp

@@ -1885 +1885 @@
     std::vector<int> pieceRanges;
     std::vector<float> params;
-    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
+    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     setSpectrum(res, whichrow);
     //
@@ -1954 +1954 @@
     std::vector<int> pieceRanges;
     std::vector<float> params;
-    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
+    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     setSpectrum(res, whichrow);
     //
@@ -1966 +1966 @@
 }
 
-std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& sectionRanges, std::vector<float>& params, int nPiece, float thresClip, int nIterClip) {
+std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& sectionRanges, std::vector<float>& params, int nPiece, float thresClip, int nIterClip, bool getResidual) {
   if (nPiece < 1) {
     throw(AipsError("wrong number of the sections for fitting"));
@@ -2004 +2004 @@
   sectionRanges.push_back(x[x.size()-1]);
 
-  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
+  std::vector<double> x1, z1, r1;
   for (int i = 0; i < nChan; ++i) {
     x1.push_back((double)i);
-    x2.push_back((double)(i*i));
-    x3.push_back((double)(i*i*i));
     z1.push_back((double)data[i]);
-    x1z1.push_back(((double)i)*(double)data[i]);
-    x2z1.push_back(((double)(i*i))*(double)data[i]);
-    x3z1.push_back(((double)(i*i*i))*(double)data[i]);
     r1.push_back(0.0);
   }
@@ -2036 +2031 @@
       for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
         if (maskArray[i] == 0) continue;
+        double xp1 = x1[i];
+        double xp2 = xp1*xp1;
+        double xp3 = xp2*xp1;
+        double xp4 = xp3*xp1;
+        double xp5 = xp4*xp1;
+        double xp6 = xp5*xp1;
         xMatrix[0][0] += 1.0;
-        xMatrix[0][1] += x1[i];
-        xMatrix[0][2] += x2[i];
-        xMatrix[0][3] += x3[i];
-        xMatrix[1][1] += x2[i];
-        xMatrix[1][2] += x3[i];
-        xMatrix[1][3] += x2[i]*x2[i];
-        xMatrix[2][2] += x2[i]*x2[i];
-        xMatrix[2][3] += x3[i]*x2[i];
-        xMatrix[3][3] += x3[i]*x3[i];
+        xMatrix[0][1] += xp1;
+        xMatrix[0][2] += xp2;
+        xMatrix[0][3] += xp3;
+        xMatrix[1][1] += xp2;
+        xMatrix[1][2] += xp3;
+        xMatrix[1][3] += xp4;
+        xMatrix[2][2] += xp4;
+        xMatrix[2][3] += xp5;
+        xMatrix[3][3] += xp6;
         zMatrix[0] += z1[i];
-        zMatrix[1] += x1z1[i];
-        zMatrix[2] += x2z1[i];
-        zMatrix[3] += x3z1[i];
+        zMatrix[1] += xp1*z1[i];
+        zMatrix[2] += xp2*z1[i];
+        zMatrix[3] += xp3*z1[i];
         for (int j = 0; j < n; ++j) {
-          double q = 1.0 - x1[i]*sectionListR[j];
+          double q = 1.0 - xp1*sectionListR[j];
           q = q*q*q;
           xMatrix[0][j+4] += q;
-          xMatrix[1][j+4] += q*x1[i];
-          xMatrix[2][j+4] += q*x2[i];
-          xMatrix[3][j+4] += q*x3[i];
+          xMatrix[1][j+4] += q*xp1;
+          xMatrix[2][j+4] += q*xp2;
+          xMatrix[3][j+4] += q*xp3;
           for (int k = 0; k < j; ++k) {
-            double r = 1.0 - x1[i]*sectionListR[k];
+            double r = 1.0 - xp1*sectionListR[k];
             r = r*r*r;
             xMatrix[k+4][j+4] += r*q;
@@ -2109 +2110 @@
     //cubic terms for the other pieces (in case nPiece>1).
     std::vector<double> y;
+    y.clear();
     for (int i = 0; i < nDOF; ++i) {
       y.push_back(0.0);
@@ -2124 +2126 @@
     for (int n = 0; n < nPiece; ++n) {
       for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
-        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
+        r1[i] = a0 + a1*x1[i] + a2*x1[i]*x1[i] + a3*x1[i]*x1[i]*x1[i];
       }
       params.push_back(a0);
@@ -2170 +2172 @@
   }
 
-  std::vector<float> residual;
-  for (int i = 0; i < nChan; ++i) {
-    residual.push_back((float)(z1[i] - r1[i]));
-  }
-  return residual;
-
+  std::vector<float> result;
+  if (getResidual) {
+    for (int i = 0; i < nChan; ++i) {
+      result.push_back((float)(z1[i] - r1[i]));
+    }
+  } else {
+    for (int i = 0; i < nChan; ++i) {
+      result.push_back((float)r1[i]);
+    }
+  }
+
+  return result;
 }
 
@@ -2207 +2215 @@
     std::vector<int> pieceRanges;
     std::vector<float> params;
-    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip);
+    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     setSpectrum(res, whichrow);
     //
@@ -2276 +2284 @@
     std::vector<int> pieceRanges;
     std::vector<float> params;
-    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip);
+    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     setSpectrum(res, whichrow);
     //
@@ -2288 +2296 @@
 }
 
-std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, std::vector<float>& params, float thresClip, int nIterClip) {
-  if (nMaxWavesInSW < 1) {
-    throw(AipsError("maximum wave number must be a positive value"));
-  }
+std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual) {
   if (data.size() != mask.size()) {
     throw(AipsError("data and mask have different size"));
+  }
+  if ((nMinWavesInSW < 0) || (nMaxWavesInSW < 0)) {
+    throw(AipsError("wave number must be positive or zero (i.e. constant)"));
+  } else {
+    if (nMaxWavesInSW < nMinWavesInSW) {
+      int temp = nMaxWavesInSW;
+      nMaxWavesInSW = nMinWavesInSW;
+      nMinWavesInSW = temp;
+    }
   }
 
@@ -2448 +2462 @@
   }
 
-  std::vector<float> residual;
-  for (int i = 0; i < nChan; ++i) {
-    residual.push_back((float)(z1[i] - r1[i]));
-  }
-  return residual;
-
+  std::vector<float> result;
+  if (getResidual) {
+    for (int i = 0; i < nChan; ++i) {
+      result.push_back((float)(z1[i] - r1[i]));
+    }
+  } else {
+    for (int i = 0; i < nChan; ++i) {
+      result.push_back((float)r1[i]);
+    }
+  }
+
+  return result;
 }
 

trunk/src/Scantable.h

@@ -673 +673 @@
                                           int nPiece=2,
                                           float thresClip=3.0, 
-                                          int nIterClip=1);
+                                          int nIterClip=1,
+                                          bool getResidual=false);
   std::vector<float> doSinusoidFitting(const std::vector<float>& data, 
                                        const std::vector<bool>& mask,
@@ -680 +681 @@
                                        std::vector<float>& params,
                                        float thresClip=3.0, 
-                                       int nIterClip=1);
+                                       int nIterClip=1,
+                                       bool getResidual=false);
   bool hasSameNchanOverIFs();
   std::string getMaskRangeList(const std::vector<bool>& mask,