source: trunk/src/STFitter.cpp@ 2369

Last change on this file since 2369 was 2163, checked in by Malte Marquarding, 14 years ago

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