|
| | LSQFit (uInt nUnknowns, uInt nConstraints=0) |
| | Construct an object with the number of unknowns and constraints, using the default collinearity factor and the default Levenberg-Marquardt adjustment factor. More...
|
| |
| | LSQFit (uInt nUnknowns, const LSQReal &, uInt nConstraints=0) |
| | Allow explicit Real specification. More...
|
| |
| | LSQFit (uInt nUnknowns, const LSQComplex &, uInt nConstraints=0) |
| | Allow explicit Complex specification. More...
|
| |
| | LSQFit () |
| | Default constructor (empty, only usable after a set(nUnknowns)) More...
|
| |
| | LSQFit (const LSQFit &other) |
| | Copy constructor (deep copy) More...
|
| |
| LSQFit & | operator= (const LSQFit &other) |
| | Assignment (deep copy) More...
|
| |
| | ~LSQFit () |
| |
| Bool | invert (uInt &nRank, Bool doSVD=False) |
| | Triangularize the normal equations and determine the rank nRank of the normal equations and, in the case of an SVD solution, the constraint equations. More...
|
| |
| template<class U > |
| void | copy (const Double *beg, const Double *end, U &sol, LSQReal) |
| | Copy date from beg to end; converting if necessary to complex data. More...
|
| |
| template<class U > |
| void | copy (const Double *beg, const Double *end, U &sol, LSQComplex) |
| |
| template<class U > |
| void | copy (const Double *beg, const Double *end, U *sol, LSQReal) |
| |
| template<class U > |
| void | copy (const Double *beg, const Double *end, U *sol, LSQComplex) |
| |
| template<class U > |
| void | uncopy (Double *beg, const Double *end, U &sol, LSQReal) |
| |
| template<class U > |
| void | uncopy (Double *beg, const Double *end, U &sol, LSQComplex) |
| |
| template<class U > |
| void | uncopy (Double *beg, const Double *end, U *sol, LSQReal) |
| |
| template<class U > |
| void | uncopy (Double *beg, const Double *end, U *sol, LSQComplex) |
| |
| template<class U > |
| void | copyDiagonal (U &errors, LSQReal) |
| |
| template<class U > |
| void | copyDiagonal (U &errors, LSQComplex) |
| |
| template<class U > |
| void | solve (U *sol) |
| | Solve normal equations. More...
|
| |
| template<class U > |
| void | solve (std::complex< U > *sol) |
| |
| template<class U > |
| void | solve (U &sol) |
| |
| template<class U > |
| Bool | solveLoop (uInt &nRank, U *sol, Bool doSVD=False) |
| | Solve a loop in a non-linear set. More...
|
| |
| template<class U > |
| Bool | solveLoop (uInt &nRank, std::complex< U > *sol, Bool doSVD=False) |
| |
| template<class U > |
| Bool | solveLoop (uInt &nRank, U &sol, Bool doSVD=False) |
| |
| template<class U > |
| Bool | solveLoop (Double &fit, uInt &nRank, U *sol, Bool doSVD=False) |
| |
| template<class U > |
| Bool | solveLoop (Double &fit, uInt &nRank, std::complex< U > *sol, Bool doSVD=False) |
| |
| template<class U > |
| Bool | solveLoop (Double &fit, uInt &nRank, U &sol, Bool doSVD=False) |
| |
| template<class U , class V > |
| void | makeNorm (const V &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True) |
| | Make normal equations using the cEq condition equation (cArray) (with nUnknowns elements) and a weight weight, given the known observed value obs. More...
|
| |
| template<class U , class V > |
| void | makeNorm (const V &cEq, const U &weight, const U &obs, LSQFit::Real, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const V &cEq, const U &weight, const std::complex< U > &obs, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Complex, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Separable, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::AsReal, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Conjugate, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNorm (uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNorm (uInt nIndex, const W &cEqIndex, const V &cEq, const V &cEq2, const U &weight, const U &obs, const U &obs2, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNorm (uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const U &obs, LSQFit::Real, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNorm (uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex< U > &obs, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNorm (uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Complex, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNorm (uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Separable, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNorm (uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::AsReal, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNorm (uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Conjugate, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const U &obs, LSQFit::Real, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const std::complex< U > &obs, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Complex, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Separable, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const std::complex< U > &obs, LSQFit::AsReal, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V > |
| void | makeNorm (const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Conjugate, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNormSorted (uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U , class V , class W > |
| void | makeNormSorted (uInt nIndex, const W &cEqIndex, const V &cEq, const V &cEq2, const U &weight, const U &obs, const U &obs2, Bool doNorm=True, Bool doKnown=True) |
| |
| template<class U > |
| Bool | getConstraint (uInt n, U *cEq) const |
| | Get the n-th (from 0 to the rank deficiency, or missing rank, see e.g. More...
|
| |
| template<class U > |
| Bool | getConstraint (uInt n, std::complex< U > *cEq) const |
| |
| template<class U > |
| Bool | getConstraint (uInt n, U &cEq) const |
| |
| template<class U , class V > |
| Bool | setConstraint (uInt n, const V &cEq, const U &obs) |
| | Add a new constraint equation (updating nConstraints); or set a numbered constraint equation (0..nConstraints-1). More...
|
| |
| template<class U , class V > |
| Bool | setConstraint (uInt n, const V &cEq, const std::complex< U > &obs) |
| |
| template<class U , class V , class W > |
| Bool | setConstraint (uInt n, uInt nIndex, const W &cEqIndex, const V &cEq, const U &obs) |
| |
| template<class U , class V , class W > |
| Bool | setConstraint (uInt n, uInt nIndex, const W &cEqIndex, const V &cEq, const std::complex< U > &obs) |
| |
| template<class U , class V > |
| Bool | addConstraint (const V &cEq, const U &obs) |
| |
| template<class U , class V > |
| Bool | addConstraint (const V &cEq, const std::complex< U > &obs) |
| |
| template<class U , class V , class W > |
| Bool | addConstraint (uInt nIndex, const W &cEqIndex, const V &cEq, const U &obs) |
| |
| template<class U , class V , class W > |
| Bool | addConstraint (uInt nIndex, const W &cEqIndex, const V &cEq, const std::complex< U > &obs) |
| |
| Bool | merge (const LSQFit &other) |
| | Merge other LSQFit object (i.e. More...
|
| |
| Bool | merge (const LSQFit &other, uInt nIndex, const uInt *nEqIndex) |
| |
| Bool | merge (const LSQFit &other, uInt nIndex, const std::vector< uInt > &nEqIndex) |
| |
| template<class W > |
| Bool | merge (const LSQFit &other, uInt nIndex, const W &nEqIndex) |
| |
| void | reset () |
| | Reset status to empty. More...
|
| |
| void | set (uInt nUnknowns, uInt nConstraints=0) |
| | Set new sizes (default is for Real) More...
|
| |
| void | set (Int nUnknowns, Int nConstraints=0) |
| |
| void | set (uInt nUnknowns, const LSQReal &, uInt nConstraints=0) |
| |
| void | set (Int nUnknowns, const LSQReal &, Int nConstraints=0) |
| |
| void | set (uInt nUnknowns, const LSQComplex &, uInt nConstraints=0) |
| |
| void | set (Int nUnknowns, const LSQComplex &, Int nConstraints=0) |
| |
| void | set (Double factor=1e-6, Double LMFactor=1e-3) |
| | Set new factors (collinearity factor, and Levenberg-Marquardt LMFactor) More...
|
| |
| void | setEpsValue (Double epsval=1e-8) |
| | Set new value solution test. More...
|
| |
| void | setEpsDerivative (Double epsder=1e-8) |
| | Set new derivative test. More...
|
| |
| void | setMaxIter (uInt maxiter=0) |
| | Set maximum number of iterations. More...
|
| |
| uInt | nIterations () const |
| | Get number of iterations done. More...
|
| |
| void | setBalanced (Bool balanced=False) |
| | Set the expected form of the normal equations. More...
|
| |
| LSQFit::ReadyCode | isReady () const |
| | Ask the state of the non-linear solutions. More...
|
| |
| const std::string & | readyText () const |
| |
| template<class U > |
| Bool | getCovariance (U *covar) |
| | Get the covariance matrix (of size nUnknowns * nUnknowns) More...
|
| |
| template<class U > |
| Bool | getCovariance (std::complex< U > *covar) |
| |
| template<class U > |
| Bool | getErrors (U *errors) |
| | Get main diagonal of covariance function (of size nUnknowns) More...
|
| |
| template<class U > |
| Bool | getErrors (std::complex< U > *errors) |
| |
| template<class U > |
| Bool | getErrors (U &errors) |
| |
| uInt | nUnknowns () const |
| | Get the number of unknowns. More...
|
| |
| uInt | nConstraints () const |
| | Get the number of constraints. More...
|
| |
| uInt | getDeficiency () const |
| | Get the rank deficiency
Warning: Note that the number is returned assuming real values; For complex values it has to be halved
More...
|
| |
| Double | getChi () const |
| | Get chi^2 (both are identical); the standard deviation (per observation) and the standard deviation per weight unit. More...
|
| |
| Double | getChi2 () const |
| |
| Double | getSD () const |
| |
| Double | getWeightedSD () const |
| |
| void | debugIt (uInt &nun, uInt &np, uInt &ncon, uInt &ner, uInt &rank, Double *&nEq, Double *&known, Double *&constr, Double *&er, uInt *&piv, Double *&sEq, Double *&sol, Double &prec, Double &nonlin) const |
| | Debug: More...
|
| |
| Bool | fromRecord (String &error, const RecordInterface &in) |
| | Create an LSQFit object from a record. More...
|
| |
| Bool | toRecord (String &error, RecordInterface &out) const |
| | Create a record from an LSQFit object. More...
|
| |
| const String & | ident () const |
| | Get identification of record. More...
|
| |
| void | toAipsIO (AipsIO &) const |
| | Save or restore using AipsIO. More...
|
| |
| void | fromAipsIO (AipsIO &) |
| |
Basic class for the least squares fitting.
Review Status
- Reviewed By:
- Neil Killeen
- Date Reviewed:
- 2000/06/01
- Test programs:
- tLSQFit
Prerequisite
-
Some knowledge of Matrix operations
-
The background information provided in Note 224.
-
Fitting module
Etymology
From Least SQuares and Fitting
Synopsis
The LSQFit class contains the basic functions to do all the fitting described in the Note about fitting. It handles real, and complex equations;
linear and non-linear (Levenberg-Marquardt) solutions;
regular (with optional constraints) or Singular Value Decomposition (SVD).
In essence they are a set of routines to generate normal equations (makeNorm()) in triangular form from a set of condition equations;
to do a Cholesky-type decomposition of the normal equations (either regular or SVD) and test its rank (invert());
to do a quasi inversion of the decomposed equations (solve()) to obtain the solution and/or the errors.
All calculations are done in place. Methods to obtain additional information about the fitting process are available.
This class can be used as a stand-alone class outside of the Casacore environment. In that case the aips.h include file can be replaced if necessary by appropriate typedefs for Double, Float and uInt.
The interface to the methods have standard data or standard STL iterator arguments only. They can be used with any container having an STL random-access iterator interface. Especially they can be used with carrays (necessary templates provided), Casacore Vectors (necessary templates provided in LSQaips), standard random access STL containers (like std::vector).
The normal operation of the class consists of the following steps:
-
Create an LSQFit object. The information that can be provided in the constructor of the object, either directly, or indirectly using the set() commands, is (see Note 224):
-
The number of unknowns that have to be solved for (mandatory)
-
The number of constraint equations you want to use explicitly (defaults to 0, but can be changed on-the-fly)
Separately settable are:
-
A collinearity test factor (defaults to 1e-8)
Warning: The collinearity factor is the square of the sine of the angle between a column in the normal equations, and the hyper-plane through all the other columns; In special cases (e;g; fitting a polynomial in a very narrow bandwidth window, it could be advisable to set this factor to zero if you want a solution (whatever the truth of it maybe);
-
A Levenberg-Marquardt adjustment factor (if appropriate, defaults to 1e-3)
-
Create the normal equations used in solving the set of condition equations of the user, by using the makeNorm() methods. Separate makenorm() methods are provided for sparse condition equations (e.g. if data for 3 antennas are provided, rather than for all 64)
-
If there are user provided constraints, either limiting constraints like the sum of the angles in a triangle is 180 degrees, or constraints to add missing information if e.g. only differences between parameters have been measured, they can be added to the normal equations with the setConstraint() or the addConstraint() methods. Lagrange multipliers will be used to solve the extended normal equations.
-
The normal equations are triangu;arised (using the collinearity factor as a check for solvability) with the invert() method. If the normal equations are non-solvable an error is returned, or a switch to an SVD solution is made if indicated in the invert call.
-
The solutions and adjustment errors are obtained with the solve() method. A non-linear loop in a Levenberg-Marquardt adjustment can be obtained (together with convergence information), with the solveLoop() method (see below) replacing the combination of invert and solve.
-
Non-linear loops are done by looping through the data using makeNorm() calls, and upgrade the solution with the solveLoop() method. The normal equations are upgraded by changing LM factors. Upgrade depends on the 'balanced' factor. The LM factor is either added in some way to all diagonal elements (if balanced) or all diagonal elements are multiplied by (1+factor) After each loop convergence can be tested by the isReady() call; which will return False or a non-zero code indicating the ready reason. Reasons for stopping can be:
-
SOLINCREMENT: the relative change in the norm of the parameter solutions is less than (a settable,
setEpsValue(), default 1e-8) value.
-
DERIVLEVEL: the inf-norm of the known vector of the equations to be solved is less than the settable,
setEpsDerivative(), default 1e-8, value.
-
MAXITER: maximum number of iterations reached (only possible if a number is explicitly set)
-
NOREDUCTION: if the Levenberg-Marquardt correction factor goes towards infinity. I.e. if no Chi2 improvement seems possible. Have to redo the solution with a different start condition for the unknowns.
-
SINGULAR: can only happen due to numeric rounding, since the LM equations are always positive-definite. Best solution is to indicate SVD needed in the
solveLoop call, which is cost-free
-
Covariance information in various forms can be obtained with the
getCovariance(), getErrors(), getChi() (or getChi2), getSD and getWeightedSD methods after a solve() or after the final loop in a non-linear solution (of course, when necessary only).
An LSQFit object can be re-used by issuing the reset() command, or set() of new values. If an unknown has not been used in the condition equations at all, the doDiagonal() will make sure a proper solution is obtained, with missing unknowns zeroed.
Most of the calculations are done in place; however, enough data is saved that it is possible to continue with the same (partial) normal equations after e.g. an interim solution.
If the normal equations are produced in separate partial sets (e.g. in a multi-processor environment) a merge() method can combine them.
Tip: It is suggested to add any possible constraint equations after the merge;
A debugIt() method provides read access to all internal information.
The member definitions are split over three files. The second one contains the templated member function definitions, to bypass the problem of duplicate definitions of non-templated members when pre-compiling them. The third contains methods for saving objects as Records or through AipsIO.
Warning: No boundary checks on input and output containers is done for faster execution; In general these tests should be done at the higher level routines, like the LinearFit and NonLinearFit classes which should be checked for usage of LSQFit;
The contents can be saved in a record (toRecord), and an object can be created from a record (fromRecord). The record identifier is 'lfit'.
The object can also be saved or restored using AipsIO.
Example
See the tLSQFit.cc and tLSQaips.cc program for extensive examples.
The following example will first create 2 condition equations for 3 unknowns (the third is degenerate). It will first create normal equations for a 2 unknown solution and solve; then it will create normal equations for a 3 unknown solution, and solve (note that the degenerate will be set to 0. The last one will use SVD and one condition equation.r
#include <iostream>
int main() {
Double ce[2][3] = {{1, 1, 0}, {1, -1, 0}};
Double m[2] = {2, 4};
for (
uInt i=0; i<2; i++) fit.makeNorm(ce[i], 1.0, m[i]);
ok = fit.invert(rank);
cout << "ok? " << ok << "; rank: " << rank << endl;
if (ok) {
fit.solve(sol, &sd, &mu);
for (
uInt i=0; i<2; i++) cout <<
"Sol" << i <<
": " << sol[i] << endl;
cout << "sd: "<< sd << "; mu: " << mu << endl;
};
cout << "----------" << endl;
for (
uInt i=0; i<2; i++) fit.makeNorm(ce[i], 1.0, m[i]);
ok = fit.invert(rank);
cout << "ok? " << ok << "; rank: " << rank << endl;
if (ok) {
fit.solve(sol, &sd, &mu);
for (
uInt i=0; i<3; i++) cout <<
"Sol" << i <<
": " << sol[i] << endl;
cout << "sd: "<< sd << "; mu: " << mu << endl;
};
cout << "----------" << endl;
fit.reset();
for (
uInt i=0; i<1; i++) fit.makeNorm(ce[i], 1.0, m[i]);
ok = fit.invert(rank, True);
cout << "ok? " << ok << "; rank: " << rank << endl;
if (ok) {
fit.solve(sol, &sd, &mu);
for (
uInt i=0; i<3; i++) cout <<
"Sol" << i <<
": " << sol[i] << endl;
cout << "sd: "<< sd << "; mu: " << mu << endl;
};
cout << "----------" << endl;
fit.reset();
for (
uInt i=0; i<1; i++) fit.makeNorm(ce[i], 1.0, m[i]);
ok = fit.invert(rank);
cout << "ok? " << ok << "; rank: " << rank << endl;
if (ok) {
fit.solve(sol, &sd, &mu);
for (
uInt i=0; i<3; i++) cout <<
"Sol" << i <<
": " << sol[i] << endl;
cout << "sd: "<< sd << "; mu: " << mu << endl;
};
cout << "----------" << endl;
exit(0);
}
Which will produce the output:
ok? 1; rank: 2
Sol0: 3
Sol1: -1
sd: 0; mu: 0
----------
ok? 1; rank: 3
Sol0: 3
Sol1: -1
Sol2: 0
sd: 0; mu: 0
----------
ok? 1; rank: 2
Sol0: 1
Sol1: 1
Sol2: 0
sd: 0; mu: 0
----------
ok? 0; rank: 2
----------
Motivation
The class was written to be able to do complex, real standard and SVD solutions in a simple and fast way.
To Do
-
a thorough check if all loops are optimal in the makeNorm() methods
-
input of condition equations with cross covariance
Definition at line 330 of file LSQFit.h.
1.8.5