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 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()); 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 aips++
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), aips++ Vectors (necessary templates
provided, standard random access STL containers (like std::vector) (templates
not provided but just a matter of replacing names).
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 (mandatory)
- The number of constraint equations you want to use explicitly (defaults to 0)
- A collinearity test factor (defaults to 1e-8)
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 scarce condition
equations (like data for 3 antennas are provided, rather than 64)
- If there are user provided constraints, they can be added to the normal
equations with the setConstraint() or
the addConstraint() methods.
- The normal equations are decomposed (using the collinearity factor as a
check) in either SVD or standard mode with the invert() method.
- 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.
- Covariance information in various forms can be obtained with the getCovariance(), getErrors(), getChi() (or getChi2), getSD and getWeightedSD methods (of course, when necessary only).
A class 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.
A debugIt() method provides access to all internal information.
The member definitions are split over two 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.
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 <aips/aips.h>
#include <aips/Fitting/LSQ.h>
#include <iostream>
int main() {
// Condition equations for x+y=2; x-y=4;
Double ce[2][3] = {{1, 1, 0}, {1, -1, 0}};
Double m[2] = {2, 4};
// Solution and error area
Double sol[3];
Double sd, mu;
uInt rank;
Bool ok;
// LSQ area
LSQ fit(2);
// Make normal equation
for (uInt i=0; i<2; i++) fit.makeNorm(ce[i], 1.0, m[i]);
// Invert(decompose) and show
ok = fit.invert(rank);
cout << "ok? " << ok << "; rank: " << rank << endl;
// Solve and show
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;
// Retry with 3 unknowns: note auto fill of unmentioned one
fit.set(uInt(3));
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;
// Retry with 3 unknowns; but 1 condition equation and use SVD
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;
// Without SVD it would be:
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 thorogh check if all loops are optimal in the makeNorm() methods
- method to load normal equations (and other information) directly from external source.
- input of condition equations with cross covariance
Member Description
explicit 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.
Assume real
LSQFit(uInt nUnknowns, const LSQReal &, 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.
Allow explicit Real specification
LSQFit(uInt nUnknowns, const LSQComplex &, 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.
Allow explicit Complex specification
LSQFit()
Default constructor (empty, only usable after a set(nUnknowns))LSQFit(const LSQFit &other)
Copy constructor (deep copy)LSQFit &operator=(const LSQFit &other)
Assignment (deep copy)~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. The collinearity factor is used to determine if the system can be solved (in essence it is the square of the sine of the angle between a column in the normal equations and the plane suspended by the other columns: if too parallel, the equations are degenerate). If doSVD is given as False, False is returned if rank not maximal, else an SVD solution is done.
template <class U> void copy(const Double *beg, const Double *end, U &sol, U)
template <class U> void copy(const Double *beg, const Double *end, U &sol, U)
template <class U> void copy(const Double *beg, const Double *end, U *sol, U)
template <class U> void copy(const Double *beg, const Double *end, U *sol, U)
template <class U> void uncopy(Double *beg, const Double *end, U &sol, U)
template <class U> void uncopy(Double *beg, const Double *end, U &sol, U)
template <class U> void uncopy(Double *beg, const Double *end, U *sol, U)
template <class U> void uncopy(Double *beg, const Double *end, U *sol, U)
template <class U> void copyDiagonal(U &errors, U)
template <class U> void copyDiagonal(U &errors, U)
Copy date from beg to end; converting if necessary to complex data
template <class U> void solve(U *sol)
template <class U> void solve(std::complex<U> *sol)
template <class U> void solve(U &sol)
Solve normal equations.
The solution will be given in sol.
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)
Solve a loop in a non-linear set. The fit argument returns
for each loop a goodness
of fit indicator. If it is >0; more loops are necessary.
If it is negative,
and has an absolute value of say less than .001, it is probably ok, and
the iterations can be stopped.
Other arguments are as for solve() and invert().
The sol is used for both input (parameter guess) and output.
This non-linear solution will only work for 1 known.
template <class U, class V> void makeNorm(const V &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True)
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 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)
Make normal equations using the cEq condition equation
(with nUnknowns elements) and a weight weight,
given the known observed value obs (with
nKnowns elements).
doNorm and doKnown can be used to e.g. re-use existing normal equations, i.e. the condition equations, but make a new known side (i.e. new observations).
The versions with cEqIndex[] indicate which of the
nUnknowns are actually present in the condition equation
(starting indexing at 0); the other terms are supposed to be zero. E.g.
if a 12-telescope array has an equation only using telescopes 2 and 4,
the lengths of
Note that the
use of const U & is due to a Float->Double conversion problem
on Solaris. Linux was ok.
Get pointer in rectangular array
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. getDeficiency())
constraint equation as determined by invert() in SVD-mode in
cEq[nUnknown]. False returned for illegal n. Note
that nMissing will be equal to the number of unknowns
(nUnknowns, or double that for the complex case) minus the
rank as returned from the invert() method.
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). False if illegal
number n. The constraints are equations with nUnknowns terms,
and a constant value. E.g. measuring three angles of a triangle
could lead to equation [1,1,1] with obs as
3.1415. Note that each complex constraint will be
converted into two real constraints (see
Note 224).
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)
void reset()
Reset status to empty
void set(uInt nUnknowns, uInt nConstraints=0)
Set new sizes (default is for Real)
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-8, Double LMFactor=1e-3)
Set new factors (collinearity factor, and Levenberg-Marquardt
LMFactor)
template <class U> Bool getCovariance(U *covar)
Get the covariance matrix (of size nUnknowns * nUnknowns)
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)
template <class U> Bool getErrors(std::complex<U> *errors)
template <class U> Bool getErrors(U &errors)
uInt nUnknowns() const
Get the number of unknowns
uInt nConstraints() const
Get the number of constraints
uInt getDeficiency() const
Get the rank deficiency Note that the numberf is
eturned assuming real values. For complex values it has to be halved
Double getChi() const
Get chi^2 (both are identical); the standard deviation (per observation)
and the standard deviation per weight unit.
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:
Note that all pointers may be 0.
enum ErrorField
Offset of fields in error_p data area.
enum StateBit
Bits that can be set/referenced
Double *rowrt(uInt i) const
Double *rowru(uInt i) const
static Double realMC(const std::complex<Double> &x, const std::complex<Double> &y)
Calculate the real or imag part of x*conj(y)
static Double imagMC(const std::complex<Double> &x, const std::complex<Double> &y)
static Float realMC(const std::complex<Float> &x, const std::complex<Float> &y)
static Float imagMC(const std::complex<Float> &x, const std::complex<Float> &y)
void init()
Initialise areas
void clear()
Clear areas
void deinit()
De-initialise area
void solveIt()
Solve normal equations
Bool solveItLoop(Double &fit, uInt &nRank, Bool doSVD=False)
One non-linear LM loop
void solveMR(uInt nin)
Solve missing rank part
Bool invertRect()
Invert rectangular matrix (i.e. when constraints present)
void save(Bool all=True)
Save current status (or part)
void restore(Bool all=True)
Restore current status
void copy(const LSQFit &other, Bool all=True)
Copy data. If all False, only the relevant data for non-linear
solution are copied (normal equations, knows and errors).
void extendConstraints(uInt n)
Extend the constraint equation area to the specify number of
equations.
void createNCEQ()
Create the solution equation area nceq_p and fill it.
void getWorkSOL()
Get work areas for solutions, covariance
void getWorkCOV()