Changeset 2064 for trunk


Ignore:
Timestamp:
03/24/11 18:04:56 (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: minor clean-up for Scantable::doCubicSplineFitting()


Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/Scantable.cpp

    r2058 r2064  
    19661966}
    19671967
    1968 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) {
     1968std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& idxEdge, std::vector<float>& params, int nPiece, float thresClip, int nIterClip, bool getResidual) {
    19691969  if (nPiece < 1) {
    19701970    throw(AipsError("wrong number of the sections for fitting"));
    … …  
    19861986  int nData = x.size();
    19871987  int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
    1988 
    1989   std::vector<int> sectionList0, sectionList1;
    1990   std::vector<double> sectionListR;
    1991   sectionList0.push_back(x[0]);
    1992   sectionRanges.clear();
    1993   sectionRanges.push_back(x[0]);
     1988  std::vector<double> invEdge;
     1989
     1990  idxEdge.clear();
     1991  idxEdge.push_back(x[0]);
     1992
    19941993  for (int i = 1; i < nPiece; ++i) {
    19951994    int valX = x[nElement*i];
    1996     sectionList0.push_back(valX);
    1997     sectionList1.push_back(valX);
    1998     sectionListR.push_back(1.0/(double)valX);
    1999 
    2000     sectionRanges.push_back(valX-1);
    2001     sectionRanges.push_back(valX);
    2002   }
    2003   sectionList1.push_back(x[x.size()-1]+1);
    2004   sectionRanges.push_back(x[x.size()-1]);
    2005 
    2006   std::vector<double> x1, z1, r1;
     1995    idxEdge.push_back(valX);
     1996    invEdge.push_back(1.0/(double)valX);
     1997  }
     1998
     1999  idxEdge.push_back(x[x.size()-1]+1);
     2000
     2001  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
    20072002  for (int i = 0; i < nChan; ++i) {
    2008     x1.push_back((double)i);
    2009     z1.push_back((double)data[i]);
     2003    double di = (double)i;
     2004    double dD = (double)data[i];
     2005    x1.push_back(di);
     2006    x2.push_back(di*di);
     2007    x3.push_back(di*di*di);
     2008    z1.push_back(dD);
     2009    x1z1.push_back(dD*di);
     2010    x2z1.push_back(dD*di*di);
     2011    x3z1.push_back(dD*di*di*di);
    20102012    r1.push_back(0.0);
    20112013  }
    … …  
    20152017
    20162018  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    2017     //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an
    2018     //identity matrix (right).
    2019     //the right part is used to calculate the inverse matrix of the left part.
     2019    // xMatrix : horizontal concatenation of
     2020    //           the least-sq. matrix (left) and an
     2021    //           identity matrix (right).
     2022    // the right part is used to calculate the inverse matrix of the left part.
    20202023    double xMatrix[nDOF][2*nDOF];
    20212024    double zMatrix[nDOF];
    … …  
    20292032
    20302033    for (int n = 0; n < nPiece; ++n) {
    2031       for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
     2034      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
     2035
    20322036        if (maskArray[i] == 0) continue;
    2033         double xp1 = x1[i];
    2034         double xp2 = xp1*xp1;
    2035         double xp3 = xp2*xp1;
    2036         double xp4 = xp3*xp1;
    2037         double xp5 = xp4*xp1;
    2038         double xp6 = xp5*xp1;
     2037
    20392038        xMatrix[0][0] += 1.0;
    2040         xMatrix[0][1] += xp1;
    2041         xMatrix[0][2] += xp2;
    2042         xMatrix[0][3] += xp3;
    2043         xMatrix[1][1] += xp2;
    2044         xMatrix[1][2] += xp3;
    2045         xMatrix[1][3] += xp4;
    2046         xMatrix[2][2] += xp4;
    2047         xMatrix[2][3] += xp5;
    2048         xMatrix[3][3] += xp6;
     2039        xMatrix[0][1] += x1[i];
     2040        xMatrix[0][2] += x2[i];
     2041        xMatrix[0][3] += x3[i];
     2042        xMatrix[1][1] += x2[i];
     2043        xMatrix[1][2] += x3[i];
     2044        xMatrix[1][3] += x2[i]*x2[i];
     2045        xMatrix[2][2] += x2[i]*x2[i];
     2046        xMatrix[2][3] += x3[i]*x2[i];
     2047        xMatrix[3][3] += x3[i]*x3[i];
    20492048        zMatrix[0] += z1[i];
    2050         zMatrix[1] += xp1*z1[i];
    2051         zMatrix[2] += xp2*z1[i];
    2052         zMatrix[3] += xp3*z1[i];
     2049        zMatrix[1] += x1z1[i];
     2050        zMatrix[2] += x2z1[i];
     2051        zMatrix[3] += x3z1[i];
     2052
    20532053        for (int j = 0; j < n; ++j) {
    2054           double q = 1.0 - xp1*sectionListR[j];
     2054          double q = 1.0 - x1[i]*invEdge[j];
    20552055          q = q*q*q;
    20562056          xMatrix[0][j+4] += q;
    2057           xMatrix[1][j+4] += q*xp1;
    2058           xMatrix[2][j+4] += q*xp2;
    2059           xMatrix[3][j+4] += q*xp3;
     2057          xMatrix[1][j+4] += q*x1[i];
     2058          xMatrix[2][j+4] += q*x2[i];
     2059          xMatrix[3][j+4] += q*x3[i];
    20602060          for (int k = 0; k < j; ++k) {
    2061             double r = 1.0 - xp1*sectionListR[k];
     2061            double r = 1.0 - x1[i]*invEdge[k];
    20622062            r = r*r*r;
    20632063            xMatrix[k+4][j+4] += r*q;
    … …  
    20662066          zMatrix[j+4] += q*z1[i];
    20672067        }
     2068
    20682069      }
    20692070    }
    … …  
    21252126
    21262127    for (int n = 0; n < nPiece; ++n) {
    2127       for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
    2128         r1[i] = a0 + a1*x1[i] + a2*x1[i]*x1[i] + a3*x1[i]*x1[i]*x1[i];
     2128      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
     2129        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
    21292130      }
    21302131      params.push_back(a0);
    … …  
    21362137
    21372138      double d = y[4+n];
    2138       a0 += d;
    2139       a1 -= 3.0*d*sectionListR[n];
    2140       a2 += 3.0*d*sectionListR[n]*sectionListR[n];
    2141       a3 -= d*sectionListR[n]*sectionListR[n]*sectionListR[n];
     2139      double iE = invEdge[n];
     2140      a0 +=     d;
     2141      a1 -= 3.0*d*iE;
     2142      a2 += 3.0*d*iE*iE;
     2143      a3 -=     d*iE*iE*iE;
    21422144    }
    21432145
    … …  
    21652167      }
    21662168      if (newNData == currentNData) {
    2167         break;                   //no additional flag. finish iteration.
     2169        break; //no more flag to add. iteration stops.
    21682170      } else {
    21692171        currentNData = newNData;
    … …  
    21862188}
    21872189
    2188   void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     2190void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
    21892191  ofstream ofs;
    21902192  String coordInfo;
    … …  
    23242326  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
    23252327  for (int i = 0; i < nChan; ++i) {
    2326     x1.push_back((double)i);
    2327     x2.push_back((double)(i*i));
    2328     x3.push_back((double)(i*i*i));
    2329     z1.push_back((double)data[i]);
    2330     x1z1.push_back(((double)i)*(double)data[i]);
    2331     x2z1.push_back(((double)(i*i))*(double)data[i]);
    2332     x3z1.push_back(((double)(i*i*i))*(double)data[i]);
     2328    double di = (double)i;
     2329    double dD = (double)data[i];
     2330    x1.push_back(di);
     2331    x2.push_back(di*di);
     2332    x3.push_back(di*di*di);
     2333    z1.push_back(dD);
     2334    x1z1.push_back(di*dD);
     2335    x2z1.push_back(di*di*dD);
     2336    x3z1.push_back(di*di*di*dD);
    23332337    r1.push_back(0.0);
    23342338  }
    … …  
    25482552
    25492553/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
    2550 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& pieceEdges, const std::vector<float>& params, const std::vector<bool>& fixed) {
     2554void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params, const std::vector<bool>& fixed) {
    25512555  if (outLogger || outTextFile) {
    25522556    float rms = getRms(chanMask, whichrow);
    … …  
    25552559    if (outLogger) {
    25562560      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
    2557       ols << formatPiecewiseBaselineParams(pieceEdges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2561      ols << formatPiecewiseBaselineParams(edge, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
    25582562    }
    25592563    if (outTextFile) {
    2560       ofs << formatPiecewiseBaselineParams(pieceEdges, params, fixed, rms, masklist, whichrow, true) << flush;
     2564      ofs << formatPiecewiseBaselineParams(edge, params, fixed, rms, masklist, whichrow, true) << flush;
    25612565    }
    25622566  }
    … …  
    26352639}
    26362640
    2637 std::string Scantable::formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
    2638 {
    2639   ostringstream oss;
    2640   oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
    2641 
    2642   if (params.size() > 0) {
    2643     for (uInt i = 0; i < params.size(); ++i) {
    2644       if (i > 0) {
     2641  std::string Scantable::formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose, int start, int count, bool resetparamid) const
     2642{
     2643  int nParam = (int)(params.size());
     2644
     2645  if (nParam < 1) {
     2646    return("  Not fitted");
     2647  } else {
     2648
     2649    ostringstream oss;
     2650    oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
     2651
     2652    if (start < 0) start = 0;
     2653    if (count < 0) count = nParam;
     2654    int end = start + count;
     2655    if (end > nParam) end = nParam;
     2656    int paramidoffset = (resetparamid) ? (-start) : 0;
     2657
     2658    for (int i = start; i < end; ++i) {
     2659      if (i > start) {
    26452660        oss << ",";
    26462661      }
    2647       std::string fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
    2648       float vParam = params[i];
    2649       std::string equals = "= " + (String)((vParam < 0.0) ? "" : " ");
    2650       oss << "  p" << i << fix << equals << setprecision(6) << vParam;
    2651     }
     2662      std::string sFix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
     2663      oss << "  p" << (i+paramidoffset) << sFix << "= " << right << setw(13) << setprecision(6) << params[i];
     2664    }
     2665
     2666    oss << endl;
     2667    oss << formatBaselineParamsFooter(rms, verbose);
     2668
     2669    return String(oss);
     2670  }
     2671
     2672}
     2673
     2674std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
     2675{
     2676  int nOutParam = (int)(params.size());
     2677  int nPiece = (int)(ranges.size()) - 1;
     2678
     2679  if (nOutParam < 1) {
     2680    return("  Not fitted");
     2681  } else if (nPiece < 0) {
     2682    return formatBaselineParams(params, fixed, rms, masklist, whichrow, verbose);
     2683  } else if (nPiece < 1) {
     2684    return("  Bad count of the piece edge info");
     2685  } else if (nOutParam % nPiece != 0) {
     2686    return("  Bad count of the output baseline parameters");
    26522687  } else {
    2653     oss << "  Not fitted";
    2654   }
    2655   oss << endl;
    2656 
    2657   oss << formatBaselineParamsFooter(rms, verbose);
    2658   oss << flush;
    2659 
    2660   return String(oss);
    2661 }
    2662 
    2663 std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
    2664 {
    2665   ostringstream oss;
    2666   oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
    2667 
    2668   if ((params.size() > 0) && (ranges.size() > 0)) {
    2669     if (((ranges.size() % 2) == 0) && ((params.size() % (ranges.size() / 2)) == 0)) {
    2670       uInt nParam = params.size() / (ranges.size() / 2);
    2671       stringstream ss;
    2672       ss << ranges[ranges.size()-1] << flush;
    2673       int wRange = ss.str().size()*2+5;
     2688
     2689    int nParam = nOutParam / nPiece;
     2690
     2691    ostringstream oss;
     2692    oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
     2693
     2694    stringstream ss;
     2695    ss << ranges[nPiece] << flush;
     2696    int wRange = ss.str().size() * 2 + 5;
     2697
     2698    for (int i = 0; i < nPiece; ++i) {
    26742699      ss.str("");
    2675       for (uInt j = 0; j < ranges.size(); j+=2) {
    2676         ss << "  [" << ranges[j] << "," << ranges[j+1] << "]" << flush;
    2677         oss << setw(wRange) << left << ss.str();
    2678         ss.str("");
    2679         for (uInt i = 0; i < nParam; ++i) {
    2680           if (i > 0) {
    2681             oss << ",";
    2682           }
    2683           std::string fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
    2684           float vParam = params[j/2*nParam+i];
    2685           std::string equals = "= " + (String)((vParam < 0.0) ? "" : " ");
    2686           oss << "  p" << i << fix << equals << setprecision(6) << vParam;
    2687         }
    2688         oss << endl;
    2689       }
    2690     } else {
    2691       oss << "  " << endl;
    2692     }
    2693   } else {
    2694     oss << "  Not fitted" << endl;
    2695   }
    2696 
    2697   oss << formatBaselineParamsFooter(rms, verbose);
    2698   oss << flush;
    2699 
    2700   return String(oss);
     2700      ss << "  [" << ranges[i] << "," << (ranges[i+1]-1) << "]";
     2701      oss << left << setw(wRange) << ss.str();
     2702      oss << formatBaselineParams(params, fixed, rms, masklist, whichrow, false, i*nParam, nParam, true);
     2703    }
     2704
     2705    oss << formatBaselineParamsFooter(rms, verbose);
     2706
     2707    return String(oss);
     2708  }
     2709
    27012710}
    27022711

  • trunk/src/Scantable.h

    r2058 r2064  
    544544                                   const std::string& masklist,
    545545                                   int whichrow,
    546                                    bool verbose=false) const;
     546                                   bool verbose=false,
     547                                   int start=-1,
     548                                   int count=-1,
     549                                   bool resetparamid=false) const;
    547550  std::string formatPiecewiseBaselineParams(const std::vector<int>& ranges,
    548551                                            const std::vector<float>& params,
    … …  
    669672  std::vector<float> doCubicSplineFitting(const std::vector<float>& data,
    670673                                          const std::vector<bool>& mask,
    671                                           std::vector<int>& sectionRanges,
     674                                          std::vector<int>& idxEdge,
    672675                                          std::vector<float>& params,
    673676                                          int nPiece=2,
    … …  
    695698  //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder);
    696699  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, Fitter& fitter);
    697   void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<int>& pieceEdges, const std::vector<float>& params, const std::vector<bool>& fixed);
     700  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params, const std::vector<bool>& fixed);
    698701  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<float>& params, const std::vector<bool>& fixed);
    699702
Note: See TracChangeset for help on using the changeset viewer.