ArrayMath.h

Go to the documentation of this file.
00001 //# ArrayMath.h: ArrayMath: Simple mathematics done on an entire array.
00002 //# Copyright (C) 1993,1994,1995,1996,1998,1999,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_ARRAYMATH_H
00029 #define CASA_ARRAYMATH_H
00030 
00031 #include <casacore/casa/aips.h>
00032 #include <casacore/casa/BasicMath/Math.h>
00033 #include <casacore/casa/BasicMath/Functors.h>
00034 #include <casacore/casa/Arrays/Array.h>
00035 //# Needed to get the proper Complex typedef's
00036 #include <casacore/casa/BasicSL/Complex.h>
00037 #include <casacore/casa/Utilities/Assert.h>
00038 #include <casacore/casa/Exceptions/Error.h>
00039 #include <numeric>
00040 #include <functional>
00041 
00042 namespace casacore { //# NAMESPACE CASACORE - BEGIN
00043 
00044 template<class T> class Matrix;
00045 
00046 // <summary>
00047 //    Mathematical operations for Arrays.
00048 // </summary>
00049 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="tArray">
00050 //
00051 // <prerequisite>
00052 //   <li> <linkto class=Array>Array</linkto>
00053 // </prerequisite>
00054 //
00055 // <etymology>
00056 // This file contains global functions which perform element by element
00057 // mathematical operations on arrays.
00058 // </etymology>
00059 //
00060 // <synopsis>
00061 // These functions perform element by element mathematical operations on
00062 // arrays.  The two arrays must conform.
00063 //
00064 // Furthermore it defines functions a la std::transform to transform one or
00065 // two arrays by means of a unary or binary operator. All math and logical
00066 // operations on arrays can be expressed by means of these transform functions.
00067 // <br>It also defines an in-place transform function because for non-trivial
00068 // iterators it works faster than a transform where the result is an iterator
00069 // on the same data object as the left operand.
00070 // <br>The transform functions distinguish between contiguous and non-contiguous
00071 // arrays because iterating through a contiguous array can be done in a faster
00072 // way.
00073 // <br> Similar to the standard transform function these functions do not check
00074 // if the shapes match. The user is responsible for that.
00075 // </synopsis>
00076 //
00077 // <example>
00078 // <srcblock>
00079 //   Vector<Int> a(10);
00080 //   Vector<Int> b(10);
00081 //   Vector<Int> c(10);
00082 //      . . .
00083 //   c = a + b;
00084 // </srcblock>
00085 // This example sets the elements of c to (a+b). It checks if a and b have the
00086 // same shape.
00087 // The result of this operation is an Array.
00088 // </example>
00089 //
00090 // <example>
00091 // <srcblock>
00092 //   c = arrayTransformResult (a, b, std::plus<Double>());
00093 // </srcblock>
00094 // This example does the same as the previous example, but expressed using
00095 // the transform function (which, in fact, is used by the + operator above).
00096 // However, it is not checked if the shapes match.
00097 // </example>
00098 
00099 // <example>
00100 // <srcblock>
00101 //   arrayContTransform (a, b, c, std::plus<Double>());
00102 // </srcblock>
00103 // This example does the same as the previous example, but is faster because
00104 // the result array already exists and does not need to be allocated.
00105 // Note that the caller must be sure that c is contiguous.
00106 // </example>
00107 
00108 // <example>
00109 // <srcblock>
00110 //   Vector<Double> a(10);
00111 //   Vector<Double> b(10);
00112 //   Vector<Double> c(10);
00113 //      . . .
00114 //   c = atan2 (a, b);
00115 // </srcblock>
00116 // This example sets the elements of c to atan2 (a,b).
00117 // The result of this operation is an Array.
00118 // </example>
00119 //
00120 // <example>
00121 // <srcblock>
00122 //   Vector<Int> a(10);
00123 //   Int result;
00124 //      . . .
00125 //   result = sum (a);
00126 // </srcblock>
00127 // This example sums a.
00128 // </example>
00129 //
00130 // <motivation>
00131 // One wants to be able to perform mathematical operations on arrays.
00132 // </motivation>
00133 //
00134 // <linkfrom anchor="Array mathematical operations" classes="Array Vector Matrix Cube">
00135 //    <here>Array mathematical operations</here> -- Mathematical operations for
00136 //    Arrays.
00137 // </linkfrom>
00138 //
00139 // <group name="Array mathematical operations">
00140 
00141 
00142   // The myxtransform functions are defined to avoid a bug in g++-4.3.
00143   // That compiler generates incorrect code when only -g is used for
00144   // a std::transform with a bind1st or bind2nd for a complex<float>.
00145   // So, for example, the multiplication of a Complex array and Complex scalar
00146   // would fail (see g++ bug 39678).
00147   // <group>
00148   // sequence = scalar OP sequence
00149   template<typename _InputIterator1, typename T,
00150            typename _OutputIterator, typename _BinaryOperation>
00151     void
00152     myltransform(_InputIterator1 __first1, _InputIterator1 __last1,
00153                  _OutputIterator __result, T left,
00154                  _BinaryOperation __binary_op)
00155     {
00156       for ( ; __first1 != __last1; ++__first1, ++__result)
00157         *__result = __binary_op(left, *__first1);
00158     }
00159   // sequence = sequence OP scalar
00160   template<typename _InputIterator1, typename T,
00161            typename _OutputIterator, typename _BinaryOperation>
00162     void
00163     myrtransform(_InputIterator1 __first1, _InputIterator1 __last1,
00164                  _OutputIterator __result, T right,
00165                  _BinaryOperation __binary_op)
00166     {
00167       for ( ; __first1 != __last1; ++__first1, ++__result)
00168         *__result = __binary_op(*__first1, right);
00169     }
00170   // sequence OP= scalar
00171   template<typename _InputIterator1, typename T,
00172            typename _BinaryOperation>
00173     void
00174     myiptransform(_InputIterator1 __first1, _InputIterator1 __last1,
00175                   T right,
00176                   _BinaryOperation __binary_op)
00177     {
00178       for ( ; __first1 != __last1; ++__first1)
00179         *__first1 = __binary_op(*__first1, right);
00180     }
00181   // </group>
00182 
00183 
00184 // Functions to apply a binary or unary operator to arrays.
00185 // They are modeled after std::transform.
00186 // They do not check if the shapes conform; as in std::transform the
00187 // user must take care that the operands conform.
00188 // <group>
00189 // Transform left and right to a result using the binary operator.
00190 // Result MUST be a contiguous array.
00191 template<typename L, typename R, typename RES, typename BinaryOperator>
00192 inline void arrayContTransform (const Array<L>& left, const Array<R>& right,
00193                                 Array<RES>& result, BinaryOperator op)
00194 {
00195   DebugAssert (result.contiguousStorage(), AipsError);
00196   if (left.contiguousStorage()  &&  right.contiguousStorage()) {
00197     std::transform (left.cbegin(), left.cend(), right.cbegin(),
00198                     result.cbegin(), op);
00199   } else {
00200     std::transform (left.begin(), left.end(), right.begin(),
00201                     result.cbegin(), op);
00202   }
00203 }
00204 
00205 // Transform left and right to a result using the binary operator.
00206 // Result MUST be a contiguous array.
00207 template<typename L, typename R, typename RES, typename BinaryOperator>
00208 inline void arrayContTransform (const Array<L>& left, R right,
00209                                 Array<RES>& result, BinaryOperator op)
00210 {
00211   DebugAssert (result.contiguousStorage(), AipsError);
00212   if (left.contiguousStorage()) {
00213     myrtransform (left.cbegin(), left.cend(),
00214                  result.cbegin(), right, op);
00217   } else {
00218     myrtransform (left.begin(), left.end(),
00219                  result.cbegin(), right, op);
00222   }
00223 }
00224 
00225 // Transform left and right to a result using the binary operator.
00226 // Result MUST be a contiguous array.
00227 template<typename L, typename R, typename RES, typename BinaryOperator>
00228 inline void arrayContTransform (L left, const Array<R>& right,
00229                                 Array<RES>& result, BinaryOperator op)
00230 {
00231   DebugAssert (result.contiguousStorage(), AipsError);
00232   if (right.contiguousStorage()) {
00233     myltransform (right.cbegin(), right.cend(),
00234                   result.cbegin(), left, op);
00237   } else {
00238     myltransform (right.begin(), right.end(),
00239                   result.cbegin(), left, op);
00242   }
00243 }
00244 
00245 // Transform array to a result using the unary operator.
00246 // Result MUST be a contiguous array.
00247 template<typename T, typename RES, typename UnaryOperator>
00248 inline void arrayContTransform (const Array<T>& arr,
00249                                 Array<RES>& result, UnaryOperator op)
00250 {
00251   DebugAssert (result.contiguousStorage(), AipsError);
00252   if (arr.contiguousStorage()) {
00253     std::transform (arr.cbegin(), arr.cend(), result.cbegin(), op);
00254   } else {
00255     std::transform (arr.begin(), arr.end(), result.cbegin(), op);
00256   }
00257 }
00258 
00259 // Transform left and right to a result using the binary operator.
00260 // Result need not be a contiguous array.
00261 template<typename L, typename R, typename RES, typename BinaryOperator>
00262 void arrayTransform (const Array<L>& left, const Array<R>& right,
00263                      Array<RES>& result, BinaryOperator op);
00264 
00265 // Transform left and right to a result using the binary operator.
00266 // Result need not be a contiguous array.
00267 template<typename L, typename R, typename RES, typename BinaryOperator>
00268 void arrayTransform (const Array<L>& left, R right,
00269                      Array<RES>& result, BinaryOperator op);
00270 
00271 // Transform left and right to a result using the binary operator.
00272 // Result need not be a contiguous array.
00273 template<typename L, typename R, typename RES, typename BinaryOperator>
00274 void arrayTransform (L left, const Array<R>& right,
00275                      Array<RES>& result, BinaryOperator op);
00276 
00277 // Transform array to a result using the unary operator.
00278 // Result need not be a contiguous array.
00279 template<typename T, typename RES, typename UnaryOperator>
00280 void arrayTransform (const Array<T>& arr,
00281                      Array<RES>& result, UnaryOperator op);
00282 
00283 // Transform left and right to a result using the binary operator.
00284 // The created and returned result array is contiguous.
00285 template<typename T, typename BinaryOperator>
00286 Array<T> arrayTransformResult (const Array<T>& left, const Array<T>& right,
00287                                BinaryOperator op);
00288 
00289 // Transform left and right to a result using the binary operator.
00290 // The created and returned result array is contiguous.
00291 template<typename T, typename BinaryOperator>
00292 Array<T> arrayTransformResult (const Array<T>& left, T right, BinaryOperator op);
00293 
00294 // Transform left and right to a result using the binary operator.
00295 // The created and returned result array is contiguous.
00296 template<typename T, typename BinaryOperator>
00297 Array<T> arrayTransformResult (T left, const Array<T>& right, BinaryOperator op);
00298 
00299 // Transform array to a result using the unary operator.
00300 // The created and returned result array is contiguous.
00301 template<typename T, typename UnaryOperator>
00302 Array<T> arrayTransformResult (const Array<T>& arr, UnaryOperator op);
00303 
00304 // Transform left and right in place using the binary operator.
00305 // The result is stored in the left array (useful for e.g. the += operation).
00306 template<typename L, typename R, typename BinaryOperator>
00307 inline void arrayTransformInPlace (Array<L>& left, const Array<R>& right,
00308                                    BinaryOperator op)
00309 {
00310   if (left.contiguousStorage()  &&  right.contiguousStorage()) {
00311     transformInPlace (left.cbegin(), left.cend(), right.cbegin(), op);
00312   } else {
00313     transformInPlace (left.begin(), left.end(), right.begin(), op);
00314   }
00315 }
00316 
00317 // Transform left and right in place using the binary operator.
00318 // The result is stored in the left array (useful for e.g. the += operation).
00319 template<typename L, typename R, typename BinaryOperator>
00320 inline void arrayTransformInPlace (Array<L>& left, R right, BinaryOperator op)
00321 {
00322   if (left.contiguousStorage()) {
00323     myiptransform (left.cbegin(), left.cend(), right, op);
00325   } else {
00326     myiptransform (left.begin(), left.end(), right, op);
00328   }
00329 }
00330 
00331 // Transform the array in place using the unary operator.
00332 // E.g. doing <src>arrayTransformInPlace(array, Sin<T>())</src> is faster than
00333 // <src>array=sin(array)</src> as it does not need to create a temporary array.
00334 template<typename T, typename UnaryOperator>
00335 inline void arrayTransformInPlace (Array<T>& arr, UnaryOperator op)
00336 {
00337   if (arr.contiguousStorage()) {
00338     transformInPlace (arr.cbegin(), arr.cend(), op);
00339   } else {
00340     transformInPlace (arr.begin(), arr.end(), op);
00341   }
00342 }
00343 // </group>
00344 
00345 // 
00346 // Element by element arithmetic modifying left in-place. left and other
00347 // must be conformant.
00348 // <group>
00349 template<class T> void operator+= (Array<T> &left, const Array<T> &other);
00350 template<class T> void operator-= (Array<T> &left, const Array<T> &other);
00351 template<class T> void operator*= (Array<T> &left, const Array<T> &other);
00352 template<class T> void operator/= (Array<T> &left, const Array<T> &other);
00353 template<class T> void operator%= (Array<T> &left, const Array<T> &other);
00354 template<class T> void operator&= (Array<T> &left, const Array<T> &other);
00355 template<class T> void operator|= (Array<T> &left, const Array<T> &other);
00356 template<class T> void operator^= (Array<T> &left, const Array<T> &other);
00357 // </group>
00358 
00359 // 
00360 // Element by element arithmetic modifying left in-place. The scalar "other"
00361 // behaves as if it were a conformant Array to left filled with constant values.
00362 // <group>
00363 template<class T> void operator+= (Array<T> &left, const T &other);
00364 template<class T> void operator-= (Array<T> &left, const T &other);
00365 template<class T> void operator*= (Array<T> &left, const T &other);
00366 template<class T> void operator/= (Array<T> &left, const T &other);
00367 template<class T> void operator%= (Array<T> &left, const T &other);
00368 template<class T> void operator&= (Array<T> &left, const T &other);
00369 template<class T> void operator|= (Array<T> &left, const T &other);
00370 template<class T> void operator^= (Array<T> &left, const T &other);
00371 // </group>
00372 
00373 // Unary arithmetic operation.
00374 // 
00375 // <group>
00376 template<class T> Array<T> operator+(const Array<T> &a);
00377 template<class T> Array<T> operator-(const Array<T> &a);
00378 template<class T> Array<T> operator~(const Array<T> &a);
00379 // </group>
00380 
00381 // 
00382 // Element by element arithmetic on two arrays, returning an array.
00383 // <group>
00384 template<class T> 
00385   Array<T> operator+ (const Array<T> &left, const Array<T> &right);
00386 template<class T> 
00387   Array<T> operator- (const Array<T> &left, const Array<T> &right);
00388 template<class T> 
00389   Array<T> operator* (const Array<T> &left, const Array<T> &right);
00390 template<class T> 
00391   Array<T> operator/ (const Array<T> &left, const Array<T> &right);
00392 template<class T> 
00393   Array<T> operator% (const Array<T> &left, const Array<T> &right);
00394 template<class T> 
00395   Array<T> operator| (const Array<T> &left, const Array<T> &right);
00396 template<class T> 
00397   Array<T> operator& (const Array<T> &left, const Array<T> &right);
00398 template<class T> 
00399   Array<T> operator^ (const Array<T> &left, const Array<T> &right);
00400 // </group>
00401 
00402 // 
00403 // Element by element arithmetic between an array and a scalar, returning
00404 // an array.
00405 // <group>
00406 template<class T> 
00407     Array<T> operator+ (const Array<T> &left, const T &right);
00408 template<class T> 
00409     Array<T> operator- (const Array<T> &left, const T &right);
00410 template<class T> 
00411     Array<T> operator* (const Array<T> &left, const T &right);
00412 template<class T> 
00413     Array<T> operator/ (const Array<T> &left, const T &right);
00414 template<class T> 
00415     Array<T> operator% (const Array<T> &left, const T &right);
00416 template<class T> 
00417     Array<T> operator| (const Array<T> &left, const T &right);
00418 template<class T> 
00419     Array<T> operator& (const Array<T> &left, const T &right);
00420 template<class T> 
00421     Array<T> operator^ (const Array<T> &left, const T &right);
00422 // </group>
00423 
00424 // 
00425 // Element by element arithmetic between a scalar and an array, returning
00426 // an array.
00427 // <group>
00428 template<class T>  
00429     Array<T> operator+ (const T &left, const Array<T> &right);
00430 template<class T>  
00431     Array<T> operator- (const T &left, const Array<T> &right);
00432 template<class T>  
00433     Array<T> operator* (const T &left, const Array<T> &right);
00434 template<class T>  
00435     Array<T> operator/ (const T &left, const Array<T> &right);
00436 template<class T>  
00437     Array<T> operator% (const T &left, const Array<T> &right);
00438 template<class T>  
00439     Array<T> operator| (const T &left, const Array<T> &right);
00440 template<class T>  
00441     Array<T> operator& (const T &left, const Array<T> &right);
00442 template<class T>  
00443     Array<T> operator^ (const T &left, const Array<T> &right);
00444 // </group>
00445 
00446 // 
00447 // Transcendental function that can be applied to essentially all numeric
00448 // types. Works on an element-by-element basis.
00449 // <group>
00450 template<class T> Array<T> cos(const Array<T> &a);
00451 template<class T> Array<T> cosh(const Array<T> &a);
00452 template<class T> Array<T> exp(const Array<T> &a);
00453 template<class T> Array<T> log(const Array<T> &a);
00454 template<class T> Array<T> log10(const Array<T> &a);
00455 template<class T> Array<T> pow(const Array<T> &a, const Array<T> &b);
00456 template<class T> Array<T> pow(const T &a, const Array<T> &b);
00457 template<class T> Array<T> sin(const Array<T> &a);
00458 template<class T> Array<T> sinh(const Array<T> &a);
00459 template<class T> Array<T> sqrt(const Array<T> &a);
00460 // </group>
00461 
00462 // 
00463 // Transcendental function applied to the array on an element-by-element
00464 // basis. Although a template function, this does not make sense for all
00465 // numeric types.
00466 // <group>
00467 template<class T> Array<T> acos(const Array<T> &a);
00468 template<class T> Array<T> asin(const Array<T> &a);
00469 template<class T> Array<T> atan(const Array<T> &a);
00470 template<class T> Array<T> atan2(const Array<T> &y, const Array<T> &x);
00471 template<class T> Array<T> atan2(const T &y, const Array<T> &x);
00472 template<class T> Array<T> atan2(const Array<T> &y, const T &x);
00473 template<class T> Array<T> ceil(const Array<T> &a);
00474 template<class T> Array<T> fabs(const Array<T> &a);
00475 template<class T> Array<T> abs(const Array<T> &a);
00476 template<class T> Array<T> floor(const Array<T> &a);
00477 template<class T> Array<T> round(const Array<T> &a);
00478 template<class T> Array<T> sign(const Array<T> &a);
00479 template<class T> Array<T> fmod(const Array<T> &a, const Array<T> &b);
00480 template<class T> Array<T> fmod(const T &a, const Array<T> &b);
00481 template<class T> Array<T> fmod(const Array<T> &a, const T &b);
00482 template<class T> Array<T> floormod(const Array<T> &a, const Array<T> &b);
00483 template<class T> Array<T> floormod(const T &a, const Array<T> &b);
00484 template<class T> Array<T> floormod(const Array<T> &a, const T &b);
00485 template<class T> Array<T> pow(const Array<T> &a, const Double &b);
00486 template<class T> Array<T> tan(const Array<T> &a);
00487 template<class T> Array<T> tanh(const Array<T> &a);
00488 // N.B. fabs is deprecated. Use abs.
00489 template<class T> Array<T> fabs(const Array<T> &a);
00490 // </group>
00491 
00492 // 
00493 // <group>
00494 // Find the minimum and maximum values of an array, including their locations.
00495 template<class ScalarType>
00496 void minMax(ScalarType &minVal, ScalarType &maxVal, IPosition &minPos, 
00497             IPosition &maxPos, const Array<ScalarType> &array);
00498 // The array is searched at locations where the mask equals <src>valid</src>.
00499 // (at least one such position must exist or an exception will be thrown).
00500 // MaskType should be an Array of Bool.
00501 template<class ScalarType>
00502 void minMax(ScalarType &minVal, ScalarType &maxVal, IPosition &minPos,
00503             IPosition &maxPos, const Array<ScalarType> &array, 
00504             const Array<Bool> &mask, Bool valid=True);
00505 // The array * weight is searched 
00506 template<class ScalarType>
00507 void minMaxMasked(ScalarType &minVal, ScalarType &maxVal, IPosition &minPos,
00508                   IPosition &maxPos, const Array<ScalarType> &array, 
00509                   const Array<ScalarType> &weight);
00510 // </group>
00511 
00512 // 
00513 // The "min" and "max" functions require that the type "T" have comparison 
00514 // operators.
00515 // <group>
00516 //
00517 // This sets min and max to the minimum and maximum of the array to 
00518 // avoid having to do two passes with max() and min() separately.
00519 template<class T> void minMax(T &min, T &max, const Array<T> &a);
00520 //
00521 // The minimum element of the array.
00522 // Requires that the type "T" has comparison operators.
00523 template<class T>  T min(const Array<T> &a);
00524 // The maximum element of the array.
00525 // Requires that the type "T" has comparison operators.
00526 template<class T>  T max(const Array<T> &a);
00527 
00528 // "result" contains the maximum of "a" and "b" at each position. "result",
00529 // "a", and "b" must be conformant.
00530 template<class T> void max(Array<T> &result, const Array<T> &a, 
00531                            const Array<T> &b);
00532 // "result" contains the minimum of "a" and "b" at each position. "result",
00533 // "a", and "b" must be conformant.
00534 template<class T> void min(Array<T> &result, const Array<T> &a, 
00535                            const Array<T> &b);
00536 // Return an array that contains the maximum of "a" and "b" at each position.
00537 // "a" and "b" must be conformant.
00538 template<class T> Array<T> max(const Array<T> &a, const Array<T> &b);
00539 template<class T> Array<T> max(const T &a, const Array<T> &b);
00540 // Return an array that contains the minimum of "a" and "b" at each position.
00541 // "a" and "b" must be conformant.
00542 template<class T> Array<T> min(const Array<T> &a, const Array<T> &b);
00543 
00544 // "result" contains the maximum of "a" and "b" at each position. "result",
00545 // and "a" must be conformant.
00546 template<class T> void max(Array<T> &result, const Array<T> &a, 
00547                            const T &b);
00548 template<class T> inline void max(Array<T> &result, const T &a, 
00549                                   const Array<T> &b)
00550   { max (result, b, a); }
00551 // "result" contains the minimum of "a" and "b" at each position. "result",
00552 // and "a" must be conformant.
00553 template<class T> void min(Array<T> &result, const Array<T> &a, 
00554                            const T &b);
00555 template<class T> inline void min(Array<T> &result, const T &a, 
00556                                   const Array<T> &b)
00557   { min (result, b, a); }
00558 // Return an array that contains the maximum of "a" and "b" at each position.
00559 template<class T> Array<T> max(const Array<T> &a, const T &b);
00560 template<class T> inline Array<T> max(const T &a, const Array<T> &b)
00561   { return max(b, a); }
00562 // Return an array that contains the minimum of "a" and "b" at each position.
00563 template<class T> Array<T> min(const Array<T> &a, const T &b);
00564 template<class T> inline Array<T> min(const T &a, const Array<T> &b)
00565   { return min(b, a); }
00566 // </group>
00567 
00568 // 
00569 // Fills all elements of "array" with a sequence starting with "start"
00570 // and incrementing by "inc" for each element. The first axis varies
00571 // most rapidly.
00572 template<class T> void indgen(Array<T> &a, T start, T inc);
00573 // 
00574 // Fills all elements of "array" with a sequence starting with 0
00575 // and ending with nelements() - 1. The first axis varies
00576 // most rapidly.
00577 template<class T> inline void indgen(Array<T> &a)
00578   { indgen(a, T(0), T(1)); }
00579 // 
00580 // Fills all elements of "array" with a sequence starting with start
00581 // incremented by one for each position in the array. The first axis varies
00582 // most rapidly.
00583 template<class T> inline void indgen(Array<T> &a, T start)
00584   { indgen(a, start, T(1)); }
00585 
00586 // Create a Vector of the given length and fill it with the start value
00587 // incremented with <code>inc</code> for each element.
00588 template<class T> inline Vector<T> indgen(uInt length, T start, T inc)
00589 {
00590   Vector<T> x(length);
00591   indgen(x, start, inc);
00592   return x;
00593 }
00594 
00595 
00596 // Sum of every element of the array.
00597 template<class T> T sum(const Array<T> &a);
00598 // 
00599 // Sum the square of every element of the array.
00600 template<class T> T sumsqr(const Array<T> &a);
00601 // 
00602 // Product of every element of the array. This could of course easily
00603 // overflow.
00604 template<class T> T product(const Array<T> &a);
00605 
00606 // 
00607 // The mean of "a" is the sum of all elements of "a" divided by the number
00608 // of elements of "a".
00609 template<class T> T mean(const Array<T> &a);
00610 
00611 // 
00612 // The variance of "a" is the sum of (a(i) - mean(a))**2/(a.nelements() - 1).
00613 // N.B. N-1, not N in the denominator).
00614 template<class T> T variance(const Array<T> &a);
00615 // 
00616 // The variance of "a" is the sum of (a(i) - mean(a))**2/(a.nelements() - 1).
00617 // N.B. N-1, not N in the denominator).
00618 // Rather than using a computed mean, use the supplied value.
00619 template<class T> T variance(const Array<T> &a, T mean);
00620 
00621 // 
00622 // The standard deviation of "a" is the square root of its variance.
00623 template<class T> T stddev(const Array<T> &a);
00624 // 
00625 // The standard deviation of "a" is the square root of its variance.
00626 // Rather than using a computed mean, use the supplied value.
00627 template<class T> T stddev(const Array<T> &a, T mean);
00628 
00629 // 
00630 // The average deviation of "a" is the sum of abs(a(i) - mean(a))/N. (N.B.
00631 // N, not N-1 in the denominator).
00632 template<class T> T avdev(const Array<T> &a);
00633 // 
00634 // The average deviation of "a" is the sum of abs(a(i) - mean(a))/N. (N.B.
00635 // N, not N-1 in the denominator).
00636 // Rather than using a computed mean, use the supplied value.
00637 template<class T> T avdev(const Array<T> &a,T mean);
00638 
00639 //
00640 // The root-mean-square of "a" is the sqrt of sum(a*a)/N.
00641 template<class T> T rms(const Array<T> &a);
00642 
00643 
00644 // The median of "a" is a(n/2).
00645 // If a has an even number of elements and the switch takeEvenMean is set,
00646 // the median is 0.5*(a(n/2) + a((n+1)/2)).
00647 // According to Numerical Recipes (2nd edition) it makes little sense to take
00648 // the mean if the array is large enough (> 100 elements). Therefore
00649 // the default for takeEvenMean is False if the array has > 100 elements,
00650 // otherwise it is True.
00651 // <br>If "sorted"==True we assume the data is already sorted and we
00652 // compute the median directly. Otherwise the function GenSort::kthLargest
00653 // is used to find the median (kthLargest is about 6 times faster
00654 // than a full quicksort).
00655 // <br>Finding the median means that the array has to be (partially)
00656 // sorted. By default a copy will be made, but if "inPlace" is in effect,
00657 // the data themselves will be sorted. That should only be used if the
00658 // data are used not thereafter.
00659 // <note>The function kthLargest in class GenSortIndirect can be used to
00660 // obtain the index of the median in an array. </note>
00661 // <group>
00662 template<class T> T median(const Array<T> &a, Block<T> &tmp, Bool sorted,
00663                            Bool takeEvenMean, Bool inPlace=False);
00664 template<class T> T median(const Array<T> &a, Bool sorted, Bool takeEvenMean,
00665                            Bool inPlace=False)
00666     { Block<T> tmp; return median (a, tmp, sorted, takeEvenMean, inPlace); }
00667 template<class T> inline T median(const Array<T> &a, Bool sorted)
00668     { return median (a, sorted, (a.nelements() <= 100), False); }
00669 template<class T> inline T median(const Array<T> &a)
00670     { return median (a, False, (a.nelements() <= 100), False); }
00671 template<class T> inline T medianInPlace(const Array<T> &a, Bool sorted=False)
00672     { return median (a, sorted, (a.nelements() <= 100), True); }
00673 // </group>
00674 
00675 // The median absolute deviation from the median. Interface is as for
00676 // the median functions
00677 // <group>
00678 template<class T> T madfm(const Array<T> &a, Block<T> &tmp, Bool sorted, 
00679                          Bool takeEvenMean, Bool inPlace = False);
00680 template<class T> T madfm(const Array<T> &a, Bool sorted, Bool takeEvenMean,
00681                           Bool inPlace=False)
00682     { Block<T> tmp; return madfm(a, tmp, sorted, takeEvenMean, inPlace); }
00683 template<class T> inline T madfm(const Array<T> &a, Bool sorted)
00684     { return madfm(a, sorted, (a.nelements() <= 100), False); }
00685 template<class T> inline T madfm(const Array<T> &a)
00686     { return madfm(a, False, (a.nelements() <= 100), False); }
00687 template<class T> inline T madfmInPlace(const Array<T> &a, Bool sorted=False)
00688     { return madfm(a, sorted, (a.nelements() <= 100), True); }
00689 // </group>
00690 
00691 // Return the fractile of an array.
00692 // It returns the value at the given fraction of the array.
00693 // A fraction of 0.5 is the same as the median, be it that no mean of
00694 // the two middle elements is taken if the array has an even nr of elements.
00695 // It uses kthLargest if the array is not sorted yet.
00696 // <note>The function kthLargest in class GenSortIndirect can be used to
00697 // obtain the index of the fractile in an array. </note>
00698 template<class T> T fractile(const Array<T> &a, Block<T> &tmp, Float fraction,
00699                              Bool sorted=False, Bool inPlace=False);
00700 template<class T> T fractile(const Array<T> &a, Float fraction,
00701                              Bool sorted=False, Bool inPlace=False)
00702   { Block<T> tmp; return fractile (a, tmp, fraction, sorted, inPlace); }
00703 
00704 // Return the inter-fractile range of an array.  
00705 // This is the full range between the bottom and the top fraction.
00706 // <group>
00707 template<class T> T interFractileRange(const Array<T> &a, Block<T> &tmp,
00708                                        Float fraction,
00709                                        Bool sorted=False, Bool inPlace=False);
00710 template<class T> T interFractileRange(const Array<T> &a, Float fraction,
00711                                        Bool sorted=False, Bool inPlace=False)
00712   { Block<T> tmp; return interFractileRange(a, tmp, fraction, sorted, inPlace); }
00713 // </group>
00714 
00715 // Return the inter-hexile range of an array.  
00716 // This is the full range between the bottom sixth and the top sixth
00717 // of ordered array values. "The semi-interhexile range is very nearly
00718 // equal to the rms for a Gaussian distribution, but it is much less
00719 // sensitive to the tails of extended distributions." (Condon et al
00720 // 1998)
00721 // <group>
00722 template<class T> T interHexileRange(const Array<T> &a, Block<T> &tmp,
00723                                      Bool sorted=False, Bool inPlace=False)
00724   { return interFractileRange(a, tmp, 1./6., sorted, inPlace); }
00725 template<class T> T interHexileRange(const Array<T> &a, Bool sorted=False,
00726                                      Bool inPlace=False)
00727   { return interFractileRange(a, 1./6., sorted, inPlace); }
00728 // </group>
00729 
00730 // Return the inter-quartile range of an array.  
00731 // This is the full range between the bottom quarter and the top
00732 // quarter of ordered array values.
00733 // <group>
00734 template<class T> T interQuartileRange(const Array<T> &a, Block<T> &tmp,
00735                                        Bool sorted=False, Bool inPlace=False)
00736   { return interFractileRange(a, tmp, 0.25, sorted, inPlace); }
00737 template<class T> T interQuartileRange(const Array<T> &a, Bool sorted=False,
00738                                        Bool inPlace=False)
00739   { return interFractileRange(a, 0.25, sorted, inPlace); }
00740 // </group>
00741 
00742 
00743 // Methods for element-by-element scaling of complex and real.
00744 // Note that Complex and DComplex are typedefs for std::complex.
00745 //<group>
00746 template<typename T>
00747 void operator*= (Array<std::complex<T> > &left, const Array<T> &other);
00748 template<typename T>
00749 void operator*= (Array<std::complex<T> > &left, const T &other);
00750 template<typename T>
00751 void operator/= (Array<std::complex<T> > &left, const Array<T> &other);
00752 template<typename T>
00753 void operator/= (Array<std::complex<T> > &left, const T &other);
00754 template<typename T>
00755 Array<std::complex<T> > operator* (const Array<std::complex<T> > &left,
00756                                    const Array<T> &right);
00757 template<typename T>
00758 Array<std::complex<T> > operator* (const Array<std::complex<T> > &left,
00759                                    const T &right);
00760 template<typename T>
00761 Array<std::complex<T> > operator* (const std::complex<T> &left,
00762                                    const Array<T> &right);
00763 template<typename T>
00764 Array<std::complex<T> > operator/ (const Array<std::complex<T> > &left,
00765                                    const Array<T> &right);
00766 template<typename T>
00767 Array<std::complex<T> > operator/ (const Array<std::complex<T> > &left,
00768                                    const T &right);
00769 template<typename T>
00770 Array<std::complex<T> > operator/ (const std::complex<T> &left,
00771                                    const Array<T> &right);
00772 // </group>
00773 
00774 // Returns the complex conjugate of a complex array.
00775 //<group>
00776 Array<Complex> conj(const Array<Complex> &carray);
00777 Array<DComplex> conj(const Array<DComplex> &carray);
00778 // Modifies rarray in place. rarray must be conformant.
00779 void         conj(Array<Complex> &rarray, const Array<Complex> &carray);
00780 void         conj(Array<DComplex> &rarray, const Array<DComplex> &carray);
00781 //# The following are implemented to make the compiler find the right conversion
00782 //# more often.
00783 Matrix<Complex> conj(const Matrix<Complex> &carray);
00784 Matrix<DComplex> conj(const Matrix<DComplex> &carray);
00785 //</group>
00786 
00787 // Form an array of complex numbers from the given real arrays.
00788 // Note that Complex and DComplex are simply typedefs for std::complex<float>
00789 // and std::complex<double>, so the result is in fact one of these types.
00790 // <group>
00791 template<typename T>
00792 Array<std::complex<T> > makeComplex(const Array<T> &real, const Array<T>& imag);
00793 template<typename T>
00794 Array<std::complex<T> > makeComplex(const T &real, const Array<T>& imag);
00795 template<typename T>
00796 Array<std::complex<T> > makeComplex(const Array<T> &real, const T& imag);
00797 // </group>
00798 
00799 // Set the real part of the left complex array to the right real array.
00800 template<typename L, typename R>
00801 void setReal(Array<L> &carray, const Array<R> &rarray);
00802 
00803 // Set the imaginary part of the left complex array to right real array.
00804 template<typename R, typename L>
00805 void setImag(Array<R> &carray, const Array<L> &rarray);
00806 
00807 // Extracts the real part of a complex array into an array of floats.
00808 // <group>
00809 Array<Float>  real(const Array<Complex> &carray);
00810 Array<Double> real(const Array<DComplex> &carray);
00811 // Modifies rarray in place. rarray must be conformant.
00812 void         real(Array<Float> &rarray, const Array<Complex> &carray);
00813 void         real(Array<Double> &rarray, const Array<DComplex> &carray);
00814 // </group>
00815 
00816 // 
00817 // Extracts the imaginary part of a complex array into an array of floats.
00818 // <group>
00819 Array<Float>  imag(const Array<Complex> &carray);
00820 Array<Double> imag(const Array<DComplex> &carray);
00821 // Modifies rarray in place. rarray must be conformant.
00822 void         imag(Array<Float> &rarray, const Array<Complex> &carray);
00823 void         imag(Array<Double> &rarray, const Array<DComplex> &carray);
00824 // </group>
00825 
00826 // 
00827 // Extracts the amplitude (i.e. sqrt(re*re + im*im)) from an array
00828 // of complex numbers. N.B. this is presently called "fabs" for a single
00829 // complex number.
00830 // <group>
00831 Array<Float>  amplitude(const Array<Complex> &carray);
00832 Array<Double> amplitude(const Array<DComplex> &carray);
00833 // Modifies rarray in place. rarray must be conformant.
00834 void         amplitude(Array<Float> &rarray, const Array<Complex> &carray);
00835 void         amplitude(Array<Double> &rarray, const Array<DComplex> &carray);
00836 // </group>
00837 
00838 // 
00839 // Extracts the phase (i.e. atan2(im, re)) from an array
00840 // of complex numbers. N.B. this is presently called "arg"
00841 // for a single complex number.
00842 // <group>
00843 Array<Float>  phase(const Array<Complex> &carray);
00844 Array<Double> phase(const Array<DComplex> &carray);
00845 // Modifies rarray in place. rarray must be conformant.
00846 void         phase(Array<Float> &rarray, const Array<Complex> &carray);
00847 void         phase(Array<Double> &rarray, const Array<DComplex> &carray);
00848 // </group>
00849 
00850 // Copy an array of complex into an array of real,imaginary pairs. The
00851 // first axis of the real array becomes twice as long as the complex array.
00852 // In the future versions which work by reference will be available; presently
00853 // a copy is made.
00854 Array<Float> ComplexToReal(const Array<Complex> &carray);
00855 Array<Double> ComplexToReal(const Array<DComplex> &carray);
00856 // Modify the array "rarray" in place. "rarray" must be the correct shape.
00857 // <group>
00858 void ComplexToReal(Array<Float> &rarray, const Array<Complex> &carray);
00859 void ComplexToReal(Array<Double> &rarray, const Array<DComplex> &carray);
00860 // </group>
00861 
00862 // Copy an array of real,imaginary pairs into a complex array. The first axis
00863 // must have an even length.
00864 // In the future versions which work by reference will be available; presently
00865 // a copy is made.
00866 Array<Complex>  RealToComplex(const Array<Float> &rarray);
00867 Array<DComplex> RealToComplex(const Array<Double> &rarray);
00868 // Modify the array "carray" in place. "carray" must be the correct shape.
00869 // <group>
00870 void  RealToComplex(Array<Complex> &carray, const Array<Float> &rarray);
00871 void  RealToComplex(Array<DComplex> &carray, const Array<Double> &rarray);
00872 // </group>
00873 
00874 // Make a copy of an array of a different type; for example make an array
00875 // of doubles from an array of floats. Arrays to and from must be conformant
00876 // (same shape). Also, it must be possible to convert a scalar of type U 
00877 // to type T.
00878 template<class T, class U> void convertArray(Array<T> &to,
00879                                              const Array<U> &from);
00880 
00881 
00882 // Returns an array where every element is squared.
00883 template<class T> Array<T> square(const Array<T> &val);
00884 
00885 // Returns an array where every element is cubed.
00886 template<class T> Array<T> cube(const Array<T> &val);
00887 
00888 // Expand the values of an array. The arrays must have the same dimensionality.
00889 // The length of each axis in the input array must be <= the length of the
00890 // corresponding axis in the output array and divide evenly.
00891 // For each axis <src>mult</src> is set to output/input.
00892 // <br>The <src>alternate</src> argument determines how the values are expanded.
00893 // If a row contains values '1 2 3', they can be expanded "linearly"
00894 // as '1 1 2 2 3 3'  or  alternately as '1 2 3 1 2 3'
00895 // This choice can be made for each axis; a value 0 means linearly,
00896 // another value means alternately. If the length of alternate is less than
00897 // the dimensionality of the arrays, the missing ones default to 0.
00898 template<class T> void expandArray (Array<T>& out, const Array<T>& in,
00899                                     const IPosition& alternate=IPosition());
00900 // Helper function for expandArray using recursion for each axis.
00901 template<class T>
00902 T* expandRecursive (int axis, const IPosition& shp, const IPosition& mult,
00903                     const IPosition& inSteps,
00904                     const T* in, T* out, const IPosition& alternate);
00905 // Check array shapes for expandArray. It returns the alternate argument,
00906 // where possibly missing values are appended (as 0).
00907 IPosition checkExpandArray (IPosition& mult,
00908                             const IPosition& inShape,
00909                             const IPosition& outShape,
00910                             const IPosition& alternate);
00911 
00912 
00913 // </group>
00914 
00915 
00916 } //# NAMESPACE CASACORE - END
00917 
00918 #ifndef CASACORE_NO_AUTO_TEMPLATES
00919 #include <casacore/casa/Arrays/ArrayMath.tcc>
00920 #endif //# CASACORE_NO_AUTO_TEMPLATES
00921 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1