casacore
Loading...
Searching...
No Matches
LinearFit.h
Go to the documentation of this file.
1//# LinearFit.h: Class for linear least-squares fit.
2//#
3//# Copyright (C) 1995,1999,2000,2001,2002,2004
4//# Associated Universities, Inc. Washington DC, USA.
5//#
6//# This library is free software; you can redistribute it and/or modify it
7//# under the terms of the GNU Library General Public License as published by
8//# the Free Software Foundation; either version 2 of the License, or (at your
9//# option) any later version.
10//#
11//# This library is distributed in the hope that it will be useful, but WITHOUT
12//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14//# License for more details.
15//#
16//# You should have received a copy of the GNU Library General Public License
17//# along with this library; if not, write to the Free Software Foundation,
18//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19//#
20//# Correspondence concerning AIPS++ should be addressed as follows:
21//# Internet email: aips2-request@nrao.edu.
22//# Postal address: AIPS++ Project Office
23//# National Radio Astronomy Observatory
24//# 520 Edgemont Road
25//# Charlottesville, VA 22903-2475 USA
26//#
27//# $Id$
28
29#ifndef SCIMATH_LINEARFIT_H
30#define SCIMATH_LINEARFIT_H
31
32//# Includes
33#include <casacore/casa/aips.h>
34#include <casacore/scimath/Fitting/GenericL2Fit.h>
35
36namespace casacore { //# NAMESPACE CASACORE - BEGIN
37
38//# Forward declarations
39
40// <summary> Class for linear least-squares fit.
41// </summary>
42//
43// <reviewed reviewer="wbrouw" date="2004/06/15" tests="tLinearFitSVD.cc"
44// demos="">
45// </reviewed>
46//
47// <prerequisite>
48// <li> <linkto class="Functional">Functional</linkto>
49// <li> <linkto class="Function">Function</linkto>
50// <li> <linkto module="Fitting">Fitting</linkto>
51// </prerequisite>
52//
53// <etymology>
54// A set of data point is fit with some functional equation.
55// The equations solved are linear equations. The functions
56// themselves however can be wildly nonlinear.
57// </etymology>
58//
59// <synopsis>
60// NOTE: Constraints added. Documentation out of date at moment, check
61// the tLinearFitSVD and tNonLinearFirLM programs for examples.
62//
63// The following is a brief summary of the linear least-squares fit problem.
64// See module header, <linkto module="Fitting">Fitting</linkto>,
65// for a more complete description.
66//
67// Given a set of N data points (measurements), (x(i), y(i)) i = 0,...,N-1,
68// along with a set of standard deviations, sigma(i), for the data points,
69// and M specified functions, f(j)(x) j = 0,...,M-1, we form a linear
70// combination of the functions:
71// <srcblock>
72// z(i) = a(0)f(0)(x(i)) + a(1)f(1)(x(i)) + ... + a(M-1)f(M-1)(x(i)),
73// </srcblock>
74// where a(j) j = 0,...,M-1 are a set of parameters to be determined.
75// The linear least-squares fit tries to minimize
76// <srcblock>
77// chi-square = [(y(0)-z(0))/sigma(0)]^2 + [(y(1)-z(1))/sigma(1)]^2 + ...
78// + [(y(N-1)-z(N-1))/sigma(N-1)]^2.
79// </srcblock>
80// by adjusting {a(j)} in the equation.
81//
82// For complex numbers, <code>[(y(i)-z(i))/sigma(i)]^2</code> in chi-square
83// is replaced by
84// <code>[(y(i)-z(i))/sigma(i)]*conjugate([(y(i)-z(i))/sigma(i)])</code>
85//
86// For multidimensional functions, x(i) is a vector, and
87// <srcblock>
88// f(j)(x(i)) = f(j)(x(i,0), x(i,1), x(i,2), ...)
89// </srcblock>
90//
91// Normally, it is necessary that N > M for the solutions to be valid, since
92// there must be more data points than model parameters to be solved.
93//
94// If the measurement errors (standard deviation sigma) are not known
95// at all, they can all be set to one initially. In this case, we assume all
96// measurements have the same standard deviation, after minimizing
97// chi-square, we recompute
98// <srcblock>
99// sigma^2 = {(y(0)-z(0))^2 + (y(1)-z(1))^2 + ...
100// + (y(N-1)-z(N-1))^2}/(N-M) = chi-square/(N-M).
101// </srcblock>
102//
103// A statistic weight can also be assigned to each measurement if the
104// standard deviation is not available. sigma can be calculated from
105// <srcblock>
106// sigma = 1/ sqrt(weight)
107// </srcblock>
108// Alternatively a 'weight' switch can be set with <src>asWeight()</src>.
109// For best arithmetic performance, weight should be normalized to a maximum
110// value of one. Having a large weight value can sometimes lead to overflow
111// problems.
112//
113// The function to be fitted to the data can be given as an instance of the
114// <linkto class="Function">Function</linkto> class.
115// One can also form a sum of functions using the
116// <linkto class="CompoundFunction">CompoundFunction</linkto>.
117//
118// For small datasets the usage of the calls is:
119// <ul>
120// <li> Create a functional description of the parameters
121// <li> Create a fitter: LinearFit<T> fitter();
122// <li> Set the functional representation: fitter.setFunction()
123// <li> Do the fit to the data: fitter.fit(x, data, sigma)
124// (or do a number of calls to buildNormalMatrix(x, data, sigma)
125// and finish of with fitter.fit() or fitter.sol())
126// <li> if needed the covariance; residuals; chiSquared, parameter errors
127// can all be obtained
128// </ul>
129// Note that the fitter is reusable. An example is given in the following.
130//
131// The solution of a fit always produces the total number of parameters given
132// to the fitter. I.e. including any parameters that were fixed. In the
133// latter case the solution returned will be the fixed value.
134//
135// <templating arg=T>
136// <li> Float
137// <li> Double
138// <li> Complex
139// <li> DComplex
140// </templating>
141//
142// If there are a large number of unknowns or a large number of data points
143// machine memory limits (or timing reasons) may not allow a complete
144// in-core fitting to be performed. In this case one can incrementally
145// build the normal equation (see buildNormalMatrix()).
146//
147// The normal operation of the class tests for real inversion problems
148// only. If tests are needed for almost collinear columns in the
149// solution matrix, the collinearity can be set as the square of the sine of
150// the minimum angle allowed.
151//
152// Singular Value Decomposition is supported by the
153// <linkto class=LinearFitSVD>LinearFitSVD</linkto> class,
154// which has a behaviour completely identical to this class (apart from a
155// default collinearity of 1e-8).
156//
157// Other information (see a.o. <linkto class=LSQFit>LSQFit</linkto>) can
158// be set and obtained as well.
159// </synopsis>
160//
161// <motivation>
162// The creation of this class was driven by the need to write code
163// to perform baseline fitting or continuum subtraction.
164// </motivation>
165
166// <example>
167//# /// redo example
168// In the following a polynomial is fitted through the first 20 prime numbers.
169// The data is given in the x vector (1 to 20) and in the primesTable
170// (2, 3, ..., 71) (see tLinearFitSVD test program). In the following
171// all four methods to calculate a polynomial through the data is used
172// <srcblock>
173// // The list of coordinate x-values
174// Vector<Double> x(nPrimes);
175// indgen((Array<Double>&)x, 1.0); // 1, 2, ...
176// Vector<Double> primesTable(nPrimes);
177// for (uInt i=1; i < nPrimes; i++) {
178// primesTable(i) =
179// Primes::nextLargerPrimeThan(Int(primesTable(i-1)+0.01));
180// };
181// Vector<Double> sigma(nPrimes);
182// sigma = 1.0;
183// // The fitter
184// LinearFit<Double> fitter;
185// Polynomial<AutoDiff<Double> > combination(2);
186// // Get the solution
187// fitter.setFunction(combination);
188// Vector<Double> solution = fitter.fit(x, primesTable, sigma);
189// // create a special function (should probably at beginning)
190// static void myfnc(Vector<Double> &y, const Double x) {
191// y(0) = 1; for (uInt i=1; i<y.nelements(); i++) y(i) = y(i-1)*x; };
192// fitter.setFunction(3, &myfnc);
193// solution = fitter.fit(x, primesTable, sigma);
194// // Create the direct coefficients table
195// fitter.setFunction(3);
196// Matrix<Double> xx(nPrimes, 3);
197// for (uInt i=0; i<nPrimes; i++) {
198// xx(i,0) = 1;
199// for (uInt j=1; j<3; j++) xx(i,j) = xx(i,j-1)*Double(i+1);
200// };
201// solution = fitter.fit(xx, primesTable, sigma);
202// </srcblock>
203// In the test program examples are given on how to get the other
204// information, and other examples.
205// </example>
206
207template<class T> class LinearFit : public GenericL2Fit<T>
208{
209public:
210 //# Constructors
211 // Create a fitter: the normal way to generate a fitter object. Necessary
212 // data will be deduced from the Functional provided with
213 // <src>setFunction()</src>
215 // Copy constructor (deep copy)
216 LinearFit(const LinearFit &other);
217 // Assignment (deep copy)
219
220 // Destructor
221 virtual ~LinearFit();
222
223 //# Member functions
224
225protected:
226 //#Data
227
228 //# Member functions
229 // Generalised fitter
230 virtual Bool fitIt
232 const Array<typename FunctionTraits<T>::BaseType> &x,
233 const Vector<typename FunctionTraits<T>::BaseType> &y,
234 const Vector<typename FunctionTraits<T>::BaseType> *const sigma,
235 const Vector<Bool> *const mask=0);
236
237private:
238 //# Data
239
240 //# Member functions
241
242protected:
243 //# Make members of parent classes known.
244 using GenericL2Fit<T>::pCount_p;
246 using GenericL2Fit<T>::sol_p;
247 using GenericL2Fit<T>::solved_p;
248 using GenericL2Fit<T>::nr_p;
249 using GenericL2Fit<T>::svd_p;
250 using GenericL2Fit<T>::condEq_p;
251 using GenericL2Fit<T>::err_p;
252 using GenericL2Fit<T>::errors_p;
256};
257
258
259} //# NAMESPACE CASACORE - END
260
261#ifndef CASACORE_NO_AUTO_TEMPLATES
262#include <casacore/scimath/Fitting/LinearFit.tcc>
263#endif //# CASACORE_NO_AUTO_TEMPLATES
264#endif
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
T BaseType
Template base type.
Bool solved_p
Have solution.
void fillSVDConstraints()
Get the SVD constraints.
uInt nr_p
The rank of the solution.
Vector< typename FunctionTraits< T >::BaseType > condEq_p
Condition equation parameters (for number of adjustable parameters)
const Double COLLINEARITY
Default collinearity test for SVD.
Vector< typename FunctionTraits< T >::BaseType > err_p
Local error area.
void buildConstraint()
Build the constraint equations.
uInt pCount_p
Number of available parameters.
Vector< typename FunctionTraits< T >::BaseType > sol_p
Local solution area.
Bool svd_p
SVD indicator.
Bool errors_p
Have errors.
Function< typename FunctionTraits< T >::DiffType, typename FunctionTraits< T >::DiffType > * ptr_derive_p
Function to use in evaluating condition equation.
static const String sol
Definition LSQFit.h:857
virtual ~LinearFit()
Destructor.
LinearFit & operator=(const LinearFit &other)
Assignment (deep copy)
LinearFit(const LinearFit &other)
Copy constructor (deep copy)
LinearFit()
Create a fitter: the normal way to generate a fitter object.
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)
Generalised fitter.
this file contains all the compiler specific defines
Definition mainpage.dox:28
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