LinearFit.h

Go to the documentation of this file.
00001 //# LinearFit.h: Class for linear least-squares fit.
00002 //#
00003 //# Copyright (C) 1995,1999,2000,2001,2002,2004
00004 //# Associated Universities, Inc. Washington DC, USA.
00005 //#
00006 //# This library is free software; you can redistribute it and/or modify it
00007 //# under the terms of the GNU Library General Public License as published by
00008 //# the Free Software Foundation; either version 2 of the License, or (at your
00009 //# option) any later version.
00010 //#
00011 //# This library is distributed in the hope that it will be useful, but WITHOUT
00012 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00013 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00014 //# License for more details.
00015 //#
00016 //# You should have received a copy of the GNU Library General Public License
00017 //# along with this library; if not, write to the Free Software Foundation,
00018 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00019 //#
00020 //# Correspondence concerning AIPS++ should be addressed as follows:
00021 //#        Internet email: aips2-request@nrao.edu.
00022 //#        Postal address: AIPS++ Project Office
00023 //#                        National Radio Astronomy Observatory
00024 //#                        520 Edgemont Road
00025 //#                        Charlottesville, VA 22903-2475 USA
00026 //#
00027 //# $Id$
00028 
00029 #ifndef SCIMATH_LINEARFIT_H
00030 #define SCIMATH_LINEARFIT_H
00031 
00032 //# Includes
00033 #include <casacore/casa/aips.h>
00034 #include <casacore/scimath/Fitting/GenericL2Fit.h>
00035 
00036 namespace casacore { //# NAMESPACE CASACORE - BEGIN
00037 
00038 //# Forward declarations
00039 
00040 // <summary> Class for linear least-squares fit.
00041 // </summary>
00042 //
00043 // <reviewed reviewer="wbrouw" date="2004/06/15" tests="tLinearFitSVD.cc"
00044 //       demos="">
00045 // </reviewed>
00046 //
00047 // <prerequisite>
00048 //   <li> <linkto class="Functional">Functional</linkto> 
00049 //   <li> <linkto class="Function">Function</linkto> 
00050 //   <li> <linkto module="Fitting">Fitting</linkto>
00051 // </prerequisite>
00052 //
00053 // <etymology>
00054 // A set of data point is fit with some functional equation.
00055 // The equations solved are linear equations.  The functions 
00056 // themselves however can be wildly nonlinear.
00057 // </etymology>
00058 //
00059 // <synopsis>
00060 // NOTE: Constraints added. Documentation out of date at moment, check
00061 // the tLinearFitSVD and tNonLinearFirLM programs for examples.
00062 //
00063 // The following is a brief summary of the linear least-squares fit problem.
00064 // See module header, <linkto module="Fitting">Fitting</linkto>,
00065 // for a more complete description.  
00066 //
00067 // Given a set of N data points (measurements), (x(i), y(i)) i = 0,...,N-1, 
00068 // along with a set of standard deviations, sigma(i), for the data points, 
00069 // and M specified functions, f(j)(x) j = 0,...,M-1, we form a linear 
00070 // combination of the functions: 
00071 // <srcblock>
00072 // z(i) = a(0)f(0)(x(i)) + a(1)f(1)(x(i)) + ... + a(M-1)f(M-1)(x(i)),
00073 // </srcblock>
00074 // where a(j) j = 0,...,M-1 are a set of parameters to be determined.
00075 // The linear least-squares fit tries to minimize
00076 // <srcblock>
00077 // chi-square = [(y(0)-z(0))/sigma(0)]^2 + [(y(1)-z(1))/sigma(1)]^2 + ... 
00078 //              + [(y(N-1)-z(N-1))/sigma(N-1)]^2.
00079 // </srcblock>
00080 // by adjusting {a(j)} in the equation. 
00081 //
00082 // For complex numbers, <code>[(y(i)-z(i))/sigma(i)]^2</code> in chi-square 
00083 // is replaced by
00084 // <code>[(y(i)-z(i))/sigma(i)]*conjugate([(y(i)-z(i))/sigma(i)])</code>
00085 //
00086 // For multidimensional functions, x(i) is a vector, and
00087 // <srcblock> 
00088 // f(j)(x(i)) = f(j)(x(i,0), x(i,1), x(i,2), ...)
00089 // </srcblock>
00090 //
00091 // Normally, it is necessary that N > M for the solutions to be valid, since 
00092 // there must be more data points than model parameters to be solved.
00093 //
00094 // If the measurement errors (standard deviation sigma) are not known 
00095 // at all, they can all be set to one initially.  In this case, we assume all 
00096 // measurements have the same standard deviation, after minimizing
00097 // chi-square, we recompute
00098 // <srcblock>  
00099 // sigma^2 = {(y(0)-z(0))^2 + (y(1)-z(1))^2 + ... 
00100 //           + (y(N-1)-z(N-1))^2}/(N-M) = chi-square/(N-M).
00101 // </srcblock> 
00102 //
00103 // A statistic weight can also be assigned to each measurement if the 
00104 // standard deviation is not available.  sigma can be calculated from
00105 // <srcblock>
00106 // sigma = 1/ sqrt(weight)
00107 // </srcblock>
00108 // Alternatively a 'weight' switch can be set with <src>asWeight()</src>.
00109 // For best arithmetic performance, weight should be normalized to a maximum
00110 // value of one. Having a large weight value can sometimes lead to overflow
00111 // problems.
00112 //
00113 // The function to be fitted to the data can be given as an instance of the
00114 // <linkto class="Function">Function</linkto> class.
00115 // One can also form a sum of functions using the
00116 // <linkto class="CompoundFunction">CompoundFunction</linkto>.  
00117 //
00118 // For small datasets the usage of the calls is:
00119 // <ul>
00120 //  <li> Create a functional description of the parameters
00121 //  <li> Create a fitter: LinearFit<T> fitter();
00122 //  <li> Set the functional representation: fitter.setFunction()
00123 //  <li> Do the fit to the data: fitter.fit(x, data, sigma)
00124 //      (or do a number of calls to buildNormalMatrix(x, data, sigma)
00125 //      and finish of with fitter.fit() or fitter.sol())
00126 //  <li> if needed the covariance; residuals; chiSquared, parameter errors
00127 //       can all be obtained
00128 // </ul>
00129 // Note that the fitter is reusable. An example is given in the following.
00130 //
00131 // The solution of a fit always produces the total number of parameters given 
00132 // to the fitter. I.e. including any parameters that were fixed. In the
00133 // latter case the solution returned will be the fixed value.
00134 // 
00135 // <templating arg=T>
00136 // <li> Float
00137 // <li> Double
00138 // <li> Complex
00139 // <li> DComplex   
00140 // </templating>
00141 //
00142 // If there are a large number of unknowns or a large number of data points
00143 // machine memory limits (or timing reasons) may not allow a complete
00144 // in-core fitting to be performed.  In this case one can incrementally
00145 // build the normal equation (see buildNormalMatrix()).
00146 //
00147 // The normal operation of the class tests for real inversion problems
00148 // only. If tests are needed for almost collinear columns in the
00149 // solution matrix, the collinearity can be set as the square of the sine of
00150 // the minimum angle allowed.
00151 //
00152 // Singular Value Decomposition is supported by the
00153 // <linkto class=LinearFitSVD>LinearFitSVD</linkto> class,
00154 // which has a behaviour completely identical to this class (apart from a
00155 // default collinearity of 1e-8). 
00156 //
00157 // Other information (see a.o. <linkto class=LSQFit>LSQFit</linkto>) can
00158 // be set and obtained as well.
00159 // </synopsis>
00160 //
00161 // <motivation>
00162 // The creation of this class was driven by the need to write code
00163 // to perform baseline fitting or continuum subtraction.
00164 // </motivation>
00165 
00166 // <example>
00167 //# /// redo example
00168 // In the following a polynomial is fitted through the first 20 prime numbers.
00169 // The data is given in the x vector (1 to 20) and in the primesTable
00170 // (2, 3, ..., 71) (see tLinearFitSVD test program). In the following
00171 // all four methods to calculate a polynomial through the data is used
00172 // <srcblock>
00173 //      // The list of coordinate x-values
00174 //      Vector<Double> x(nPrimes);
00175 //      indgen((Array<Double>&)x, 1.0);  // 1, 2, ...
00176 //      Vector<Double> primesTable(nPrimes);
00177 //      for (uInt i=1; i < nPrimes; i++) {
00178 //        primesTable(i) =
00179 //         Primes::nextLargerPrimeThan(Int(primesTable(i-1)+0.01));
00180 //      };   
00181 //      Vector<Double> sigma(nPrimes);
00182 //      sigma = 1.0;
00183 //      // The fitter
00184 //      LinearFit<Double> fitter;
00185 //      Polynomial<AutoDiff<Double> > combination(2);
00186 //      // Get the solution
00187 //      fitter.setFunction(combination);
00188 //      Vector<Double> solution = fitter.fit(x, primesTable, sigma);
00189 //      // create a special function (should probably at beginning)
00190 //      static void myfnc(Vector<Double> &y, const Double x) {
00191 //      y(0) = 1; for (uInt i=1; i<y.nelements(); i++) y(i) = y(i-1)*x; };
00192 //      fitter.setFunction(3, &myfnc);
00193 //      solution = fitter.fit(x, primesTable, sigma);
00194 //      // Create the direct coefficients table
00195 //      fitter.setFunction(3);
00196 //      Matrix<Double> xx(nPrimes, 3);
00197 //      for (uInt i=0; i<nPrimes; i++) {
00198 //        xx(i,0) = 1;
00199 //        for (uInt j=1; j<3; j++) xx(i,j) = xx(i,j-1)*Double(i+1);
00200 //      };
00201 //      solution = fitter.fit(xx, primesTable, sigma);
00202 // </srcblock>
00203 // In the test program examples are given on how to get the other
00204 // information, and other examples.
00205 // </example>
00206 
00207 template<class T> class LinearFit : public GenericL2Fit<T>
00208 {
00209 public: 
00210   //# Constructors
00211   // Create a fitter: the normal way to generate a fitter object. Necessary
00212   // data will be deduced from the Functional provided with
00213   // <src>setFunction()</src>
00214   LinearFit();
00215   // Copy constructor (deep copy)
00216   LinearFit(const LinearFit &other);
00217   // Assignment (deep copy)
00218   LinearFit &operator=(const LinearFit &other);
00219 
00220   // Destructor
00221   virtual ~LinearFit();
00222   
00223   //# Member functions
00224 
00225 protected:
00226   //#Data
00227 
00228   //# Member functions
00229   // Generalised fitter
00230   virtual Bool fitIt
00231     (Vector<typename FunctionTraits<T>::BaseType> &sol,
00232      const Array<typename FunctionTraits<T>::BaseType> &x, 
00233      const Vector<typename FunctionTraits<T>::BaseType> &y,
00234      const Vector<typename FunctionTraits<T>::BaseType> *const sigma,
00235      const Vector<Bool> *const mask=0);
00236 
00237 private:
00238   //# Data
00239 
00240   //# Member functions
00241 
00242 protected:
00243   //# Make members of parent classes known.
00244   using GenericL2Fit<T>::pCount_p;
00245   using GenericL2Fit<T>::ptr_derive_p;
00246   using GenericL2Fit<T>::sol_p;
00247   using GenericL2Fit<T>::solved_p;
00248   using GenericL2Fit<T>::nr_p;
00249   using GenericL2Fit<T>::svd_p;
00250   using GenericL2Fit<T>::condEq_p;
00251   using GenericL2Fit<T>::err_p;
00252   using GenericL2Fit<T>::errors_p;
00253   using GenericL2Fit<T>::COLLINEARITY;
00254   using GenericL2Fit<T>::buildConstraint;
00255   using GenericL2Fit<T>::fillSVDConstraints;
00256 };
00257 
00258 
00259 } //# NAMESPACE CASACORE - END
00260 
00261 #ifndef CASACORE_NO_AUTO_TEMPLATES
00262 #include <casacore/scimath/Fitting/LinearFit.tcc>
00263 #endif //# CASACORE_NO_AUTO_TEMPLATES
00264 #endif
00265 
00266 
00267 
00268 
00269 
00270 
00271 
00272 
00273 
00274 
00275 
00276 
00277 
00278 
00279 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1