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: aips2-request@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             : }
      58             : 
      59             : template <class MT>
      60      545637 : 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     1091274 :         casacore::NonLinearFitLM<MT> fitter;
      65      545637 :         iter_p = 0;
      66             :         // The functions to fit
      67     1091274 :         casacore::CompoundFunction<casacore::AutoDiff<MT> > func;
      68      545637 :         casacore::uInt ncomps = slist_p.nelements();
      69      545637 :         std::unique_ptr<casacore::Function<casacore::AutoDiff<MT> > > autodiff;
      70     1107561 :         for (casacore::uInt i=0; i<ncomps; i++) {
      71      561924 :                 SpectralElement *elem = slist_p[i];
      72      561924 :                 casacore::uInt nparms = elem->getOrder();
      73      561924 :                 SpectralElement::Types type = slist_p[i]->getType();
      74      561924 :                 switch(type) {
      75       32428 :                 case SpectralElement::GAUSSIAN: {
      76       32428 :                         autodiff.reset(new casacore::Gaussian1D<casacore::AutoDiff<MT> >());
      77             :                 }
      78       32428 :                 break;
      79      529222 :                 case SpectralElement::POLYNOMIAL: {
      80      529222 :                         PolynomialSpectralElement *x = dynamic_cast<PolynomialSpectralElement *>(elem);
      81      529222 :                         autodiff.reset(new casacore::Polynomial<casacore::AutoDiff<MT> >(x->getDegree()));
      82             :                 }
      83      529222 :                 break;
      84          26 :                 case SpectralElement::COMPILED:
      85             :                         // Allow fall through; these use the same code
      86             :                 case SpectralElement::GMULTIPLET: {
      87          26 :                         CompiledSpectralElement *x = dynamic_cast<CompiledSpectralElement *>(elem);
      88          26 :                         autodiff.reset(new casacore::CompiledFunction<casacore::AutoDiff<MT> >());
      89          26 :                         dynamic_cast<casacore::CompiledFunction<casacore::AutoDiff<MT> > *>(
      90             :                                 autodiff.get()
      91             :                         )->setFunction(x->getFunction());
      92             :                 }
      93          26 :                 break;
      94          50 :                 case SpectralElement::LORENTZIAN: {
      95          50 :                         autodiff.reset(new casacore::Lorentzian1D<casacore::AutoDiff<MT> >());
      96             :                 }
      97          50 :                 break;
      98         155 :                 case SpectralElement::POWERLOGPOLY: {
      99         310 :                         casacore::Vector<casacore::Double> parms = elem->get();
     100         155 :                         autodiff.reset(new casacore::PowerLogarithmicPolynomial<casacore::AutoDiff<MT> > (nparms));
     101             :                 }
     102         155 :                 break;
     103          43 :                 case SpectralElement::LOGTRANSPOLY: {
     104          43 :                         LogTransformedPolynomialSpectralElement *x = dynamic_cast<
     105             :                                 LogTransformedPolynomialSpectralElement*
     106             :                         >(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          43 :                         autodiff.reset(new casacore::Polynomial<casacore::AutoDiff<MT> > (x->getDegree()));
     110             :                 }
     111          43 :                 break;
     112           0 :                 default:
     113           0 :                         throw casacore::AipsError("SpectralFit::fit(): Logic Error: Unhandled SpectralElement type");
     114             :                 }
     115     1123848 :                 casacore::Vector<casacore::Double> parms = elem->get();
     116     1123848 :                 casacore::Vector<casacore::Bool> fixed = elem->fixed();
     117     1452105 :                 for (casacore::uInt j=0; j<nparms; j++) {
     118      890181 :                         (*autodiff)[j] = casacore::AutoDiff<MT>(parms[j], nparms, j);
     119      890181 :                         if (j == PCFSpectralElement::WIDTH && type == SpectralElement::GAUSSIAN) {
     120       32428 :                                 (*autodiff)[j] *= GaussianSpectralElement::SigmaToFWHM;
     121             :                         }
     122      890181 :                         autodiff->mask(j) = ! fixed[j];
     123             :                 }
     124      561924 :                 func.addFunction(*autodiff);
     125             :         }
     126      545637 :         fitter.setFunction(func);
     127             :         // Max. number of iterations
     128      545637 :         fitter.setMaxIter(50+ ncomps*10);
     129             :         // Convergence criterium
     130      545637 :         fitter.setCriteria(0.001);
     131             :         // Fit
     132     1091274 :         casacore::Vector<MT> sol;
     133     1091274 :         casacore::Vector<MT> err;
     134      545637 :         sol = fitter.fit(x, y, sigma, mask);
     135      545546 :         err = fitter.errors();
     136             :         // Number of iterations
     137      545546 :         iter_p = fitter.currentIteration();
     138      545546 :         chiSq_p = fitter.getChi2();
     139      545546 :         casacore::uInt j = 0;
     140     1091095 :         casacore::Vector<casacore::Double> tmp, terr;
     141     1107283 :         for (casacore::uInt i=0; i<ncomps; i++) {
     142      561740 :                 SpectralElement *element = slist_p[i];
     143      561740 :                 casacore::uInt nparms = element->getOrder();
     144      561740 :                 tmp.resize(nparms);
     145      561740 :                 terr.resize(nparms);
     146      561740 :                 SpectralElement::Types type = element->getType();
     147     1451369 :                 for (casacore::uInt k=0; k<nparms; k++) {
     148      889629 :                         casacore::Bool convertGaussWidth = k==PCFSpectralElement::WIDTH && type == SpectralElement::GAUSSIAN;
     149      889629 :                         tmp(k) = convertGaussWidth
     150      889629 :                                 ? sol(j) / GaussianSpectralElement::SigmaToFWHM
     151      857385 :                                 : sol(j);
     152      889629 :                         terr(k) = convertGaussWidth
     153      889629 :                                 ? err(j) / GaussianSpectralElement::SigmaToFWHM
     154      857385 :                                 : err(j);
     155      889629 :                         j++;
     156             :                 };
     157      561740 :                 element->set(tmp);
     158      561737 :                 element->setError(terr);
     159             :         }
     160     1091086 :         return fitter.converged();
     161             : }
     162             : 
     163             : 
     164             : } //# NAMESPACE CASA - END
     165             :