source: branches/alma/src/STFitter.cpp@ 1669

Last change on this file since 1669 was 1616, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-729, CAS-1147

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Use LogIO instead of std::cout and std::cerr.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 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 1616 2009-08-07 05:49:27Z TakeshiNakazato $
30//#---------------------------------------------------------------------------
31#include <casa/aips.h>
32#include <casa/Arrays/ArrayMath.h>
33#include <casa/Arrays/ArrayLogical.h>
34#include <casa/Logging/LogIO.h>
35#include <scimath/Fitting.h>
36#include <scimath/Fitting/LinearFit.h>
37#include <scimath/Functionals/CompiledFunction.h>
38#include <scimath/Functionals/CompoundFunction.h>
39#include <scimath/Functionals/Gaussian1D.h>
40#include <scimath/Functionals/Polynomial.h>
41#include <scimath/Mathematics/AutoDiff.h>
42#include <scimath/Mathematics/AutoDiffMath.h>
43#include <scimath/Fitting/NonLinearFitLM.h>
44#include <components/SpectralComponents/SpectralEstimate.h>
45
46#include "STFitter.h"
47
48using namespace asap;
49using namespace casa;
50
51Fitter::Fitter()
52{
53}
54
55Fitter::~Fitter()
56{
57 reset();
58}
59
60void Fitter::clear()
61{
62 for (uInt i=0;i< funcs_.nelements();++i) {
63 delete funcs_[i]; funcs_[i] = 0;
64 }
65 funcs_.resize(0,True);
66 parameters_.resize();
67 fixedpar_.resize();
68 error_.resize();
69 thefit_.resize();
70 estimate_.resize();
71 chisquared_ = 0.0;
72}
73
74void Fitter::reset()
75{
76 clear();
77 x_.resize();
78 y_.resize();
79 m_.resize();
80}
81
82
83bool Fitter::computeEstimate() {
84 if (x_.nelements() == 0 || y_.nelements() == 0)
85 throw (AipsError("No x/y data specified."));
86
87 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) == 0)
88 return false;
89 uInt n = funcs_.nelements();
90 SpectralEstimate estimator(n);
91 estimator.setQ(5);
92 Int mn,mx;
93 mn = 0;
94 mx = m_.nelements()-1;
95 for (uInt i=0; i<m_.nelements();++i) {
96 if (m_[i]) {
97 mn = i;
98 break;
99 }
100 }
101 for (uInt j=m_.nelements()-1; j>=0;--j) {
102 if (m_[j]) {
103 mx = j;
104 break;
105 }
106 }
107 //mn = 0+x_.nelements()/10;
108 //mx = x_.nelements()-x_.nelements()/10;
109 estimator.setRegion(mn,mx);
110 //estimator.setWindowing(True);
111 SpectralList listGauss = estimator.estimate(x_, y_);
112 parameters_.resize(n*3);
113 Gaussian1D<Float>* g = 0;
114 for (uInt i=0; i<n;i++) {
115 g = dynamic_cast<Gaussian1D<Float>* >(funcs_[i]);
116 if (g) {
117 (*g)[0] = listGauss[i].getAmpl();
118 (*g)[1] = listGauss[i].getCenter();
119 (*g)[2] = listGauss[i].getFWHM();
120 }
121 }
122 estimate_.resize();
123 listGauss.evaluate(estimate_,x_);
124 return true;
125}
126
127std::vector<float> Fitter::getEstimate() const
128{
129 if (estimate_.nelements() == 0)
130 throw (AipsError("No estimate set."));
131 std::vector<float> stlout;
132 estimate_.tovector(stlout);
133 return stlout;
134}
135
136
137bool Fitter::setExpression(const std::string& expr, int ncomp)
138{
139 clear();
140 if (expr == "gauss") {
141 if (ncomp < 1) throw (AipsError("Need at least one gaussian to fit."));
142 funcs_.resize(ncomp);
143 for (Int k=0; k<ncomp; ++k) {
144 funcs_[k] = new Gaussian1D<Float>();
145 }
146 } else if (expr == "poly") {
147 funcs_.resize(1);
148 funcs_[0] = new Polynomial<Float>(ncomp);
149 } else {
150 //cerr << " compiled functions not yet implemented" << endl;
151 LogIO os( LogOrigin( "Fitter", "setExpression()", WHERE ) ) ;
152 os << LogIO::WARN << " compiled functions not yet implemented" << LogIO::POST;
153 //funcs_.resize(1);
154 //funcs_[0] = new CompiledFunction<Float>();
155 //funcs_[0]->setFunction(String(expr));
156 return false;
157 }
158 return true;
159}
160
161bool Fitter::setData(std::vector<float> absc, std::vector<float> spec,
162 std::vector<bool> mask)
163{
164 x_.resize();
165 y_.resize();
166 m_.resize();
167 // convert std::vector to casa Vector
168 Vector<Float> tmpx(absc);
169 Vector<Float> tmpy(spec);
170 Vector<Bool> tmpm(mask);
171 AlwaysAssert(tmpx.nelements() == tmpy.nelements(), AipsError);
172 x_ = tmpx;
173 y_ = tmpy;
174 m_ = tmpm;
175 return true;
176}
177
178std::vector<float> Fitter::getResidual() const
179{
180 if (residual_.nelements() == 0)
181 throw (AipsError("Function not yet fitted."));
182 std::vector<float> stlout;
183 residual_.tovector(stlout);
184 return stlout;
185}
186
187std::vector<float> Fitter::getFit() const
188{
189 Vector<Float> out = thefit_;
190 std::vector<float> stlout;
191 out.tovector(stlout);
192 return stlout;
193
194}
195
196std::vector<float> Fitter::getErrors() const
197{
198 Vector<Float> out = error_;
199 std::vector<float> stlout;
200 out.tovector(stlout);
201 return stlout;
202}
203
204bool Fitter::setParameters(std::vector<float> params)
205{
206 Vector<Float> tmppar(params);
207 if (funcs_.nelements() == 0)
208 throw (AipsError("Function not yet set."));
209 if (parameters_.nelements() > 0 && tmppar.nelements() != parameters_.nelements())
210 throw (AipsError("Number of parameters inconsistent with function."));
211 if (parameters_.nelements() == 0) {
212 parameters_.resize(tmppar.nelements());
213 if (tmppar.nelements() != fixedpar_.nelements()) {
214 fixedpar_.resize(tmppar.nelements());
215 fixedpar_ = False;
216 }
217 }
218 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) != 0) {
219 uInt count = 0;
220 for (uInt j=0; j < funcs_.nelements(); ++j) {
221 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
222 (funcs_[j]->parameters())[i] = tmppar[count];
223 parameters_[count] = tmppar[count];
224 ++count;
225 }
226 }
227 } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
228 for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
229 parameters_[i] = tmppar[i];
230 (funcs_[0]->parameters())[i] = tmppar[i];
231 }
232 }
233 // reset
234 if (params.size() == 0) {
235 parameters_.resize();
236 fixedpar_.resize();
237 }
238 return true;
239}
240
241bool Fitter::setFixedParameters(std::vector<bool> fixed)
242{
243 if (funcs_.nelements() == 0)
244 throw (AipsError("Function not yet set."));
245 if (fixedpar_.nelements() > 0 && fixed.size() != fixedpar_.nelements())
246 throw (AipsError("Number of mask elements inconsistent with function."));
247 if (fixedpar_.nelements() == 0) {
248 fixedpar_.resize(parameters_.nelements());
249 fixedpar_ = False;
250 }
251 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) != 0) {
252 uInt count = 0;
253 for (uInt j=0; j < funcs_.nelements(); ++j) {
254 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
255 funcs_[j]->mask(i) = !fixed[count];
256 fixedpar_[count] = fixed[count];
257 ++count;
258 }
259 }
260 } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
261 for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
262 fixedpar_[i] = fixed[i];
263 funcs_[0]->mask(i) = !fixed[i];
264 }
265 }
266 return true;
267}
268
269std::vector<float> Fitter::getParameters() const {
270 Vector<Float> out = parameters_;
271 std::vector<float> stlout;
272 out.tovector(stlout);
273 return stlout;
274}
275
276std::vector<bool> Fitter::getFixedParameters() const {
277 Vector<Bool> out(parameters_.nelements());
278 if (fixedpar_.nelements() == 0) {
279 return std::vector<bool>();
280 //throw (AipsError("No parameter mask set."));
281 } else {
282 out = fixedpar_;
283 }
284 std::vector<bool> stlout;
285 out.tovector(stlout);
286 return stlout;
287}
288
289float Fitter::getChisquared() const {
290 return chisquared_;
291}
292
293bool Fitter::fit() {
294 NonLinearFitLM<Float> fitter;
295 CompoundFunction<Float> func;
296
297 uInt n = funcs_.nelements();
298 for (uInt i=0; i<n; ++i) {
299 func.addFunction(*funcs_[i]);
300 }
301
302 fitter.setFunction(func);
303 fitter.setMaxIter(50+n*10);
304 // Convergence criterium
305 fitter.setCriteria(0.001);
306
307 // Fit
308 Vector<Float> sigma(x_.nelements());
309 sigma = 1.0;
310
311 parameters_.resize();
312 parameters_ = fitter.fit(x_, y_, sigma, &m_);
313 if ( !fitter.converged() ) {
314 return false;
315 }
316 std::vector<float> ps;
317 parameters_.tovector(ps);
318 setParameters(ps);
319
320 error_.resize();
321 error_ = fitter.errors();
322
323 chisquared_ = fitter.getChi2();
324
325 residual_.resize();
326 residual_ = y_;
327 fitter.residual(residual_,x_);
328 // use fitter.residual(model=True) to get the model
329 thefit_.resize(x_.nelements());
330 fitter.residual(thefit_,x_,True);
331 return true;
332}
333
334bool Fitter::lfit() {
335 LinearFit<Float> fitter;
336 CompoundFunction<Float> func;
337
338 uInt n = funcs_.nelements();
339 for (uInt i=0; i<n; ++i) {
340 func.addFunction(*funcs_[i]);
341 }
342
343 fitter.setFunction(func);
344 //fitter.setMaxIter(50+n*10);
345 // Convergence criterium
346 //fitter.setCriteria(0.001);
347
348 // Fit
349 Vector<Float> sigma(x_.nelements());
350 sigma = 1.0;
351
352 parameters_.resize();
353 parameters_ = fitter.fit(x_, y_, sigma, &m_);
354 std::vector<float> ps;
355 parameters_.tovector(ps);
356 setParameters(ps);
357
358 error_.resize();
359 error_ = fitter.errors();
360
361 chisquared_ = fitter.getChi2();
362
363 residual_.resize();
364 residual_ = y_;
365 fitter.residual(residual_,x_);
366 // use fitter.residual(model=True) to get the model
367 thefit_.resize(x_.nelements());
368 fitter.residual(thefit_,x_,True);
369 return true;
370}
371
372std::vector<float> Fitter::evaluate(int whichComp) const
373{
374 std::vector<float> stlout;
375 uInt idx = uInt(whichComp);
376 Float y;
377 if ( idx < funcs_.nelements() ) {
378 for (uInt i=0; i<x_.nelements(); ++i) {
379 y = (*funcs_[idx])(x_[i]);
380 stlout.push_back(float(y));
381 }
382 }
383 return stlout;
384}
385
Note: See TracBrowser for help on using the repository browser.