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 |
|
---|
51 | using namespace asap;
|
---|
52 | using namespace casa;
|
---|
53 |
|
---|
54 | Fitter::Fitter()
|
---|
55 | {
|
---|
56 | }
|
---|
57 |
|
---|
58 | Fitter::~Fitter()
|
---|
59 | {
|
---|
60 | reset();
|
---|
61 | }
|
---|
62 |
|
---|
63 | void 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 |
|
---|
77 | void Fitter::reset()
|
---|
78 | {
|
---|
79 | clear();
|
---|
80 | x_.resize();
|
---|
81 | y_.resize();
|
---|
82 | m_.resize();
|
---|
83 | constraints_.clear();
|
---|
84 | }
|
---|
85 |
|
---|
86 |
|
---|
87 | bool 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 |
|
---|
139 | std::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 |
|
---|
149 | bool 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 |
|
---|
200 | bool 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 |
|
---|
217 | std::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 |
|
---|
226 | std::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 |
|
---|
235 | std::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 |
|
---|
243 | bool 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 |
|
---|
298 | void 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 |
|
---|
306 | void 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 |
|
---|
316 | bool 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 |
|
---|
362 | std::vector<float> Fitter::getParameters() const {
|
---|
363 | Vector<Float> out = parameters_;
|
---|
364 | std::vector<float> stlout;
|
---|
365 | out.tovector(stlout);
|
---|
366 | return stlout;
|
---|
367 | }
|
---|
368 |
|
---|
369 | std::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 |
|
---|
382 | float Fitter::getChisquared() const {
|
---|
383 | return chisquared_;
|
---|
384 | }
|
---|
385 |
|
---|
386 | bool 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 |
|
---|
428 | bool 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 |
|
---|
458 | std::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 |
|
---|
472 | STFitEntry 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 | }
|
---|