source: trunk/src/SDFitter.cc @ 108

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

added fixed parameters
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.2 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    Int mn,mx;
87    mn = 0;
88    mx = m_.nelements()-1;
89    for (uInt i=0; i<m_.nelements();++i) {
90      if (m_[i]) {
91        mn = i;
92        break;
93      }
94    }
95    for (uInt j=m_.nelements()-1; j>=0;--j) {
96      if (m_[j]) {
97        mx = j;
98        break;
99      }
100    }
101    //mn = 0+x_.nelements()/10;
102    //mx = x_.nelements()-x_.nelements()/10;
103    estimator.setRegion(mn,mx);
104    //estimator.setWindowing(True);
105    SpectralList listGauss = estimator.estimate(x_, y_);
106    Gaussian1D<Float>* g;
107    parameters_.resize(n*3);
108    uInt count = 0;
109    for (uInt i=0; i<n;i++) {
110        g = dynamic_cast<Gaussian1D<Float>* >(funcs_[i]);
111        if (g) {
112            (*g)[0] = listGauss[i].getAmpl();
113            (*g)[1] = listGauss[i].getCenter();
114            (*g)[2] = listGauss[i].getFWHM();
115            ++count;
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    //CompiledFunction<AutoDiff<Float> > comp;
280    //Polynomial<AutoDiff<Float> > poly;
281    CompoundFunction<AutoDiff<Float> > func;
282    if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) != 0) {
283        //computeEstimates();
284        for (uInt i=0; i<funcs_.nelements(); i++) {
285            Gaussian1D<AutoDiff<Float> > gauss;
286            for (uInt j=0; j<funcs_[i]->nparameters(); j++) {
287                gauss[j] = AutoDiff<Float>((*funcs_[i])[j], gauss.nparameters(), j);
288                gauss.mask(j) = funcs_[i]->mask(j);
289            }
290            func.addFunction(gauss);
291        }
292    } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
293        Polynomial<AutoDiff<Float> > poly(funcs_[0]->nparameters()-1);
294        for (uInt j=0; j<funcs_[0]->nparameters(); j++) {
295            poly[j] = AutoDiff<Float>(0, poly.nparameters(), j);
296            poly.mask(j) = funcs_[0]->mask(j);
297        }
298        func.addFunction(poly);
299    } else if (dynamic_cast<CompiledFunction<Float>* >(funcs_[0]) != 0) {
300
301//         CompiledFunction<AutoDiff<Float> > comp;
302//         for (uInt j=0; j<funcs_[0]->nparameters(); j++) {
303//             comp[j] = AutoDiff<Float>(0, comp.nparameters(), j);
304//             comp.mask(j) = funcs_[0]->mask(j);
305//         }
306//         func.addFunction(comp);
307
308        cout << "NYI." << endl;
309    } else {
310        throw (AipsError("Fitter not set up correctly."));
311    }
312    fitter.setFunction(func);
313    fitter.setMaxIter(50+funcs_.nelements()*10);
314    // Convergence criterium
315    fitter.setCriteria(0.001);
316    // Fit
317    Vector<Float> sigma(x_.nelements());
318    sigma = 1.0;
319    //Vector<Float> sol;
320    parameters_.resize();
321    //Vector<Float> err;
322    error_.resize();
323    parameters_ = fitter.fit(x_, y_, sigma, &m_);
324    /*
325    CompoundFunction<Float> f;
326    for (uInt i=0; i<funcs_.nelements(); i++) {
327        f.addFunction(*funcs_[i]);
328    }
329    f.parameters().setParameters(parameters_);
330    */
331    error_ = fitter.errors();
332    chisquared_ = fitter.getChi2();
333    residual_.resize();
334    residual_ =  y_;
335    fitter.residual(residual_,x_);
336    // use fitter.residual(model=True) to get the model
337    thefit_.resize(x_.nelements());
338    fitter.residual(thefit_,x_,True);
339    return true;
340}
Note: See TracBrowser for help on using the repository browser.