casacore
ArrayMath.h
Go to the documentation of this file.
1 //# ArrayMath.h: ArrayMath: Simple mathematics done on an entire array.
2 //# Copyright (C) 1993,1994,1995,1996,1998,1999,2001,2003
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 CASA_ARRAYMATH_H
29 #define CASA_ARRAYMATH_H
30 
31 #include <casacore/casa/aips.h>
32 #include <casacore/casa/BasicMath/Math.h>
33 #include <casacore/casa/BasicMath/Functors.h>
34 #include <casacore/casa/Arrays/Array.h>
35 //# Needed to get the proper Complex typedef's
36 #include <casacore/casa/BasicSL/Complex.h>
37 #include <casacore/casa/Utilities/Assert.h>
38 #include <casacore/casa/Exceptions/Error.h>
39 #include <numeric>
40 #include <functional>
41 
42 namespace casacore { //# NAMESPACE CASACORE - BEGIN
43 
44 template<class T> class Matrix;
45 
46 // <summary>
47 // Mathematical operations for Arrays.
48 // </summary>
49 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="tArray">
50 //
51 // <prerequisite>
52 // <li> <linkto class=Array>Array</linkto>
53 // </prerequisite>
54 //
55 // <etymology>
56 // This file contains global functions which perform element by element
57 // mathematical operations on arrays.
58 // </etymology>
59 //
60 // <synopsis>
61 // These functions perform element by element mathematical operations on
62 // arrays. The two arrays must conform.
63 //
64 // Furthermore it defines functions a la std::transform to transform one or
65 // two arrays by means of a unary or binary operator. All math and logical
66 // operations on arrays can be expressed by means of these transform functions.
67 // <br>It also defines an in-place transform function because for non-trivial
68 // iterators it works faster than a transform where the result is an iterator
69 // on the same data object as the left operand.
70 // <br>The transform functions distinguish between contiguous and non-contiguous
71 // arrays because iterating through a contiguous array can be done in a faster
72 // way.
73 // <br> Similar to the standard transform function these functions do not check
74 // if the shapes match. The user is responsible for that.
75 // </synopsis>
76 //
77 // <example>
78 // <srcblock>
79 // Vector<Int> a(10);
80 // Vector<Int> b(10);
81 // Vector<Int> c(10);
82 // . . .
83 // c = a + b;
84 // </srcblock>
85 // This example sets the elements of c to (a+b). It checks if a and b have the
86 // same shape.
87 // The result of this operation is an Array.
88 // </example>
89 //
90 // <example>
91 // <srcblock>
92 // c = arrayTransformResult (a, b, std::plus<Double>());
93 // </srcblock>
94 // This example does the same as the previous example, but expressed using
95 // the transform function (which, in fact, is used by the + operator above).
96 // However, it is not checked if the shapes match.
97 // </example>
98 
99 // <example>
100 // <srcblock>
101 // arrayContTransform (a, b, c, std::plus<Double>());
102 // </srcblock>
103 // This example does the same as the previous example, but is faster because
104 // the result array already exists and does not need to be allocated.
105 // Note that the caller must be sure that c is contiguous.
106 // </example>
107 
108 // <example>
109 // <srcblock>
110 // Vector<Double> a(10);
111 // Vector<Double> b(10);
112 // Vector<Double> c(10);
113 // . . .
114 // c = atan2 (a, b);
115 // </srcblock>
116 // This example sets the elements of c to atan2 (a,b).
117 // The result of this operation is an Array.
118 // </example>
119 //
120 // <example>
121 // <srcblock>
122 // Vector<Int> a(10);
123 // Int result;
124 // . . .
125 // result = sum (a);
126 // </srcblock>
127 // This example sums a.
128 // </example>
129 //
130 // <motivation>
131 // One wants to be able to perform mathematical operations on arrays.
132 // </motivation>
133 //
134 // <linkfrom anchor="Array mathematical operations" classes="Array Vector Matrix Cube">
135 // <here>Array mathematical operations</here> -- Mathematical operations for
136 // Arrays.
137 // </linkfrom>
138 //
139 // <group name="Array mathematical operations">
141 
142  // The myxtransform functions are defined to avoid a bug in g++-4.3.
143  // That compiler generates incorrect code when only -g is used for
144  // a std::transform with a bind1st or bind2nd for a complex<float>.
145  // So, for example, the multiplication of a Complex array and Complex scalar
146  // would fail (see g++ bug 39678).
147  // <group>
148  // sequence = scalar OP sequence
149  template<typename _InputIterator1, typename T,
150  typename _OutputIterator, typename _BinaryOperation>
151  void
152  myltransform(_InputIterator1 __first1, _InputIterator1 __last1,
153  _OutputIterator __result, T left,
154  _BinaryOperation __binary_op)
155  {
156  for ( ; __first1 != __last1; ++__first1, ++__result)
157  *__result = __binary_op(left, *__first1);
158  }
159  // sequence = sequence OP scalar
160  template<typename _InputIterator1, typename T,
161  typename _OutputIterator, typename _BinaryOperation>
162  void
163  myrtransform(_InputIterator1 __first1, _InputIterator1 __last1,
164  _OutputIterator __result, T right,
165  _BinaryOperation __binary_op)
166  {
167  for ( ; __first1 != __last1; ++__first1, ++__result)
168  *__result = __binary_op(*__first1, right);
169  }
170  // sequence OP= scalar
171  template<typename _InputIterator1, typename T,
172  typename _BinaryOperation>
173  void
174  myiptransform(_InputIterator1 __first1, _InputIterator1 __last1,
175  T right,
176  _BinaryOperation __binary_op)
177  {
178  for ( ; __first1 != __last1; ++__first1)
179  *__first1 = __binary_op(*__first1, right);
180  }
181  // </group>
182 
183 
184 // Function to check the shapes. It throws an exception if not equal.
185 // <group>
186 void throwArrayShapes (const char* name);
187 inline void checkArrayShapes (const ArrayBase& left, const ArrayBase& right,
188  const char* name)
189 {
190  if (! left.shape().isEqual (right.shape())) {
191  throwArrayShapes (name);
192  }
193 }
194 // </group>
195 
196 
197 // Functions to apply a binary or unary operator to arrays.
198 // They are modeled after std::transform.
199 // They do not check if the shapes conform; as in std::transform the
200 // user must take care that the operands conform.
201 // <group>
202 // Transform left and right to a result using the binary operator.
203 // Result MUST be a contiguous array.
204 template<typename L, typename R, typename RES, typename BinaryOperator>
205 inline void arrayContTransform (const Array<L>& left, const Array<R>& right,
206  Array<RES>& result, BinaryOperator op)
207 {
209  if (left.contiguousStorage() && right.contiguousStorage()) {
210  std::transform (left.cbegin(), left.cend(), right.cbegin(),
211  result.cbegin(), op);
212  } else {
213  std::transform (left.begin(), left.end(), right.begin(),
214  result.cbegin(), op);
215  }
216 }
217 
218 // Transform left and right to a result using the binary operator.
219 // Result MUST be a contiguous array.
220 template<typename L, typename R, typename RES, typename BinaryOperator>
221 inline void arrayContTransform (const Array<L>& left, R right,
222  Array<RES>& result, BinaryOperator op)
223 {
225  if (left.contiguousStorage()) {
226  myrtransform (left.cbegin(), left.cend(),
227  result.cbegin(), right, op);
230  } else {
231  myrtransform (left.begin(), left.end(),
232  result.cbegin(), right, op);
235  }
236 }
237 
238 // Transform left and right to a result using the binary operator.
239 // Result MUST be a contiguous array.
240 template<typename L, typename R, typename RES, typename BinaryOperator>
241 inline void arrayContTransform (L left, const Array<R>& right,
242  Array<RES>& result, BinaryOperator op)
243 {
245  if (right.contiguousStorage()) {
246  myltransform (right.cbegin(), right.cend(),
247  result.cbegin(), left, op);
250  } else {
251  myltransform (right.begin(), right.end(),
252  result.cbegin(), left, op);
255  }
256 }
257 
258 // Transform array to a result using the unary operator.
259 // Result MUST be a contiguous array.
260 template<typename T, typename RES, typename UnaryOperator>
261 inline void arrayContTransform (const Array<T>& arr,
262  Array<RES>& result, UnaryOperator op)
263 {
265  if (arr.contiguousStorage()) {
266  std::transform (arr.cbegin(), arr.cend(), result.cbegin(), op);
267  } else {
268  std::transform (arr.begin(), arr.end(), result.cbegin(), op);
269  }
270 }
271 
272 // Transform left and right to a result using the binary operator.
273 // Result need not be a contiguous array.
274 template<typename L, typename R, typename RES, typename BinaryOperator>
275 void arrayTransform (const Array<L>& left, const Array<R>& right,
276  Array<RES>& result, BinaryOperator op);
277 
278 // Transform left and right to a result using the binary operator.
279 // Result need not be a contiguous array.
280 template<typename L, typename R, typename RES, typename BinaryOperator>
281 void arrayTransform (const Array<L>& left, R right,
282  Array<RES>& result, BinaryOperator op);
283 
284 // Transform left and right to a result using the binary operator.
285 // Result need not be a contiguous array.
286 template<typename L, typename R, typename RES, typename BinaryOperator>
287 void arrayTransform (L left, const Array<R>& right,
288  Array<RES>& result, BinaryOperator op);
289 
290 // Transform array to a result using the unary operator.
291 // Result need not be a contiguous array.
292 template<typename T, typename RES, typename UnaryOperator>
293 void arrayTransform (const Array<T>& arr,
294  Array<RES>& result, UnaryOperator op);
295 
296 // Transform left and right to a result using the binary operator.
297 // The created and returned result array is contiguous.
298 template<typename T, typename BinaryOperator>
299 Array<T> arrayTransformResult (const Array<T>& left, const Array<T>& right,
300  BinaryOperator op);
301 
302 // Transform left and right to a result using the binary operator.
303 // The created and returned result array is contiguous.
304 template<typename T, typename BinaryOperator>
305 Array<T> arrayTransformResult (const Array<T>& left, T right, BinaryOperator op);
306 
307 // Transform left and right to a result using the binary operator.
308 // The created and returned result array is contiguous.
309 template<typename T, typename BinaryOperator>
310 Array<T> arrayTransformResult (T left, const Array<T>& right, BinaryOperator op);
311 
312 // Transform array to a result using the unary operator.
313 // The created and returned result array is contiguous.
314 template<typename T, typename UnaryOperator>
315 Array<T> arrayTransformResult (const Array<T>& arr, UnaryOperator op);
316 
317 // Transform left and right in place using the binary operator.
318 // The result is stored in the left array (useful for e.g. the += operation).
319 template<typename L, typename R, typename BinaryOperator>
320 inline void arrayTransformInPlace (Array<L>& left, const Array<R>& right,
321  BinaryOperator op)
322 {
323  if (left.contiguousStorage() && right.contiguousStorage()) {
324  transformInPlace (left.cbegin(), left.cend(), right.cbegin(), op);
325  } else {
326  transformInPlace (left.begin(), left.end(), right.begin(), op);
327  }
328 }
329 
330 // Transform left and right in place using the binary operator.
331 // The result is stored in the left array (useful for e.g. the += operation).
332 template<typename L, typename R, typename BinaryOperator>
333 inline void arrayTransformInPlace (Array<L>& left, R right, BinaryOperator op)
334 {
335  if (left.contiguousStorage()) {
336  myiptransform (left.cbegin(), left.cend(), right, op);
338  } else {
339  myiptransform (left.begin(), left.end(), right, op);
341  }
342 }
343 
344 // Transform the array in place using the unary operator.
345 // E.g. doing <src>arrayTransformInPlace(array, Sin<T>())</src> is faster than
346 // <src>array=sin(array)</src> as it does not need to create a temporary array.
347 template<typename T, typename UnaryOperator>
348 inline void arrayTransformInPlace (Array<T>& arr, UnaryOperator op)
349 {
350  if (arr.contiguousStorage()) {
351  transformInPlace (arr.cbegin(), arr.cend(), op);
352  } else {
353  transformInPlace (arr.begin(), arr.end(), op);
354  }
355 }
356 // </group>
357 
358 //
359 // Element by element arithmetic modifying left in-place. left and other
360 // must be conformant.
361 // <group>
362 template<class T> void operator+= (Array<T> &left, const Array<T> &other);
363 template<class T> void operator-= (Array<T> &left, const Array<T> &other);
364 template<class T> void operator*= (Array<T> &left, const Array<T> &other);
365 template<class T> void operator/= (Array<T> &left, const Array<T> &other);
366 template<class T> void operator%= (Array<T> &left, const Array<T> &other);
367 template<class T> void operator&= (Array<T> &left, const Array<T> &other);
368 template<class T> void operator|= (Array<T> &left, const Array<T> &other);
369 template<class T> void operator^= (Array<T> &left, const Array<T> &other);
370 // </group>
371 
372 //
373 // Element by element arithmetic modifying left in-place. The scalar "other"
374 // behaves as if it were a conformant Array to left filled with constant values.
375 // <group>
376 template<class T> void operator+= (Array<T> &left, const T &other);
377 template<class T> void operator-= (Array<T> &left, const T &other);
378 template<class T> void operator*= (Array<T> &left, const T &other);
379 template<class T> void operator/= (Array<T> &left, const T &other);
380 template<class T> void operator%= (Array<T> &left, const T &other);
381 template<class T> void operator&= (Array<T> &left, const T &other);
382 template<class T> void operator|= (Array<T> &left, const T &other);
383 template<class T> void operator^= (Array<T> &left, const T &other);
384 // </group>
385 
386 // Unary arithmetic operation.
387 //
388 // <group>
389 template<class T> Array<T> operator+(const Array<T> &a);
390 template<class T> Array<T> operator-(const Array<T> &a);
391 template<class T> Array<T> operator~(const Array<T> &a);
392 // </group>
393 
394 //
395 // Element by element arithmetic on two arrays, returning an array.
396 // <group>
397 template<class T>
398  Array<T> operator+ (const Array<T> &left, const Array<T> &right);
399 template<class T>
400  Array<T> operator- (const Array<T> &left, const Array<T> &right);
401 template<class T>
402  Array<T> operator* (const Array<T> &left, const Array<T> &right);
403 template<class T>
404  Array<T> operator/ (const Array<T> &left, const Array<T> &right);
405 template<class T>
406  Array<T> operator% (const Array<T> &left, const Array<T> &right);
407 template<class T>
408  Array<T> operator| (const Array<T> &left, const Array<T> &right);
409 template<class T>
410  Array<T> operator& (const Array<T> &left, const Array<T> &right);
411 template<class T>
412  Array<T> operator^ (const Array<T> &left, const Array<T> &right);
413 // </group>
414 
415 //
416 // Element by element arithmetic between an array and a scalar, returning
417 // an array.
418 // <group>
419 template<class T>
420  Array<T> operator+ (const Array<T> &left, const T &right);
421 template<class T>
422  Array<T> operator- (const Array<T> &left, const T &right);
423 template<class T>
424  Array<T> operator* (const Array<T> &left, const T &right);
425 template<class T>
426  Array<T> operator/ (const Array<T> &left, const T &right);
427 template<class T>
428  Array<T> operator% (const Array<T> &left, const T &right);
429 template<class T>
430  Array<T> operator| (const Array<T> &left, const T &right);
431 template<class T>
432  Array<T> operator& (const Array<T> &left, const T &right);
433 template<class T>
434  Array<T> operator^ (const Array<T> &left, const T &right);
435 // </group>
436 
437 //
438 // Element by element arithmetic between a scalar and an array, returning
439 // an array.
440 // <group>
441 template<class T>
442  Array<T> operator+ (const T &left, const Array<T> &right);
443 template<class T>
444  Array<T> operator- (const T &left, const Array<T> &right);
445 template<class T>
446  Array<T> operator* (const T &left, const Array<T> &right);
447 template<class T>
448  Array<T> operator/ (const T &left, const Array<T> &right);
449 template<class T>
450  Array<T> operator% (const T &left, const Array<T> &right);
451 template<class T>
452  Array<T> operator| (const T &left, const Array<T> &right);
453 template<class T>
454  Array<T> operator& (const T &left, const Array<T> &right);
455 template<class T>
456  Array<T> operator^ (const T &left, const Array<T> &right);
457 // </group>
458 
459 //
460 // Transcendental function that can be applied to essentially all numeric
461 // types. Works on an element-by-element basis.
462 // <group>
463 template<class T> Array<T> cos(const Array<T> &a);
464 template<class T> Array<T> cosh(const Array<T> &a);
465 template<class T> Array<T> exp(const Array<T> &a);
466 template<class T> Array<T> log(const Array<T> &a);
467 template<class T> Array<T> log10(const Array<T> &a);
468 template<class T> Array<T> pow(const Array<T> &a, const Array<T> &b);
469 template<class T> Array<T> pow(const T &a, const Array<T> &b);
470 template<class T> Array<T> sin(const Array<T> &a);
471 template<class T> Array<T> sinh(const Array<T> &a);
472 template<class T> Array<T> sqrt(const Array<T> &a);
473 // </group>
474 
475 //
476 // Transcendental function applied to the array on an element-by-element
477 // basis. Although a template function, this does not make sense for all
478 // numeric types.
479 // <group>
480 template<class T> Array<T> acos(const Array<T> &a);
481 template<class T> Array<T> asin(const Array<T> &a);
482 template<class T> Array<T> atan(const Array<T> &a);
483 template<class T> Array<T> atan2(const Array<T> &y, const Array<T> &x);
484 template<class T> Array<T> atan2(const T &y, const Array<T> &x);
485 template<class T> Array<T> atan2(const Array<T> &y, const T &x);
486 template<class T> Array<T> ceil(const Array<T> &a);
487 template<class T> Array<T> fabs(const Array<T> &a);
488 template<class T> Array<T> abs(const Array<T> &a);
489 template<class T> Array<T> floor(const Array<T> &a);
490 template<class T> Array<T> round(const Array<T> &a);
491 template<class T> Array<T> sign(const Array<T> &a);
492 template<class T> Array<T> fmod(const Array<T> &a, const Array<T> &b);
493 template<class T> Array<T> fmod(const T &a, const Array<T> &b);
494 template<class T> Array<T> fmod(const Array<T> &a, const T &b);
495 template<class T> Array<T> pow(const Array<T> &a, const Double &b);
496 template<class T> Array<T> tan(const Array<T> &a);
497 template<class T> Array<T> tanh(const Array<T> &a);
498 // N.B. fabs is deprecated. Use abs.
499 template<class T> Array<T> fabs(const Array<T> &a);
500 // </group>
501 
502 //
503 // <group>
504 // Find the minimum and maximum values of an array, including their locations.
505 template<class ScalarType>
506 void minMax(ScalarType &minVal, ScalarType &maxVal, IPosition &minPos,
507  IPosition &maxPos, const Array<ScalarType> &array);
508 // The array is searched at locations where the mask equals <src>valid</src>.
509 // (at least one such position must exist or an exception will be thrown).
510 // MaskType should be an Array of Bool.
511 template<class ScalarType>
512 void minMax(ScalarType &minVal, ScalarType &maxVal, IPosition &minPos,
513  IPosition &maxPos, const Array<ScalarType> &array,
514  const Array<Bool> &mask, Bool valid=True);
515 // The array * weight is searched
516 template<class ScalarType>
517 void minMaxMasked(ScalarType &minVal, ScalarType &maxVal, IPosition &minPos,
518  IPosition &maxPos, const Array<ScalarType> &array,
519  const Array<ScalarType> &weight);
520 // </group>
521 
522 //
523 // The "min" and "max" functions require that the type "T" have comparison
524 // operators.
525 // <group>
526 //
527 // This sets min and max to the minimum and maximum of the array to
528 // avoid having to do two passes with max() and min() separately.
529 template<class T> void minMax(T &min, T &max, const Array<T> &a);
530 //
531 // The minimum element of the array.
532 // Requires that the type "T" has comparison operators.
533 template<class T> T min(const Array<T> &a);
534 // The maximum element of the array.
535 // Requires that the type "T" has comparison operators.
536 template<class T> T max(const Array<T> &a);
537 
538 // "result" contains the maximum of "a" and "b" at each position. "result",
539 // "a", and "b" must be conformant.
540 template<class T> void max(Array<T> &result, const Array<T> &a,
541  const Array<T> &b);
542 // "result" contains the minimum of "a" and "b" at each position. "result",
543 // "a", and "b" must be conformant.
544 template<class T> void min(Array<T> &result, const Array<T> &a,
545  const Array<T> &b);
546 // Return an array that contains the maximum of "a" and "b" at each position.
547 // "a" and "b" must be conformant.
548 template<class T> Array<T> max(const Array<T> &a, const Array<T> &b);
549 template<class T> Array<T> max(const T &a, const Array<T> &b);
550 // Return an array that contains the minimum of "a" and "b" at each position.
551 // "a" and "b" must be conformant.
552 template<class T> Array<T> min(const Array<T> &a, const Array<T> &b);
553 
554 // "result" contains the maximum of "a" and "b" at each position. "result",
555 // and "a" must be conformant.
556 template<class T> void max(Array<T> &result, const Array<T> &a,
557  const T &b);
558 template<class T> inline void max(Array<T> &result, const T &a,
559  const Array<T> &b)
560  { max (result, b, a); }
561 // "result" contains the minimum of "a" and "b" at each position. "result",
562 // and "a" must be conformant.
563 template<class T> void min(Array<T> &result, const Array<T> &a,
564  const T &b);
565 template<class T> inline void min(Array<T> &result, const T &a,
566  const Array<T> &b)
567  { min (result, b, a); }
568 // Return an array that contains the maximum of "a" and "b" at each position.
569 template<class T> Array<T> max(const Array<T> &a, const T &b);
570 template<class T> inline Array<T> max(const T &a, const Array<T> &b)
571  { return max(b, a); }
572 // Return an array that contains the minimum of "a" and "b" at each position.
573 template<class T> Array<T> min(const Array<T> &a, const T &b);
574 template<class T> inline Array<T> min(const T &a, const Array<T> &b)
575  { return min(b, a); }
576 // </group>
577 
578 //
579 // Fills all elements of "array" with a sequence starting with "start"
580 // and incrementing by "inc" for each element. The first axis varies
581 // most rapidly.
582 template<class T> void indgen(Array<T> &a, T start, T inc);
583 //
584 // Fills all elements of "array" with a sequence starting with 0
585 // and ending with nelements() - 1. The first axis varies
586 // most rapidly.
587 template<class T> inline void indgen(Array<T> &a)
588  { indgen(a, T(0), T(1)); }
589 //
590 // Fills all elements of "array" with a sequence starting with start
591 // incremented by one for each position in the array. The first axis varies
592 // most rapidly.
593 template<class T> inline void indgen(Array<T> &a, T start)
594  { indgen(a, start, T(1)); }
595 
596 // Create a Vector of the given length and fill it with the start value
597 // incremented with <code>inc</code> for each element.
598 template<class T> inline Vector<T> indgen(uInt length, T start, T inc)
599 {
600  Vector<T> x(length);
601  indgen(x, start, inc);
602  return x;
603 }
604 
605 
606 // Sum of every element of the array.
607 template<class T> T sum(const Array<T> &a);
608 //
609 // Product of every element of the array. This could of course easily
610 // overflow.
611 template<class T> T product(const Array<T> &a);
612 
613 //
614 // The mean of "a" is the sum of all elements of "a" divided by the number
615 // of elements of "a".
616 template<class T> T mean(const Array<T> &a);
617 
618 //
619 // The variance of "a" is the sum of (a(i) - mean(a))**2/(a.nelements() - 1).
620 // N.B. N-1, not N in the denominator).
621 template<class T> T variance(const Array<T> &a);
622 //
623 // The variance of "a" is the sum of (a(i) - mean(a))**2/(a.nelements() - 1).
624 // N.B. N-1, not N in the denominator).
625 // Rather than using a computed mean, use the supplied value.
626 template<class T> T variance(const Array<T> &a, T mean);
627 
628 //
629 // The standard deviation of "a" is the square root of its variance.
630 template<class T> T stddev(const Array<T> &a);
631 //
632 // The standard deviation of "a" is the square root of its variance.
633 // Rather than using a computed mean, use the supplied value.
634 template<class T> T stddev(const Array<T> &a, T mean);
635 
636 //
637 // The average deviation of "a" is the sum of abs(a(i) - mean(a))/N. (N.B.
638 // N, not N-1 in the denominator).
639 template<class T> T avdev(const Array<T> &a);
640 //
641 // The average deviation of "a" is the sum of abs(a(i) - mean(a))/N. (N.B.
642 // N, not N-1 in the denominator).
643 // Rather than using a computed mean, use the supplied value.
644 template<class T> T avdev(const Array<T> &a,T mean);
645 
646 //
647 // The root-mean-square of "a" is the sqrt of sum(a*a)/N.
648 template<class T> T rms(const Array<T> &a);
649 
650 
651 // The median of "a" is a(n/2).
652 // If a has an even number of elements and the switch takeEvenMean is set,
653 // the median is 0.5*(a(n/2) + a((n+1)/2)).
654 // According to Numerical Recipes (2nd edition) it makes little sense to take
655 // the mean if the array is large enough (> 100 elements). Therefore
656 // the default for takeEvenMean is False if the array has > 100 elements,
657 // otherwise it is True.
658 // <br>If "sorted"==True we assume the data is already sorted and we
659 // compute the median directly. Otherwise the function GenSort::kthLargest
660 // is used to find the median (kthLargest is about 6 times faster
661 // than a full quicksort).
662 // <br>Finding the median means that the array has to be (partially)
663 // sorted. By default a copy will be made, but if "inPlace" is in effect,
664 // the data themselves will be sorted. That should only be used if the
665 // data are used not thereafter.
666 // <note>The function kthLargest in class GenSortIndirect can be used to
667 // obtain the index of the median in an array. </note>
668 // <group>
669 template<class T> T median(const Array<T> &a, Block<T> &tmp, Bool sorted,
670  Bool takeEvenMean, Bool inPlace=False);
671 template<class T> T median(const Array<T> &a, Bool sorted, Bool takeEvenMean,
672  Bool inPlace=False)
673  { Block<T> tmp; return median (a, tmp, sorted, takeEvenMean, inPlace); }
674 template<class T> inline T median(const Array<T> &a, Bool sorted)
675  { return median (a, sorted, (a.nelements() <= 100), False); }
676 template<class T> inline T median(const Array<T> &a)
677  { return median (a, False, (a.nelements() <= 100), False); }
678 template<class T> inline T medianInPlace(const Array<T> &a, Bool sorted=False)
679  { return median (a, sorted, (a.nelements() <= 100), True); }
680 // </group>
681 
682 // The median absolute deviation from the median. Interface is as for
683 // the median functions
684 // <group>
685 template<class T> T madfm(const Array<T> &a, Block<T> &tmp, Bool sorted,
686  Bool takeEvenMean, Bool inPlace = False);
687 template<class T> T madfm(const Array<T> &a, Bool sorted, Bool takeEvenMean,
688  Bool inPlace=False)
689  { Block<T> tmp; return madfm(a, tmp, sorted, takeEvenMean, inPlace); }
690 template<class T> inline T madfm(const Array<T> &a, Bool sorted)
691  { return madfm(a, sorted, (a.nelements() <= 100), False); }
692 template<class T> inline T madfm(const Array<T> &a)
693  { return madfm(a, False, (a.nelements() <= 100), False); }
694 template<class T> inline T madfmInPlace(const Array<T> &a, Bool sorted=False)
695  { return madfm(a, sorted, (a.nelements() <= 100), True); }
696 // </group>
697 
698 // Return the fractile of an array.
699 // It returns the value at the given fraction of the array.
700 // A fraction of 0.5 is the same as the median, be it that no mean of
701 // the two middle elements is taken if the array has an even nr of elements.
702 // It uses kthLargest if the array is not sorted yet.
703 // <note>The function kthLargest in class GenSortIndirect can be used to
704 // obtain the index of the fractile in an array. </note>
705 template<class T> T fractile(const Array<T> &a, Block<T> &tmp, Float fraction,
706  Bool sorted=False, Bool inPlace=False);
707 template<class T> T fractile(const Array<T> &a, Float fraction,
708  Bool sorted=False, Bool inPlace=False)
709  { Block<T> tmp; return fractile (a, tmp, fraction, sorted, inPlace); }
710 
711 // Return the inter-fractile range of an array.
712 // This is the full range between the bottom and the top fraction.
713 // <group>
714 template<class T> T interFractileRange(const Array<T> &a, Block<T> &tmp,
715  Float fraction,
716  Bool sorted=False, Bool inPlace=False);
717 template<class T> T interFractileRange(const Array<T> &a, Float fraction,
718  Bool sorted=False, Bool inPlace=False)
719  { Block<T> tmp; return interFractileRange(a, tmp, fraction, sorted, inPlace); }
720 // </group>
721 
722 // Return the inter-hexile range of an array.
723 // This is the full range between the bottom sixth and the top sixth
724 // of ordered array values. "The semi-interhexile range is very nearly
725 // equal to the rms for a Gaussian distribution, but it is much less
726 // sensitive to the tails of extended distributions." (Condon et al
727 // 1998)
728 // <group>
729 template<class T> T interHexileRange(const Array<T> &a, Block<T> &tmp,
730  Bool sorted=False, Bool inPlace=False)
731  { return interFractileRange(a, tmp, 1./6., sorted, inPlace); }
732 template<class T> T interHexileRange(const Array<T> &a, Bool sorted=False,
733  Bool inPlace=False)
734  { return interFractileRange(a, 1./6., sorted, inPlace); }
735 // </group>
736 
737 // Return the inter-quartile range of an array.
738 // This is the full range between the bottom quarter and the top
739 // quarter of ordered array values.
740 // <group>
741 template<class T> T interQuartileRange(const Array<T> &a, Block<T> &tmp,
742  Bool sorted=False, Bool inPlace=False)
743  { return interFractileRange(a, tmp, 0.25, sorted, inPlace); }
744 template<class T> T interQuartileRange(const Array<T> &a, Bool sorted=False,
745  Bool inPlace=False)
746  { return interFractileRange(a, 0.25, sorted, inPlace); }
747 // </group>
748 
749 
750 // Methods for element-by-element scaling of complex and real.
751 // Note that Complex and DComplex are typedefs for std::complex.
752 //<group>
753 template<typename T>
754 void operator*= (Array<std::complex<T> > &left, const Array<T> &other);
755 template<typename T>
756 void operator*= (Array<std::complex<T> > &left, const T &other);
757 template<typename T>
758 void operator/= (Array<std::complex<T> > &left, const Array<T> &other);
759 template<typename T>
760 void operator/= (Array<std::complex<T> > &left, const T &other);
761 template<typename T>
762 Array<std::complex<T> > operator* (const Array<std::complex<T> > &left,
763  const Array<T> &right);
764 template<typename T>
765 Array<std::complex<T> > operator* (const Array<std::complex<T> > &left,
766  const T &right);
767 template<typename T>
768 Array<std::complex<T> > operator* (const std::complex<T> &left,
769  const Array<T> &right);
770 template<typename T>
771 Array<std::complex<T> > operator/ (const Array<std::complex<T> > &left,
772  const Array<T> &right);
773 template<typename T>
774 Array<std::complex<T> > operator/ (const Array<std::complex<T> > &left,
775  const T &right);
776 template<typename T>
777 Array<std::complex<T> > operator/ (const std::complex<T> &left,
778  const Array<T> &right);
779 // </group>
780 
781 // Returns the complex conjugate of a complex array.
782 //<group>
783 Array<Complex> conj(const Array<Complex> &carray);
784 Array<DComplex> conj(const Array<DComplex> &carray);
785 // Modifies rarray in place. rarray must be conformant.
786 void conj(Array<Complex> &rarray, const Array<Complex> &carray);
787 void conj(Array<DComplex> &rarray, const Array<DComplex> &carray);
788 //# The following are implemented to make the compiler find the right conversion
789 //# more often.
790 Matrix<Complex> conj(const Matrix<Complex> &carray);
791 Matrix<DComplex> conj(const Matrix<DComplex> &carray);
792 //</group>
793 
794 // Form an array of complex numbers from the given real arrays.
795 // Note that Complex and DComplex are simply typedefs for std::complex<float>
796 // and std::complex<double>, so the result is in fact one of these types.
797 template<typename T>
798 Array<std::complex<T> > makeComplex(const Array<T> &real, const Array<T>& imag);
799 
800 // Set the real part of the left complex array to the right real array.
801 template<typename L, typename R>
802 void setReal(Array<L> &carray, const Array<R> &rarray);
803 
804 // Set the imaginary part of the left complex array to right real array.
805 template<typename R, typename L>
806 void setImag(Array<R> &carray, const Array<L> &rarray);
807 
808 // Extracts the real part of a complex array into an array of floats.
809 // <group>
810 Array<Float> real(const Array<Complex> &carray);
811 Array<Double> real(const Array<DComplex> &carray);
812 // Modifies rarray in place. rarray must be conformant.
813 void real(Array<Float> &rarray, const Array<Complex> &carray);
814 void real(Array<Double> &rarray, const Array<DComplex> &carray);
815 // </group>
816 
817 //
818 // Extracts the imaginary part of a complex array into an array of floats.
819 // <group>
820 Array<Float> imag(const Array<Complex> &carray);
821 Array<Double> imag(const Array<DComplex> &carray);
822 // Modifies rarray in place. rarray must be conformant.
823 void imag(Array<Float> &rarray, const Array<Complex> &carray);
824 void imag(Array<Double> &rarray, const Array<DComplex> &carray);
825 // </group>
826 
827 //
828 // Extracts the amplitude (i.e. sqrt(re*re + im*im)) from an array
829 // of complex numbers. N.B. this is presently called "fabs" for a single
830 // complex number.
831 // <group>
832 Array<Float> amplitude(const Array<Complex> &carray);
834 // Modifies rarray in place. rarray must be conformant.
835 void amplitude(Array<Float> &rarray, const Array<Complex> &carray);
836 void amplitude(Array<Double> &rarray, const Array<DComplex> &carray);
837 // </group>
838 
839 //
840 // Extracts the phase (i.e. atan2(im, re)) from an array
841 // of complex numbers. N.B. this is presently called "arg"
842 // for a single complex number.
843 // <group>
844 Array<Float> phase(const Array<Complex> &carray);
845 Array<Double> phase(const Array<DComplex> &carray);
846 // Modifies rarray in place. rarray must be conformant.
847 void phase(Array<Float> &rarray, const Array<Complex> &carray);
848 void phase(Array<Double> &rarray, const Array<DComplex> &carray);
849 // </group>
850 
851 // Copy an array of complex into an array of real,imaginary pairs. The
852 // first axis of the real array becomes twice as long as the complex array.
853 // In the future versions which work by reference will be available; presently
854 // a copy is made.
855 Array<Float> ComplexToReal(const Array<Complex> &carray);
856 Array<Double> ComplexToReal(const Array<DComplex> &carray);
857 // Modify the array "rarray" in place. "rarray" must be the correct shape.
858 // <group>
859 void ComplexToReal(Array<Float> &rarray, const Array<Complex> &carray);
860 void ComplexToReal(Array<Double> &rarray, const Array<DComplex> &carray);
861 // </group>
862 
863 // Copy an array of real,imaginary pairs into a complex array. The first axis
864 // must have an even length.
865 // In the future versions which work by reference will be available; presently
866 // a copy is made.
867 Array<Complex> RealToComplex(const Array<Float> &rarray);
868 Array<DComplex> RealToComplex(const Array<Double> &rarray);
869 // Modify the array "carray" in place. "carray" must be the correct shape.
870 // <group>
871 void RealToComplex(Array<Complex> &carray, const Array<Float> &rarray);
872 void RealToComplex(Array<DComplex> &carray, const Array<Double> &rarray);
873 // </group>
874 
875 // Make a copy of an array of a different type; for example make an array
876 // of doubles from an array of floats. Arrays to and from must be conformant
877 // (same shape). Also, it must be possible to convert a scalar of type U
878 // to type T.
879 template<class T, class U> void convertArray(Array<T> &to,
880  const Array<U> &from);
881 
882 
883 // Returns an array where every element is squared.
884 template<class T> Array<T> square(const Array<T> &val);
885 
886 // Returns an array where every element is cubed.
887 template<class T> Array<T> cube(const Array<T> &val);
888 
889 // Expand the values of an array. The arrays must have the same dimensionality.
890 // The length of each axis in the input array must be <= the length of the
891 // corresponding axis in the output array and divide evenly.
892 // For each axis <src>mult</src> is set to output/input.
893 // <br>The <src>alternate</src> argument determines how the values are expanded.
894 // If a row contains values '1 2 3', they can be expanded "linearly"
895 // as '1 1 2 2 3 3' or alternately as '1 2 3 1 2 3'
896 // This choice can be made for each axis; a value 0 means linearly,
897 // another value means alternately. If the length of alternate is less than
898 // the dimensionality of the arrays, the missing ones default to 0.
899 template<class T> void expandArray (Array<T>& out, const Array<T>& in,
900  const IPosition& alternate=IPosition());
901 // Helper function for expandArray using recursion for each axis.
902 template<class T>
903 T* expandRecursive (int axis, const IPosition& shp, const IPosition& mult,
904  const IPosition& inSteps,
905  const T* in, T* out, const IPosition& alternate);
906 // Check array shapes for expandArray. It returns the alternate argument,
907 // where possibly missing values are appended (as 0).
908 IPosition checkExpandArray (IPosition& mult,
909  const IPosition& inShape,
910  const IPosition& outShape,
911  const IPosition& alternate);
912 
913 
914 // </group>
915 
916 
917 } //# NAMESPACE CASACORE - END
918 
919 #ifndef CASACORE_NO_AUTO_TEMPLATES
920 #include <casacore/casa/Arrays/ArrayMath.tcc>
921 #endif //# CASACORE_NO_AUTO_TEMPLATES
922 #endif
Array< T > max(const T &a, const Array< T > &b)
Definition: ArrayMath.h:570
A Vector of integers, for indexing into Array<T> objects.
Definition: IPosition.h:119
LatticeExprNode log10(const LatticeExprNode &expr)
A 1-D Specialization of the Array class.
Definition: ArrayIO.h:45
T interHexileRange(const Array< T > &a, Bool sorted=False, Bool inPlace=False)
Definition: ArrayMath.h:732
LatticeExprNode log(const LatticeExprNode &expr)
Non-templated base class for templated Array class.
Definition: ArrayBase.h:74
T interHexileRange(const Array< T > &a, Block< T > &tmp, Bool sorted=False, Bool inPlace=False)
Return the inter-hexile range of an array.
Definition: ArrayMath.h:729
T median(const Array< T > &a, Bool sorted, Bool takeEvenMean, Bool inPlace=False)
Definition: ArrayMath.h:671
LatticeExprNode median(const LatticeExprNode &expr)
LatticeExprNode operator/(const LatticeExprNode &left, const LatticeExprNode &right)
LatticeExprNode mask(const LatticeExprNode &expr)
This function returns the mask of the given expression.
T product(const TableVector< T > &tv)
Definition: TabVecMath.h:385
const IPosition & shape() const
The length of each axis.
Definition: ArrayBase.h:121
LatticeExprNode imag(const LatticeExprNode &expr)
void myiptransform(_InputIterator1 __first1, _InputIterator1 __last1, T right, _BinaryOperation __binary_op)
sequence OP= scalar
Definition: ArrayMath.h:174
T medianInPlace(const Array< T > &a, Bool sorted=False)
Definition: ArrayMath.h:678
void myrtransform(_InputIterator1 __first1, _InputIterator1 __last1, _OutputIterator __result, T right, _BinaryOperation __binary_op)
sequence = sequence OP scalar
Definition: ArrayMath.h:163
TableExprNode array(const TableExprNode &values, const TableExprNodeSet &shape)
Create an array of the given shape and fill it with the values.
Definition: ExprNode.h:1986
LatticeExprNode sum(const LatticeExprNode &expr)
LatticeExprNode max(const LatticeExprNode &left, const LatticeExprNode &right)
LatticeExprNode operator%(const LatticeExprNode &left, const LatticeExprNode &right)
T fractile(const Array< T > &a, Float fraction, Bool sorted=False, Bool inPlace=False)
Definition: ArrayMath.h:707
TableExprNode phase(const TableExprNode &node)
The phase (i.e.
Definition: ExprNode.h:1510
Array< T > min(const T &a, const Array< T > &b)
Definition: ArrayMath.h:574
void max(Array< T > &result, const T &a, const Array< T > &b)
Definition: ArrayMath.h:558
void arrayTransformInPlace(Array< L > &left, const Array< R > &right, BinaryOperator op)
Transform left and right in place using the binary operator.
Definition: ArrayMath.h:320
LatticeExprNode fractile(const LatticeExprNode &expr, const LatticeExprNode &fraction)
Determine the value of the element at the part fraction from the beginning of the given lattice...
LatticeExprNode exp(const LatticeExprNode &expr)
iterator begin()
Get the begin iterator object for any array.
Definition: Array.h:855
LatticeExprNode floor(const LatticeExprNode &expr)
LatticeExprNode cos(const LatticeExprNode &expr)
contiter cbegin()
Get the begin iterator object for a contiguous array.
Definition: Array.h:867
T interFractileRange(const Array< T > &a, Float fraction, Bool sorted=False, Bool inPlace=False)
Definition: ArrayMath.h:717
contiter cend()
Definition: Array.h:871
void arrayTransformInPlace(Array< L > &left, R right, BinaryOperator op)
Transform left and right in place using the binary operator.
Definition: ArrayMath.h:333
LatticeExprNode conj(const LatticeExprNode &expr)
void arrayContTransform(const Array< L > &left, R right, Array< RES > &result, BinaryOperator op)
Transform left and right to a result using the binary operator.
Definition: ArrayMath.h:221
Float pow(Float f1, Float f2)
Definition: math.h:90
LatticeExprNode tanh(const LatticeExprNode &expr)
LatticeExprNode min(const LatticeExprNode &left, const LatticeExprNode &right)
void transformInPlace(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, BinaryOperator op)
Define a function to do a binary transform in place.
Definition: Functors.h:44
LatticeExprNode avdev(const LatticeExprNode &expr)
size_t nelements() const
How many elements does this array have? Product of all axis lengths.
Definition: ArrayBase.h:99
Bool contiguousStorage() const
Are the array data contiguous? If they are not contiguous, getStorage (see below) needs to make a cop...
Definition: ArrayBase.h:112
Bool isEqual(const IPosition &other) const
Element-by-element comparison for equality.
void arrayContTransform(L left, const Array< R > &right, Array< RES > &result, BinaryOperator op)
Transform left and right to a result using the binary operator.
Definition: ArrayMath.h:241
double Double
Definition: aipstype.h:52
LatticeExprNode abs(const LatticeExprNode &expr)
Numerical 1-argument functions which result in a real number regardless of input expression type...
LatticeExprNode length(const LatticeExprNode &expr, const LatticeExprNode &axis)
2-argument function to get the length of an axis.
T interQuartileRange(const Array< T > &a, Bool sorted=False, Bool inPlace=False)
Definition: ArrayMath.h:744
LatticeExprNode sign(const LatticeExprNode &expr)
LatticeExprNode sqrt(const LatticeExprNode &expr)
#define DebugAssert(expr, exception)
Definition: Assert.h:185
LatticeExprNode tan(const LatticeExprNode &expr)
void indgen(TableVector< T > &tv, Int start, Int inc)
Definition: TabVecMath.h:400
LatticeExprNode atan(const LatticeExprNode &expr)
void arrayContTransform(const Array< T > &arr, Array< RES > &result, UnaryOperator op)
Transform array to a result using the unary operator.
Definition: ArrayMath.h:261
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:39
void minMax(T &min, T &max, const TableVector< T > &tv)
Definition: TabVecMath.h:390
void min(Array< T > &result, const T &a, const Array< T > &b)
Definition: ArrayMath.h:565
TableExprNode cube(const TableExprNode &node)
Definition: ExprNode.h:1417
LatticeExprNode stddev(const LatticeExprNode &expr)
LatticeExprNode round(const LatticeExprNode &expr)
float Float
Definition: aipstype.h:51
Vector< T > indgen(uInt length, T start, T inc)
Create a Vector of the given length and fill it with the start value incremented with inc for each el...
Definition: ArrayMath.h:598
const Bool False
Definition: aipstype.h:41
template <class T, class U> class vector;
Definition: Array.h:169
TableExprNode operator|(const TableExprNode &left, const TableExprNode &right)
Definition: ExprNode.h:1252
LatticeExprNode atan2(const LatticeExprNode &left, const LatticeExprNode &right)
Numerical 2-argument functions.
simple 1-D array
Definition: ArrayIO.h:47
LatticeExprNode operator+(const LatticeExprNode &expr)
Global functions operating on a LatticeExprNode.
LatticeExprNode fmod(const LatticeExprNode &left, const LatticeExprNode &right)
iterator end()
Definition: Array.h:859
T madfmInPlace(const Array< T > &a, Bool sorted=False)
Definition: ArrayMath.h:694
Base class for all Casacore library errors.
Definition: Error.h:135
LatticeExprNode asin(const LatticeExprNode &expr)
LatticeExprNode mean(const LatticeExprNode &expr)
void arrayTransformInPlace(Array< T > &arr, UnaryOperator op)
Transform the array in place using the unary operator.
Definition: ArrayMath.h:348
void myltransform(_InputIterator1 __first1, _InputIterator1 __last1, _OutputIterator __result, T left, _BinaryOperation __binary_op)
The myxtransform functions are defined to avoid a bug in g++-4.3.
Definition: ArrayMath.h:152
TableExprNode rms(const TableExprNode &array)
Definition: ExprNode.h:1737
LatticeExprNode sinh(const LatticeExprNode &expr)
LatticeExprNode acos(const LatticeExprNode &expr)
TableExprNode square(const TableExprNode &node)
Definition: ExprNode.h:1412
LatticeExprNode operator-(const LatticeExprNode &expr)
LatticeExprNode operator^(const LatticeExprNode &left, const LatticeExprNode &right)
TableExprNode operator &(const TableExprNode &left, const TableExprNode &right)
Definition: ExprNode.h:1247
void arrayContTransform(const Array< L > &left, const Array< R > &right, Array< RES > &result, BinaryOperator op)
Functions to apply a binary or unary operator to arrays.
Definition: ArrayMath.h:205
void indgen(Array< T > &a)
Fills all elements of "array" with a sequence starting with 0 and ending with nelements() - 1...
Definition: ArrayMath.h:587
void checkArrayShapes(const ArrayBase &left, const ArrayBase &right, const char *name)
Definition: ArrayMath.h:187
LatticeExprNode variance(const LatticeExprNode &expr)
T madfm(const Array< T > &a, Bool sorted, Bool takeEvenMean, Bool inPlace=False)
Definition: ArrayMath.h:687
LatticeExprNode ceil(const LatticeExprNode &expr)
T interQuartileRange(const Array< T > &a, Block< T > &tmp, Bool sorted=False, Bool inPlace=False)
Return the inter-quartile range of an array.
Definition: ArrayMath.h:741
void indgen(Array< T > &a, T start)
Fills all elements of "array" with a sequence starting with start incremented by one for each positio...
Definition: ArrayMath.h:593
const Bool True
Definition: aipstype.h:40
this file contains all the compiler specific defines
Definition: mainpage.dox:28
MVBaseline operator*(const RotMatrix &left, const MVBaseline &right)
Rotate a Baseline vector with rotation matrix and other multiplications.
LatticeExprNode cosh(const LatticeExprNode &expr)
LatticeExprNode real(const LatticeExprNode &expr)
LatticeExprNode sin(const LatticeExprNode &expr)
Numerical 1-argument functions.
unsigned int uInt
Definition: aipstype.h:48
TableExprNode amplitude(const TableExprNode &node)
The amplitude (i.e.
Definition: ExprNode.h:1502