LSQBase.h

Classes

LSQBase -- Basic class for the least squares fitting (full description)

class LSQBase

Types

enum normType

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)...

enum ErrorField

NC

Number of condition equations

SUMWEIGHT

Sum weights of condition equations

SUMLL

Sum known terms squared

CHI2

Calculated chi^2

N_ErrorField

Number of error fields

enum StateBit

INVERTED = 1

Inverted matrix present

TRIANGLE = 2*INVERTED,

Triangularised

NONLIN = 2*TRIANGLE,

Non-linear solution Non-linear solution

N_StateBit

Filler for cxx2html Filler for cxx2html

Interface

Public Members

LSQBase(uInt nUnknowns, uInt nKnowns=1, uInt nConstraints=0)

LSQBase(uInt nUnknowns, LSQBase::normType type, uInt nKnowns=1, uInt nConstraints=0)

LSQBase()

LSQBase(const LSQBase &other)

LSQBase &operator=(const LSQBase &other)

~LSQBase()

Bool invert(uInt &nRank, Bool doSVD=False)

void solve(Double sol[], Double sd[], Double mu[])

void solve(Float sol[], Double sd[], Double mu[])

Bool solveLoop(Double &fit, uInt &nRank, Double sol[], Double sd[], Double mu[], Bool doSVD=False)

Bool solveLoop(Double &fit, uInt &nRank, Float sol[], Double sd[], Double mu[], Bool doSVD=False)

void doDiagonal()

void makeNorm(const Double cEq[], Double weight, const Double obs[], Bool doNorm=True, Bool doKnown=True)

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 getConstraint(uInt &nMissing, Double *cEq) const

void makeConstraint(const Double cEq[])

void makeConstraint(const Float cEq[])

void reset()

void set(LSQBase::normType type)

void set(uInt nUnknowns, uInt nKnowns=1, uInt nConstraints=0)

void set(Int nUnknowns, Int nKnowns=1, Int nConstraints=0)

void set(Double factor=1e-8, Double LMFactor=1e-3)

Bool getCovariance(Double *covar)

Bool getCovariance(Float *covar)

Bool getErrors(Double *errors)

Bool getErrors(Float *errors)

Double getChi(uInt m=0) const

void debugIt(Double *&nEq, Double *&er, Double *&known, uInt *&piv, uInt &rank) const

Protected Members

Double *rowne(uInt i) const

Double *rowrt(uInt i) const

Double *rowru(uInt i) const

uInt lentri() const

static void swap(Double &x, Double &y)

static void swap(uInt &x, uInt &y)

void init()

void clear()

void deinit()

void solveMR(uInt nin)

Bool invertRect()

void save(Bool all=True)

void restore(Bool all=True)

void copy(const LSQBase &other, Bool all=True)

void mulDiagonal(Double fac)

void getWorkCEQ()

void getWorkSOL()

void getWorkCOV()

Description

Review Status

Reviewed By:

Neil Killeen

Date Reviewed:

2000/06/01

Programs:

Demos:

• dLSQBase

Tests:

• tLSQ

Prerequisite

• Some knowledge of Matrix operations
• The background information provided in Note 224
• Fitting module

Etymology

Synopsis

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 LSQ class is derived from this class, and contains the Complex routines (and can be used outside the aips++ environment). As soon as aips++ understands the standard complex classes, the LSQBase and LSQ classes will be combined into the one LSQ class. Hence the examples given below are including LSQ.h, and creating LSQ objects.
The FitLSQ class is derived from LSQ, and provides an interface using the (aips++) Array classes.
All three types of objects are identical (LSQBase, LSQ, FitLSQ). They only differ in the allowable interface to generate the normal equations, and the interface to obtain information (like solution) from them.

The normal operation of the class consists of the following steps:

• Create an LSQ (or LSQBase; or FitLSQ) 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 knowns that have to be solved against simultaneously (defaults to 1)
  • The number of constraint equations you want to use explicitly (defaults to 0)
  Separately settable are:
  • 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.

• If there are user provided constraints, they can be added to the normal equations with the makeConstraint() method.

• 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() and getChi() 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, the copy constructor and assignment copy all information, and make it possible to continue with the same (partial) normal equations after e.g. an interim solution.

A debugIt() method will provide access to all internal information.

Example

      #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);
      }
    

   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

To Do

• method to load normal equations (and other information) directly from external source.
• auto scaling of input (imagine the result of a 6-degree polynomial as function of frequency in Hz, and with a standard deviation of 1e-6 Jy per equation).
• indirect addressing for sparse condition equations
• templated in/out partial interface
• input of condition equations with cross covariance

Member Description

enum normType

LSQBase(uInt nUnknowns, uInt nKnowns=1, uInt nConstraints=0)

Assume real

LSQBase(uInt nUnknowns, LSQBase::normType type, uInt nKnowns=1, uInt nConstraints=0)

Allow complex/real specification

LSQBase()

LSQBase(const LSQBase &other)

LSQBase &operator=(const LSQBase &other)

~LSQBase()

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.

void solve(Double sol[], Double sd[], Double mu[])
void solve(Float sol[], Double sd[], Double mu[])