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