00001 //# Matrix.h: A 2-D Specialization of the Array Class 00002 //# Copyright (C) 1993,1994,1995,1996,1999,2000,2001,2003 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$ 00027 00028 #ifndef CASA_MATRIX_H 00029 #define CASA_MATRIX_H 00030 00031 00032 //# Includes 00033 #include <casacore/casa/aips.h> 00034 #include <casacore/casa/Arrays/Array.h> 00035 00036 namespace casacore { //#Begin casa namespace 00037 00038 //# Forward Declarations 00039 template<class T> class Vector; 00040 00041 00042 // <summary> A 2-D Specialization of the Array class </summary> 00043 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="" demos=""> 00044 // </reviewed> 00045 // 00046 // Matrix objects are two-dimensional specializations (e.g., more convenient 00047 // and efficient indexing) of the general Array class. You might also want 00048 // to look at the Array documentation to see inherited functionality. A 00049 // tutorial on using the array classes in general is available in the 00050 // "AIPS++ Programming Manual". 00051 // 00052 // Generally the member functions of Array are also available in 00053 // Matrix versions which take a pair of integers where the array 00054 // needs an IPosition. Since the Matrix 00055 // is two-dimensional, the IPositions are overkill, although you may 00056 // use those versions if you want to. 00057 // <srcblock> 00058 // Matrix<Int> mi(100,100); // Shape is 100x100 00059 // mi.resize(50,50); // Shape now 50x50 00060 // </srcblock> 00061 // 00062 // Slices may be taken with the Slice class. To take a slice, one "indexes" 00063 // with one Slice(start, length, inc) for each axis, 00064 // where end and inc are optional. 00065 // Additionally, there are row(), column() and diagonal() 00066 // member functions which return Vector's which refer to the storage back 00067 // in the Matrix: 00068 // <srcblock> 00069 // Matrix<Float> mf(100, 100); 00070 // mf.diagonal() = 1; 00071 // </srcblock> 00072 // 00073 // Correct indexing order of a matrix is: 00074 // <srcblock> 00075 // Matrix<Int> mi(n1,n2) // [nrow, ncolumn] 00076 // for (uInt j=0; j<mi.ncolumn(); j++) { 00077 // for (uInt i=0; i<mi.nrow(); i++) { 00078 // mi(i,j) = i*j; 00079 // } 00080 // } 00081 // </srcblock> 00082 // 00083 // 00084 // Element-by-element arithmetic and logical operations are available (in 00085 // aips/ArrayMath.h and aips/ArrayLogical.h). Other Matrix operations (e.g. 00086 // LU decomposition) are available, and more appear periodically. 00087 // 00088 // As with the Arrays, if the preprocessor symbol AIPS_DEBUG is 00089 // defined at compile time invariants will be checked on entry to most 00090 // member functions. Additionally, if AIPS_ARRAY_INDEX_CHECK is defined 00091 // index operations will be bounds-checked. Neither of these should 00092 // be defined for production code. 00093 00094 template<class T> class Matrix : public Array<T> 00095 { 00096 public: 00097 // A Matrix of length zero in each dimension; zero origin. 00098 Matrix(); 00099 00100 // A Matrix with "l1" rows and "l2" columns. 00101 Matrix(size_t l1, size_t l2); 00102 00103 // A Matrix with "l1" rows and "l2" columns. 00104 Matrix(size_t l1, size_t l2, ArrayInitPolicy initPolicy); 00105 00106 // A Matrix with "l1" rows and "l2" columns. 00107 // Fill it with the initial value. 00108 Matrix(size_t l1, size_t l2, const T &initialValue); 00109 00110 // A matrix of shape with shape "len". 00111 Matrix(const IPosition &len); 00112 00113 // A matrix of shape with shape "len". 00114 Matrix(const IPosition &len, ArrayInitPolicy initPolicy); 00115 00116 // A matrix of shape with shape "len". 00117 // Fill it with the initial value. 00118 Matrix(const IPosition &len, const T &initialValue); 00119 00120 // The copy constructor uses reference semantics. 00121 Matrix(const Matrix<T> &other); 00122 00123 // Construct a Matrix by reference from "other". "other must have 00124 // ndim() of 2 or less. 00125 Matrix(const Array<T> &other); 00126 00127 // Create an Matrix of a given shape from a pointer. 00128 Matrix(const IPosition &shape, T *storage, StorageInitPolicy policy = COPY); 00129 // Create an Matrix of a given shape from a pointer. 00130 Matrix(const IPosition &shape, T *storage, StorageInitPolicy policy, AbstractAllocator<T> const &allocator); 00131 // Create an Matrix of a given shape from a pointer. Because the pointer 00132 // is const, a copy is always made. 00133 Matrix(const IPosition &shape, const T *storage); 00134 00135 // Define a destructor, otherwise the (SUN) compiler makes a static one. 00136 virtual ~Matrix(); 00137 00138 // Create an identity matrix of side length n. (Could not do this as a constructor 00139 // because of ambiguities with other constructors). 00140 static Matrix<T> identity (size_t n); 00141 00142 // Assign the other array (which must be dimension 2) to this matrix. 00143 // If the shapes mismatch, this array is resized. 00144 virtual void assign (const Array<T>& other); 00145 00146 // Make this matrix a reference to other. Other must be of dimensionality 00147 // 2 or less. 00148 virtual void reference(const Array<T> &other); 00149 00150 // Resize to the given shape (must be 2-dimensional). 00151 // Resize without argument is equal to resize(0,0). 00152 // <group> 00153 using Array<T>::resize; 00154 void resize(size_t nx, size_t ny, Bool copyValues=False) { 00155 Matrix<T>::resize(nx, ny, copyValues, Array<T>::defaultArrayInitPolicy()); 00156 } 00157 void resize(size_t nx, size_t ny, Bool copyValues, ArrayInitPolicy policy); 00158 virtual void resize(); 00159 virtual void resize(const IPosition &newShape, Bool copyValues, ArrayInitPolicy policy); 00160 // </group> 00161 00162 // Copy the values from other to this Matrix. If this matrix has zero 00163 // elements then it will resize to be the same shape as other; otherwise 00164 // other must conform to this. 00165 // Note that the assign function can be used to assign a 00166 // non-conforming matrix. 00167 // <group> 00168 Matrix<T> &operator=(const Matrix<T> &other); 00169 virtual Array<T> &operator=(const Array<T> &other); 00170 // </group> 00171 00172 // Copy val into every element of this Matrix; i.e. behaves as if 00173 // val were a constant conformant matrix. 00174 Array<T> &operator=(const T &val) 00175 { return Array<T>::operator=(val); } 00176 00177 // Copy to this those values in marray whose corresponding elements 00178 // in marray's mask are True. 00179 Matrix<T> &operator= (const MaskedArray<T> &marray) 00180 { Array<T> (*this) = marray; return *this; } 00181 00182 00183 // Single-pixel addressing. If AIPS_ARRAY_INDEX_CHECK is defined, 00184 // bounds checking is performed. 00185 // <group> 00186 T &operator()(const IPosition &i) 00187 { return Array<T>::operator()(i); } 00188 const T &operator()(const IPosition &i) const 00189 { return Array<T>::operator()(i); } 00190 T &operator()(size_t i1, size_t i2) 00191 { 00192 #if defined(AIPS_ARRAY_INDEX_CHECK) 00193 this->validateIndex(i1, i2); // Throws an exception on failure 00194 #endif 00195 return this->contiguous_p ? this->begin_p[i1 + i2*yinc_p] : 00196 this->begin_p[i1*xinc_p + i2*yinc_p]; 00197 } 00198 00199 const T &operator()(size_t i1, size_t i2) const 00200 { 00201 #if defined(AIPS_ARRAY_INDEX_CHECK) 00202 this->validateIndex(i1, i2); // Throws an exception on failure 00203 #endif 00204 return this->contiguous_p ? this->begin_p[i1 + i2*yinc_p] : 00205 this->begin_p[i1*xinc_p + i2*yinc_p]; 00206 } 00207 // </group> 00208 00209 00210 // The array is masked by the input LogicalArray. 00211 // This mask must conform to the array. 00212 // <group> 00213 00214 // Return a MaskedArray. 00215 MaskedArray<T> operator() (const LogicalArray &mask) const 00216 { return Array<T>::operator() (mask); } 00217 00218 // Return a MaskedArray. 00219 MaskedArray<T> operator() (const LogicalArray &mask) 00220 { return Array<T>::operator() (mask); } 00221 00222 // </group> 00223 00224 00225 // The array is masked by the input MaskedLogicalArray. 00226 // The mask is effectively the AND of the internal LogicalArray 00227 // and the internal mask of the MaskedLogicalArray. 00228 // The MaskedLogicalArray must conform to the array. 00229 // <group> 00230 00231 // Return a MaskedArray. 00232 MaskedArray<T> operator() (const MaskedLogicalArray &mask) const 00233 { return Array<T>::operator() (mask); } 00234 00235 // Return a MaskedArray. 00236 MaskedArray<T> operator() (const MaskedLogicalArray &mask) 00237 { return Array<T>::operator() (mask); } 00238 00239 // </group> 00240 00241 00242 // Returns a reference to the i'th row. 00243 // <group> 00244 Vector<T> row(size_t i); 00245 const Vector<T> row(size_t i) const; 00246 // </group> 00247 00248 // Returns a reference to the j'th column 00249 // <group> 00250 Vector<T> column(size_t j); 00251 const Vector<T> column(size_t j) const; 00252 // </group> 00253 00254 // Returns a diagonal from the Matrix. The Matrix must be square. 00255 // n==0 is the main diagonal. n>0 is above the main diagonal, n<0 00256 // is below it. 00257 // <group> 00258 Vector<T> diagonal(Int64 n=0); 00259 const Vector<T> diagonal(Int64 n=0) const; 00260 // </group> 00261 00262 // Take a slice of this matrix. Slices are always indexed starting 00263 // at zero. This uses reference semantics, i.e. changing a value 00264 // in the slice changes the original. 00265 // <srcblock> 00266 // Matrix<Double> vd(100,100); 00267 // //... 00268 // vd(Slice(0,10),Slice(10,10)) = -1.0; // 10x10 sub-matrix set to -1.0 00269 // </srcblock> 00270 // <group> 00271 Matrix<T> operator()(const Slice &sliceX, const Slice &sliceY); 00272 const Matrix<T> operator()(const Slice &sliceX, const Slice &sliceY) const; 00273 // </group> 00274 00275 // Slice using IPositions. Required to be defined, otherwise the base 00276 // class versions are hidden. 00277 // <group> 00278 Array<T> operator()(const IPosition &blc, const IPosition &trc, 00279 const IPosition &incr) 00280 { return Array<T>::operator()(blc,trc,incr); } 00281 const Array<T> operator()(const IPosition &blc, const IPosition &trc, 00282 const IPosition &incr) const 00283 { return Array<T>::operator()(blc,trc,incr); } 00284 Array<T> operator()(const IPosition &blc, const IPosition &trc) 00285 { return Array<T>::operator()(blc,trc); } 00286 const Array<T> operator()(const IPosition &blc, const IPosition &trc) const 00287 { return Array<T>::operator()(blc,trc); } 00288 Array<T> operator()(const Slicer& slicer) 00289 { return Array<T>::operator()(slicer); } 00290 const Array<T> operator()(const Slicer& slicer) const 00291 { return Array<T>::operator()(slicer); } 00292 // </group> 00293 00294 // The length of each axis of the Matrix. 00295 const IPosition &shape() const 00296 { return this->length_p; } 00297 void shape(Int &s1, Int &s2) const 00298 { s1 = this->length_p(0); s2=this->length_p(1); } 00299 00300 // The number of rows in the Matrix, i.e. the length of the first axis. 00301 size_t nrow() const 00302 { return this->length_p(0); } 00303 00304 // The number of columns in the Matrix, i.e. the length of the 2nd axis. 00305 size_t ncolumn() const 00306 { return this->length_p(1); } 00307 00308 // Checks that the Matrix is consistent (invariants check out). 00309 virtual Bool ok() const; 00310 00311 protected: 00312 virtual void preTakeStorage(const IPosition &shape); 00313 virtual void postTakeStorage(); 00314 // Remove the degenerate axes from other and store result in this matrix. 00315 // An exception is thrown if removing degenerate axes does not result 00316 // in a matrix. 00317 virtual void doNonDegenerate(const Array<T> &other, 00318 const IPosition &ignoreAxes); 00319 00320 private: 00321 // Cached constants to improve indexing. 00322 size_t xinc_p, yinc_p; 00323 00324 // Helper fn to calculate the indexing constants. 00325 void makeIndexingConstants(); 00326 }; 00327 00328 //# Declare extern templates for often used types. 00329 #ifdef AIPS_CXX11 00330 extern template class Matrix<Bool>; 00331 extern template class Matrix<Float>; 00332 extern template class Matrix<Double>; 00333 extern template class Matrix<Complex>; 00334 extern template class Matrix<DComplex>; 00335 #endif 00336 00337 } //#End casa namespace 00338 #ifndef CASACORE_NO_AUTO_TEMPLATES 00339 #include <casacore/casa/Arrays/Matrix.tcc> 00340 #endif //# CASACORE_NO_AUTO_TEMPLATES 00341 #endif