source: branches/Release-1-fixes/src/SDFitter.cc@ 573

Last change on this file since 573 was 517, checked in by mar637, 20 years ago
  • updated to reflect Wim's changes to Functionals
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.8 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDFitter.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:
30//#---------------------------------------------------------------------------
31#include <casa/aips.h>
32#include <casa/Arrays/ArrayMath.h>
33#include <casa/Arrays/ArrayLogical.h>
34#include <scimath/Fitting.h>
35#include <scimath/Fitting/LinearFit.h>
36#include <scimath/Functionals/CompiledFunction.h>
37#include <scimath/Functionals/CompoundFunction.h>
38#include <scimath/Functionals/Gaussian1D.h>
39#include <scimath/Functionals/Polynomial.h>
40#include <scimath/Mathematics/AutoDiff.h>
41#include <scimath/Mathematics/AutoDiffMath.h>
42#include <scimath/Fitting/NonLinearFitLM.h>
43#include <components/SpectralComponents/SpectralEstimate.h>
44
45#include "SDFitter.h"
46using namespace asap;
47using namespace casa;
48
49SDFitter::SDFitter()
50{
51}
52
53SDFitter::~SDFitter()
54{
55 reset();
56}
57
58void SDFitter::clear()
59{
60 for (uInt i=0;i< funcs_.nelements();++i) {
61 delete funcs_[i]; funcs_[i] = 0;
62 }
63 parameters_.resize();
64 error_.resize();
65 thefit_.resize();
66 estimate_.resize();
67 chisquared_ = 0.0;
68}
69
70void SDFitter::reset()
71{
72 clear();
73 x_.resize();
74 y_.resize();
75 m_.resize();
76}
77
78
79bool SDFitter::computeEstimate() {
80 if (x_.nelements() == 0 || y_.nelements() == 0)
81 throw (AipsError("No x/y data specified."));
82
83 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) == 0)
84 return false;
85 uInt n = funcs_.nelements();
86 SpectralEstimate estimator(n);
87 estimator.setQ(5);
88 Int mn,mx;
89 mn = 0;
90 mx = m_.nelements()-1;
91 for (uInt i=0; i<m_.nelements();++i) {
92 if (m_[i]) {
93 mn = i;
94 break;
95 }
96 }
97 for (uInt j=m_.nelements()-1; j>=0;--j) {
98 if (m_[j]) {
99 mx = j;
100 break;
101 }
102 }
103 mn = 0+x_.nelements()/10;
104 mx = x_.nelements()-x_.nelements()/10;
105 estimator.setRegion(mn,mx);
106 //estimator.setWindowing(True);
107 SpectralList listGauss = estimator.estimate(x_, y_);
108 parameters_.resize(n*3);
109 Gaussian1D<Float>* g = 0;
110 for (uInt i=0; i<n;i++) {
111 g = dynamic_cast<Gaussian1D<Float>* >(funcs_[i]);
112 if (g) {
113 (*g)[0] = listGauss[i].getAmpl();
114 (*g)[1] = listGauss[i].getCenter();
115 (*g)[2] = listGauss[i].getFWHM();
116 }
117 }
118 estimate_.resize();
119 listGauss.evaluate(estimate_,x_);
120 return true;
121}
122
123std::vector<float> SDFitter::getEstimate() const
124{
125 if (estimate_.nelements() == 0)
126 throw (AipsError("No estimate set."));
127 std::vector<float> stlout;
128 estimate_.tovector(stlout);
129 return stlout;
130}
131
132
133bool SDFitter::setExpression(const std::string& expr, int ncomp)
134{
135 clear();
136 if (expr == "gauss") {
137 if (ncomp < 1) throw (AipsError("Need at least one gaussian to fit."));
138 funcs_.resize(ncomp);
139 for (Int k=0; k<ncomp; ++k) {
140 funcs_[k] = new Gaussian1D<Float>();
141 }
142 } else if (expr == "poly") {
143 funcs_.resize(1);
144 funcs_[0] = new Polynomial<Float>(ncomp);
145 } else {
146 cerr << " compiled functions not yet implemented" << endl;
147 //funcs_.resize(1);
148 //funcs_[0] = new CompiledFunction<Float>();
149 //funcs_[0]->setFunction(String(expr));
150 return false;
151 }
152 return true;
153}
154
155bool SDFitter::setData(std::vector<float> absc, std::vector<float> spec,
156 std::vector<bool> mask)
157{
158 x_.resize();
159 y_.resize();
160 m_.resize();
161 // convert std::vector to casa Vector
162 Vector<Float> tmpx(absc);
163 Vector<Float> tmpy(spec);
164 Vector<Bool> tmpm(mask);
165 AlwaysAssert(tmpx.nelements() == tmpy.nelements(), AipsError);
166 x_ = tmpx;
167 y_ = tmpy;
168 m_ = tmpm;
169 return true;
170}
171
172std::vector<float> SDFitter::getResidual() const
173{
174 if (residual_.nelements() == 0)
175 throw (AipsError("Function not yet fitted."));
176 std::vector<float> stlout;
177 residual_.tovector(stlout);
178 return stlout;
179}
180
181std::vector<float> SDFitter::getFit() const
182{
183 Vector<Float> out = thefit_;
184 std::vector<float> stlout;
185 out.tovector(stlout);
186 return stlout;
187
188}
189
190std::vector<float> SDFitter::getErrors() const
191{
192 Vector<Float> out = error_;
193 std::vector<float> stlout;
194 out.tovector(stlout);
195 return stlout;
196}
197
198bool SDFitter::setParameters(std::vector<float> params)
199{
200 Vector<Float> tmppar(params);
201 if (funcs_.nelements() == 0)
202 throw (AipsError("Function not yet set."));
203 if (parameters_.nelements() > 0 && tmppar.nelements() != parameters_.nelements())
204 throw (AipsError("Number of parameters inconsistent with function."));
205 if (parameters_.nelements() == 0)
206 parameters_.resize(tmppar.nelements());
207 fixedpar_.resize(tmppar.nelements());
208 fixedpar_ = False;
209 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) != 0) {
210 uInt count = 0;
211 for (uInt j=0; j < funcs_.nelements(); ++j) {
212 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
213 (funcs_[j]->parameters())[i] = tmppar[count];
214 parameters_[count] = tmppar[count];
215 ++count;
216 }
217 }
218 } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
219 for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
220 parameters_[i] = tmppar[i];
221 (funcs_[0]->parameters())[i] = tmppar[i];
222 }
223 }
224 return true;
225}
226
227bool SDFitter::setFixedParameters(std::vector<bool> fixed)
228{
229 Vector<Bool> tmp(fixed);
230 if (funcs_.nelements() == 0)
231 throw (AipsError("Function not yet set."));
232 if (fixedpar_.nelements() > 0 && tmp.nelements() != fixedpar_.nelements())
233 throw (AipsError("Number of mask elements inconsistent with function."));
234 if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) != 0) {
235 uInt count = 0;
236 for (uInt j=0; j < funcs_.nelements(); ++j) {
237 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
238 funcs_[j]->mask(i) = !tmp[count];
239 fixedpar_[count] = !tmp[count];
240 ++count;
241 }
242 }
243 } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
244 for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
245 fixedpar_[i] = tmp[i];
246 funcs_[0]->mask(i) = tmp[i];
247 }
248 }
249 //fixedpar_ = !tmpmsk;
250 return true;
251}
252
253std::vector<float> SDFitter::getParameters() const {
254 Vector<Float> out = parameters_;
255 std::vector<float> stlout;
256 out.tovector(stlout);
257 return stlout;
258}
259
260std::vector<bool> SDFitter::getFixedParameters() const {
261 Vector<Bool> out(parameters_.nelements());
262 if (fixedpar_.nelements() == 0) {
263 out = False;
264 //throw (AipsError("No parameter mask set."));
265 } else {
266 out = fixedpar_;
267 }
268 std::vector<bool> stlout;
269 out.tovector(stlout);
270 return stlout;
271}
272
273float SDFitter::getChisquared() const {
274 return chisquared_;
275}
276
277bool SDFitter::fit() {
278 NonLinearFitLM<Float> fitter;
279 CompoundFunction<Float> func;
280 const uInt n = funcs_.nelements();
281 for (uInt i=0; i<n; ++i) {
282 func.addFunction(*funcs_[i]);
283 }
284 fitter.setFunction(func);
285 fitter.setMaxIter(50+n*10);
286 // Convergence criterium
287 fitter.setCriteria(0.001);
288
289 // Fit
290 Vector<Float> sigma(x_.nelements());
291 sigma = 1.0;
292
293 parameters_.resize();
294 parameters_ = fitter.fit(x_, y_, sigma, &m_);
295 std::vector<float> ps;
296 parameters_.tovector(ps);
297 setParameters(ps);
298 error_.resize();
299 error_ = fitter.errors();
300
301 chisquared_ = fitter.getChi2();
302
303 residual_.resize();
304 residual_ = y_;
305 fitter.residual(residual_,x_);
306
307 // use fitter.residual(model=True) to get the model
308 thefit_.resize(x_.nelements());
309 fitter.residual(thefit_,x_,True);
310 return true;
311}
312
313
314std::vector<float> SDFitter::evaluate(int whichComp) const
315{
316 std::vector<float> stlout;
317 uInt idx = uInt(whichComp);
318 Float y;
319 if ( idx < funcs_.nelements() ) {
320 for (uInt i=0; i<x_.nelements(); ++i) {
321 y = (*funcs_[idx])(x_[i]);
322 stlout.push_back(float(y));
323 }
324 }
325 return stlout;
326}
327
Note: See TracBrowser for help on using the repository browser.