source: trunk/external-alma/components/SpectralComponents/SpectralFit2.tcc@ 3016

Last change on this file since 3016 was 2980, checked in by Malte Marquarding, 10 years ago

Add a copy of casacore/components/SpectralComponents to external-alma directory to prepare for its removal from casacore-trunk. DOn;t activate in SConscript yet.

File size: 5.8 KB
Line 
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 21465 2014-06-19 05:56:56Z gervandiepen $
27
28//# Includes
29#include <components/SpectralComponents/SpectralFit.h>
30
31#include <casa/Utilities/PtrHolder.h>
32#include <components/SpectralComponents/CompiledSpectralElement.h>
33#include <components/SpectralComponents/GaussianSpectralElement.h>
34#include <components/SpectralComponents/LogTransformedPolynomialSpectralElement.h>
35#include <components/SpectralComponents/LorentzianSpectralElement.h>
36#include <components/SpectralComponents/PolynomialSpectralElement.h>
37#include <components/SpectralComponents/PowerLogPolynomialSpectralElement.h>
38#include <scimath/Fitting/NonLinearFitLM.h>
39#include <scimath/Functionals/CompiledFunction.h>
40#include <scimath/Functionals/CompoundFunction.h>
41#include <scimath/Functionals/CompoundParam.h>
42#include <scimath/Functionals/Gaussian1D.h>
43#include <scimath/Functionals/Lorentzian1D.h>
44#include <scimath/Functionals/Polynomial.h>
45#include <scimath/Functionals/PowerLogarithmicPolynomial.h>
46
47namespace casa { //# NAMESPACE CASA - BEGIN
48
49//# Templated member functions
50
51template <class MT>
52Bool SpectralFit::fit(const Vector<MT> &y,
53 const Vector<MT> &x,
54 const Vector<Bool> *mask) {
55 Vector<MT> sigma(x.nelements());
56 sigma = 1.0;
57 return fit(sigma, y, x, mask);
58}
59
60template <class MT>
61Bool SpectralFit::fit(
62 const Vector<MT> &sigma, const Vector<MT> &y,
63 const Vector<MT> &x, const Vector<Bool> *mask
64) {
65 NonLinearFitLM<MT> fitter;
66 iter_p = 0;
67 // The functions to fit
68 CompoundFunction<AutoDiff<MT> > func;
69 uInt ncomps = slist_p.nelements();
70 PtrHolder<Function<AutoDiff<MT> > > autodiff;
71 for (uInt i=0; i<ncomps; i++) {
72 SpectralElement *elem = slist_p[i];
73 uInt nparms = elem->getOrder();
74 SpectralElement::Types type = slist_p[i]->getType();
75 switch(type) {
76 case SpectralElement::GAUSSIAN: {
77 autodiff.set(new Gaussian1D<AutoDiff<MT> >());
78 }
79 break;
80 case SpectralElement::POLYNOMIAL: {
81 PolynomialSpectralElement *x = dynamic_cast<PolynomialSpectralElement *>(elem);
82 autodiff.set(new Polynomial<AutoDiff<MT> >(x->getDegree()));
83 }
84 break;
85 case SpectralElement::COMPILED:
86 // Allow fall through; these use the same code
87 case SpectralElement::GMULTIPLET: {
88 CompiledSpectralElement *x = dynamic_cast<CompiledSpectralElement *>(elem);
89 autodiff.set(new CompiledFunction<AutoDiff<MT> >());
90 dynamic_cast<CompiledFunction<AutoDiff<MT> > *>(
91 autodiff.ptr()
92 )->setFunction(x->getFunction());
93 }
94 break;
95 case SpectralElement::LORENTZIAN: {
96 autodiff.set(new Lorentzian1D<AutoDiff<MT> >());
97 }
98 break;
99 case SpectralElement::POWERLOGPOLY: {
100 Vector<Double> parms = elem->get();
101 autodiff.set(new PowerLogarithmicPolynomial<AutoDiff<MT> > (nparms));
102 }
103 break;
104 case SpectralElement::LOGTRANSPOLY: {
105 LogTransformedPolynomialSpectralElement *x = dynamic_cast<
106 LogTransformedPolynomialSpectralElement*
107 >(elem);
108 // treated as a polynomial for fitting purposes. The caller is responsible for passing the ln's of
109 // the ordinate and obscissa values to the fitter.
110 autodiff.set(new Polynomial<AutoDiff<MT> > (x->getDegree()));
111 }
112 break;
113 default:
114 throw AipsError("SpectralFit::fit(): Logic Error: Unhandled SpectralElement type");
115 }
116 Vector<Double> parms = elem->get();
117 Vector<Bool> fixed = elem->fixed();
118 for (uInt j=0; j<nparms; j++) {
119 (*autodiff)[j] = AutoDiff<MT>(parms[j], nparms, j);
120 if (j == PCFSpectralElement::WIDTH && type == SpectralElement::GAUSSIAN) {
121 (*autodiff)[j] *= GaussianSpectralElement::SigmaToFWHM;
122 }
123 autodiff->mask(j) = ! fixed[j];
124 }
125 func.addFunction(*autodiff);
126 }
127 fitter.setFunction(func);
128 // Max. number of iterations
129 fitter.setMaxIter(50+ ncomps*10);
130 // Convergence criterium
131 fitter.setCriteria(0.001);
132 // Fit
133 Vector<MT> sol;
134 Vector<MT> err;
135 sol = fitter.fit(x, y, sigma, mask);
136 err = fitter.errors();
137 // Number of iterations
138 iter_p = fitter.currentIteration();
139 chiSq_p = fitter.getChi2();
140 uInt j = 0;
141 Vector<Double> tmp, terr;
142 for (uInt i=0; i<ncomps; i++) {
143 SpectralElement *element = slist_p[i];
144 uInt nparms = element->getOrder();
145 tmp.resize(nparms);
146 terr.resize(nparms);
147 SpectralElement::Types type = element->getType();
148 for (uInt k=0; k<nparms; k++) {
149 Bool convertGaussWidth = k==PCFSpectralElement::WIDTH && type == SpectralElement::GAUSSIAN;
150 tmp(k) = convertGaussWidth
151 ? sol(j) / GaussianSpectralElement::SigmaToFWHM
152 : sol(j);
153 terr(k) = convertGaussWidth
154 ? err(j) / GaussianSpectralElement::SigmaToFWHM
155 : err(j);
156 j++;
157 };
158 element->set(tmp);
159 element->setError(terr);
160 }
161 return fitter.converged();
162}
163
164
165} //# NAMESPACE CASA - END
166
Note: See TracBrowser for help on using the repository browser.