source: trunk/src/SDFitter.cc @ 91

Last change on this file since 91 was 91, checked in by mar637, 20 years ago

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