Changeset 2344


Ignore:
Timestamp:
11/08/11 17:02:32 (12 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): xmlcasa

Description: The cubic spline fitting was returning zero values both for the fitted results and the residual within masked regions touching the edge of the given spetrum. For such masked regions, the value of the fitted result at the nearest unmasked channel is now adopted, and residual is correctly computed also.


File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/Scantable.cpp

    r2332 r2344  
    1111//
    1212#include <map>
     13#include <mcheck.h>
    1314
    1415#include <atnf/PKSIO/SrcType.h>
    … …  
    9091  initFactories();
    9192  setupMainTable();
    92   freqTable_ = STFrequencies(*this);
    93   table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
     93  freqTable_ = new STFrequencies(*this);
     94  table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_->table());
    9495  weatherTable_ = STWeather(*this);
    9596  table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
    … …  
    194195void Scantable::copySubtables(const Scantable& other) {
    195196  Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
    196   TableCopy::copyRows(t, other.freqTable_.table());
     197  TableCopy::copyRows(t, other.freqTable_->table());
    197198  t = table_.rwKeywordSet().asTable("FOCUS");
    198199  TableCopy::copyRows(t, other.focusTable_.table());
    … …  
    211212void Scantable::attachSubtables()
    212213{
    213   freqTable_ = STFrequencies(table_);
     214  freqTable_ = new STFrequencies(table_);
    214215  focusTable_ = STFocus(table_);
    215216  weatherTable_ = STWeather(table_);
    … …  
    222223Scantable::~Scantable()
    223224{
    224   //cout << "~Scantable() " << this << endl;
     225  delete freqTable_;
    225226}
    226227
    … …  
    15091510  //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
    15101511  Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
    1511   return freqTable_.getSpectralCoordinate(md, mp, me, rf,
     1512  return freqTable_->getSpectralCoordinate(md, mp, me, rf,
    15121513                                          mfreqidCol_(whichrow));
    15131514}
    … …  
    15181519  std::vector<double> stlout;
    15191520  int nchan = specCol_(whichrow).nelements();
    1520   String us = freqTable_.getUnitString();
     1521  String us = freqTable_->getUnitString();
    15211522  if ( us == "" || us == "pixel" || us == "channel" ) {
    15221523    for (int i=0; i<nchan; ++i) {
    … …  
    15851586  const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
    15861587  SpectralCoordinate spc =
    1587     freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
     1588    freqTable_->getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
    15881589
    15891590  String s = "Channel";
    1590   Unit u = Unit(freqTable_.getUnitString());
     1591  Unit u = Unit(freqTable_->getUnitString());
    15911592  if (u == Unit("km/s")) {
    15921593    s = CoordinateUtil::axisLabel(spc, 0, True,True,  True);
    … …  
    18491850  Double refval ;
    18501851  Double increment ;
    1851   int freqnrow = freqTable_.table().nrow() ;
     1852  int freqnrow = freqTable_->table().nrow() ;
    18521853  Vector<uInt> oldId( freqnrow ) ;
    18531854  Vector<uInt> newId( freqnrow ) ;
    18541855  for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
    1855     freqTable_.getEntry( refpix, refval, increment, irow ) ;
     1856    freqTable_->getEntry( refpix, refval, increment, irow ) ;
    18561857    /***
    18571858     * need to shift refpix to nmin
    … …  
    18601861    refval = refval - ( refpix - nmin ) * increment ;
    18611862    refpix = 0 ;
    1862     freqTable_.setEntry( refpix, refval, increment, irow ) ;
     1863    freqTable_->setEntry( refpix, refval, increment, irow ) ;
    18631864  }
    18641865
    … …  
    23052306    //fitter.setIterClipping(thresClip, nIterClip);
    23062307
    2307     int nRow = nrow();
    2308     std::vector<bool> chanMask;
    23092308    bool showProgress;
    23102309    int minNRow;
    23112310    parseProgressInfo(progressInfo, showProgress, minNRow);
    23122311
     2312    int nRow = nrow();
     2313    std::vector<bool> chanMask;
     2314
     2315    //--------------------------------
    23132316    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    23142317      chanMask = getCompositeChanMask(whichrow, mask);
    23152318      //fitBaseline(chanMask, whichrow, fitter);
    23162319      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2317       std::vector<int> pieceEdges;
    2318       std::vector<float> params;
     2320      std::vector<int> pieceEdges(nPiece+1);
     2321      std::vector<float> params(nPiece*4);
    23192322      int nClipped = 0;
    23202323      std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
    … …  
    23252328      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    23262329    }
    2327 
     2330    //--------------------------------
     2331   
    23282332    if (outTextFile) ofs.close();
    23292333
    … …  
    23942398      //fitBaseline(chanMask, whichrow, fitter);
    23952399      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2396       std::vector<int> pieceEdges;
    2397       std::vector<float> params;
     2400      std::vector<int> pieceEdges(nPiece+1);
     2401      std::vector<float> params(nPiece*4);
    23982402      int nClipped = 0;
    23992403      std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
    … …  
    24222426
    24232427  int nChan = data.size();
    2424   std::vector<int> maskArray;
    2425   std::vector<int> x;
     2428  std::vector<int> maskArray(nChan);
     2429  std::vector<int> x(nChan);
     2430  int j = 0;
    24262431  for (int i = 0; i < nChan; ++i) {
    2427     maskArray.push_back(mask[i] ? 1 : 0);
     2432    maskArray[i] = mask[i] ? 1 : 0;
    24282433    if (mask[i]) {
    2429       x.push_back(i);
    2430     }
    2431   }
    2432 
    2433   int initNData = x.size();
     2434      x[j] = i;
     2435      j++;
     2436    }
     2437  }
     2438  int initNData = j;
     2439
    24342440  if (initNData < nPiece) {
    24352441    throw(AipsError("too few non-flagged channels"));
    … …  
    24372443
    24382444  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
    2439   std::vector<double> invEdge;
    2440   idxEdge.clear();
    2441   idxEdge.push_back(x[0]);
     2445  std::vector<double> invEdge(nPiece-1);
     2446  idxEdge[0] = x[0];
    24422447  for (int i = 1; i < nPiece; ++i) {
    24432448    int valX = x[nElement*i];
    2444     idxEdge.push_back(valX);
    2445     invEdge.push_back(1.0/(double)valX);
    2446   }
    2447   idxEdge.push_back(x[x.size()-1]+1);
     2449    idxEdge[i] = valX;
     2450    invEdge[i-1] = 1.0/(double)valX;
     2451  }
     2452  idxEdge[nPiece] = x[initNData-1]+1;
    24482453
    24492454  int nData = initNData;
    24502455  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
    24512456
    2452   std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1, residual;
     2457  std::vector<double> x1(nChan), x2(nChan), x3(nChan);
     2458  std::vector<double> z1(nChan), x1z1(nChan), x2z1(nChan), x3z1(nChan);
     2459  std::vector<double> r1(nChan), residual(nChan);
    24532460  for (int i = 0; i < nChan; ++i) {
    24542461    double di = (double)i;
    24552462    double dD = (double)data[i];
    2456     x1.push_back(di);
    2457     x2.push_back(di*di);
    2458     x3.push_back(di*di*di);
    2459     z1.push_back(dD);
    2460     x1z1.push_back(dD*di);
    2461     x2z1.push_back(dD*di*di);
    2462     x3z1.push_back(dD*di*di*di);
    2463     r1.push_back(0.0);
    2464     residual.push_back(0.0);
     2463    x1[i]   = di;
     2464    x2[i]   = di*di;
     2465    x3[i]   = di*di*di;
     2466    z1[i]   = dD;
     2467    x1z1[i] = dD*di;
     2468    x2z1[i] = dD*di*di;
     2469    x3z1[i] = dD*di*di*di;
     2470    r1[i]   = 0.0;
     2471    residual[i] = 0.0;
    24652472  }
    24662473
    … …  
    25402547    }
    25412548
    2542     std::vector<double> invDiag;
     2549    std::vector<double> invDiag(nDOF);
    25432550    for (int i = 0; i < nDOF; ++i) {
    2544       invDiag.push_back(1.0/xMatrix[i][i]);
     2551      invDiag[i] = 1.0/xMatrix[i][i];
    25452552      for (int j = 0; j < nDOF; ++j) {
    25462553        xMatrix[i][j] *= invDiag[i];
    … …  
    25742581    //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
    25752582    //cubic terms for the other pieces (in case nPiece>1).
    2576     std::vector<double> y;
    2577     y.clear();
     2583    std::vector<double> y(nDOF);
    25782584    for (int i = 0; i < nDOF; ++i) {
    2579       y.push_back(0.0);
     2585      y[i] = 0.0;
    25802586      for (int j = 0; j < nDOF; ++j) {
    25812587        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
    … …  
    25872593    double a2 = y[2];
    25882594    double a3 = y[3];
    2589     params.clear();
    2590 
     2595
     2596    int j = 0;
    25912597    for (int n = 0; n < nPiece; ++n) {
    25922598      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    25932599        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
    2594         residual[i] = z1[i] - r1[i];
    2595       }
    2596       params.push_back(a0);
    2597       params.push_back(a1);
    2598       params.push_back(a2);
    2599       params.push_back(a3);
     2600      }
     2601      params[j]   = a0;
     2602      params[j+1] = a1;
     2603      params[j+2] = a2;
     2604      params[j+3] = a3;
     2605      j += 4;
    26002606
    26012607      if (n == nPiece-1) break;
    … …  
    26072613      a2 += 3.0*d*iE*iE;
    26082614      a3 -=     d*iE*iE*iE;
     2615    }
     2616
     2617    //subtract constant value for masked regions at the edge of spectrum
     2618    if (idxEdge[0] > 0) {
     2619      int n = idxEdge[0];
     2620      for (int i = 0; i < idxEdge[0]; ++i) {
     2621        //--cubic extrapolate--
     2622        //r1[i] = params[0] + params[1]*x1[i] + params[2]*x2[i] + params[3]*x3[i];
     2623        //--linear extrapolate--
     2624        //r1[i] = (r1[n+1] - r1[n])/(x1[n+1] - x1[n])*(x1[i] - x1[n]) + r1[n];
     2625        //--constant--
     2626        r1[i] = r1[n];
     2627      }
     2628    }
     2629    if (idxEdge[nPiece] < nChan) {
     2630      int n = idxEdge[nPiece]-1;
     2631      for (int i = idxEdge[nPiece]; i < nChan; ++i) {
     2632        //--cubic extrapolate--
     2633        //int m = 4*(nPiece-1);
     2634        //r1[i] = params[m] + params[m+1]*x1[i] + params[m+2]*x2[i] + params[m+3]*x3[i];
     2635        //--linear extrapolate--
     2636        //r1[i] = (r1[n-1] - r1[n])/(x1[n-1] - x1[n])*(x1[i] - x1[n]) + r1[n];
     2637        //--constant--
     2638        r1[i] = r1[n];
     2639      }
     2640    }
     2641
     2642    for (int i = 0; i < nChan; ++i) {
     2643      residual[i] = z1[i] - r1[i];
    26092644    }
    26102645
    … …  
    26382673  nClipped = initNData - nData;
    26392674
    2640   std::vector<float> result;
     2675  std::vector<float> result(nChan);
    26412676  if (getResidual) {
    26422677    for (int i = 0; i < nChan; ++i) {
    2643       result.push_back((float)residual[i]);
     2678      result[i] = (float)residual[i];
    26442679    }
    26452680  } else {
    26462681    for (int i = 0; i < nChan; ++i) {
    2647       result.push_back((float)r1[i]);
     2682      result[i] = (float)r1[i];
    26482683    }
    26492684  }
    … …  
    26522687}
    26532688
    2654   void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
     2689void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
    26552690{
    26562691  nWaves.clear();
Note: See TracChangeset for help on using the changeset viewer.