casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
LSQMatrix.h
Go to the documentation of this file.
00001 //# LSQMatrix.h: Support class for the LSQ package
00002 //# Copyright (C) 2004,2005,2006
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //# $Id: LSQMatrix.h 20652 2009-07-06 05:04:32Z Malte.Marquarding $
00027 
00028 #ifndef SCIMATH_LSQMATRIX_H
00029 #define SCIMATH_LSQMATRIX_H
00030 
00031 //# Includes
00032 #include <casa/aips.h>
00033 #include <algorithm>
00034 #include <casa/Utilities/RecordTransformable.h>
00035 
00036 namespace casa { //# NAMESPACE CASA - BEGIN
00037 
00038 //# Forward Declarations
00039 class AipsIO;
00040 
00041 // <summary> Support class for the LSQ package </summary>
00042 // <reviewed reviewer="Wim Brouw" date="2004/03/20" tests="tLSQFit"
00043 //       demos="">
00044 // </reviewed>
00045 
00046 // <prerequisite>
00047 //   <li> Some knowledge of Matrix operations
00048 // </prerequisite>
00049 //
00050 // <etymology>
00051 // From Least SQuares and Matrix
00052 // </etymology>
00053 //
00054 // <synopsis>
00055 // The LSQMatrix class contains the handling of the basic container used
00056 // in the <linkto class="LSQFit">LSQFit</linkto> class and its derivatives.
00057 // This basic container is a triangular matrix.
00058 //
00059 // The basic operations provided are referencing and indexing of cells,
00060 // rows, columns and diagonal of the triangular matrix.
00061 // The class is a private structure, with explicit friends.
00062 //
00063 // The class contains a number of public methods (with _pub in name) that
00064 // can be used anywhere, and which perform index range checking.
00065 //
00066 // The contents can be saved in a record (<src>toRecord</src>), 
00067 // and an object can be created from a record (<src>fromRecord</src>).
00068 // The record identifier is 'tmat'.
00069 // </synopsis>
00070 //
00071 // <example>
00072 // See the <linkto class="LSQFit">LSQFit</linkto> class for its use.
00073 // </example>
00074 //
00075 // <motivation>
00076 // The class was written to isolate the handling of the normal equations
00077 // used in the <src>LSQ</src> classes.
00078 // </motivation>
00079 //
00080 // <todo asof="2004/03/20">
00081 //      <li> Look in possibility of an STL iterator along row, column and
00082 //         diagonal
00083 // </todo>
00084 
00085 class LSQMatrix : public RecordTransformable {
00086   //# Friends
00087   friend class LSQFit;
00088 
00089  public:
00090   // A set of public interface functions. Checks for index ranges are made, 
00091   // and zero or null returned if error.
00092   // <group>
00093   // Get row pointer in normal equation (points to element <src>[i][0]</src>)
00094   Double *row_pub(uInt i) const { return (i<n_p) ? row(i) : 0; };
00095   // Get next row or previous row pointer in normal equation if the pointer
00096   // <src>row</src> is at row <src>i</src>.
00097   // <group>
00098   void incRow_pub(Double *&row, uInt i) const { if (i<n_p-1) incRow(row,i); };
00099   void decRow_pub(Double *&row, uInt i) const { if (i>0) decRow(row,i); };
00100   // </group>
00101   // Get diagonal element pointer <src>[i][i]</src>
00102   Double *diag_pub(uInt i) const { return ((i<n_p) ? diag(i) : 0); };
00103   // Get length of triangular array
00104   uInt nelements_pub() const { return (len_p); };
00105   // Get number of rows
00106   uInt nrows_pub() const { return n_p; };
00107   // Make diagonal element 1 if zero (Note that this is always called when
00108   // <src>invert()</src> is called). Only n-length sub-matrix is done.
00109   void doDiagonal_pub(uInt n) { if (n<n_p) doDiagonal(n); }
00110   // Multiply n-length of diagonal with <src>1+fac</src>
00111   void mulDiagonal_pub(uInt n, Double fac) { if (n<n_p) mulDiagonal(n,fac); };
00112   // Add <src>fac</src> to n-length of diagonal
00113   void addDiagonal_pub(uInt n, Double fac) { if (n<n_p) addDiagonal(n,fac); };
00114   // Determine max of abs values of n-length of diagonal
00115   Double maxDiagonal_pub(uInt n){ return ((n<n_p) ? maxDiagonal(n) : 0); };
00116   // </group>
00117 
00118  private:
00119   //# Constructors
00120   // Default constructor (empty, only usable after a <src>set(n)</src>)
00121   LSQMatrix();
00122   // Construct an object with the number of rows and columns indicated.
00123   // If a <src>Bool</src> argument is present, the number
00124   // will be taken as double the number given (assumes complex). 
00125   // <group>
00126   explicit LSQMatrix(uInt n);
00127   LSQMatrix(uInt n, Bool);
00128   // </group>
00129   // Copy constructor (deep copy)
00130   LSQMatrix(const LSQMatrix &other);
00131   // Assignment (deep copy)
00132   LSQMatrix &operator=(const LSQMatrix &other);
00133 
00134   //# Destructor
00135   ~LSQMatrix();
00136 
00137   //# Operators
00138   // Index an element in the triangularised matrix
00139   // <group>
00140   Double &operator[](uInt index) { return (trian_p[index]); };
00141   Double operator[](uInt index) const { return (trian_p[index]); };
00142   // </group>
00143 
00144   //# General Member Functions
00145   // Reset all data to zero
00146   void reset() { clear(); };
00147   // Set new sizes (default is for Real, a Bool argument will make it complex)
00148   // <group>
00149   void set(uInt n);
00150   void set(uInt n, Bool);
00151   // </group>
00152   // Get row pointer in normal equation (points to element <src>[i][0]</src>)
00153   Double *row(uInt i) const { return &trian_p[((n2m1_p-i)*i)/2]; };
00154   // Get next row or previous row pointer in normal equation if the pointer
00155   // <src>row</src> is at row <src>i</src>.
00156   // <group>
00157   void incRow(Double *&row, uInt i) const { row += nm1_p-i; };
00158   void decRow(Double *&row, uInt i) const { row -= n_p-i; };
00159   // </group>
00160   // Get diagonal element pointer <src>[i][i]</src>
00161   Double *diag(uInt i) const { return &trian_p[((n2p1_p-i)*i)/2]; };
00162   // Get length of triangular array
00163   uInt nelements() const { return (len_p); };
00164   // Get number of rows
00165   uInt nrows() const { return n_p; };
00166   // Copy data.
00167   void copy(const LSQMatrix &other);
00168   // Initialise matrix
00169   void init();
00170   // Clear matrix
00171   void clear();
00172   // De-initialise matrix
00173   void deinit();
00174   // Make diagonal element 1 if zero (Note that this is always called when
00175   // <src>invert()</src> is called). Only n-length sub-matrix is done.
00176   void doDiagonal(uInt n);
00177   // Multiply n-length of diagonal with <src>1+fac</src>
00178   void mulDiagonal(uInt n, Double fac);
00179   // Add <src>fac</src> to n-length of diagonal
00180   void addDiagonal(uInt n, Double fac);
00181   // Determine max of abs values of n-length of diagonal
00182   Double maxDiagonal(uInt n);
00183   // Create a Matrix from a record. An error message is generated, and False
00184   // returned if an invalid record is given. A valid record will return True.
00185   // Error messages are postfixed to error.
00186   // <group>
00187   Bool fromRecord(String &error, const RecordInterface &in);
00188   // </group>
00189   // Create a record from an LSQMatrix. The return will be False and an error
00190   // message generated only if the object does not contain a valid Matrix.
00191   // Error messages are postfixed to error.
00192   Bool toRecord(String &error, RecordInterface &out) const;
00193   // Get identification of record
00194   const String &ident() const;
00195   // Convert a <src>carray</src> to/from a record. Field only written if 
00196   // non-zero length. No carray created if field does not exist on input.
00197   // False returned if unexpectedly no data available for non-zero length
00198   // (put), or a field has zero length vector(get).
00199   // <group>
00200   static Bool putCArray(String &error, RecordInterface &out,
00201                         const String &fname,
00202                         uInt len, const Double * const in);
00203   static Bool getCArray(String &error, const RecordInterface &in,
00204                         const String &fname,
00205                         uInt len, Double *&out);
00206   static Bool putCArray(String &error, RecordInterface &out,
00207                         const String &fname,
00208                         uInt len, const uInt * const in);
00209   static Bool getCArray(String &error, const RecordInterface &in,
00210                         const String &fname,
00211                         uInt len, uInt *&out);
00212   // </group>
00213 
00214   // Save or restore using AipsIO.
00215   void fromAipsIO (AipsIO& in);
00216   void toAipsIO (AipsIO& out) const;
00217   static void putCArray (AipsIO& out, uInt len, const Double* const in);
00218   static void getCArray (AipsIO& in, uInt len, Double*& out);
00219   static void putCArray (AipsIO& out, uInt len, const uInt* const in);
00220   static void getCArray (AipsIO& in, uInt len, uInt*& out);
00221 
00222   //# Data
00223   // Matrix size (linear size)
00224   uInt n_p;
00225   // Derived sizes (all 0 if n_p equals 0)
00226   // <group>
00227   // Total size
00228   uInt len_p;
00229   // <src>n-1</src>
00230   uInt nm1_p;
00231   // <src>2n-1</src>
00232   Int n2m1_p;
00233   // <src>2n+1</src>
00234   Int n2p1_p;
00235   // </group>
00236   // Matrix (triangular n_p * n_p)
00237   Double *trian_p;
00238   // Record field names
00239   static const String tmatsiz;
00240   static const String tmatdat;
00241   // <group>
00242   // </group>
00243   //
00244 };
00245 
00246 
00247 } //# NAMESPACE CASA - END
00248 
00249 #endif