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
38namespace 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">
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.
185template<typename L, typename AllocL, typename R, typename AllocR, typename RES, typename AllocRES, typename BinaryOperator>
186inline 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.
201template<typename L, typename AllocL, typename R, typename RES, typename AllocRES, typename BinaryOperator>
202inline 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.
221template<typename L, typename R, typename AllocR, typename RES, typename AllocRES, typename BinaryOperator>
222inline 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.
241template<typename T, typename Alloc, typename RES, typename AllocRES, typename UnaryOperator>
242inline 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.
255template<typename L, typename R, typename RES, typename BinaryOperator, typename AllocL, typename AllocR, typename AllocRES>
256void 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.
261template<typename L, typename R, typename RES, typename BinaryOperator, typename Alloc, typename AllocRES>
262void 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.
267template<typename L, typename R, typename RES, typename BinaryOperator, typename Alloc, typename AllocRES>
268void 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.
273template<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.
279template<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.
285template<typename T, typename BinaryOperator, typename Alloc>
286Array<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.
290template<typename T, typename BinaryOperator, typename Alloc>
291Array<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.
295template<typename T, typename UnaryOperator, typename Alloc>
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).
300template<typename L, typename R, typename BinaryOperator, typename AllocL, typename AllocR>
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).
313template<typename L, typename R, typename BinaryOperator, typename Alloc>
314inline 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.
328template<typename T, typename UnaryOperator, typename Alloc>
329inline 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>
343template<typename T, typename Alloc> void operator+= (Array<T, Alloc> &left, const Array<T, Alloc> &other);
344template<typename T, typename Alloc> void operator-= (Array<T, Alloc> &left, const Array<T, Alloc> &other);
345template<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
351template<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}
356template<typename T, typename Alloc> void operator%= (Array<T, Alloc> &left, const Array<T, Alloc> &other);
357template<typename T, typename Alloc> void operator&= (Array<T, Alloc> &left, const Array<T, Alloc> &other);
358template<typename T, typename Alloc> void operator|= (Array<T, Alloc> &left, const Array<T, Alloc> &other);
359template<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>
366template<typename T, typename Alloc> void operator+= (Array<T, Alloc> &left, const T &other);
367template<typename T, typename Alloc> void operator-= (Array<T, Alloc> &left, const T &other);
368template<typename T, typename Alloc> void operator*= (Array<T, Alloc> &left, const T &other)
369{
370 arrayTransformInPlace (left, other, std::multiplies<T>());
371}
372template<typename T, typename Alloc> void operator/= (Array<T, Alloc> &left, const T &other)
373{
374 arrayTransformInPlace (left, other, std::divides<T>());
375}
376template<typename T, typename Alloc> void operator%= (Array<T, Alloc> &left, const T &other);
377template<typename T, typename Alloc> void operator&= (Array<T, Alloc> &left, const T &other);
378template<typename T, typename Alloc> void operator|= (Array<T, Alloc> &left, const T &other);
379template<typename T, typename Alloc> void operator^= (Array<T, Alloc> &left, const T &other);
380// </group>
381
382// Unary arithmetic operation.
383//
384// <group>
385template<typename T, typename Alloc> Array<T, Alloc> operator+(const Array<T, Alloc> &a);
386template<typename T, typename Alloc> Array<T, Alloc> operator-(const Array<T, Alloc> &a);
387template<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>
393template<typename T, typename Alloc>
395template<typename T, typename Alloc>
397template<typename T, typename Alloc>
399{
400 checkArrayShapes (left, right, "*");
401 return arrayTransformResult (left, right, std::multiplies<T>());
402}
403template<typename T, typename Alloc>
405template<typename T, typename Alloc>
407template<typename T, typename Alloc>
409template<typename T, typename Alloc>
411template<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>
419template<typename T, typename Alloc>
420 Array<T, Alloc> operator+ (const Array<T, Alloc> &left, const T &right);
421template<typename T, typename Alloc>
422 Array<T, Alloc> operator- (const Array<T, Alloc> &left, const T &right);
423template<class T, typename Alloc>
424Array<T, Alloc> operator* (const Array<T, Alloc> &left, const T &right)
425{
426 return arrayTransformResult (left, right, std::multiplies<T>());
427}
428template<typename T, typename Alloc>
429 Array<T, Alloc> operator/ (const Array<T, Alloc> &left, const T &right);
430template<typename T, typename Alloc>
431 Array<T, Alloc> operator% (const Array<T, Alloc> &left, const T &right);
432template<typename T, typename Alloc>
433 Array<T, Alloc> operator| (const Array<T, Alloc> &left, const T &right);
434template<typename T, typename Alloc>
435 Array<T, Alloc> operator& (const Array<T, Alloc> &left, const T &right);
436template<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>
444template<typename T, typename Alloc>
445 Array<T, Alloc> operator+ (const T &left, const Array<T, Alloc> &right);
446template<typename T, typename Alloc>
447 Array<T, Alloc> operator- (const T &left, const Array<T, Alloc> &right);
448template<class T, typename Alloc>
449Array<T, Alloc> operator* (const T &left, const Array<T, Alloc> &right)
450{
451 return arrayTransformResult (left, right, std::multiplies<T>());
452}
453
454template<typename T, typename Alloc>
455 Array<T, Alloc> operator/ (const T &left, const Array<T, Alloc> &right);
456template<typename T, typename Alloc>
457 Array<T, Alloc> operator% (const T &left, const Array<T, Alloc> &right);
458template<typename T, typename Alloc>
459 Array<T, Alloc> operator| (const T &left, const Array<T, Alloc> &right);
460template<typename T, typename Alloc>
461 Array<T, Alloc> operator& (const T &left, const Array<T, Alloc> &right);
462template<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>
470template<typename T, typename Alloc> Array<T, Alloc> cos(const Array<T, Alloc> &a);
471template<typename T, typename Alloc> Array<T, Alloc> cosh(const Array<T, Alloc> &a);
472template<typename T, typename Alloc> Array<T, Alloc> exp(const Array<T, Alloc> &a);
473template<typename T, typename Alloc> Array<T, Alloc> log(const Array<T, Alloc> &a);
474template<typename T, typename Alloc> Array<T, Alloc> log10(const Array<T, Alloc> &a);
475template<typename T, typename Alloc> Array<T, Alloc> pow(const Array<T, Alloc> &a, const Array<T, Alloc> &b);
476template<typename T, typename Alloc> Array<T, Alloc> pow(const T &a, const Array<T, Alloc> &b);
477template<typename T, typename Alloc> Array<T, Alloc> sin(const Array<T, Alloc> &a);
478template<typename T, typename Alloc> Array<T, Alloc> sinh(const Array<T, Alloc> &a);
479template<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>
487template<typename T, typename Alloc> Array<T, Alloc> acos(const Array<T, Alloc> &a);
488template<typename T, typename Alloc> Array<T, Alloc> asin(const Array<T, Alloc> &a);
489template<typename T, typename Alloc> Array<T, Alloc> atan(const Array<T, Alloc> &a);
490template<typename T, typename Alloc> Array<T, Alloc> atan2(const Array<T, Alloc> &y, const Array<T, Alloc> &x);
491template<typename T, typename Alloc> Array<T, Alloc> atan2(const T &y, const Array<T, Alloc> &x);
492template<typename T, typename Alloc> Array<T, Alloc> atan2(const Array<T, Alloc> &y, const T &x);
493template<typename T, typename Alloc> Array<T, Alloc> ceil(const Array<T, Alloc> &a);
494template<typename T, typename Alloc> Array<T, Alloc> fabs(const Array<T, Alloc> &a);
495template<typename T, typename Alloc> Array<T, Alloc> abs(const Array<T, Alloc> &a);
496template<typename T, typename Alloc> Array<T, Alloc> floor(const Array<T, Alloc> &a);
497template<typename T, typename Alloc> Array<T, Alloc> round(const Array<T, Alloc> &a);
498template<typename T, typename Alloc> Array<T, Alloc> sign(const Array<T, Alloc> &a);
499template<typename T, typename Alloc> Array<T, Alloc> fmod(const Array<T, Alloc> &a, const Array<T, Alloc> &b);
500template<typename T, typename Alloc> Array<T, Alloc> fmod(const T &a, const Array<T, Alloc> &b);
501template<typename T, typename Alloc> Array<T, Alloc> fmod(const Array<T, Alloc> &a, const T &b);
502template<typename T, typename Alloc> Array<T, Alloc> floormod(const Array<T, Alloc> &a, const Array<T, Alloc> &b);
503template<typename T, typename Alloc> Array<T, Alloc> floormod(const T &a, const Array<T, Alloc> &b);
504template<typename T, typename Alloc> Array<T, Alloc> floormod(const Array<T, Alloc> &a, const T &b);
505template<typename T, typename Alloc> Array<T, Alloc> pow(const Array<T, Alloc> &a, const double &b);
506template<typename T, typename Alloc> Array<T, Alloc> tan(const Array<T, Alloc> &a);
507template<typename T, typename Alloc> Array<T, Alloc> tanh(const Array<T, Alloc> &a);
508// N.B. fabs is deprecated. Use abs.
509template<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.
515template<typename ScalarType, typename Alloc>
516void minMax(ScalarType &minVal, ScalarType &maxVal, IPosition &minPos,
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.
521template<typename ScalarType, typename Alloc>
522void minMax(ScalarType &minVal, ScalarType &maxVal, IPosition &minPos,
524 const Array<bool> &mask, bool valid=true);
525// The array * weight is searched
526template<typename ScalarType, typename Alloc>
527void minMaxMasked(ScalarType &minVal, ScalarType &maxVal, IPosition &minPos,
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.
539template<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.
543template<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.
546template<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.
550template<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.
554template<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.
558template<typename T, typename Alloc> Array<T, Alloc> max(const Array<T, Alloc> &a, const Array<T, Alloc> &b);
559template<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.
562template<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.
566template<typename T, typename Alloc> void max(Array<T, Alloc> &result, const Array<T, Alloc> &a,
567 const T &b);
568template<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.
573template<typename T, typename Alloc> void min(Array<T, Alloc> &result, const Array<T, Alloc> &a,
574 const T &b);
575template<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.
579template<typename T, typename Alloc> Array<T, Alloc> max(const Array<T, Alloc> &a, const T &b);
580template<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.
583template<typename T, typename Alloc> Array<T, Alloc> min(const Array<T, Alloc> &a, const T &b);
584template<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.
592template<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.
597template<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.
603template<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.
608template<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.
617template<typename T, typename Alloc> T sum(const Array<T, Alloc> &a);
618//
619// Sum the square of every element of the array.
620template<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.
624template<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".
629template<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.
638template<typename T, typename Alloc> T variance(const Array<T, Alloc> &a);
639template<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.
641template<typename T, typename Alloc> T variance(const Array<T, Alloc> &a, T mean);
642template<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.
645template<typename T, typename Alloc> T stddev(const Array<T, Alloc> &a);
646template<typename T, typename Alloc> T pstddev(const Array<T, Alloc> &a, size_t ddof=0);
647template<typename T, typename Alloc> T stddev(const Array<T, Alloc> &a, T mean);
648template<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).
653template<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.
658template<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.
662template<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
684template<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
687template<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); }
690template<typename T, typename Alloc> inline T median(const Array<T, Alloc> &a, bool sorted)
691 { return median (a, sorted, (a.nelements() <= 100), false); }
692template<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
695template<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
703template<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
706template<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); }
709template<typename T, typename Alloc> inline T madfm(const Array<T, Alloc> &a, bool sorted)
710 { return madfm(a, sorted, (a.nelements() <= 100), false); }
711template<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
714template<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
726template<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
729template<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
737template<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
741template<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
754template<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
758template<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
768template<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
772template<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>
781template<typename T, typename AllocC, typename AllocR>
782void 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
789template<typename T, typename Alloc>
790void 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
796template<typename T, typename AllocC, typename AllocR>
797void 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
804template<typename T, typename Alloc>
805void 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
811template<typename T, typename AllocC, typename AllocR>
812Array<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}
821template<typename T, typename Alloc>
822Array<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}
830template<typename T, typename Alloc>
831Array<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
840template<typename T, typename AllocC, typename AllocR>
841Array<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}
850template<typename T, typename Alloc>
851Array<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}
859template<typename T, typename Alloc>
860Array<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>
872Array<std::complex<float>> conj(const Array<std::complex<float>> &carray);
873Array<std::complex<double>> conj(const Array<std::complex<double>> &carray);
874// Modifies rarray in place. rarray must be conformant.
875void conj(Array<std::complex<float>> &rarray, const Array<std::complex<float>> &carray);
876void 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.
879Matrix<std::complex<float>> conj(const Matrix<std::complex<float>> &carray);
880Matrix<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>
887template<typename T, typename Alloc>
889template<typename T, typename Alloc>
891template<typename T, typename Alloc>
893// </group>
894
895// Set the real part of the left complex array to the right real array.
896template<typename C, typename R, typename AllocC, typename AllocR>
897void 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.
900template<typename C, typename R, typename AllocC, typename AllocR>
901void 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>
905Array<float> real(const Array<std::complex<float>> &carray);
906Array<double> real(const Array<std::complex<double>> &carray);
907// Modifies rarray in place. rarray must be conformant.
908void real(Array<float> &rarray, const Array<std::complex<float>> &carray);
909void 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>
915Array<float> imag(const Array<std::complex<float>> &carray);
916Array<double> imag(const Array<std::complex<double>> &carray);
917// Modifies rarray in place. rarray must be conformant.
918void imag(Array<float> &rarray, const Array<std::complex<float>> &carray);
919void 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>
927Array<float> amplitude(const Array<std::complex<float>> &carray);
928Array<double> amplitude(const Array<std::complex<double>> &carray);
929// Modifies rarray in place. rarray must be conformant.
930void amplitude(Array<float> &rarray, const Array<std::complex<float>> &carray);
931void 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>
939Array<float> phase(const Array<std::complex<float>> &carray);
940Array<double> phase(const Array<std::complex<double>> &carray);
941// Modifies rarray in place. rarray must be conformant.
942void phase(Array<float> &rarray, const Array<std::complex<float>> &carray);
943void 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.
950Array<float> ComplexToReal(const Array<std::complex<float>> &carray);
951Array<double> ComplexToReal(const Array<std::complex<double>> &carray);
952// Modify the array "rarray" in place. "rarray" must be the correct shape.
953// <group>
954void ComplexToReal(Array<float> &rarray, const Array<std::complex<float>> &carray);
955void 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>
966void RealToComplex(Array<std::complex<float>> &carray, const Array<float> &rarray);
967void 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.
974template<typename T, typename U, typename AllocT, typename AllocU>
976
977// Returns an array where every element is squared.
978template<typename T, typename Alloc> Array<T, Alloc> square(const Array<T, Alloc> &val);
979
980// Returns an array where every element is cubed.
981template<typename T, typename Alloc> Array<T, Alloc> cube(const Array<T, Alloc> &val);
982
983// Helper function for expandArray using recursion for each axis.
984template<typename T>
985T* 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.
1042template<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 & steps() const
Return steps to be made if stepping one element in a dimension.
Definition: ArrayBase.h:136
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 & shape() const
The length of each axis.
Definition: ArrayBase.h:125
virtual void reference(const Array< T, Alloc > &other)
After invocation, this array and other reference the same storage.
void putStorage(T *&storage, bool deleteAndCopy)
putStorage() is normally called after a call to getStorage() (cf).
T * data()
Get a pointer to the beginning of the array.
Definition: Array.h:604
contiter cbegin()
Get the begin iterator object for a contiguous array.
Definition: Array.h:871
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...
T * getStorage(bool &deleteIt)
Generally use of this should be shunned, except to use a FORTRAN routine or something similar.
contiter cend()
Definition: Array.h:875
iterator begin()
Get the begin iterator object for any array.
Definition: Array.h:859
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)
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 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
void operator|=(Array< T, Alloc > &left, const Array< T, Alloc > &other)
Array< T, Alloc > floormod(const Array< T, Alloc > &a, const T &b)
Array< double > amplitude(const Array< std::complex< double > > &carray)
Array< double > imag(const Array< std::complex< double > > &carray)
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.
Array< T, Alloc > ceil(const Array< T, Alloc > &a)
Array< T, Alloc > asin(const Array< T, Alloc > &a)
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.
void ComplexToReal(Array< float > &rarray, const Array< std::complex< float > > &carray)
Modify the array "rarray" in place.
void imag(Array< float > &rarray, const Array< std::complex< float > > &carray)
Modifies rarray in place.
Array< T, Alloc > sign(const Array< T, Alloc > &a)
Array< T, Alloc > atan(const Array< T, Alloc > &a)
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
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.
Array< float > phase(const Array< std::complex< float > > &carray)
Extracts the phase (i.e.
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
Array< T, Alloc > operator/(const Array< T, Alloc > &left, const Array< T, Alloc > &right)
T product(const Array< T, Alloc > &a)
Product of every element of the array.
Array< float > real(const Array< std::complex< float > > &carray)
Extracts the real part of a complex array into an array of floats.
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 RealToComplex(Array< std::complex< double > > &carray, const Array< double > &rarray)
Array< T, Alloc > fmod(const T &a, const Array< T, Alloc > &b)
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
Array< T, Alloc > operator&(const Array< T, Alloc > &left, const Array< T, Alloc > &right)
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 operator%=(Array< T, Alloc > &left, const Array< T, Alloc > &other)
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.
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.
T min(const Array< T, Alloc > &a)
The minimum element of the array.
Array< double > ComplexToReal(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 > acos(const Array< T, Alloc > &a)
Transcendental function applied to the array on an element-by-element basis.
void conj(Array< std::complex< float > > &rarray, const Array< std::complex< float > > &carray)
Modifies rarray in place.
T pstddev(const Array< T, Alloc > &a, T mean, size_t ddof=0)
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< T, Alloc > floormod(const T &a, const Array< T, Alloc > &b)
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< T, Alloc > abs(const Array< T, Alloc > &a)
Array< double > phase(const Array< std::complex< double > > &carray)
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 > tan(const Array< T, Alloc > &a)
void operator*=(Array< T, Alloc > &left, const Array< T, Alloc > &other)
Definition: ArrayMath.h:345
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 > round(const Array< T, Alloc > &a)
Array< T, Alloc > operator~(const Array< T, Alloc > &a)
Array< T, Alloc > pow(const T &a, const Array< T, Alloc > &b)
void amplitude(Array< double > &rarray, const Array< std::complex< double > > &carray)
Array< T, Alloc > arrayTransformResult(T left, const Array< T, Alloc > &right, BinaryOperator op)
Transform left and right to a result using the binary operator.
Array< T, Alloc > fmod(const Array< T, Alloc > &a, const Array< T, Alloc > &b)
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< float > > conj(const Array< std::complex< float > > &carray)
Returns the complex conjugate of a complex array.
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< 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.
void ComplexToReal(Array< double > &rarray, const Array< std::complex< double > > &carray)
void real(Array< float > &rarray, const Array< std::complex< float > > &carray)
Modifies rarray in place.
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 > arrayTransformResult(const Array< T, Alloc > &arr, UnaryOperator op)
Transform array to a result using the unary operator.
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
Array< std::complex< T > > makeComplex(const T &real, const Array< T, Alloc > &imag)
Array< float > amplitude(const Array< std::complex< float > > &carray)
Extracts the amplitude (i.e.
Array< std::complex< T > > makeComplex(const Array< T, Alloc > &real, const T &imag)
Array< T, Alloc > tanh(const Array< T, Alloc > &a)
Array< T, Alloc > floormod(const Array< T, Alloc > &a, const Array< T, Alloc > &b)
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< float > imag(const Array< std::complex< float > > &carray)
Extracts the imaginary part of a complex array into an array of floats.
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
void phase(Array< double > &rarray, const Array< std::complex< double > > &carray)
Array< double > real(const Array< std::complex< double > > &carray)
Matrix< std::complex< float > > conj(const Matrix< std::complex< float > > &carray)
Array< std::complex< T > > operator*(const std::complex< T > &left, const Array< T, Alloc > &other)
Definition: ArrayMath.h:831
Matrix< std::complex< double > > conj(const Matrix< std::complex< double > > &carray)
Array< T, Alloc > fmod(const Array< T, Alloc > &a, const T &b)
Array< T, Alloc > atan2(const Array< T, Alloc > &y, const Array< T, Alloc > &x)
void operator+=(Array< T, Alloc > &left, const Array< T, Alloc > &other)
Element by element arithmetic modifying left in-place.
Array< T, Alloc > atan2(const Array< T, Alloc > &y, const T &x)
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
Array< std::complex< float > > RealToComplex(const Array< float > &rarray)
Copy an array of real,imaginary pairs into a complex array.
Array< T, Alloc > exp(const Array< T, Alloc > &a)
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< std::complex< T > > operator/(const std::complex< T > &left, const Array< T, Alloc > &other)
Definition: ArrayMath.h:860
Array< T, Alloc > operator-(const Array< T, Alloc > &a)
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
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.
Array< T, Alloc > operator|(const Array< T, Alloc > &left, const Array< T, Alloc > &right)
Array< T, Alloc > cosh(const Array< T, Alloc > &a)
T rms(const Array< T, Alloc > &a)
The root-mean-square of "a" is the sqrt of sum(a*a)/N.
Array< T, Alloc > pow(const Array< T, Alloc > &a, const Array< T, Alloc > &b)
Array< T, Alloc > min(const T &a, const Array< T, Alloc > &b)
Definition: ArrayMath.h:584
void imag(Array< double > &rarray, 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 > floor(const Array< T, Alloc > &a)
Array< T, Alloc > max(const T &a, const Array< T, Alloc > &b)
Array< T, Alloc > operator%(const Array< T, Alloc > &left, const Array< T, Alloc > &right)
void operator-=(Array< T, Alloc > &left, const Array< T, Alloc > &other)
void myrtransform(_InputIterator1 __first1, _InputIterator1 __last1, _OutputIterator __result, T right, _BinaryOperation __binary_op)
sequence = sequence OP scalar
Definition: ArrayMath.h:157
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< T, Alloc > square(const Array< T, Alloc > &val)
Returns an array where every element is squared.
void phase(Array< float > &rarray, const Array< std::complex< float > > &carray)
Modifies rarray in place.
Array< T, Alloc > operator^(const Array< T, Alloc > &left, const Array< T, Alloc > &right)
Array< std::complex< double > > RealToComplex(const Array< double > &rarray)
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...
Array< T, Alloc > cos(const Array< T, Alloc > &a)
Transcendental function that can be applied to essentially all numeric types.
T median(const Array< T, Alloc > &a, bool sorted)
Definition: ArrayMath.h:690
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.
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".
void RealToComplex(Array< std::complex< float > > &carray, const Array< float > &rarray)
Modify the array "carray" in place.
Array< T, Alloc > sin(const Array< T, Alloc > &a)
void operator^=(Array< T, Alloc > &left, const Array< T, Alloc > &other)
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.
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 > arrayTransformResult(const Array< T, Alloc > &left, T right, BinaryOperator op)
Transform left and right to a result using the binary operator.
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.
Array< T, Alloc > log10(const Array< T, Alloc > &a)
Array< T, Alloc > operator+(const Array< T, Alloc > &a)
Unary arithmetic operation.
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
Array< T, Alloc > fabs(const Array< T, Alloc > &a)
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
Array< T, Alloc > operator*(const Array< T, Alloc > &left, const Array< T, Alloc > &right)
Definition: ArrayMath.h:398
T sum(const Array< T, Alloc > &a)
Sum of every element of the array.
Array< T, Alloc > atan2(const T &y, const Array< T, Alloc > &x)
Array< T, Alloc > log(const Array< T, Alloc > &a)
void amplitude(Array< float > &rarray, const Array< std::complex< float > > &carray)
Modifies rarray in place.
void real(Array< double > &rarray, const Array< std::complex< double > > &carray)
Array< std::complex< double > > conj(const Array< std::complex< double > > &carray)
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
T variance(const Array< T, Alloc > &a)
The variance of "a" is the sum of (a(i) - mean(a))**2/(a.nelements() - ddof).
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
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)
Array< T, Alloc > sinh(const Array< T, Alloc > &a)
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.
void operator&=(Array< T, Alloc > &left, const Array< T, Alloc > &other)
void conj(Array< std::complex< double > > &rarray, const Array< std::complex< double > > &carray)
Array< T, Alloc > sqrt(const Array< T, Alloc > &a)
void operator/=(Array< T, Alloc > &left, const Array< T, Alloc > &other)
Definition: ArrayMath.h:351
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.
void arrayTransform(const Array< T, Alloc > &arr, Array< RES, AllocRES > &result, UnaryOperator op)
Transform array to a result using the unary operator.