casacore
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 
36 namespace 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 
207 template<class T> class LinearFit : public GenericL2Fit<T>
208 {
209 public:
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)
218  LinearFit &operator=(const LinearFit &other);
219 
220  // Destructor
221  virtual ~LinearFit();
222 
223  //# Member functions
224 
225 protected:
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 
237 private:
238  //# Data
239 
240  //# Member functions
241 
242 protected:
243  //# Make members of parent classes known.
248  using GenericL2Fit<T>::nr_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.
static const String sol
Definition: LSQFit.h:857
LinearFit & operator=(const LinearFit &other)
Assignment (deep copy)
virtual ~LinearFit()
Destructor.
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