source: branches/hpc33/src/STFitter.cpp@ 2429

Last change on this file since 2429 was 2394, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

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...

Change according to an update of SpectralElement class.


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