Changeset 2713 for trunk

Timestamp:

12/28/12 15:52:45 (12 years ago)

Author:

WataruKawasaki

Message:

New Development: Yes

JIRA Issue: Yes CAS-3618

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sd

Description: new method for scantable to calculate AIC, AICc, BIC and GCV.

Location:

trunk

Files:

• python/scantable.py (modified) (1 diff)
• src/Scantable.cpp (modified) (3 diffs)
• src/Scantable.h (modified) (2 diffs)
• src/ScantableWrapper.h (modified) (1 diff)
• src/python_Scantable.cpp (modified) (1 diff)

trunk/python/scantable.py

@@ -2539 +2539 @@
             raise RuntimeError(msg)
 
+    @asaplog_post_dec
+    def calc_aic(self, value=None, blfunc=None, order=None, mask=None,
+                 whichrow=None, uselinefinder=None, edge=None,
+                 threshold=None, chan_avg_limit=None):
+        """\
+        Calculates and returns model selection criteria for a specified
+        baseline model and a given spectrum data. 
+        Available values include Akaike Information Criterion (AIC), the
+        corrected Akaike Information Criterion (AICc) by Sugiura(1978),
+        Bayesian Information Criterion (BIC) and the Generalised Cross
+        Validation (GCV). 
+
+        Parameters:
+            value:         name of model selection criteria to calculate.
+                           available ones include 'aic', 'aicc', 'bic' and
+                           'gcv'. default is 'aicc'.
+            blfunc:        baseline function name. available ones include 
+                           'chebyshev', 'cspline' and 'sinusoid'.
+                           default is 'chebyshev'.
+            order:         parameter for basline function. actually stands for
+                           order of polynomial (order) for 'chebyshev',
+                           number of spline pieces (npiece) for 'cspline' and
+                           maximum wave number for 'sinusoid', respectively.
+                           default is 5 (which is also the default order value
+                           for [auto_]chebyshev_baseline()).
+            mask:          an optional mask. default is [].
+            whichrow:      row number. default is 0 (the first row)
+            uselinefinder: use sd.linefinder() to flag out line regions
+                           default is True.
+            edge:           an optional number of channel to drop at
+                            the edge of spectrum. If only one value is
+                            specified, the same number will be dropped
+                            from both sides of the spectrum. Default
+                            is to keep all channels. Nested tuples
+                            represent individual edge selection for
+                            different IFs (a number of spectral channels
+                            can be different)
+                            default is (0, 0).
+            threshold:      the threshold used by line finder. It is
+                            better to keep it large as only strong lines
+                            affect the baseline solution.
+                            default is 3.
+            chan_avg_limit: a maximum number of consequtive spectral
+                            channels to average during the search of
+                            weak and broad lines. The default is no
+                            averaging (and no search for weak lines).
+                            If such lines can affect the fitted baseline
+                            (e.g. a high order polynomial is fitted),
+                            increase this parameter (usually values up
+                            to 8 are reasonable). Most users of this
+                            method should find the default value sufficient.
+                            default is 1.
+
+        Example:
+            aic = scan.calc_aic(blfunc='chebyshev', order=5, whichrow=0)
+        """
+
+        try:
+            varlist = vars()
+
+            if value          is None: value          = 'aicc'
+            if blfunc         is None: blfunc         = 'chebyshev'
+            if order          is None: order          = 5
+            if mask           is None: mask           = []
+            if whichrow       is None: whichrow       = 0
+            if uselinefinder  is None: uselinefinder  = True
+            if edge           is None: edge           = (0, 0)
+            if threshold      is None: threshold      = 3
+            if chan_avg_limit is None: chan_avg_limit = 1
+
+            return self._calc_aic(value, blfunc, order, mask,
+                                  whichrow, uselinefinder, edge,
+                                  threshold, chan_avg_limit)
+            
+        except RuntimeError, e:
+            raise_fitting_failure_exception(e)
 
     @asaplog_post_dec

trunk/src/Scantable.cpp

@@ -2656 +2656 @@
 }
 
+double Scantable::calculateModelSelectionCriteria(const std::string& valname, const std::string& blfunc, int order, const std::vector<bool>& inMask, int whichrow, bool useLineFinder, const std::vector<int>& edge, float threshold, int chanAvgLimit)
+{
+  if (useLineFinder) {
+    int minEdgeSize = getIFNos().size()*2;
+
+
+    int edgeSize = edge.size();
+    std::vector<int> currentEdge;
+    if (edgeSize >= 2) {
+      int idx = 0;
+      if (edgeSize > 2) {
+        if (edgeSize < minEdgeSize) {
+          throw(AipsError("Length of edge element info is less than that of IFs"));
+        }
+        idx = 2 * getIF(whichrow);
+      }
+      currentEdge.push_back(edge[idx]);
+      currentEdge.push_back(edge[idx+1]);
+    } else {
+      throw(AipsError("Wrong length of edge element"));
+    }
+
+    STLineFinder lineFinder = STLineFinder();
+    std::vector<float> sp = getSpectrum(whichrow);
+    lineFinder.setData(sp);
+    lineFinder.setOptions(threshold, 3, chanAvgLimit);
+    lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
+    std::vector<bool> chanMask = lineFinder.getMask();
+    
+    return doCalculateModelSelectionCriteria(valname, sp, chanMask, blfunc, order);
+
+  } else {
+    return doCalculateModelSelectionCriteria(valname, getSpectrum(whichrow), getCompositeChanMask(whichrow, inMask), blfunc, order);
+  }
+
+}
+
+double Scantable::doCalculateModelSelectionCriteria(const std::string& valname, const std::vector<float>& spec, const std::vector<bool>& mask, const std::string& blfunc, int order)
+{
+  int nparam;
+  std::vector<float> params;
+  int nClipped = 0;
+  float thresClip = 0.0;
+  int nIterClip = 0;
+  std::vector<float> res;
+  if (blfunc == "chebyshev") {
+    nparam = order + 1;
+    res = doChebyshevFitting(spec, mask, order, params, nClipped, thresClip, nIterClip, true);
+  } else if (blfunc == "sinusoid") {
+    std::vector<int> nWaves;
+    nWaves.clear();
+    for (int i = 0; i <= order; ++i) {
+      nWaves.push_back(i);
+    }
+    nparam = 2*order + 1;  // order = nwave
+    res = doSinusoidFitting(spec, mask, nWaves, params, nClipped, thresClip, nIterClip, true);
+  } else if (blfunc == "cspline") {
+    std::vector<int> pieceEdges(order+1);  //order = npiece
+    nparam = order + 3;
+    params.resize(4*order);
+    res = doCubicSplineFitting(spec, mask, order, pieceEdges, params, nClipped, thresClip, nIterClip, true);
+  } else {
+    throw(AipsError("blfunc must be chebyshev, cspline or sinusoid."));
+  }
+
+  double msq = 0.0;
+  int nusedchan = 0;
+  int nChan = res.size();
+  for (int i = 0; i < nChan; ++i) {
+    if (mask[i]) {
+      msq += (double)res[i]*(double)res[i];
+      nusedchan++;
+    }
+  }
+  if (nusedchan == 0) {
+    throw(AipsError("all channels masked."));
+  }
+  msq /= (double)nusedchan;
+
+  nparam++;  //add 1 for sigma of Gaussian distribution
+  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
+
+  if (valname.find("aic") == 0) {
+    // Original Akaike Information Criterion (AIC)
+    double aic = nusedchan * (log(2.0 * PI * msq) + 1.0) + 2.0 * nparam;
+
+    // Corrected AIC by Sugiura(1978)
+    if (valname == "aicc") {
+      if (nusedchan - nparam - 1 <= 0) {
+        throw(AipsError("channel size is too small to calculate AICc."));
+      }
+      aic += 2.0*nparam*(nparam + 1)/(double)(nusedchan - nparam - 1);
+    }
+
+    return aic;
+
+  } else if (valname == "bic") {
+    // Bayesian Information Criterion (BIC)
+    double bic = nusedchan * log(msq) + nparam * log((double)nusedchan);
+    return bic;
+
+  } else if (valname == "gcv") {
+    // Generalised Cross Validation
+    double x = 1.0 - (double)nparam / (double)nusedchan;
+    double gcv = msq / (x * x);
+    return gcv;
+
+  } else {
+    throw(AipsError("valname must be aic, aicc, bic or gcv."));
+  }
+}
+
 void Scantable::autoChebyshevBaseline(const std::vector<bool>& mask, int order, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
 {
@@ -2751 +2863 @@
   if ((x < -1.0)||(x > 1.0)) {
     throw(AipsError("out of definition range (-1 <= x <= 1)."));
+  } else if (x == 1.0) {
+    return 1.0;
+  } else if (x == 0.0) {
+    double res;
+    if (n%2 == 0) {
+      if (n%4 == 0) {
+        res = 1.0;
+      } else {
+        res = -1.0;
+      }
+    } else {
+      res = 0.0;
+    }
+    return res;
+  } else if (x == -1.0) {
+    double res = (n%2 == 0 ? 1.0 : -1.0);
+    return res;
   } else if (n < 0) {
     throw(AipsError("the order must be zero or positive."));
@@ -2766 +2895 @@
         }
       }
-      res += (m%2 == 0 ? 1.0 : -1.0)*(double)n/(double)(n-m)*pow(2.0*x, (double)(n-2*m-1))/2.0*c;
+      res += (m%2 == 0 ? 1.0 : -1.0)*(double)n/(double)(n-m)*pow(2.0*x, (double)(n-2*m))/2.0*c;
     }
     return res;

trunk/src/Scantable.h

@@ -619 +619 @@
   std::vector<uint> getMoleculeIdColumnData() const;
   void setMoleculeIdColumnData(const std::vector<uint>& molids);
+  double calculateModelSelectionCriteria(const std::string& valname,
+                                         const std::string& blfunc,
+                                         int order,
+                                         const std::vector<bool>& inMask,
+                                         int whichrow,
+                                         bool useLineFinder,
+                                         const std::vector<int>& edge,
+                                         float threshold,
+                                         int chanAvgLimit);
 
 
@@ -804 +813 @@
   void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);
 
+  double doCalculateModelSelectionCriteria(const std::string& valname, 
+                                           const std::vector<float>& spec, 
+                                           const std::vector<bool>& mask, 
+                                           const std::string& blfunc, 
+                                           int order);
+
 };
 

trunk/src/ScantableWrapper.h

@@ -321 +321 @@
   void setMoleculeIdColumnData(const std::vector<uint>& molids)
   { table_->setMoleculeIdColumnData(molids); }
-
+  double calculateModelSelectionCriteria(const std::string& valname, const std::string& blfunc, int order, const std::vector<bool>& inMask, int whichrow, bool useLineFinder, const std::vector<int>& edge, float threshold, int chanAvgLimit)
+  { return table_->calculateModelSelectionCriteria(valname, blfunc, order, inMask, whichrow, useLineFinder, edge, threshold, chanAvgLimit); }
 
 private:

trunk/src/python_Scantable.cpp

@@ -163 +163 @@
     .def("_getmolidcol_list", &ScantableWrapper::getMoleculeIdColumnData)
     .def("_setmolidcol_list", &ScantableWrapper::setMoleculeIdColumnData)
+    .def("_calc_aic", &ScantableWrapper::calculateModelSelectionCriteria)
   ;
 };