source: trunk/src/STFitter.cpp@ 2674

Last change on this file since 2674 was 2666, checked in by Malte Marquarding, 12 years ago

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