Line data    Source code

       1             : //# SpectralFit2.cc: Least Squares fitting of spectral elements: templated part
       2             : //# Copyright (C) 2001,2002,2004
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: casa-feedback@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //# $Id: SpectralFit2.tcc 19935 2007-02-27 05:07:40Z Malte.Marquarding $
      27             : 
      28             : //# Includes
      29             : #include <components/SpectralComponents/SpectralFit.h>
      30             : 
      31             : #include <components/SpectralComponents/CompiledSpectralElement.h>
      32             : #include <components/SpectralComponents/GaussianSpectralElement.h>
      33             : #include <components/SpectralComponents/LogTransformedPolynomialSpectralElement.h>
      34             : #include <components/SpectralComponents/LorentzianSpectralElement.h>
      35             : #include <components/SpectralComponents/PolynomialSpectralElement.h>
      36             : #include <components/SpectralComponents/PowerLogPolynomialSpectralElement.h>
      37             : #include <casacore/scimath/Fitting/NonLinearFitLM.h>
      38             : #include <casacore/scimath/Functionals/CompiledFunction.h>
      39             : #include <casacore/scimath/Functionals/CompoundFunction.h>
      40             : #include <casacore/scimath/Functionals/CompoundParam.h>
      41             : #include <casacore/scimath/Functionals/Gaussian1D.h>
      42             : #include <casacore/scimath/Functionals/Lorentzian1D.h>
      43             : #include <casacore/scimath/Functionals/Polynomial.h>
      44             : #include <casacore/scimath/Functionals/PowerLogarithmicPolynomial.h>
      45             : 
      46             : namespace casa { //# NAMESPACE CASA - BEGIN
      47             : 
      48             : //# Templated member functions
      49             : 
      50             : template <class MT>
      51           0 : bool SpectralFit::fit(const casacore::Vector<MT> &y,
      52             :                       const casacore::Vector<MT> &x,
      53             :                       const casacore::Vector<casacore::Bool> *mask) {
      54           0 :   casacore::Vector<MT> sigma(x.nelements());
      55           0 :   sigma = 1.0;
      56           0 :   return fit(sigma, y, x, mask);
      57           0 : }
      58             : 
      59             : template <class MT>
      60           0 : bool SpectralFit::fit(
      61             :         const casacore::Vector<MT> &sigma, const casacore::Vector<MT> &y,
      62             :         const casacore::Vector<MT> &x, const casacore::Vector<casacore::Bool> *mask
      63             : ) {
      64           0 :         casacore::NonLinearFitLM<MT> fitter;
      65           0 :         iter_p = 0;
      66             :         // The functions to fit
      67           0 :         casacore::CompoundFunction<casacore::AutoDiff<MT> > func;
      68           0 :         casacore::uInt ncomps = slist_p.nelements();
      69           0 :         std::unique_ptr<casacore::Function<casacore::AutoDiff<MT> > > autodiff;
      70           0 :         for (casacore::uInt i=0; i<ncomps; i++) {
      71           0 :                 SpectralElement *elem = slist_p[i];
      72           0 :                 casacore::uInt nparms = elem->getOrder();
      73           0 :                 SpectralElement::Types type = slist_p[i]->getType();
      74           0 :                 switch(type) {
      75           0 :                 case SpectralElement::GAUSSIAN: {
      76           0 :                         autodiff.reset(new casacore::Gaussian1D<casacore::AutoDiff<MT> >());
      77             :                 }
      78           0 :                 break;
      79           0 :                 case SpectralElement::POLYNOMIAL: {
      80           0 :                         PolynomialSpectralElement *x = dynamic_cast<PolynomialSpectralElement *>(elem);
      81           0 :                         autodiff.reset(new casacore::Polynomial<casacore::AutoDiff<MT> >(x->getDegree()));
      82             :                 }
      83           0 :                 break;
      84           0 :                 case SpectralElement::COMPILED:
      85             :                         // Allow fall through; these use the same code
      86             :                 case SpectralElement::GMULTIPLET: {
      87           0 :                         CompiledSpectralElement *x = dynamic_cast<CompiledSpectralElement *>(elem);
      88           0 :                         autodiff.reset(new casacore::CompiledFunction<casacore::AutoDiff<MT> >());
      89           0 :                         dynamic_cast<casacore::CompiledFunction<casacore::AutoDiff<MT> > *>(
      90           0 :                                 autodiff.get()
      91           0 :                         )->setFunction(x->getFunction());
      92             :                 }
      93           0 :                 break;
      94           0 :                 case SpectralElement::LORENTZIAN: {
      95           0 :                         autodiff.reset(new casacore::Lorentzian1D<casacore::AutoDiff<MT> >());
      96             :                 }
      97           0 :                 break;
      98           0 :                 case SpectralElement::POWERLOGPOLY: {
      99           0 :                         casacore::Vector<casacore::Double> parms = elem->get();
     100           0 :                         autodiff.reset(new casacore::PowerLogarithmicPolynomial<casacore::AutoDiff<MT> > (nparms));
     101           0 :                 }
     102           0 :                 break;
     103           0 :                 case SpectralElement::LOGTRANSPOLY: {
     104           0 :                         LogTransformedPolynomialSpectralElement *x = dynamic_cast<
     105             :                                 LogTransformedPolynomialSpectralElement*
     106           0 :                         >(elem);
     107             :                         // treated as a polynomial for fitting purposes. The caller is responsible for passing the ln's of
     108             :                         // the ordinate and obscissa values to the fitter.
     109           0 :                         autodiff.reset(new casacore::Polynomial<casacore::AutoDiff<MT> > (x->getDegree()));
     110             :                 }
     111           0 :                 break;
     112           0 :                 default:
     113           0 :                         throw casacore::AipsError("SpectralFit::fit(): Logic Error: Unhandled SpectralElement type");
     114             :                 }
     115           0 :                 casacore::Vector<casacore::Double> parms = elem->get();
     116           0 :                 casacore::Vector<casacore::Bool> fixed = elem->fixed();
     117           0 :                 for (casacore::uInt j=0; j<nparms; j++) {
     118           0 :                         (*autodiff)[j] = casacore::AutoDiff<MT>(parms[j], nparms, j);
     119           0 :                         if (j == PCFSpectralElement::WIDTH && type == SpectralElement::GAUSSIAN) {
     120           0 :                                 (*autodiff)[j] *= GaussianSpectralElement::SigmaToFWHM;
     121             :                         }
     122           0 :                         autodiff->mask(j) = ! fixed[j];
     123             :                 }
     124           0 :                 func.addFunction(*autodiff);
     125             :         }
     126           0 :         fitter.setFunction(func);
     127             :         // Max. number of iterations
     128           0 :         fitter.setMaxIter(50+ ncomps*10);
     129             :         // Convergence criterium
     130           0 :         fitter.setCriteria(0.001);
     131             :         // Fit
     132           0 :         casacore::Vector<MT> sol;
     133           0 :         casacore::Vector<MT> err;
     134           0 :         sol = fitter.fit(x, y, sigma, mask);
     135           0 :         err = fitter.errors();
     136             :         // Number of iterations
     137           0 :         iter_p = fitter.currentIteration();
     138           0 :         chiSq_p = fitter.getChi2();
     139           0 :         casacore::uInt j = 0;
     140           0 :         casacore::Vector<casacore::Double> tmp, terr;
     141           0 :         for (casacore::uInt i=0; i<ncomps; i++) {
     142           0 :                 SpectralElement *element = slist_p[i];
     143           0 :                 casacore::uInt nparms = element->getOrder();
     144           0 :                 tmp.resize(nparms);
     145           0 :                 terr.resize(nparms);
     146           0 :                 SpectralElement::Types type = element->getType();
     147           0 :                 for (casacore::uInt k=0; k<nparms; k++) {
     148           0 :                         casacore::Bool convertGaussWidth = k==PCFSpectralElement::WIDTH && type == SpectralElement::GAUSSIAN;
     149           0 :                         tmp(k) = convertGaussWidth
     150           0 :                                 ? sol(j) / GaussianSpectralElement::SigmaToFWHM
     151           0 :                                 : sol(j);
     152           0 :                         terr(k) = convertGaussWidth
     153           0 :                                 ? err(j) / GaussianSpectralElement::SigmaToFWHM
     154           0 :                                 : err(j);
     155           0 :                         j++;
     156             :                 };
     157           0 :                 element->set(tmp);
     158           0 :                 element->setError(terr);
     159             :         }
     160           0 :         return fitter.converged();
     161           0 : }
     162             : 
     163             : 
     164             : } //# NAMESPACE CASA - END
     165             :