casacore
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 
41 namespace 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 
330 class 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;
340  static Complex COMPLEX;
342  static AsReal ASREAL;
344 
345  //# Public enums
346  // State of the non-linear solution
347  enum ReadyCode {
355  };
356  // Offset of fields in error_p data area.
357  enum ErrorField {
358  // Number of condition equations
359  NC,
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,
518  LSQFit::Real,
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,
557  LSQFit::Real,
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,
596  LSQFit::Real,
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  };
725  void set(Int nUnknowns, const LSQReal &, Int nConstraints=0) {
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>
755  Bool getCovariance(U *covar);
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>
812  Bool fromRecord(String &error, const RecordInterface &in);
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
833  INVERTED = 1,
834  // Triangularised
836  // Non-linear solution
838  // Filler for cxx2html
839  N_StateBit
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
966  Double normSolution(const Double *sol) const;
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)
972  void save(Bool all=True);
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.
982  void createNCEQ();
983  // Get work areas for solutions, covariance
984  // <group>
985  void getWorkSOL();
986  void getWorkCOV();
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)
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 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)
LSQFit & operator=(const LSQFit &other)
Assignment (deep copy)
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
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
Double * rowru(uInt i) const
Definition: LSQFit.h:934
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)
const std::string & readyText() const
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.
const String & ident() const
Get identification of record.
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.
Double * rowrt(uInt i) const
Get pointer in rectangular array.
Definition: LSQFit.h:933
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
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
const Double e
e and functions thereof:
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