- REAL = 0
- Real solution (default)
- COMPLEX = 1
- Normal complex variables and condition equations
- SEPARABLE = 3
- Separable complex. Solve for n complex variables, but have condition equations with 2n terms. Alternating for the real and imaginary part of the unknowns
- ASREAL = 7
- Separable as if complex are just two reals. Will give the same result as having 2n real variables to solve for, with independent condition equations for the two sets of n unknowns.
- CONJUGATE = 11
- Separable with consecutive unknowns are each other's conjugate. Solve for n complex unknowns, with 2n condition equation arguments, alternatively valid for x0, conj(x0), x1, conj(x1)...
The only difference with the LSQBase class is the set of methods makeNorm(), solve(), solveLoop(), which can accept complex values.
The constructors can accept a complex code, indicating the type of complex solution wanted (see Note 224).
#include <aips/aips.h>
#include <aips/Fitting/LSQ.h>
#include <aips/Mathematics/Complex.h>
#include <iostream>
int main() {
// Condition equations for x+y=2,3i; x-y=4,1i;
DComplex ce[2][3] = {{DComplex(1,0), DComplex(1,0), DComplex(0,0)},
{DComplex(1,0), DComplex(-1,0), DComplex(0,0)}};
DComplex m[2] = {DComplex(2,3), DComplex(4,1)};
// Solution and error area
DComplex sol[3];
Double sd, mu;
uInt rank;
Bool ok;
// LSQ area
LSQ fit(2, LSQ::COMPLEX);
// Make normal equation
for (uInt i=0; i<2; i++) fit.makeNorm(ce[i], 1.0, m[i]);
// Invert 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 ASREAL type
fit.set(LSQ::ASREAL);
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;
// Retry with CONJUGATE type: note # of unknowns!
fit.set(LSQ::CONJUGATE);
fit.set(1);
m[0] = DComplex(2,0); m[1] = DComplex(0,1);
for (uInt i=0; i<2; 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<1; i++) cout << "Sol" << i << ": " << sol[i] << endl;
cout << "sd: "<< sd << "; mu: " << mu << endl;
};
cout << "----------" << endl;
// Do a real non-linear one for solving x^2 and y^2
{
// Condition equations
Double ce[4][2] = {{1, 1}, {1, -1}, {1, 2}, {3,2}};
Double m[4] = {13, -5, 22.5, 29.5};
Double wce[2], wm;
// Solution and error area (and guess)
Double sol[2] = {1.6, 3.3};
Double sd, mu;
Double covar[2][2];
// How is fit?
Double howfit(1);
uInt rank;
Bool ok;
// Number of iterations allowed
uInt nIter(20);
cout << "Iter " << nIter << ": sol: (" << sol[0] << ", " <<
sol[1] << ")";
cout << endl;
// LSQ area
LSQ fit(2, LSQ::REAL);
// Iterate till fitting param ok or iterations used up
while (nIter > 0 && (howfit > 0 || howfit < -0.001)) {
// Make normal equation
for (uInt i=0; i<4; i++) {
// Adjust
wm = m[i];
for (uInt j=0; j<2; j++) {
wce[j] = 2*ce[i][j]*sol[j];
wm -= ce[i][j]*sol[j]*sol[j];
};
fit.makeNorm(wce, 1.0, &wm);
};
ok = fit.solveLoop(howfit, rank, sol, &sd, &mu);
cout << "Iter " << --nIter << ": sol: (" <<
sol[0] << ", " << sol[1] << ")";
cout << "; how: " << howfit << ", rank: " << rank << endl;
};
// Show covariance -- but first to solve last time:
// Make normal equation
for (uInt i=0; i<4; i++) {
// Adjust
wm = m[i];
for (uInt j=0; j<2; j++) {
wce[j] = 2*ce[i][j]*sol[j];
wm -= ce[i][j]*sol[j]*sol[j];
};
fit.makeNorm(wce, 1.0, &wm);
};
ok = fit.invert(rank);
fit.solve(sol, &sd, &mu);
cout << "Iter " << --nIter << ": sol: (" <<
sol[0] << ", " << sol[1] << ")";
cout << endl;
fit.getCovariance(covar[0]);
for (uInt i=0; i<2; i++) {
for (uInt j=0; j<2; j++) {
cout << "Covar(" << i << ", " << j << ") = " <<
covar[i][j] << endl;
};
};
}
exit(0);
}
Which will produce:
ok? 1; rank: 4 Sol0: (3, 2) Sol1: (-1, 1) sd: 0; mu: 0 ---------- ok? 1; rank: 4 Sol0: (3, 0) Sol1: (-1, 0) sd: 3.16228; mu: 3.16228 ---------- ok? 1; rank: 2 Sol0: (1, 0.5) sd: 0; mu: 0 ---------- Iter 20: sol: (1.6, 3.3) Iter 19: sol: (1.99258, 3.03617); how: 1, rank: 2 Iter 18: sol: (1.95521, 3.02374); how: -1.78727, rank: 2 Iter 17: sol: (1.95485, 3.02372); how: -0.877943, rank: 2 Iter 16: sol: (1.95485, 3.02372); how: -8.59124e-05, rank: 2 Iter 15: sol: (-9.54194e-14, 5.1996e-14) Covar(0, 0) = 0.0116822 Covar(0, 1) = -0.0060421 Covar(1, 0) = -0.0060421 Covar(1, 1) = 0.00585938
Construct an object with the number of unknowns, knowns and constraints, and optionally the type, specified, using the default collinearity factor and the default Levenberg-Marquardt adjustment factor
Assume real
Construct an object with the number of unknowns, knowns and constraints, and optionally the type, specified, using the default collinearity factor and the default Levenberg-Marquardt adjustment factor
Allow complex/real specification
Default constructor (empty, Complex, only usable after a set(nUnknowns))
Copy constructor (deep copy)
Assignment (deep copy)
Solve normal equations.
The solution will be given in sol, with the
adjustment error mu, and the standard
deviation
Bool solveLoop(Double &fit, uInt &nRank, Double sol[], Double sd[], Double mu[], 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().
This non-linear solution wil only work for 1 known.
Bool solveLoop(Double &fit, uInt &nRank, Float sol[], Double sd[], Double mu[], Bool doSVD=False)
Bool solveLoop(Double &fit, uInt &nRank, Complex sol[], Double sd[], Double mu[], Bool doSVD=False)
Bool solveLoop(Double &fit, uInt &nRank, DComplex sol[], Double sd[], Double mu[], Bool doSVD=False)
void makeNorm(const Double cEq[], Double weight, const Double obs[], Bool doNorm=True, Bool doKnown=True)
Make normal equations using the cEq condition equations
(with nUnknowns elements) and a weight weight,
given the known observed values obs (with
nKnowns elements).
doNorm and doKnown can be used
to e.g. re-use existing normal equations, but make a new known side.
void makeNorm(const Double cEq[], Double weight, Double obs, Bool doNorm=True, Bool doKnown=True)
void makeNorm(const Float cEq[], Float weight, const Float obs[], Bool doNorm=True, Bool doKnown=True)
void makeNorm(const Float cEq[], Float weight, Float obs, Bool doNorm=True, Bool doKnown=True)
void makeNorm(const DComplex cEq[], Double weight, const DComplex obs[], Bool doNorm=True, Bool doKnown=True)
void makeNorm(const DComplex cEq[], Double weight, const DComplex &obs, Bool doNorm=True, Bool doKnown=True)
void makeNorm(const Complex cEq[], Float weight, const Complex obs[], Bool doNorm=True, Bool doKnown=True)
void makeNorm(const Complex cEq[], Float weight, const Complex &obs, Bool doNorm=True, Bool doKnown=True)
void makeConstraint(const Double cEq[])
Set constraint equations cEq (with
[nUnknowns][nConstraints] elements).
void makeConstraint(const Float cEq[])
void makeConstraint(const DComplex cEq[])
void makeConstraint(const Complex cEq[])
Bool getCovariance(Double *covar)
Get the covariance matrix (of size nUnknowns * nUnknowns)
Bool getCovariance(Float *covar)
Bool getCovariance(DComplex *covar)
Bool getCovariance(Complex *covar)
Bool getErrors(Double *errors)
Get main diagonal of covariance function (of size nUnknowns)
Bool getErrors(Float *errors)
Bool getErrors(Complex *errors)
Bool getErrors(DComplex *errors)