source: trunk/src/STFitter.cpp@ 2683

Last change on this file since 2683 was 2675, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CAS-4429

Ready for Test: 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...

convSupport for gauss and gjinc is set to ceil of truncation
radius, which was round before.


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