source: trunk/src/STFitter.cpp @ 2445

Last change on this file since 2445 was 2445, checked in by Malte Marquarding, 12 years ago

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