casacore
Loading...
Searching...
No Matches
LSQFit.h
Go to the documentation of this file.
1//# LSQFit.h: Basic class for least squares fitting
2//# Copyright (C) 1999-2001,2004-2008
3//# Associated Universities, Inc. Washington DC, USA.
4//#
5//# This library is free software; you can redistribute it and/or modify it
6//# under the terms of the GNU Library General Public License as published by
7//# the Free Software Foundation; either version 2 of the License, or (at your
8//# option) any later version.
9//#
10//# This library is distributed in the hope that it will be useful, but WITHOUT
11//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13//# License for more details.
14//#
15//# You should have received a copy of the GNU Library General Public License
16//# along with this library; if not, write to the Free Software Foundation,
17//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18//#
19//# Correspondence concerning AIPS++ should be addressed as follows:
20//# Internet email: aips2-request@nrao.edu.
21//# Postal address: AIPS++ Project Office
22//# National Radio Astronomy Observatory
23//# 520 Edgemont Road
24//# Charlottesville, VA 22903-2475 USA
25//#
26//# $Id$
27
28#ifndef SCIMATH_LSQFIT_H
29#define SCIMATH_LSQFIT_H
30
31//# Includes
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>
36#include <complex>
37#include <string>
38#include <utility>
39#include <vector>
40
41namespace casacore { //# NAMESPACE CASACORE - BEGIN
42
43//# Forward Declarations
44
45// <summary> Basic class for the least squares fitting </summary>
46// <reviewed reviewer="Neil Killeen" date="2000/06/01" tests="tLSQFit"
47// demos="">
48// </reviewed>
49
50// <prerequisite>
51// <li> Some knowledge of Matrix operations
52// <li> The background information provided in
53// <a href="../notes/224.html">Note 224</a>.
54// <li> <linkto module="Fitting">Fitting module</linkto>
55// </prerequisite>
56//
57// <etymology>
58// From Least SQuares and Fitting
59// </etymology>
60//
61// <synopsis>
62// The LSQFit class contains the basic functions to do all the fitting
63// described in the
64// <a href="../notes/224.html">Note</a>
65// about fitting.
66// It handles real, and complex equations;<br>
67// linear and non-linear (Levenberg-Marquardt) solutions;<br>
68// regular (with optional constraints) or Singular Value Decomposition
69// (<src>SVD</src>).<br>
70// In essence they are a set of routines to generate normal equations
71// (<src>makeNorm()</src>) in triangular form from a set of condition
72// equations;<br>
73// to do a Cholesky-type decomposition of the normal
74// equations (either regular or <src>SVD</src>) and test its rank
75// (<src>invert()</src>);<br>
76// to do a quasi inversion of the decomposed equations (<src>solve()</src>) to
77// obtain the solution and/or the errors.
78//
79// All calculations are done in place.
80// Methods to obtain additional information about the fitting process are
81// available.
82//
83// This class can be used as a stand-alone class outside of the Casacore
84// environment. In that case the aips.h include file
85// can be replaced if necessary by appropriate typedefs for Double, Float and
86// uInt.<br>
87// The interface to the methods have standard data or standard STL iterator
88// arguments only. They can be used with any container having an STL
89// random-access iterator interface. Especially they can be used with
90// <src>carrays</src> (necessary templates provided),
91// Casacore Vectors (necessary templates
92// provided in <src>LSQaips</src>),
93// standard random access STL containers (like <src>std::vector</src>).
94//
95// The normal operation of the class consists of the following steps:
96// <ul>
97// <li> Create an LSQFit object.
98// The information that can be provided in the constructor of the object,
99// either directly, or indirectly using the <src>set()</src> commands, is
100// (see <a href="../notes/224.html">Note 224</a>):
101// <ul>
102// <li> The number of unknowns that have to be solved for (mandatory)
103// <li> The number of constraint equations you want to use explicitly
104// (defaults to 0, but can be changed on-the-fly)
105// </ul>
106// Separately settable are:
107// <ul>
108// <li> A collinearity test factor (defaults to 1e-8)
109// <note role=warning>
110// The collinearity factor is the square of the sine of the angle between
111// a column in the normal equations, and the hyper-plane through
112// all the other columns. In special cases (e.g. fitting a polynomial
113// in a very narrow bandwidth window, it could be advisable to set this
114// factor to zero if you want a solution (whatever the truth of it maybe).
115// </note>
116// <li> A Levenberg-Marquardt adjustment factor (if appropriate,
117// defaults to 1e-3)
118// </ul>
119//
120// <li>Create the normal equations used in solving the set of condition
121// equations of the user, by using the <src>makeNorm()</src> methods.
122// Separate <src>makenorm()</src> methods are provided for sparse condition
123// equations (e.g. if data for 3 antennas are provided, rather than for all 64)
124//
125// <li>If there are user provided constraints, either limiting constraints like
126// the sum of the angles in a triangle is 180 degrees, or constraints to add
127// missing information if e.g. only differences between parameters have been
128// measured, they can be added to the normal
129// equations with the <src>setConstraint()</src> or
130// the <src>addConstraint()</src> methods. Lagrange multipliers will be used to
131// solve the extended normal equations.
132//
133// <li>The normal equations are triangu;arised (using the collinearity factor
134// as a check for solvability) with the <src>invert()</src> method. If the
135// normal equations are non-solvable an error is returned, or a switch to
136// an SVD solution is made if indicated in the <src>invert</src> call.
137//
138// <li>The solutions and adjustment errors are obtained with the
139// <src>solve()</src> method.
140// A non-linear loop in a Levenberg-Marquardt adjustment can be obtained
141// (together with convergence information), with the <src>solveLoop()</src>
142// method (see below) replacing the combination of
143// <src>invert</src> and <src>solve</src>.
144//
145// <li>Non-linear loops are done by looping through the data using
146// <src>makeNorm()</src> calls, and upgrade the solution with the
147// <src>solveLoop()</src> method.
148// The normal equations are upgraded by changing LM factors. Upgrade depends
149// on the 'balanced' factor. The LM factor is either added in some way to all
150// diagonal elements (if balanced) or all diagonal elements are multiplied by
151// <src>(1+factor)</src> After each loop convergence can be tested
152// by the <src>isReady()</src> call; which will return <src>False</src> or
153// a non-zero code indicating the ready reason. Reasons for stopping can be:
154// <ul>
155// <li> SOLINCREMENT: the relative change in the norm of the parameter
156// solutions is less than
157// (a settable, <src>setEpsValue()</src>, default 1e-8) value.
158// <li> DERIVLEVEL: the inf-norm of the known vector of the equations to be
159// solved is less than the settable, <src>setEpsDerivative()</src>, default
160// 1e-8, value.
161// <li> MAXITER: maximum number of iterations reached (only possible if a
162// number is explicitly set)
163// <li> NOREDUCTION: if the Levenberg-Marquardt correction factor goes towards
164// infinity. I.e. if no Chi2 improvement seems possible. Have to redo the
165// solution with a different start condition for the unknowns.
166// <li> SINGULAR: can only happen due to numeric rounding, since the LM
167// equations are always positive-definite. Best solution is to indicate SVD
168// needed in the <src>solveLoop</src> call, which is cost-free
169// </ul>
170//
171// <li>Covariance information in various forms can be obtained with the
172// <src>getCovariance(), getErrors()</src>, <src>getChi()</src>
173// (or <src>getChi2</src>), <src>getSD</src> and <src>getWeightedSD</src>
174// methods after a <src>solve()</src> or after the final loop in a non-linear
175// solution (of course, when necessary only).
176// </ul>
177//
178// An LSQFit object can be re-used by issuing the <src>reset()</src> command,
179// or <src>set()</src> of new
180// values. If an unknown has not been used in the condition equations at all,
181// the <src>doDiagonal()</src> will make sure a proper solution is obtained,
182// with missing unknowns zeroed.
183//
184// Most of the calculations are done in place; however, enough data is saved
185// that it is possible to continue
186// with the same (partial) normal equations after e.g. an interim solution.
187//
188// If the normal equations are produced in separate partial sets (e.g.
189// in a multi-processor environment) a <src>merge()</src> method can combine
190// them.
191// <note role=tip>
192// It is suggested to add any possible constraint equations after the merge.
193// </note>
194//
195// A <src>debugIt()</src> method provides read access to all internal
196// information.
197//
198// The member definitions are split over three files. The second
199// one contains the templated member function definitions, to bypass the
200// problem of duplicate definitions of non-templated members when
201// pre-compiling them. The third contains methods for saving objects as
202// Records or through AipsIO.
203//
204// <note role=warning> No boundary checks on input and output containers
205// is done for faster execution. In general these tests should be done at
206// the higher level routines, like the
207// <linkto class=LinearFit>LinearFit</linkto> and
208// <linkto class=NonLinearFitLM>NonLinearFit</linkto> classes which should be
209// checked for usage of LSQFit.
210// </note>
211//
212// The contents can be saved in a record (<src>toRecord</src>),
213// and an object can be created from a record (<src>fromRecord</src>).
214// The record identifier is 'lfit'.
215// <br>The object can also be saved or restored using AipsIO.
216// </synopsis>
217//
218// <example>
219// See the tLSQFit.cc and tLSQaips.cc program for extensive examples.
220//
221// The following example will first create 2 condition equations for
222// 3 unknowns (the third is degenerate). It will first create normal equations
223// for a 2 unknown solution and solve; then it will create normal equations
224// for a 3 unknown solution, and solve (note that the degenerate will be
225// set to 0. The last one will use SVD and one condition equation.r
226// <srcblock>
227// #include <casacore/casa/aips.h>
228// #include <casacore/scimath/Fitting/LSQFit.h>
229// #include <iostream>
230//
231// int main() {
232// // Condition equations for x+y=2; x-y=4;
233// Double ce[2][3] = {{1, 1, 0}, {1, -1, 0}};
234// Double m[2] = {2, 4};
235// // Solution and error area
236// Double sol[3];
237// Double sd, mu;
238// uInt rank;
239// Bool ok;
240//
241// // LSQ object
242// LSQFit fit(2);
243//
244// // Make normal equation
245// for (uInt i=0; i<2; i++) fit.makeNorm(ce[i], 1.0, m[i]);
246// // Invert(decompose) and show
247// ok = fit.invert(rank);
248// cout << "ok? " << ok << "; rank: " << rank << endl;
249// // Solve and show
250// if (ok) {
251// fit.solve(sol, &sd, &mu);
252// for (uInt i=0; i<2; i++) cout << "Sol" << i << ": " << sol[i] << endl;
253// cout << "sd: "<< sd << "; mu: " << mu << endl;
254// };
255// cout << "----------" << endl;
256//
257// // Retry with 3 unknowns: note auto fill of unmentioned one
258// fit.set(uInt(3));
259// for (uInt i=0; i<2; i++) fit.makeNorm(ce[i], 1.0, m[i]);
260// ok = fit.invert(rank);
261// cout << "ok? " << ok << "; rank: " << rank << endl;
262// if (ok) {
263// fit.solve(sol, &sd, &mu);
264// for (uInt i=0; i<3; i++) cout << "Sol" << i << ": " << sol[i] << endl;
265// cout << "sd: "<< sd << "; mu: " << mu << endl;
266// };
267// cout << "----------" << endl;
268//
269// // Retry with 3 unknowns; but 1 condition equation and use SVD
270// fit.reset();
271// for (uInt i=0; i<1; i++) fit.makeNorm(ce[i], 1.0, m[i]);
272// ok = fit.invert(rank, True);
273// cout << "ok? " << ok << "; rank: " << rank << endl;
274// if (ok) {
275// fit.solve(sol, &sd, &mu);
276// for (uInt i=0; i<3; i++) cout << "Sol" << i << ": " << sol[i] << endl;
277// cout << "sd: "<< sd << "; mu: " << mu << endl;
278// };
279// cout << "----------" << endl;
280//
281// // Without SVD it would be:
282// fit.reset();
283// for (uInt i=0; i<1; i++) fit.makeNorm(ce[i], 1.0, m[i]);
284// ok = fit.invert(rank);
285// cout << "ok? " << ok << "; rank: " << rank << endl;
286// if (ok) {
287// fit.solve(sol, &sd, &mu);
288// for (uInt i=0; i<3; i++) cout << "Sol" << i << ": " << sol[i] << endl;
289// cout << "sd: "<< sd << "; mu: " << mu << endl;
290// };
291// cout << "----------" << endl;
292//
293// exit(0);
294// }
295// </srcblock>
296// Which will produce the output:
297// <srcblock>
298// ok? 1; rank: 2
299// Sol0: 3
300// Sol1: -1
301// sd: 0; mu: 0
302// ----------
303// ok? 1; rank: 3
304// Sol0: 3
305// Sol1: -1
306// Sol2: 0
307// sd: 0; mu: 0
308// ----------
309// ok? 1; rank: 2
310// Sol0: 1
311// Sol1: 1
312// Sol2: 0
313// sd: 0; mu: 0
314// ----------
315// ok? 0; rank: 2
316// ----------
317// </srcblock>
318// </example>
319//
320// <motivation>
321// The class was written to be able to do complex, real standard and SVD
322// solutions in a simple and fast way.
323// </motivation>
324//
325// <todo asof="2006/04/02">
326// <li> a thorough check if all loops are optimal in the makeNorm() methods
327// <li> input of condition equations with cross covariance
328// </todo>
329
330class LSQFit {
331 public:
332 // Simple classes to overload templated memberfunctions
333 struct Real { enum normType { REAL }; };
334 struct Complex { enum normType { COMPLEX }; };
335 struct Separable { enum normType { SEPARABLE }; };
336 struct AsReal { enum normType { ASREAL }; };
337 struct Conjugate { enum normType { CONJUGATE }; };
338 // And values to use
339 static Real REAL;
344
345 //# Public enums
346 // State of the non-linear solution
356 // Offset of fields in error_p data area.
358 // Number of condition equations
360 // Sum weights of condition equations
362 // Sum known terms squared
364 // Calculated chi^2
366 // Number of error fields
368 };
369 //# Constructors
370 // Construct an object with the number of unknowns and
371 // constraints, using the default collinearity factor and the
372 // default Levenberg-Marquardt adjustment factor.
373 // <group>
374 // Assume real
376 // Allow explicit Real specification
378 // Allow explicit Complex specification
380 // </group>
381 // Default constructor (empty, only usable after a <src>set(nUnknowns)</src>)
383 // Copy constructor (deep copy)
384 LSQFit(const LSQFit &other);
385 // Assignment (deep copy)
386 LSQFit &operator=(const LSQFit &other);
387
388 //# Destructor
390
391 //# Operators
392
393 //# General Member Functions
394 // Triangularize the normal equations and determine
395 // the rank <src>nRank</src> of the normal equations and, in the case of
396 // an <src>SVD</src> solution, the constraint
397 // equations. The collinearity factor is used
398 // to determine if the system can be solved (in essence it is the square
399 // of the sine of the angle between a column in the normal equations and
400 // the plane suspended by the other columns: if too
401 // parallel, the equations are degenerate).
402 // If <src>doSVD</src> is given as False, False is returned if rank not
403 // maximal, else an <src>SVD</src> solution is done.
404 Bool invert(uInt &nRank, Bool doSVD=False);
405 // Copy date from beg to end; converting if necessary to complex data
406 // <group>
407 template <class U>
408 void copy(const Double *beg, const Double *end, U &sol, LSQReal);
409 template <class U>
410 void copy(const Double *beg, const Double *end, U &sol, LSQComplex);
411 template <class U>
412 void copy(const Double *beg, const Double *end, U *sol, LSQReal);
413 template <class U>
414 void copy(const Double *beg, const Double *end, U *sol, LSQComplex);
415 template <class U>
416 void uncopy(Double *beg, const Double *end, U &sol, LSQReal);
417 template <class U>
418 void uncopy(Double *beg, const Double *end, U &sol, LSQComplex);
419 template <class U>
420 void uncopy(Double *beg, const Double *end, U *sol, LSQReal);
421 template <class U>
422 void uncopy(Double *beg, const Double *end, U *sol, LSQComplex);
423 template <class U>
425 template <class U>
427 // </group>
428 // Solve normal equations.
429 // The solution will be given in <src>sol</src>.
430 // <group>
431 template <class U>
432 void solve(U *sol);
433 template <class U>
434 void solve(std::complex<U> *sol);
435 template <class U>
436 void solve(U &sol);
437 // </group>
438 // Solve a loop in a non-linear set.
439 // The methods with the <src>fit</src> argument are deprecated. Use
440 // the combination without the 'fit' parameter, and the <src>isReady()</src>
441 // call. The 'fit' parameter returns
442 // for each loop a goodness
443 // of fit indicator. If it is >0; more loops are necessary.
444 // If it is negative,
445 // and has an absolute value of say less than .001, it is probably ok, and
446 // the iterations can be stopped.
447 // Other arguments are as for <src>solve()</src> and <src>invert()</src>.
448 // The <src>sol</src> is used for both input (parameter guess) and output.
449 // <group>
450 template <class U>
452 U *sol,
453 Bool doSVD=False);
454 template <class U>
456 std::complex<U> *sol,
457 Bool doSVD=False);
458 template <class U>
460 U &sol,
461 Bool doSVD=False);
462 template <class U>
463 Bool solveLoop(Double &fit, uInt &nRank,
464 U *sol,
465 Bool doSVD=False);
466 template <class U>
467 Bool solveLoop(Double &fit, uInt &nRank,
468 std::complex<U> *sol,
469 Bool doSVD=False);
470 template <class U>
471 Bool solveLoop(Double &fit, uInt &nRank,
472 U &sol,
473 Bool doSVD=False);
474 // </group>
475 // Make normal equations using the <src>cEq</src> condition equation (cArray)
476 // (with <src>nUnknowns</src> elements) and a weight <src>weight</src>,
477 // given the known observed value <src>obs</src>.
478 //
479 // <src>doNorm</src> and <src>doKnown</src> can be used
480 // to e.g. re-use existing normal equations, i.e. the condition equations,
481 // but make a new known side (i.e. new observations).
482 //
483 // The versions with <src>cEqIndex[]</src> indicate which of the
484 // <src>nUnknowns</src> are actually present in the condition equation
485 // (starting indexing at 0); the other terms are supposed to be zero. E.g.
486 // if a 12-telescope array has an equation only using telescopes 2 and 4,
487 // the lengths of <src>cEqIndex</src> and <src>cEq</src> will be both 2,
488 // and the index will contain 1 and 3 (when telescope numbering starts at 1)
489 // or 2 and 4 (when telescope numbering starts at 0. The index is given
490 // as an iterator (and hence can be a raw pointer)
491 //
492 // The complex versions can have different interpretation of the inputs,
493 // where the complex number can be seen either as a complex number; as two
494 // real numbers, or as coefficients of equations with complex conjugates.
495 // See <a href="../notes/224.html">Note 224</a>)
496 // for the details.
497 //
498 // Versions with <em>pair</em> assume that the pairs are created by the
499 // <em>SparseDiff</em> automatic differentiation class. The pair is an index
500 // and a value. The indices are assumed to be sorted.
501 //
502 // Special (<em>makeNormSorted</em>) indexed versions exist which assume
503 // that the given indices are sorted (which is the case for the
504 // LOFAR BBS environment).
505 //
506 // Some versions exist with two sets of equations (<em>cEq2, obs2</em>).
507 // If two simultaneous equations are created they will be faster.
508 //
509 // Note that the
510 // use of <src>const U &</src> is due to a Float->Double conversion problem
511 // on Solaris. Linux was ok.
512 // <group>
513 template <class U, class V>
514 void makeNorm(const V &cEq, const U &weight, const U &obs,
515 Bool doNorm=True, Bool doKnown=True);
516 template <class U, class V>
517 void makeNorm(const V &cEq, const U &weight, const U &obs,
519 Bool doNorm=True, Bool doKnown=True);
520 template <class U, class V>
521 void makeNorm(const V &cEq, const U &weight,
522 const std::complex<U> &obs,
523 Bool doNorm=True, Bool doKnown=True);
524 template <class U, class V>
525 void makeNorm(const V &cEq, const U &weight,
526 const std::complex<U> &obs,
528 Bool doNorm=True, Bool doKnown=True);
529 template <class U, class V>
530 void makeNorm(const V &cEq, const U &weight,
531 const std::complex<U> &obs,
533 Bool doNorm=True, Bool doKnown=True);
534 template <class U, class V>
535 void makeNorm(const V &cEq, const U &weight,
536 const std::complex<U> &obs,
538 Bool doNorm=True, Bool doKnown=True);
539 template <class U, class V>
540 void makeNorm(const V &cEq, const U &weight,
541 const std::complex<U> &obs,
543 Bool doNorm=True, Bool doKnown=True);
544 //
545 template <class U, class V, class W>
546 void makeNorm(uInt nIndex, const W &cEqIndex,
547 const V &cEq, const U &weight, const U &obs,
548 Bool doNorm=True, Bool doKnown=True);
549 template <class U, class V, class W>
550 void makeNorm(uInt nIndex, const W &cEqIndex,
551 const V &cEq, const V &cEq2,
552 const U &weight, const U &obs, const U &obs2,
553 Bool doNorm=True, Bool doKnown=True);
554 template <class U, class V, class W>
555 void makeNorm(uInt nIndex, const W &cEqIndex,
556 const V &cEq, const U &weight, const U &obs,
558 Bool doNorm=True, Bool doKnown=True);
559 template <class U, class V, class W>
560 void makeNorm(uInt nIndex, const W &cEqIndex,
561 const V &cEq, const U &weight,
562 const std::complex<U> &obs,
563 Bool doNorm=True, Bool doKnown=True);
564 template <class U, class V, class W>
565 void makeNorm(uInt nIndex, const W &cEqIndex,
566 const V &cEq, const U &weight,
567 const std::complex<U> &obs,
569 Bool doNorm=True, Bool doKnown=True);
570 template <class U, class V, class W>
571 void makeNorm(uInt nIndex, const W &cEqIndex,
572 const V &cEq, const U &weight,
573 const std::complex<U> &obs,
575 Bool doNorm=True, Bool doKnown=True);
576 template <class U, class V, class W>
577 void makeNorm(uInt nIndex, const W &cEqIndex,
578 const V &cEq, const U &weight,
579 const std::complex<U> &obs,
581 Bool doNorm=True, Bool doKnown=True);
582 template <class U, class V, class W>
583 void makeNorm(uInt nIndex, const W &cEqIndex,
584 const V &cEq, const U &weight,
585 const std::complex<U> &obs,
587 Bool doNorm=True, Bool doKnown=True);
588 //
589 template <class U, class V>
590 void makeNorm(const std::vector<std::pair<uInt, V> > &cEq,
591 const U &weight, const U &obs,
592 Bool doNorm=True, Bool doKnown=True);
593 template <class U, class V>
594 void makeNorm(const std::vector<std::pair<uInt, V> > &cEq,
595 const U &weight, const U &obs,
597 Bool doNorm=True, Bool doKnown=True);
598 template <class U, class V>
599 void makeNorm(const std::vector<std::pair<uInt, V> > &cEq,
600 const U &weight,
601 const std::complex<U> &obs,
602 Bool doNorm=True, Bool doKnown=True);
603 template <class U, class V>
604 void makeNorm(const std::vector<std::pair<uInt, V> > &cEq,
605 const U &weight,
606 const std::complex<U> &obs,
608 Bool doNorm=True, Bool doKnown=True);
609 template <class U, class V>
610 void makeNorm(const std::vector<std::pair<uInt, V> > &cEq,
611 const U &weight,
612 const std::complex<U> &obs,
614 Bool doNorm=True, Bool doKnown=True);
615 template <class U, class V>
616 void makeNorm(const std::vector<std::pair<uInt, V> > &cEq,
617 const U &weight,
618 const std::complex<U> &obs,
620 Bool doNorm=True, Bool doKnown=True);
621 template <class U, class V>
622 void makeNorm(const std::vector<std::pair<uInt, V> > &cEq,
623 const U &weight,
624 const std::complex<U> &obs,
626 Bool doNorm=True, Bool doKnown=True);
627 //
628 template <class U, class V, class W>
629 void makeNormSorted(uInt nIndex, const W &cEqIndex,
630 const V &cEq, const U &weight,
631 const U &obs,
632 Bool doNorm=True, Bool doKnown=True);
633 template <class U, class V, class W>
634 void makeNormSorted(uInt nIndex, const W &cEqIndex,
635 const V &cEq, const V &cEq2, const U &weight,
636 const U &obs, const U &obs2,
637 Bool doNorm=True, Bool doKnown=True);
638 // </group>
639 // Get the <src>n-th</src> (from 0 to the rank deficiency, or missing rank,
640 // see e.g. <src>getDeficiency()</src>)
641 // constraint equation as determined by <src>invert()</src> in SVD-mode in
642 // <src> cEq[nUnknown]</src>. False returned for illegal n. Note
643 // that nMissing will be equal to the number of unknowns
644 // (<src>nUnknowns</src>, or double that for the complex case) minus the
645 // rank as returned from the <src>invert()</src> method.
646 // <group>
647 template <class U>
648 Bool getConstraint(uInt n, U *cEq) const;
649 template <class U>
650 Bool getConstraint(uInt n, std::complex<U> *cEq) const;
651 template <class U>
652 Bool getConstraint(uInt n, U &cEq) const;
653 // </group>
654 // Add a new constraint equation (updating nConstraints); or set a
655 // numbered constraint equation (0..nConstraints-1). False if illegal
656 // number n. The constraints are equations with <src>nUnknowns</src> terms,
657 // and a constant value. E.g. measuring three angles of a triangle
658 // could lead to equation <src>[1,1,1]</src> with obs as
659 // <src>3.1415</src>. Note that each complex constraint will be
660 // converted into two real constraints (see
661 // <a href="../notes/224.html">Note 224</a>).
662 // <group>
663 template <class U, class V>
664 Bool setConstraint(uInt n, const V &cEq, const U &obs);
665 template <class U, class V>
666 Bool setConstraint(uInt n, const V &cEq,
667 const std::complex<U> &obs);
668 template <class U, class V, class W>
669 Bool setConstraint(uInt n, uInt nIndex, const W &cEqIndex,
670 const V &cEq, const U &obs);
671 template <class U, class V, class W>
672 Bool setConstraint(uInt n, uInt nIndex, const W &cEqIndex,
673 const V &cEq,
674 const std::complex<U> &obs);
675 template <class U, class V>
676 Bool addConstraint(const V &cEq, const U &obs);
677 template <class U, class V>
678 Bool addConstraint(const V &cEq,
679 const std::complex<U> &obs);
680 template <class U, class V, class W>
681 Bool addConstraint(uInt nIndex, const W &cEqIndex,
682 const V &cEq, const U &obs);
683 template <class U, class V, class W>
684 Bool addConstraint(uInt nIndex, const W &cEqIndex,
685 const V &cEq,
686 const std::complex<U> &obs);
687 // </group>
688 // Merge other <src>LSQFit</src> object (i.e. the normal equation and
689 // related information) into <src>this</src>. Both objects must have the
690 // same number of unknowns, and be pure normal equations (i.e. no
691 // <src>invert(), solve(), solveLoop()</src> or statistics calls
692 // should have been made). If merging cannot be done, <src>False</src>
693 // is returned. The index case (the index is an iterator) assumes that
694 // the normal equations to be merged are a sparse subset of the complete
695 // matrix. The index 'vector' specifies which unknowns are present. An index
696 // outside the scope of the final equations will be skipped.
697 // <note role=tip> For highest numerical precision in the case of a larger
698 // number of partial normal equations to be merged, it is best to merge
699 // them in pairs (and repeat).
700 // </note>
701 // <group>
702 Bool merge(const LSQFit &other);
703 Bool merge(const LSQFit &other, uInt nIndex, const uInt *nEqIndex) {
704 return mergeIt(other, nIndex, nEqIndex); }
705 Bool merge(const LSQFit &other, uInt nIndex,
706 const std::vector<uInt> &nEqIndex) {
707 return mergeIt(other, nIndex, &nEqIndex[0]); }
708 template <class W>
709 Bool merge(const LSQFit &other, uInt nIndex, const W &nEqIndex) {
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]); }
713 // </group>
714 // Reset status to empty
715 void reset();
716 // Set new sizes (default is for Real)
717 // <group>
720 set (static_cast<uInt>(nUnknowns), static_cast<uInt>(nConstraints));
721 };
724 };
727 };
730 set (static_cast<uInt>(nUnknowns), LSQComplex(),
731 static_cast<uInt>(nConstraints));
732 };
733 // </group>
734 // Set new factors (collinearity <src>factor</src>, and Levenberg-Marquardt
735 // <src>LMFactor</src>)
736 void set(Double factor=1e-6, Double LMFactor=1e-3);
737 // Set new value solution test
738 void setEpsValue(Double epsval=1e-8) {epsval_p = epsval; };
739 // Set new derivative test
740 void setEpsDerivative(Double epsder=1e-8) {epsder_p = epsder; };
741 // Set maximum number of iterations
742 void setMaxIter(uInt maxiter=0) { maxiter_p = maxiter; };
743 // Get number of iterations done
744 uInt nIterations() const { return (maxiter_p>0 ? maxiter_p-niter_p : 0); };
745 // Set the expected form of the normal equations
746 void setBalanced(Bool balanced=False) { balanced_p = balanced; };
747 // Ask the state of the non-linear solutions
748 // <group>
749 LSQFit::ReadyCode isReady() const { return ready_p; };
750 const std::string &readyText() const;
751 // </group>
752// Get the covariance matrix (of size <src>nUnknowns * nUnknowns</src>)
753 // <group>
754 template <class U>
756 template <class U>
757 Bool getCovariance(std::complex<U> *covar);
758 // </group>
759 // Get main diagonal of covariance function (of size <src>nUnknowns</src>)
760 // <group>
761 template <class U>
763 template <class U>
764 Bool getErrors(std::complex<U> *errors);
765 template <class U>
767 // </group>
768 // Get the number of unknowns
769 uInt nUnknowns() const { return nun_p; };
770 // Get the number of constraints
771 uInt nConstraints() const { return ncon_p; };
772 // Get the rank deficiency <note role=warning>Note that the number is
773 // returned assuming real values. For complex values it has to be halved
774 // </note>
775 uInt getDeficiency() const { return n_p-r_p; };
776 // Get chi^2 (both are identical); the standard deviation (per observation)
777 // and the standard deviation per weight unit.
778 // <group>
779 Double getChi() const;
780 Double getChi2() const { return getChi(); };
781 Double getSD() const;
783 // </group>
784 // Debug:
785 // <ul>
786 // <li> <src>nun = </src> number of unknowns
787 // <li> <src>np = </src> total number of solved unknowns (nun+ncon)
788 // <li> <src>ncon = </src> number of constraint equations
789 // <li> <src>ner = </src> number of elements in chi<sup>2</sup> vector
790 // <li> <src>rank = </src> rank)
791 // <li> <src>nEq = </src> normal equation (nun*nun as triangular matrix)
792 // <li> <src>known = </src> known vector (np)
793 // <li> <src>constr = </src> constraint matrix (ncon*nun)
794 // <li> <src>er = </src> error info vector (ner)
795 // <li> <src>piv = </src> pivot vector (np)
796 // <li> <src>sEq = </src> normal solution equation (np*np triangular)
797 // <li> <src>sol = </src> internal solution vector (np)
798 // <li> <src>prec = </src> collinearity precision
799 // <li> <src>nonlin = </src> current Levenberg factor-1
800 // </ul>
801 // Note that all pointers may be 0.
802 void debugIt(uInt &nun, uInt &np, uInt &ncon, uInt &ner, uInt &rank,
803 Double *&nEq, Double *&known, Double *&constr, Double *&er,
804 uInt *&piv, Double *&sEq, Double *&sol,
805 Double &prec, Double &nonlin) const;
806 //
807 // Create an LSQFit object from a record.
808 // An error message is generated, and False
809 // returned if an invalid record is given. A valid record will return True.
810 // Error messages are postfixed to error.
811 // <group>
813 // </group>
814 // Create a record from an LSQFit object.
815 // The return will be False and an error
816 // message generated only if the object does not contain a valid object.
817 // Error messages are postfixed to error.
818 Bool toRecord(String &error, RecordInterface &out) const;
819 // Get identification of record
820 const String &ident() const;
821 //
822 // Save or restore using AipsIO.
823 // <group>
824 void toAipsIO (AipsIO&) const;
826 // </group>
827 //
828 protected:
829 //# enum
830 // Bits that can be set/referenced
831 enum StateBit {
832 // Inverted matrix present
834 // Triangularised
836 // Non-linear solution
838 // Filler for cxx2html
840 };
841
842 // Record field names
843 // <group>
844 static const String recid;
845 static const String state;
846 static const String nun;
847 static const String ncon;
848 static const String prec;
849 static const String startnon;
850 static const String nonlin;
851 static const String rank;
852 static const String nnc;
853 static const String piv;
854 static const String constr;
855 static const String known;
856 static const String errors;
857 static const String sol;
858 static const String lar;
859 static const String wsol;
860 static const String wcov;
861 static const String nceq;
862 static const String nar;
863 // </group>
864
865 //# Data
866 // Bits set to indicate state
868 // Number of unknowns
870 // Number of constraints
872 // Matrix size (will be n_p = nun_p + ncon_p)
874 // Rank of normal equations (normally n_p)
876 // Collinearity precision
878 // Levenberg start factor
880 // Levenberg current factor
882 // Levenberg step factor
884 // Test value for [incremental] solution in non-linear loop.
885 // The <src>||sol increment||/||sol||</src> is tested
887 // Test value for known vector in non-linear loop.
888 // ||known||<sub>inf</sub> is tested
890 // Indicator for a well balanced normal equation. A balanced equation is
891 // one with similar values in the main diagonal.
893 // Maximum number of iterations for non-linear solution. If a non-zero
894 // maximum number of iterations is set, the value is tested in non-linear
895 // loops
897 // Iteration count for non-linear solution
899 // Indicate the non-linear state. A non-zero code indicates that non-linear
900 // looping is ready.
902
903 // Pivot table (n_p)
905 // Normal equations (triangular nun_p * nun_p)
907 // Current length nceq_p
909 // Normal combined with constraint equations for solutions
910 // (triangular nnc_p*nnc_p)
912 // Known part equations (n_p)
914 // Counts for errors (N_ErrorField)
916 // Constraint equation area (nun_p*ncon_p))
918 // Solution area (n_p)
920 // Save area for non-linear case (size determined internally)
922 // Save area for non-symmetric (i.e. with constraints) (n_p * n_p)
924 // Work areas for interim solutions and covariance
925 // <group>
928 // </group>
929
930 //# Member functions
931 // Get pointer in rectangular array
932 // <group>
933 Double *rowrt(uInt i) const { return &lar_p[n_p*i]; };
934 Double *rowru(uInt i) const { return &lar_p[nun_p*i]; };
935 // </group>
936 // Calculate the real or imag part of <src>x*conj(y)</src>
937 // <group>
938 static Double realMC(const std::complex<Double> &x,
939 const std::complex<Double> &y) {
940 return (x.real()*y.real() + x.imag()*y.imag()); };
941 static Double imagMC(const std::complex<Double> &x,
942 const std::complex<Double> &y) {
943 return (x.imag()*y.real() - x.real()*y.imag()); };
944 static Float realMC(const std::complex<Float> &x,
945 const std::complex<Float> &y) {
946 return (x.real()*y.real() + x.imag()*y.imag()); };
947 static Float imagMC(const std::complex<Float> &x,
948 const std::complex<Float> &y) {
949 return (x.imag()*y.real() - x.real()*y.imag()); };
950 // </group>
951 // Initialise areas
952 void init();
953 // Clear areas
954 void clear();
955 // De-initialise area
956 void deinit();
957 // Solve normal equations
958 void solveIt();
959 // One non-linear LM loop
960 Bool solveItLoop(Double &fit, uInt &nRank, Bool doSVD=False);
961 // Solve missing rank part
962 void solveMR(uInt nin);
963 // Invert rectangular matrix (i.e. when constraints present)
965 // Get the norm of the current solution vector
967 // Get the infinite norm of the known vector
969 // Merge sparse normal equations
970 Bool mergeIt(const LSQFit &other, uInt nIndex, const uInt *nEqIndex);
971 // Save current status (or part)
973 // Restore current status
975 // Copy data. If all False, only the relevant data for non-linear
976 // solution are copied (normal equations, knows and errors).
977 void copy(const LSQFit &other, Bool all=True);
978 // Extend the constraint equation area to the specify number of
979 // equations.
981 // Create the solution equation area nceq_p and fill it.
983 // Get work areas for solutions, covariance
984 // <group>
987 // </group>
988 //
989};
990
991
992} //# NAMESPACE CASACORE - END
993
994#ifndef CASACORE_NO_AUTO_TEMPLATES
995#include <casacore/scimath/Fitting/LSQFit2.tcc>
996#endif //# CASACORE_NO_AUTO_TEMPLATES
997#endif
Type of complex numeric class indicator
Definition LSQTraits.h:71
Double stepfactor_p
Levenberg step factor.
Definition LSQFit.h:883
void solve(U *sol)
Solve normal equations.
Double prec_p
Collinearity precision.
Definition LSQFit.h:877
uInt nUnknowns() const
Get the number of unknowns.
Definition LSQFit.h:769
static const String nceq
Definition LSQFit.h:861
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)
void createNCEQ()
Create the solution equation area nceq_p and fill it.
Bool merge(const LSQFit &other)
Merge other LSQFit object (i.e.
void set(Int nUnknowns, Int nConstraints=0)
Definition LSQFit.h:719
Bool invert(uInt &nRank, Bool doSVD=False)
Triangularize the normal equations and determine the rank nRank of the normal equations and,...
void extendConstraints(uInt n)
Extend the constraint equation area to the specify number of equations.
Bool solveLoop(uInt &nRank, U *sol, Bool doSVD=False)
Solve a loop in a non-linear set.
static const String state
Definition LSQFit.h:845
Bool setConstraint(uInt n, const V &cEq, const std::complex< U > &obs)
void makeNorm(const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const U &obs, LSQFit::Real, Bool doNorm=True, Bool doKnown=True)
static Separable SEPARABLE
Definition LSQFit.h:341
Double * wsol_p
Work areas for interim solutions and covariance.
Definition LSQFit.h:926
Bool getCovariance(std::complex< U > *covar)
static const String recid
Record field names.
Definition LSQFit.h:844
void uncopy(Double *beg, const Double *end, U &sol, LSQComplex)
void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True)
ReadyCode ready_p
Indicate the non-linear state.
Definition LSQFit.h:901
Bool getErrors(std::complex< U > *errors)
void set(uInt nUnknowns, uInt nConstraints=0)
Set new sizes (default is for Real)
Bool getConstraint(uInt n, U *cEq) const
Get the n-th (from 0 to the rank deficiency, or missing rank, see e.g.
Bool solveLoop(Double &fit, uInt &nRank, U &sol, Bool doSVD=False)
const String & ident() const
Get identification of record.
void save(Bool all=True)
Save current status (or part)
void init()
Initialise areas.
uInt nIterations() const
Get number of iterations done.
Definition LSQFit.h:744
void makeNorm(const V &cEq, const U &weight, const std::complex< U > &obs, Bool doNorm=True, Bool doKnown=True)
static const String nonlin
Definition LSQFit.h:850
void set(Double factor=1e-6, Double LMFactor=1e-3)
Set new factors (collinearity factor, and Levenberg-Marquardt LMFactor)
Double startnon_p
Levenberg start factor.
Definition LSQFit.h:879
LSQFit(uInt nUnknowns, const LSQComplex &, uInt nConstraints=0)
Allow explicit Complex specification.
Bool invertRect()
Invert rectangular matrix (i.e.
Double normInfKnown(const Double *known) const
Get the infinite norm of the known vector.
uInt r_p
Rank of normal equations (normally n_p)
Definition LSQFit.h:875
static const String constr
Definition LSQFit.h:854
LSQFit * nar_p
Save area for non-linear case (size determined internally)
Definition LSQFit.h:921
Double normSolution(const Double *sol) const
Get the norm of the current solution vector.
Bool solveLoop(Double &fit, uInt &nRank, std::complex< U > *sol, Bool doSVD=False)
void solveIt()
Solve normal equations.
uInt nnc_p
Current length nceq_p.
Definition LSQFit.h:908
uInt state_p
Bits set to indicate state.
Definition LSQFit.h:867
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)
void makeNorm(const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const std::complex< U > &obs, Bool doNorm=True, Bool doKnown=True)
Double * rowru(uInt i) const
Definition LSQFit.h:934
Double * rowrt(uInt i) const
Get pointer in rectangular array.
Definition LSQFit.h:933
Double getChi() const
Get chi^2 (both are identical); the standard deviation (per observation) and the standard deviation p...
Bool getConstraint(uInt n, U &cEq) const
void makeNorm(const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::AsReal, Bool doNorm=True, Bool doKnown=True)
uInt nun_p
Number of unknowns.
Definition LSQFit.h:869
Double * constr_p
Constraint equation area (nun_p*ncon_p))
Definition LSQFit.h:917
static AsReal ASREAL
Definition LSQFit.h:342
void toAipsIO(AipsIO &) const
Save or restore using AipsIO.
LSQFit(const LSQFit &other)
Copy constructor (deep copy)
uInt n_p
Matrix size (will be n_p = nun_p + ncon_p)
Definition LSQFit.h:873
static Complex COMPLEX
Definition LSQFit.h:340
Bool getErrors(U &errors)
Bool solveLoop(Double &fit, uInt &nRank, U *sol, Bool doSVD=False)
uInt * piv_p
Pivot table (n_p)
Definition LSQFit.h:904
static const String wsol
Definition LSQFit.h:859
Bool addConstraint(const V &cEq, const U &obs)
void uncopy(Double *beg, const Double *end, U *sol, LSQReal)
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)
void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex< U > &obs, Bool doNorm=True, Bool doKnown=True)
Double epsval_p
Test value for [incremental] solution in non-linear loop.
Definition LSQFit.h:886
LSQMatrix * nceq_p
Normal combined with constraint equations for solutions (triangular nnc_p*nnc_p)
Definition LSQFit.h:911
Bool addConstraint(const V &cEq, const std::complex< U > &obs)
Double * wcov_p
Definition LSQFit.h:927
void makeNorm(const std::vector< std::pair< uInt, V > > &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True)
Bool merge(const LSQFit &other, uInt nIndex, const W &nEqIndex)
Definition LSQFit.h:709
Double * lar_p
Save area for non-symmetric (i.e.
Definition LSQFit.h:923
LSQFit(uInt nUnknowns, const LSQReal &, uInt nConstraints=0)
Allow explicit Real specification.
Bool solveLoop(uInt &nRank, std::complex< U > *sol, Bool doSVD=False)
uInt ncon_p
Number of constraints.
Definition LSQFit.h:871
static Double realMC(const std::complex< Double > &x, const std::complex< Double > &y)
Calculate the real or imag part of x*conj(y)
Definition LSQFit.h:938
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)
Bool getConstraint(uInt n, std::complex< U > *cEq) const
Bool addConstraint(uInt nIndex, const W &cEqIndex, const V &cEq, const U &obs)
Bool getCovariance(U *covar)
Get the covariance matrix (of size nUnknowns * nUnknowns)
void solve(U &sol)
static const String known
Definition LSQFit.h:855
Bool setConstraint(uInt n, uInt nIndex, const W &cEqIndex, const V &cEq, const std::complex< U > &obs)
Double getChi2() const
Definition LSQFit.h:780
const std::string & readyText() const
Bool solveItLoop(Double &fit, uInt &nRank, Bool doSVD=False)
One non-linear LM loop.
Bool toRecord(String &error, RecordInterface &out) const
Create a record from an LSQFit object.
Double * error_p
Counts for errors (N_ErrorField)
Definition LSQFit.h:915
void solve(std::complex< U > *sol)
void setMaxIter(uInt maxiter=0)
Set maximum number of iterations.
Definition LSQFit.h:742
Bool getErrors(U *errors)
Get main diagonal of covariance function (of size nUnknowns)
void set(uInt nUnknowns, const LSQReal &, uInt nConstraints=0)
Definition LSQFit.h:722
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)
static const String errors
Definition LSQFit.h:856
uInt nConstraints() const
Get the number of constraints.
Definition LSQFit.h:771
void copy(const Double *beg, const Double *end, U *sol, LSQComplex)
void makeNorm(const V &cEq, const U &weight, const U &obs, LSQFit::Real, Bool doNorm=True, Bool doKnown=True)
void clear()
Clear areas.
static const String ncon
Definition LSQFit.h:847
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...
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)
void makeNorm(const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Separable, Bool doNorm=True, Bool doKnown=True)
void set(Int nUnknowns, const LSQComplex &, Int nConstraints=0)
Definition LSQFit.h:729
void copy(const LSQFit &other, Bool all=True)
Copy data.
static const String rank
Definition LSQFit.h:851
LSQFit()
Default constructor (empty, only usable after a set(nUnknowns))
Double * sol_p
Solution area (n_p)
Definition LSQFit.h:919
Bool setConstraint(uInt n, uInt nIndex, const W &cEqIndex, const V &cEq, const U &obs)
static const String nun
Definition LSQFit.h:846
void setEpsValue(Double epsval=1e-8)
Set new value solution test.
Definition LSQFit.h:738
void getWorkSOL()
Get work areas for solutions, covariance.
static Double imagMC(const std::complex< Double > &x, const std::complex< Double > &y)
Definition LSQFit.h:941
static const String wcov
Definition LSQFit.h:860
LSQFit(uInt nUnknowns, uInt nConstraints=0)
Construct an object with the number of unknowns and constraints, using the default collinearity facto...
Bool merge(const LSQFit &other, uInt nIndex, const std::vector< uInt > &nEqIndex)
Definition LSQFit.h:705
uInt niter_p
Iteration count for non-linear solution.
Definition LSQFit.h:898
void makeNorm(const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Complex, Bool doNorm=True, Bool doKnown=True)
LSQFit::ReadyCode isReady() const
Ask the state of the non-linear solutions.
Definition LSQFit.h:749
Bool setConstraint(uInt n, const V &cEq, const U &obs)
Add a new constraint equation (updating nConstraints); or set a numbered constraint equation (0....
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)
static const String sol
Definition LSQFit.h:857
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)
uInt maxiter_p
Maximum number of iterations for non-linear solution.
Definition LSQFit.h:896
void setBalanced(Bool balanced=False)
Set the expected form of the normal equations.
Definition LSQFit.h:746
void fromAipsIO(AipsIO &)
Bool balanced_p
Indicator for a well balanced normal equation.
Definition LSQFit.h:892
void set(uInt nUnknowns, const LSQComplex &, uInt nConstraints=0)
void restore(Bool all=True)
Restore current status.
StateBit
Bits that can be set/referenced.
Definition LSQFit.h:831
@ NONLIN
Non-linear solution.
Definition LSQFit.h:837
@ TRIANGLE
Triangularised.
Definition LSQFit.h:835
@ N_StateBit
Filler for cxx2html.
Definition LSQFit.h:839
@ INVERTED
Inverted matrix present.
Definition LSQFit.h:833
void uncopy(Double *beg, const Double *end, U &sol, LSQReal)
static const String lar
Definition LSQFit.h:858
ErrorField
Offset of fields in error_p data area.
Definition LSQFit.h:357
@ NC
Number of condition equations.
Definition LSQFit.h:359
@ CHI2
Calculated chi^2.
Definition LSQFit.h:365
@ SUMWEIGHT
Sum weights of condition equations.
Definition LSQFit.h:361
@ N_ErrorField
Number of error fields.
Definition LSQFit.h:367
@ SUMLL
Sum known terms squared.
Definition LSQFit.h:363
Bool merge(const LSQFit &other, uInt nIndex, const uInt *nEqIndex)
Definition LSQFit.h:703
Bool solveLoop(uInt &nRank, U &sol, Bool doSVD=False)
void makeNorm(const V &cEq, const U &weight, const std::complex< U > &obs, LSQFit::Conjugate, Bool doNorm=True, Bool doKnown=True)
Double * known_p
Known part equations (n_p)
Definition LSQFit.h:913
void uncopy(Double *beg, const Double *end, U *sol, LSQComplex)
void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const U &obs, LSQFit::Real, Bool doNorm=True, Bool doKnown=True)
Double nonlin_p
Levenberg current factor.
Definition LSQFit.h:881
static Float imagMC(const std::complex< Float > &x, const std::complex< Float > &y)
Definition LSQFit.h:947
static const String startnon
Definition LSQFit.h:849
uInt getDeficiency() const
Get the rank deficiency Warning: Note that the number is returned assuming real values; For complex ...
Definition LSQFit.h:775
void copyDiagonal(U &errors, LSQComplex)
static const String piv
Definition LSQFit.h:853
static Real REAL
And values to use.
Definition LSQFit.h:339
Bool mergeIt(const LSQFit &other, uInt nIndex, const uInt *nEqIndex)
Merge sparse normal equations.
ReadyCode
State of the non-linear solution.
Definition LSQFit.h:347
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)
void reset()
Reset status to empty.
static const String nnc
Definition LSQFit.h:852
void copy(const Double *beg, const Double *end, U &sol, LSQReal)
Copy date from beg to end; converting if necessary to complex data.
Double getWeightedSD() const
Bool fromRecord(String &error, const RecordInterface &in)
Create an LSQFit object from a record.
static Conjugate CONJUGATE
Definition LSQFit.h:343
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)
static const String prec
Definition LSQFit.h:848
static Float realMC(const std::complex< Float > &x, const std::complex< Float > &y)
Definition LSQFit.h:944
void setEpsDerivative(Double epsder=1e-8)
Set new derivative test.
Definition LSQFit.h:740
void solveMR(uInt nin)
Solve missing rank part.
void set(Int nUnknowns, const LSQReal &, Int nConstraints=0)
Definition LSQFit.h:725
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:
void deinit()
De-initialise area.
LSQMatrix * norm_p
Normal equations (triangular nun_p * nun_p)
Definition LSQFit.h:906
LSQFit & operator=(const LSQFit &other)
Assignment (deep copy)
Bool addConstraint(uInt nIndex, const W &cEqIndex, const V &cEq, const std::complex< U > &obs)
Double getSD() const
void makeNormSorted(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True)
void copy(const Double *beg, const Double *end, U *sol, LSQReal)
Double epsder_p
Test value for known vector in non-linear loop.
Definition LSQFit.h:889
void copy(const Double *beg, const Double *end, U &sol, LSQComplex)
void copyDiagonal(U &errors, LSQReal)
static const String nar
Definition LSQFit.h:862
String: the storage and methods of handling collections of characters.
Definition String.h:225
this file contains all the compiler specific defines
Definition mainpage.dox:28
const Bool False
Definition aipstype.h:44
unsigned int uInt
Definition aipstype.h:51
float Float
Definition aipstype.h:54
int Int
Definition aipstype.h:50
bool Bool
Define the standard types used by Casacore.
Definition aipstype.h:42
const Bool True
Definition aipstype.h:43
double Double
Definition aipstype.h:55
LatticeExprNode all(const LatticeExprNode &expr)
Simple classes to overload templated memberfunctions.
Definition LSQFit.h:333