Go to the documentation of this file.
28 #ifndef SCIMATH_LSQFIT_H
29 #define SCIMATH_LSQFIT_H
32 #include <casacore/casa/aips.h>
33 #include <casacore/casa/Utilities/RecordTransformable.h>
34 #include <casacore/scimath/Fitting/LSQMatrix.h>
35 #include <casacore/scimath/Fitting/LSQTraits.h>
456 std::complex<U> *
sol,
468 std::complex<U> *
sol,
513 template <
class U,
class V>
514 void makeNorm(
const V &cEq,
const U &weight,
const U &obs,
516 template <
class U,
class V>
517 void makeNorm(
const V &cEq,
const U &weight,
const U &obs,
520 template <
class U,
class V>
521 void makeNorm(
const V &cEq,
const U &weight,
522 const std::complex<U> &obs,
524 template <
class U,
class V>
525 void makeNorm(
const V &cEq,
const U &weight,
526 const std::complex<U> &obs,
529 template <
class U,
class V>
530 void makeNorm(
const V &cEq,
const U &weight,
531 const std::complex<U> &obs,
534 template <
class U,
class V>
535 void makeNorm(
const V &cEq,
const U &weight,
536 const std::complex<U> &obs,
539 template <
class U,
class V>
540 void makeNorm(
const V &cEq,
const U &weight,
541 const std::complex<U> &obs,
545 template <
class U,
class V,
class W>
547 const V &cEq,
const U &weight,
const U &obs,
549 template <
class U,
class V,
class W>
551 const V &cEq,
const V &cEq2,
552 const U &weight,
const U &obs,
const U &obs2,
554 template <
class U,
class V,
class W>
556 const V &cEq,
const U &weight,
const U &obs,
559 template <
class U,
class V,
class W>
561 const V &cEq,
const U &weight,
562 const std::complex<U> &obs,
564 template <
class U,
class V,
class W>
566 const V &cEq,
const U &weight,
567 const std::complex<U> &obs,
570 template <
class U,
class V,
class W>
572 const V &cEq,
const U &weight,
573 const std::complex<U> &obs,
576 template <
class U,
class V,
class W>
578 const V &cEq,
const U &weight,
579 const std::complex<U> &obs,
582 template <
class U,
class V,
class W>
584 const V &cEq,
const U &weight,
585 const std::complex<U> &obs,
589 template <
class U,
class V>
590 void makeNorm(
const std::vector<std::pair<uInt, V> > &cEq,
591 const U &weight,
const U &obs,
593 template <
class U,
class V>
594 void makeNorm(
const std::vector<std::pair<uInt, V> > &cEq,
595 const U &weight,
const U &obs,
598 template <
class U,
class V>
599 void makeNorm(
const std::vector<std::pair<uInt, V> > &cEq,
601 const std::complex<U> &obs,
603 template <
class U,
class V>
604 void makeNorm(
const std::vector<std::pair<uInt, V> > &cEq,
606 const std::complex<U> &obs,
609 template <
class U,
class V>
610 void makeNorm(
const std::vector<std::pair<uInt, V> > &cEq,
612 const std::complex<U> &obs,
615 template <
class U,
class V>
616 void makeNorm(
const std::vector<std::pair<uInt, V> > &cEq,
618 const std::complex<U> &obs,
621 template <
class U,
class V>
622 void makeNorm(
const std::vector<std::pair<uInt, V> > &cEq,
624 const std::complex<U> &obs,
628 template <
class U,
class V,
class W>
630 const V &cEq,
const U &weight,
633 template <
class U,
class V,
class W>
635 const V &cEq,
const V &cEq2,
const U &weight,
636 const U &obs,
const U &obs2,
663 template <
class U,
class V>
665 template <
class U,
class V>
667 const std::complex<U> &obs);
668 template <
class U,
class V,
class W>
670 const V &cEq,
const U &obs);
671 template <
class U,
class V,
class W>
674 const std::complex<U> &obs);
675 template <
class U,
class V>
677 template <
class U,
class V>
679 const std::complex<U> &obs);
680 template <
class U,
class V,
class W>
682 const V &cEq,
const U &obs);
683 template <
class U,
class V,
class W>
686 const std::complex<U> &obs);
704 return mergeIt(other, nIndex, nEqIndex); }
706 const std::vector<uInt> &nEqIndex) {
707 return mergeIt(other, nIndex, &nEqIndex[0]); }
710 std::vector<uInt> ix(nIndex);
711 for (
uInt i=0; i<nIndex; ++i) ix[i] = nEqIndex[i];
712 return mergeIt(other, nIndex, &ix[0]); }
939 const std::complex<Double> &y) {
940 return (x.real()*y.real() + x.imag()*y.imag()); };
942 const std::complex<Double> &y) {
943 return (x.imag()*y.real() - x.real()*y.imag()); };
945 const std::complex<Float> &y) {
946 return (x.real()*y.real() + x.imag()*y.imag()); };
948 const std::complex<Float> &y) {
949 return (x.imag()*y.real() - x.real()*y.imag()); };
994 #ifndef CASACORE_NO_AUTO_TEMPLATES
995 #include <casacore/scimath/Fitting/LSQFit2.tcc>
996 #endif //# CASACORE_NO_AUTO_TEMPLATES
LSQFit()
Default constructor (empty, only usable after a set(nUnknowns))
Double normSolution(const Double *sol) const
Get the norm of the current solution vector.
LatticeExprNode all(const LatticeExprNode &expr)
void extendConstraints(uInt n)
Extend the constraint equation area to the specify number of equations.
uInt nUnknowns() const
Get the number of unknowns.
static Separable SEPARABLE
Type of complex numeric class indicator
Bool setConstraint(uInt n, const V &cEq, const U &obs)
Add a new constraint equation (updating nConstraints); or set a numbered constraint equation (0....
Number of condition equations.
AipsIO is the object persistency mechanism of Casacore
ReadyCode
State of the non-linear solution.
Bool addConstraint(const V &cEq, const U &obs)
void setMaxIter(uInt maxiter=0)
Set maximum number of iterations.
void toAipsIO(AipsIO &) const
Save or restore using AipsIO.
static const String state
uInt maxiter_p
Maximum number of iterations for non-linear solution.
void restore(Bool all=True)
Restore current status.
Bool merge(const LSQFit &other, uInt nIndex, const uInt *nEqIndex)
Double * lar_p
Save area for non-symmetric (i.e.
Double * constr_p
Constraint equation area (nun_p*ncon_p))
Bool merge(const LSQFit &other, uInt nIndex, const std::vector< uInt > &nEqIndex)
LSQMatrix * norm_p
Normal equations (triangular nun_p * nun_p)
Basic class for the least squares fitting.
Double normInfKnown(const Double *known) const
Get the infinite norm of the known vector.
Double * known_p
Known part equations (n_p)
Double getChi() const
Get chi^2 (both are identical); the standard deviation (per observation) and the standard deviation p...
void set(uInt nUnknowns, const LSQReal &, uInt nConstraints=0)
const String & ident() const
Get identification of record.
StateBit
Bits that can be set/referenced.
static const String errors
void setBalanced(Bool balanced=False)
Set the expected form of the normal equations.
void fromAipsIO(AipsIO &)
void solve(U *sol)
Solve normal equations.
static Conjugate CONJUGATE
Double stepfactor_p
Levenberg step factor.
void set(Int nUnknowns, const LSQComplex &, Int nConstraints=0)
ErrorField
Offset of fields in error_p data area.
void createNCEQ()
Create the solution equation area nceq_p and fill it.
LSQFit * nar_p
Save area for non-linear case (size determined internally)
void save(Bool all=True)
Save current status (or part)
void copy(const Double *beg, const Double *end, U &sol, LSQReal)
Copy date from beg to end; converting if necessary to complex data.
void setEpsValue(Double epsval=1e-8)
Set new value solution test.
void solveMR(uInt nin)
Solve missing rank part.
Sum weights of condition equations.
void init()
Initialise areas.
uInt ncon_p
Number of constraints.
Double startnon_p
Levenberg start factor.
Bool invert(uInt &nRank, Bool doSVD=False)
Triangularize the normal equations and determine the rank nRank of the normal equations and,...
static Float imagMC(const std::complex< Float > &x, const std::complex< Float > &y)
uInt nIterations() const
Get number of iterations done.
LSQFit & operator=(const LSQFit &other)
Assignment (deep copy)
void uncopy(Double *beg, const Double *end, U &sol, LSQReal)
Double * rowrt(uInt i) const
Get pointer in rectangular array.
uInt niter_p
Iteration count for non-linear solution.
Bool getCovariance(U *covar)
Get the covariance matrix (of size nUnknowns * nUnknowns)
Bool merge(const LSQFit &other, uInt nIndex, const W &nEqIndex)
Simple classes to overload templated memberfunctions.
static const String recid
Record field names.
uInt nConstraints() const
Get the number of constraints.
static const String constr
Bool getConstraint(uInt n, U *cEq) const
Get the n-th (from 0 to the rank deficiency, or missing rank, see e.g.
Double epsval_p
Test value for [incremental] solution in non-linear loop.
static const String nonlin
Double * wsol_p
Work areas for interim solutions and covariance.
Double prec_p
Collinearity precision.
void copyDiagonal(U &errors, LSQReal)
Bool mergeIt(const LSQFit &other, uInt nIndex, const uInt *nEqIndex)
Merge sparse normal equations.
uInt getDeficiency() const
Get the rank deficiency Warning: Note that the number is returned assuming real values; For complex ...
Bool invertRect()
Invert rectangular matrix (i.e.
static const String startnon
Bool getErrors(U *errors)
Get main diagonal of covariance function (of size nUnknowns)
this file contains all the compiler specific defines
Double * error_p
Counts for errors (N_ErrorField)
Double epsder_p
Test value for known vector in non-linear loop.
Support class for the LSQ package.
Bool toRecord(String &error, RecordInterface &out) const
Create a record from an LSQFit object.
static const String known
void set(Int nUnknowns, const LSQReal &, Int nConstraints=0)
uInt n_p
Matrix size (will be n_p = nun_p + ncon_p)
const Double e
e and functions thereof:
Abstract base class for Record classes.
Bool balanced_p
Indicator for a well balanced normal equation.
uInt * piv_p
Pivot table (n_p)
ReadyCode ready_p
Indicate the non-linear state.
uInt r_p
Rank of normal equations (normally n_p)
Double * rowru(uInt i) const
Bool solveItLoop(Double &fit, uInt &nRank, Bool doSVD=False)
One non-linear LM loop.
void deinit()
De-initialise area.
String: the storage and methods of handling collections of characters.
Double nonlin_p
Levenberg current factor.
static Real REAL
And values to use.
Bool solveLoop(uInt &nRank, U *sol, Bool doSVD=False)
Solve a loop in a non-linear set.
LSQMatrix * nceq_p
Normal combined with constraint equations for solutions (triangular nnc_p*nnc_p)
bool Bool
Define the standard types used by Casacore.
void set(Int nUnknowns, Int nConstraints=0)
void getWorkSOL()
Get work areas for solutions, covariance.
void setEpsDerivative(Double epsder=1e-8)
Set new derivative test.
Bool fromRecord(String &error, const RecordInterface &in)
Create an LSQFit object from a record.
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:
static Float realMC(const std::complex< Float > &x, const std::complex< Float > &y)
void solveIt()
Solve normal equations.
static Double realMC(const std::complex< Double > &x, const std::complex< Double > &y)
Calculate the real or imag part of x*conj(y)
Bool merge(const LSQFit &other)
Merge other LSQFit object (i.e.
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 weigh...
const std::string & readyText() const
Double * sol_p
Solution area (n_p)
static Double imagMC(const std::complex< Double > &x, const std::complex< Double > &y)
void makeNormSorted(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True)
uInt state_p
Bits set to indicate state.
void reset()
Reset status to empty.
Typing support classes for LSQ classes.
uInt nun_p
Number of unknowns.
void set(uInt nUnknowns, uInt nConstraints=0)
Set new sizes (default is for Real)
Double getWeightedSD() const
LSQFit::ReadyCode isReady() const
Ask the state of the non-linear solutions.
uInt nnc_p
Current length nceq_p.