source: trunk/src/STFitter.cpp@ 3101

Last change on this file since 3101 was 3087, checked in by Takeshi Nakazato, 9 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes/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...


Bug fix on STFitter::computeEstimate. Introduced size check of func_ member variable to avoid segmentation fault. Uncommented lines that seems to be commented unintentionally.

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