source: trunk/src/STFitter.cpp@ 2454

Last change on this file since 2454 was 2445, checked in by Malte Marquarding, 13 years ago

Temporary fix for divergence of casacore and casapy
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.0 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 2445 2012-03-29 01:37:19Z 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}
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#ifdef USE_CASAPY
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#else
127 (*g)[0] = listGauss[i].getAmpl();
128 (*g)[1] = listGauss[i].getCenter();
129 (*g)[2] = listGauss[i].getFWHM();
130#endif
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
297bool Fitter::setFixedParameters(std::vector<bool> fixed)
298{
299 if (funcs_.nelements() == 0)
300 throw (AipsError("Function not yet set."));
301 if (fixedpar_.nelements() > 0 && fixed.size() != fixedpar_.nelements())
302 throw (AipsError("Number of mask elements inconsistent with function."));
303 if (fixedpar_.nelements() == 0) {
304 fixedpar_.resize(parameters_.nelements());
305 fixedpar_ = False;
306 }
307 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) != 0) {
308 uInt count = 0;
309 for (uInt j=0; j < funcs_.nelements(); ++j) {
310 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
311 funcs_[j]->mask(i) = !fixed[count];
312 fixedpar_[count] = fixed[count];
313 ++count;
314 }
315 }
316 } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
317 uInt count = 0;
318 for (uInt j=0; j < funcs_.nelements(); ++j) {
319 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
320 funcs_[j]->mask(i) = !fixed[count];
321 fixedpar_[count] = fixed[count];
322 ++count;
323 }
324 }
325 } else if (dynamic_cast<Sinusoid1D<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<Polynomial<Float>* >(funcs_[0]) != 0) {
335 for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
336 fixedpar_[i] = fixed[i];
337 funcs_[0]->mask(i) = !fixed[i];
338 }
339 }
340 return true;
341}
342
343std::vector<float> Fitter::getParameters() const {
344 Vector<Float> out = parameters_;
345 std::vector<float> stlout;
346 out.tovector(stlout);
347 return stlout;
348}
349
350std::vector<bool> Fitter::getFixedParameters() const {
351 Vector<Bool> out(parameters_.nelements());
352 if (fixedpar_.nelements() == 0) {
353 return std::vector<bool>();
354 //throw (AipsError("No parameter mask set."));
355 } else {
356 out = fixedpar_;
357 }
358 std::vector<bool> stlout;
359 out.tovector(stlout);
360 return stlout;
361}
362
363float Fitter::getChisquared() const {
364 return chisquared_;
365}
366
367bool Fitter::fit() {
368 NonLinearFitLM<Float> fitter;
369 CompoundFunction<Float> func;
370
371 uInt n = funcs_.nelements();
372 for (uInt i=0; i<n; ++i) {
373 func.addFunction(*funcs_[i]);
374 }
375
376 fitter.setFunction(func);
377 fitter.setMaxIter(50+n*10);
378 // Convergence criterium
379 fitter.setCriteria(0.001);
380
381 // Fit
382 Vector<Float> sigma(x_.nelements());
383 sigma = 1.0;
384
385 parameters_.resize();
386 parameters_ = fitter.fit(x_, y_, sigma, &m_);
387 if ( !fitter.converged() ) {
388 return false;
389 }
390 std::vector<float> ps;
391 parameters_.tovector(ps);
392 setParameters(ps);
393
394 error_.resize();
395 error_ = fitter.errors();
396
397 chisquared_ = fitter.getChi2();
398
399 residual_.resize();
400 residual_ = y_;
401 fitter.residual(residual_,x_);
402 // use fitter.residual(model=True) to get the model
403 thefit_.resize(x_.nelements());
404 fitter.residual(thefit_,x_,True);
405 return true;
406}
407
408bool Fitter::lfit() {
409 LinearFit<Float> fitter;
410 CompoundFunction<Float> func;
411
412 uInt n = funcs_.nelements();
413 for (uInt i=0; i<n; ++i) {
414 func.addFunction(*funcs_[i]);
415 }
416
417 fitter.setFunction(func);
418 //fitter.setMaxIter(50+n*10);
419 // Convergence criterium
420 //fitter.setCriteria(0.001);
421
422 // Fit
423 Vector<Float> sigma(x_.nelements());
424 sigma = 1.0;
425
426 parameters_.resize();
427 parameters_ = fitter.fit(x_, y_, sigma, &m_);
428 std::vector<float> ps;
429 parameters_.tovector(ps);
430 setParameters(ps);
431
432 error_.resize();
433 error_ = fitter.errors();
434
435 chisquared_ = fitter.getChi2();
436
437 residual_.resize();
438 residual_ = y_;
439 fitter.residual(residual_,x_);
440 // use fitter.residual(model=True) to get the model
441 thefit_.resize(x_.nelements());
442 fitter.residual(thefit_,x_,True);
443 return true;
444}
445
446std::vector<float> Fitter::evaluate(int whichComp) const
447{
448 std::vector<float> stlout;
449 uInt idx = uInt(whichComp);
450 Float y;
451 if ( idx < funcs_.nelements() ) {
452 for (uInt i=0; i<x_.nelements(); ++i) {
453 y = (*funcs_[idx])(x_[i]);
454 stlout.push_back(float(y));
455 }
456 }
457 return stlout;
458}
459
460STFitEntry Fitter::getFitEntry() const
461{
462 STFitEntry fit;
463 fit.setParameters(getParameters());
464 fit.setErrors(getErrors());
465 fit.setComponents(funccomponents_);
466 fit.setFunctions(funcnames_);
467 fit.setParmasks(getFixedParameters());
468 return fit;
469}
Note: See TracBrowser for help on using the repository browser.