source: branches/Release2.1.1/src/STFitter.cpp@ 1679

Last change on this file since 1679 was 1232, checked in by mar637, 19 years ago

Merge from Release2.1.0b tag

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.2 KB
Line 
1//#---------------------------------------------------------------------------
2//# Fitter.cc: A Fitter class for spectra
3//#--------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id: STFitter.cpp 1232 2006-09-01 03:14:11Z mar637 $
30//#---------------------------------------------------------------------------
31#include <casa/aips.h>
32#include <casa/Arrays/ArrayMath.h>
33#include <casa/Arrays/ArrayLogical.h>
34#include <scimath/Fitting.h>
35#include <scimath/Fitting/LinearFit.h>
36#include <scimath/Functionals/CompiledFunction.h>
37#include <scimath/Functionals/CompoundFunction.h>
38#include <scimath/Functionals/Gaussian1D.h>
39#include <scimath/Functionals/Polynomial.h>
40#include <scimath/Mathematics/AutoDiff.h>
41#include <scimath/Mathematics/AutoDiffMath.h>
42#include <scimath/Fitting/NonLinearFitLM.h>
43#include <components/SpectralComponents/SpectralEstimate.h>
44
45#include "STFitter.h"
46
47using namespace asap;
48using namespace casa;
49
50Fitter::Fitter()
51{
52}
53
54Fitter::~Fitter()
55{
56 reset();
57}
58
59void Fitter::clear()
60{
61 for (uInt i=0;i< funcs_.nelements();++i) {
62 delete funcs_[i]; funcs_[i] = 0;
63 }
64 funcs_.resize(0,True);
65 parameters_.resize();
66 fixedpar_.resize();
67 error_.resize();
68 thefit_.resize();
69 estimate_.resize();
70 chisquared_ = 0.0;
71}
72
73void Fitter::reset()
74{
75 clear();
76 x_.resize();
77 y_.resize();
78 m_.resize();
79}
80
81
82bool Fitter::computeEstimate() {
83 if (x_.nelements() == 0 || y_.nelements() == 0)
84 throw (AipsError("No x/y data specified."));
85
86 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) == 0)
87 return false;
88 uInt n = funcs_.nelements();
89 SpectralEstimate estimator(n);
90 estimator.setQ(5);
91 Int mn,mx;
92 mn = 0;
93 mx = m_.nelements()-1;
94 for (uInt i=0; i<m_.nelements();++i) {
95 if (m_[i]) {
96 mn = i;
97 break;
98 }
99 }
100 for (uInt j=m_.nelements()-1; j>=0;--j) {
101 if (m_[j]) {
102 mx = j;
103 break;
104 }
105 }
106 //mn = 0+x_.nelements()/10;
107 //mx = x_.nelements()-x_.nelements()/10;
108 estimator.setRegion(mn,mx);
109 //estimator.setWindowing(True);
110 SpectralList listGauss = estimator.estimate(x_, y_);
111 parameters_.resize(n*3);
112 Gaussian1D<Float>* g = 0;
113 for (uInt i=0; i<n;i++) {
114 g = dynamic_cast<Gaussian1D<Float>* >(funcs_[i]);
115 if (g) {
116 (*g)[0] = listGauss[i].getAmpl();
117 (*g)[1] = listGauss[i].getCenter();
118 (*g)[2] = listGauss[i].getFWHM();
119 }
120 }
121 estimate_.resize();
122 listGauss.evaluate(estimate_,x_);
123 return true;
124}
125
126std::vector<float> Fitter::getEstimate() const
127{
128 if (estimate_.nelements() == 0)
129 throw (AipsError("No estimate set."));
130 std::vector<float> stlout;
131 estimate_.tovector(stlout);
132 return stlout;
133}
134
135
136bool Fitter::setExpression(const std::string& expr, int ncomp)
137{
138 clear();
139 if (expr == "gauss") {
140 if (ncomp < 1) throw (AipsError("Need at least one gaussian to fit."));
141 funcs_.resize(ncomp);
142 for (Int k=0; k<ncomp; ++k) {
143 funcs_[k] = new Gaussian1D<Float>();
144 }
145 } else if (expr == "poly") {
146 funcs_.resize(1);
147 funcs_[0] = new Polynomial<Float>(ncomp);
148 } else {
149 cerr << " compiled functions not yet implemented" << endl;
150 //funcs_.resize(1);
151 //funcs_[0] = new CompiledFunction<Float>();
152 //funcs_[0]->setFunction(String(expr));
153 return false;
154 }
155 return true;
156}
157
158bool Fitter::setData(std::vector<float> absc, std::vector<float> spec,
159 std::vector<bool> mask)
160{
161 x_.resize();
162 y_.resize();
163 m_.resize();
164 // convert std::vector to casa Vector
165 Vector<Float> tmpx(absc);
166 Vector<Float> tmpy(spec);
167 Vector<Bool> tmpm(mask);
168 AlwaysAssert(tmpx.nelements() == tmpy.nelements(), AipsError);
169 x_ = tmpx;
170 y_ = tmpy;
171 m_ = tmpm;
172 return true;
173}
174
175std::vector<float> Fitter::getResidual() const
176{
177 if (residual_.nelements() == 0)
178 throw (AipsError("Function not yet fitted."));
179 std::vector<float> stlout;
180 residual_.tovector(stlout);
181 return stlout;
182}
183
184std::vector<float> Fitter::getFit() const
185{
186 Vector<Float> out = thefit_;
187 std::vector<float> stlout;
188 out.tovector(stlout);
189 return stlout;
190
191}
192
193std::vector<float> Fitter::getErrors() const
194{
195 Vector<Float> out = error_;
196 std::vector<float> stlout;
197 out.tovector(stlout);
198 return stlout;
199}
200
201bool Fitter::setParameters(std::vector<float> params)
202{
203 Vector<Float> tmppar(params);
204 if (funcs_.nelements() == 0)
205 throw (AipsError("Function not yet set."));
206 if (parameters_.nelements() > 0 && tmppar.nelements() != parameters_.nelements())
207 throw (AipsError("Number of parameters inconsistent with function."));
208 if (parameters_.nelements() == 0) {
209 parameters_.resize(tmppar.nelements());
210 if (tmppar.nelements() != fixedpar_.nelements()) {
211 fixedpar_.resize(tmppar.nelements());
212 fixedpar_ = False;
213 }
214 }
215 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) != 0) {
216 uInt count = 0;
217 for (uInt j=0; j < funcs_.nelements(); ++j) {
218 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
219 (funcs_[j]->parameters())[i] = tmppar[count];
220 parameters_[count] = tmppar[count];
221 ++count;
222 }
223 }
224 } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
225 for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
226 parameters_[i] = tmppar[i];
227 (funcs_[0]->parameters())[i] = tmppar[i];
228 }
229 }
230 // reset
231 if (params.size() == 0) {
232 parameters_.resize();
233 fixedpar_.resize();
234 }
235 return true;
236}
237
238bool Fitter::setFixedParameters(std::vector<bool> fixed)
239{
240 if (funcs_.nelements() == 0)
241 throw (AipsError("Function not yet set."));
242 if (fixedpar_.nelements() > 0 && fixed.size() != fixedpar_.nelements())
243 throw (AipsError("Number of mask elements inconsistent with function."));
244 if (fixedpar_.nelements() == 0) {
245 fixedpar_.resize(parameters_.nelements());
246 fixedpar_ = False;
247 }
248 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) != 0) {
249 uInt count = 0;
250 for (uInt j=0; j < funcs_.nelements(); ++j) {
251 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
252 funcs_[j]->mask(i) = !fixed[count];
253 fixedpar_[count] = fixed[count];
254 ++count;
255 }
256 }
257 } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
258 for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
259 fixedpar_[i] = fixed[i];
260 funcs_[0]->mask(i) = !fixed[i];
261 }
262 }
263 return true;
264}
265
266std::vector<float> Fitter::getParameters() const {
267 Vector<Float> out = parameters_;
268 std::vector<float> stlout;
269 out.tovector(stlout);
270 return stlout;
271}
272
273std::vector<bool> Fitter::getFixedParameters() const {
274 Vector<Bool> out(parameters_.nelements());
275 if (fixedpar_.nelements() == 0) {
276 return std::vector<bool>();
277 //throw (AipsError("No parameter mask set."));
278 } else {
279 out = fixedpar_;
280 }
281 std::vector<bool> stlout;
282 out.tovector(stlout);
283 return stlout;
284}
285
286float Fitter::getChisquared() const {
287 return chisquared_;
288}
289
290bool Fitter::fit() {
291 NonLinearFitLM<Float> fitter;
292 CompoundFunction<Float> func;
293
294 uInt n = funcs_.nelements();
295 for (uInt i=0; i<n; ++i) {
296 func.addFunction(*funcs_[i]);
297 }
298
299 fitter.setFunction(func);
300 fitter.setMaxIter(50+n*10);
301 // Convergence criterium
302 fitter.setCriteria(0.001);
303
304 // Fit
305 Vector<Float> sigma(x_.nelements());
306 sigma = 1.0;
307
308 parameters_.resize();
309 parameters_ = fitter.fit(x_, y_, sigma, &m_);
310 if ( !fitter.converged() ) {
311 return false;
312 }
313 std::vector<float> ps;
314 parameters_.tovector(ps);
315 setParameters(ps);
316
317 error_.resize();
318 error_ = fitter.errors();
319
320 chisquared_ = fitter.getChi2();
321
322 residual_.resize();
323 residual_ = y_;
324 fitter.residual(residual_,x_);
325 // use fitter.residual(model=True) to get the model
326 thefit_.resize(x_.nelements());
327 fitter.residual(thefit_,x_,True);
328 return true;
329}
330
331
332std::vector<float> Fitter::evaluate(int whichComp) const
333{
334 std::vector<float> stlout;
335 uInt idx = uInt(whichComp);
336 Float y;
337 if ( idx < funcs_.nelements() ) {
338 for (uInt i=0; i<x_.nelements(); ++i) {
339 y = (*funcs_[idx])(x_[i]);
340 stlout.push_back(float(y));
341 }
342 }
343 return stlout;
344}
345
Note: See TracBrowser for help on using the repository browser.