casacore
NonLinearFit.h
Go to the documentation of this file.
1//# NonLinearFit.h: Class for non-linear least-squares fit.
2//# Copyright (C) 1994,1995,1996,1999,2000,2001,2002,2004
3//# Associated Universities, Inc. Washington DC, USA.
4//#
5//# This library is free software; you can redistribute it and/or modify it
6//# under the terms of the GNU Library General Public License as published by
7//# the Free Software Foundation; either version 2 of the License, or (at your
8//# option) any later version.
9//#
10//# This library is distributed in the hope that it will be useful, but WITHOUT
11//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13//# License for more details.
14//#
15//# You should have received a copy of the GNU Library General Public License
16//# along with this library; if not, write to the Free Software Foundation,
17//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18//#
19//# Correspondence concerning AIPS++ should be addressed as follows:
20//# Internet email: aips2-request@nrao.edu.
21//# Postal address: AIPS++ Project Office
22//# National Radio Astronomy Observatory
23//# 520 Edgemont Road
24//# Charlottesville, VA 22903-2475 USA
25//#
26//# $Id$
27
28#ifndef SCIMATH_NONLINEARFIT_H
29#define SCIMATH_NONLINEARFIT_H
30
31//# Includes
32#include <casacore/casa/aips.h>
33#include <casacore/scimath/Fitting/GenericL2Fit.h>
34namespace casacore { //# begin namesapce casa
35//# Forward declarations
36
37//
38// <summary> Class for non-linear least-squares fit.
39// </summary>
40//
41// <reviewed reviewer="wbrouw" date="2006/06/15" tests="tNonLinearFitLM.cc">
42// </reviewed>
43//
44// <prerequisite>
45// <li> <linkto class="Functional">Functional</linkto>
46// <li> <linkto class="Function">Function</linkto>
47// <li> <linkto module="Fitting">Fitting</linkto>
48// </prerequisite>
49//
50// <etymology>
51// A nonlinear function is used to fit a set of data points.
52// </etymology>
53//
54// <synopsis>
55// NOTE: Constraints added. Documentation out of date at moment, check
56// the tLinearFitSVD and tNonLinearFirLM programs for examples.
57//
58// The following is a brief summary of the non-linear least-squares fit
59// problem.
60// See module header, <linkto module="Fitting">Fitting</linkto>,
61// for a more complete description.
62//
63// Given a set of N data points (measurements), (x(i), y(i)) i = 0,...,N-1,
64// along with a set of standard deviations, sigma(i), for the data points,
65// and a specified non-linear function, f(x;a) where a = a(j) j = 0,...,M-1
66// are a set of parameters to be
67// determined, the non-linear least-squares fit tries to minimize
68// <srcblock>
69// chi-square = [(y(0)-f(x(0);a)/sigma(0)]^2 + [(y(1)-f(x(1);a)/sigma(1)]^2 +
70// ... + [(y(N-1)-f(x(N-1);a))/sigma(N-1)]^2.
71// </srcblock>
72// by adjusting {a(j)} in the equation.
73//
74// For multidimensional functions, x(i) is a vector, and
75// <srcblock>
76// f(x(i);a) = f(x(i,0), x(i,1), x(i,2), ...;a)
77// </srcblock>
78//
79// If the measurement errors (standard deviation sigma) are not known
80// at all, they can all be set to one initially. In this case, we assume all
81// measurements have the same standard deviation, after minimizing
82// chi-square, we recompute
83// <srcblock>
84// sigma^2 = {(y(0)-z(0))^2 + (y(1)-z(1))^2 + ...
85// + (y(N-1)-z(N-1))^2}/(N-M) = chi-square/(N-M).
86// </srcblock>
87//
88// A statistic weight can be also be assigned to each measurement if the
89// standard deviation is not available. sigma can be calculated from
90// <srcblock>
91// sigma = 1/ sqrt(weight)
92// </srcblock>
93// Alternatively a 'weight' switch can be set with <src>asWeight()</src>.
94// For best arithmetic performance, weight should be normalized to a maximum
95// value of one. Having a large weight value can sometimes lead to overflow
96// problems.
97//
98// The function to be fitted to the data can be given as an instance of the
99// <linkto class="Function">Function</linkto> class.
100// One can also form a sum of functions using the
101// <linkto class="CompoundFunction">CompoundFunction</linkto>.
102//
103// For small datasets the usage of the calls is:
104// <ul>
105// <li> Create a functional description of the parameters
106// <li> Create a fitter: NonLinearFit<T> fitter();
107// <li> Set the functional representation: fitter.setFunction()
108// <li> Do the fit to the data: fitter.fit(x, data, sigma)
109// (or do a number of calls to buildNormalMatrix(x, data, sigma)
110// and finish of with fitter.fit() or fitter.sol())
111// <li> if needed the covariance; residuals; chiSquared, parameter errors
112// can all be obtained
113// </ul>
114// Note that the fitter is reusable. An example is given in the following.
115//
116// The solution of a fit always produces the total number of parameters given
117// to the fitter. I.e. including any parameters that were fixed. In the
118// latter case the solution returned will be the fixed value.
119//
120// <templating arg=T>
121// <li> Float
122// <li> Double
123// <li> Complex
124// <li> DComplex
125// </templating>
126//
127// If there are a large number of unknowns or a large number of data points
128// machine memory limits (or timing reasons) may not allow a complete
129// in-core fitting to be performed. In this case one can incrementally
130// build the normal equation (see buildNormalMatrix()).
131//
132// The normal operation of the class tests for real inversion problems
133// only. If tests are needed for almost collinear columns in the
134// solution matrix, the collinearity can be set as the square of the sine of
135// the minimum angle allowed.
136//
137// Singular Value Decomposition is supported by setting the 'svd' switch,
138// which has a behaviour completely identical to, apart from a
139// default collinearity check of 1e-8.
140//
141// Other information (see a.o. <linkto class=LSQaips>LSQaips</linkto>) can
142// be set and obtained as well.
143// </synopsis>
144//
145// <motivation>
146// The creation of the class module was driven by the need to write code
147// to fit Gaussian functions to data points.
148// </motivation>
149//
150// <example>
151// </example>
152
153template<class T> class NonLinearFit : public GenericL2Fit<T>
154{
155public:
156 //# Constants
157 // Default maximum number of iterations (30)
158 static const uInt MAXITER = 30;
159 // Default convergence criterium (0.001)
160 static const Double CRITERIUM;
161
162 //# Constructors
163 // Create a fitter: the normal way to generate a fitter object. Necessary
164 // data will be deduced from the Functional provided with
165 // <src>setFunction()</src>.
166 // Create optionally a fitter with SVD behaviour specified.
167 explicit NonLinearFit(Bool svd=False);
168 // Copy constructor (deep copy)
170 // Assignment (deep copy)
172
173 // Destructor
174 virtual ~NonLinearFit();
175
176 // setMaxIter() sets the maximum number of iterations to do before stopping.
177 // Default value is 30.
178 void setMaxIter(uInt maxIter=MAXITER);
179
180 // getMaxIter() queries what the maximum number of iterations currently is
181 uInt getMaxIter() const { return maxiter_p; };
182
183 // currentIteration() queries what the current iteration is
185
186 // setCriteria() sets the convergence criteria. The actual value and
187 // its interpretation depends on the derived class used to do the
188 // actual iteration. Default value is 0.001.
189 void setCriteria(const Double criteria=CRITERIUM) {criterium_p = criteria; };
190
191 // getCriteria() queries the current criteria
192 Double getCriteria() const { return criterium_p; };
193
194 // Check to see if the fit has converged
195 Bool converged() const { return converge_p; };
196
197protected:
198 //#Data
199 // Maximum number of iterations
201 // Current iteration number
203 // Convergence criteria
205 // Has fit converged
207
208 //# Member functions
209 // Generalised fitter
210 virtual Bool fitIt
212 const Array<typename FunctionTraits<T>::BaseType> &x,
213 const Vector<typename FunctionTraits<T>::BaseType> &y,
214 const Vector<typename FunctionTraits<T>::BaseType> *const sigma,
215 const Vector<Bool> *const mask=0) = 0;
216
217private:
218 //# Data
219
220 //# Member functions
221
222protected:
223 //# Make members of parent classes known.
224 using GenericL2Fit<T>::pCount_p;
226 using GenericL2Fit<T>::arg_p;
227 using GenericL2Fit<T>::sol_p;
228 using GenericL2Fit<T>::fsol_p;
229 using GenericL2Fit<T>::solved_p;
230 using GenericL2Fit<T>::nr_p;
231 using GenericL2Fit<T>::svd_p;
232 using GenericL2Fit<T>::condEq_p;
233 using GenericL2Fit<T>::fullEq_p;
234 using GenericL2Fit<T>::err_p;
235 using GenericL2Fit<T>::ferr_p;
236 using GenericL2Fit<T>::errors_p;
237 using GenericL2Fit<T>::valder_p;
238 using GenericL2Fit<T>::set;
241 using GenericL2Fit<T>::isReady;
242};
243
244} //# End namespace casacore
245#ifndef CASACORE_NO_AUTO_TEMPLATES
246#include <casacore/scimath/Fitting/NonLinearFit.tcc>
247#endif //# CASACORE_NO_AUTO_TEMPLATES
248#endif
T BaseType
Template base type.
Bool solved_p
Have solution.
Definition: GenericL2Fit.h:486
void fillSVDConstraints()
Get the SVD constraints.
uInt nr_p
The rank of the solution.
Definition: GenericL2Fit.h:493
Vector< typename FunctionTraits< T >::BaseType > condEq_p
Condition equation parameters (for number of adjustable parameters)
Definition: GenericL2Fit.h:495
Vector< typename FunctionTraits< T >::BaseType > err_p
Local error area.
Definition: GenericL2Fit.h:510
void buildConstraint()
Build the constraint equations.
Vector< typename FunctionTraits< T >::BaseType > ferr_p
Definition: GenericL2Fit.h:511
uInt pCount_p
Number of available parameters.
Definition: GenericL2Fit.h:480
Vector< typename FunctionTraits< T >::BaseType > sol_p
Local solution area.
Definition: GenericL2Fit.h:505
Bool svd_p
SVD indicator.
Definition: GenericL2Fit.h:462
Bool errors_p
Have errors.
Definition: GenericL2Fit.h:488
Function< typename FunctionTraits< T >::DiffType, typename FunctionTraits< T >::DiffType > * ptr_derive_p
Function to use in evaluating condition equation.
Definition: GenericL2Fit.h:465
Vector< typename FunctionTraits< T >::BaseType > fsol_p
Definition: GenericL2Fit.h:506
Vector< typename FunctionTraits< T >::ArgType > arg_p
Contiguous argument areas.
Definition: GenericL2Fit.h:500
FunctionTraits< T >::DiffType valder_p
Local value and derivatives.
Definition: GenericL2Fit.h:514
Vector< typename FunctionTraits< T >::BaseType > fullEq_p
Equation for all available parameters.
Definition: GenericL2Fit.h:497
void set(uInt nUnknowns, uInt nConstraints=0)
Set new sizes (default is for Real)
LSQFit::ReadyCode isReady() const
Ask the state of the non-linear solutions.
Definition: LSQFit.h:749
static const String sol
Definition: LSQFit.h:857
static const Double CRITERIUM
Default convergence criterium (0.001)
Definition: NonLinearFit.h:160
virtual ~NonLinearFit()
Destructor.
uInt curiter_p
Current iteration number.
Definition: NonLinearFit.h:202
Double getCriteria() const
getCriteria() queries the current criteria
Definition: NonLinearFit.h:192
virtual Bool fitIt(Vector< typename FunctionTraits< T >::BaseType > &sol, const Array< typename FunctionTraits< T >::BaseType > &x, const Vector< typename FunctionTraits< T >::BaseType > &y, const Vector< typename FunctionTraits< T >::BaseType > *const sigma, const Vector< Bool > *const mask=0)=0
Generalised fitter.
static const uInt MAXITER
Default maximum number of iterations (30)
Definition: NonLinearFit.h:158
NonLinearFit(Bool svd=False)
Create a fitter: the normal way to generate a fitter object.
void setMaxIter(uInt maxIter=MAXITER)
setMaxIter() sets the maximum number of iterations to do before stopping.
void setCriteria(const Double criteria=CRITERIUM)
setCriteria() sets the convergence criteria.
Definition: NonLinearFit.h:189
uInt getMaxIter() const
getMaxIter() queries what the maximum number of iterations currently is
Definition: NonLinearFit.h:181
Bool converged() const
Check to see if the fit has converged.
Definition: NonLinearFit.h:195
uInt currentIteration() const
currentIteration() queries what the current iteration is
Definition: NonLinearFit.h:184
NonLinearFit & operator=(const NonLinearFit &other)
Assignment (deep copy)
uInt maxiter_p
Maximum number of iterations.
Definition: NonLinearFit.h:200
NonLinearFit(const NonLinearFit &other)
Copy constructor (deep copy)
Double criterium_p
Convergence criteria.
Definition: NonLinearFit.h:204
Bool converge_p
Has fit converged.
Definition: NonLinearFit.h:206
this file contains all the compiler specific defines
Definition: mainpage.dox:28
const Bool False
Definition: aipstype.h:44
unsigned int uInt
Definition: aipstype.h:51
LatticeExprNode mask(const LatticeExprNode &expr)
This function returns the mask of the given expression.
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:42
double Double
Definition: aipstype.h:55