casa  5.7.0-16
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
34 namespace 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 
153 template<class T> class NonLinearFit : public GenericL2Fit<T>
154 {
155 public:
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)
169  NonLinearFit(const NonLinearFit &other);
170  // Assignment (deep copy)
171  NonLinearFit &operator=(const NonLinearFit &other);
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
184  uInt currentIteration() const { return maxiter_p - curiter_p; };
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 
197 protected:
198  //#Data
199  // Maximum number of iterations
200  uInt maxiter_p;
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 
217 private:
218  //# Data
219 
220  //# Member functions
221 
222 protected:
223  //# Make members of parent classes known.
230  using GenericL2Fit<T>::nr_p;
238  using GenericL2Fit<T>::set;
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
A 1-D Specialization of the Array class.
uInt getMaxIter() const
getMaxIter() queries what the maximum number of iterations currently is
Definition: NonLinearFit.h:181
LatticeExprNode mask(const LatticeExprNode &expr)
This function returns the mask of the given expression.
T BaseType
Template base type.
Bool converge_p
Has fit converged.
Definition: NonLinearFit.h:206
void setCriteria(const Double criteria=CRITERIUM)
setCriteria() sets the convergence criteria.
Definition: NonLinearFit.h:189
NonLinearFit(Bool svd=False)
Create a fitter: the normal way to generate a fitter object.
uInt curiter_p
Current iteration number.
Definition: NonLinearFit.h:202
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.
NonLinearFit & operator=(const NonLinearFit &other)
Assignment (deep copy)
uInt maxiter_p
Maximum number of iterations.
Definition: NonLinearFit.h:195
void setMaxIter(uInt maxIter=MAXITER)
setMaxIter() sets the maximum number of iterations to do before stopping.
double Double
Definition: aipstype.h:55
Generic base class for least-squares fit.
Definition: FittingProxy.h:39
Class for non-linear least-squares fit.
Definition: NonLinearFit.h:153
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:42
const Bool False
Definition: aipstype.h:44
template &lt;class T, class U&gt; class vector;
Definition: MSFlagger.h:37
Bool converged() const
Check to see if the fit has converged.
Definition: NonLinearFit.h:195
static const uInt MAXITER
Default maximum number of iterations (30)
Definition: NonLinearFit.h:158
static const String sol
Definition: LSQFit.h:857
virtual ~NonLinearFit()
Destructor.
Double getCriteria() const
getCriteria() queries the current criteria
Definition: NonLinearFit.h:192
static const Double CRITERIUM
Default convergence criterium (0.001)
Definition: NonLinearFit.h:160
Double criterium_p
Convergence criteria.
Definition: NonLinearFit.h:204
uInt currentIteration() const
currentIteration() queries what the current iteration is
Definition: NonLinearFit.h:184
unsigned int uInt
Definition: aipstype.h:51
#define casacore
&lt;X11/Intrinsic.h&gt; #defines true, false, casacore::Bool, and String.
Definition: X11Intrinsic.h:42