//# MaskArrMath.h: Simple mathematics done with MaskedArray's.
//# Copyright (C) 1993,1994,1995,1996,1999,2001
//# Associated Universities, Inc. Washington DC, USA.
//#
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//#
//# This library is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//# Internet email: aips2-request@nrao.edu.
//# Postal address: AIPS++ Project Office
//# National Radio Astronomy Observatory
//# 520 Edgemont Road
//# Charlottesville, VA 22903-2475 USA
//#
//# $Id$
#ifndef CASA_MASKARRMATH_2_H
#define CASA_MASKARRMATH_2_H
#include "Array.h"
#include "MaskedArray.h"
#include "IPosition.h"
namespace casacore { //# NAMESPACE CASACORE - BEGIN
// Mathematical operations for MaskedArrays (and with Arrays)
//
//
//
// Array
// MaskedArray
//
//
//
// MaskArrMath is short for MaskedArrayMath, which is too long by the old
// AIPS++ file naming conventions. This file contains global functions
// which perform element by element mathematical operations on masked arrays.
//
//
//
// These functions perform element by element mathematical operations on
// masked arrays. With two arrays, they must both conform, and the result
// is done element by element, for those locations where the mask of the
// MaskedArray is true. For two MaskedArrays, the "and" of the masks is used.
//
//
//
//
// Vector a(10);
// Vector b(10);
// Vector c(10);
// . . .
// c = a(a>0) + b(b>0);
//
// This example sets those elements of c where ((a>0) && (b>0)) to (a+b).
// Elements of c where !((a>0) && (b>0)) are unchanged. The result of
// this operation is a MaskedArray. The assignment from this
// MaskedArray to the Vector c only assigns those elements
// where the mask is true.
//
//
//
//
// Vector a(10);
// Vector b(10);
// Vector c(10);
// . . .
// c = atan2 (a, b(b>0);
//
// This example sets those elements of c where (b>0) to atan2 (a,b).
// Elements of c where !(b>0) are unchanged. The result of
// this operation is a MaskedArray. The assignment from this
// MaskedArray to the Vector c only assigns those elements
// where the mask is true.
//
//
//
//
// Vector a(10);
// int result;
// . . .
// result = sum (a(a>0));
//
// This example sums a, for those elements of a which are greater than 0.
//
//
//
// One wants to be able to mask arrays and perform mathematical operations on
// those masked arrays. Since the masked arrays are only defined where
// the masks are true, the result must be a MaskedArray, or a simple number.
//
//
//
// MaskedArray mathematical operations -- Mathematical
// operations for MaskedArrays, and between MaskedArrays and Arrays.
//
//
//
// Element by element arithmetic modifying left in-place. left and other
// must be conformant.
//
//
// ArrayConformanceError
//
//
//
template const MaskedArray & operator+= (const MaskedArray &left, const Array &other);
template const MaskedArray & operator-= (const MaskedArray &left, const Array &other);
template const MaskedArray & operator*= (const MaskedArray &left, const Array &other);
template const MaskedArray & operator/= (const MaskedArray &left, const Array &other);
template Array & operator+= (Array &left, const MaskedArray &other);
template Array & operator-= (Array &left, const MaskedArray &other);
template Array & operator*= (Array &left, const MaskedArray &other);
template Array & operator/= (Array &left, const MaskedArray &other);
template const MaskedArray & operator+= (const MaskedArray &left, const MaskedArray &other);
template const MaskedArray & operator-= (const MaskedArray &left, const MaskedArray &other);
template const MaskedArray & operator*= (const MaskedArray &left, const MaskedArray &other);
template const MaskedArray & operator/= (const MaskedArray &left,const MaskedArray &other);
template const MaskedArray & operator/= (const MaskedArray &left,const MaskedArray &other);
//
//
// Element by element arithmetic modifying left in-place. The scalar "other"
// behaves as if it were a conformant Array to left filled with constant values.
//
template const MaskedArray & operator+= (const MaskedArray &left,const T &other);
template const MaskedArray & operator-= (const MaskedArray &left,const T &other);
template const MaskedArray & operator*= (const MaskedArray &left,const T &other);
template const MaskedArray & operator/= (const MaskedArray &left,const T &other);
//
// Unary arithmetic operation.
//
//
template MaskedArray operator+(const MaskedArray &a);
template MaskedArray operator-(const MaskedArray &a);
//
//
// Element by element arithmetic on MaskedArrays, returns a MaskedArray.
//
//
// ArrayConformanceError
//
//
//
template MaskedArray operator+ (const MaskedArray &left, const Array &right);
template MaskedArray operator- (const MaskedArray &left, const Array &right);
template MaskedArray operator* (const MaskedArray &left, const Array &right);
template MaskedArray operator/ (const MaskedArray &left, const Array &right);
template MaskedArray operator+ (const Array &left, const MaskedArray &right);
template MaskedArray operator- (const Array &left, const MaskedArray &right);
template MaskedArray operator* (const Array &left, const MaskedArray &right);
template MaskedArray operator/ (const Array &left, const MaskedArray &right);
template MaskedArray operator+ (const MaskedArray &left,const MaskedArray &right);
template MaskedArray operator- (const MaskedArray &left,const MaskedArray &right);
template MaskedArray operator* (const MaskedArray &left,const MaskedArray &right);
template MaskedArray operator/ (const MaskedArray &left,const MaskedArray &right);
//
//
// Element by element arithmetic between a MaskedArray and a scalar, returning
// a MaskedArray.
//
template MaskedArray operator+ (const MaskedArray &left, const T &right);
template MaskedArray operator- (const MaskedArray &left, const T &right);
template MaskedArray operator* (const MaskedArray &left, const T &right);
template MaskedArray operator/ (const MaskedArray &left, const T &right);
MaskedArray> operator* (const MaskedArray> &left, const float &right);
//
//
// Element by element arithmetic between a scalar and a MaskedArray, returning
// a MaskedArray.
//
template MaskedArray operator+ (const T &left, const MaskedArray &right);
template MaskedArray operator- (const T &left, const MaskedArray &right);
template MaskedArray operator* (const T &left, const MaskedArray &right);
template MaskedArray operator/ (const T &left, const MaskedArray &right);
MaskedArray> operator* (const float &left, const MaskedArray> &right);
//
//
// Transcendental function applied to the array on an element-by-element
// basis. Although a template function, this may not make sense for all
// numeric types.
//
template MaskedArray sin(const MaskedArray &left);
template MaskedArray cos(const MaskedArray &left);
template MaskedArray tan(const MaskedArray &left);
template MaskedArray asin(const MaskedArray &left);
template MaskedArray acos(const MaskedArray &left);
template MaskedArray atan(const MaskedArray &left);
template MaskedArray sinh(const MaskedArray &left);
template MaskedArray cosh(const MaskedArray &left);
template MaskedArray tanh(const MaskedArray &left);
template MaskedArray exp(const MaskedArray &left);
template MaskedArray log(const MaskedArray &left);
template MaskedArray log10(const MaskedArray &left);
template MaskedArray sqrt(const MaskedArray &left);
template MaskedArray abs(const MaskedArray &left);
template MaskedArray fabs(const MaskedArray &left);
template MaskedArray ceil(const MaskedArray &left);
template MaskedArray floor(const MaskedArray &left);
//
// Transcendental functions requiring two arguments applied on an element-by-element
// basis. Although a template function, this may not make sense for all
// numeric types.
//
// ArrayConformanceError
//
//
//
template MaskedArray atan2(const MaskedArray &left, const Array &right);
template MaskedArray fmod(const MaskedArray &left, const Array &right);
template MaskedArray atan2(const Array &left, const MaskedArray &right);
template MaskedArray fmod(const Array &left, const MaskedArray &right);
template MaskedArray atan2(const MaskedArray &left,const MaskedArray &right);
template MaskedArray fmod(const MaskedArray &left,const MaskedArray &right);
template MaskedArray atan2(const MaskedArray &left, const T &right);
template MaskedArray fmod(const MaskedArray &left, const T &right);
template MaskedArray atan2(const T &left, const MaskedArray &right);
template MaskedArray fmod(const T &left, const MaskedArray &right);
template MaskedArray pow(const MaskedArray &left, const Array &right);
template MaskedArray pow(const Array &left, const MaskedArray &right);
template MaskedArray pow(const MaskedArray &left,const MaskedArray &right);
template MaskedArray pow(const MaskedArray &left, const double &right);
//
// Extracts the real part of a complex array into an array of floats.
template
MaskedArray real(const MaskedArray> &carray)
{ return MaskedArray (real(carray.getArray()), carray.getMask()); }
//
// Extracts the imaginary part of a complex array into an array of floats.
template
MaskedArray imag(const MaskedArray> &carray)
{ return MaskedArray (imag(carray.getArray()), carray.getMask()); }
//
// Find the minimum and maximum values of a MaskedArray.
// Also find the IPositions of the minimum and maximum values.
//
//
// ArrayError
//
//
//
template void minMax(T &minVal, T &maxVal, IPosition &minPos, IPosition &maxPos,const MaskedArray &marray);
template void minMax(T &minVal, T &maxVal,const MaskedArray &marray);
//
//
// The "min" and "max" functions require that the type "T" have comparison
// operators.
// The minimum element of the array.
template T min(const MaskedArray &left);
// Return an array that contains the minimum of "left" and "right" at each
// position.
//
// "left" and "right" must be conformant.
//
//
// ArrayError
//
//
template MaskedArray min(const MaskedArray &left, const Array &right);
template MaskedArray min(const Array &left, const MaskedArray &right);
template MaskedArray min(const MaskedArray &left, const MaskedArray &right);
template MaskedArray min(const T &left, const MaskedArray &right);
template MaskedArray min(const MaskedArray &left, const T &right);
//
// "result" contains the minimum of "left" and "right" at each position.
// "result", "left", and "right" must be conformant.
//
//
// ArrayConformanceError
//
//
template void min(const MaskedArray &result, const Array &left, const Array &right);
// The maximum element of the array.
template T max(const MaskedArray &left);
// Return an array that contains the maximum of "left" and "right" at each
// position.
//
// "left" and "right" must be conformant.
//
// ArrayError
//
//
//
template MaskedArray max(const MaskedArray &left, const Array &right);
template MaskedArray max(const Array &left, const MaskedArray &right);
template MaskedArray max(const MaskedArray &left, const MaskedArray &right);
template MaskedArray max(const T &left, const MaskedArray &right);
template MaskedArray max(const MaskedArray &left, const T &right);
//
// "result" contains the maximum of "left" and "right" at each position.
// "result", "left", and "right" must be conformant.
//
//
// ArrayConformanceError
//
//
template void max(const MaskedArray &result,const Array &left, const Array &right);
//
// Fills all elements of "array" where the mask is true with a sequence
// starting with "start" and incrementing by "inc" for each element
// where the mask is true.
// The first axis varies most rapidly.
template void indgen(MaskedArray &a, T start, T inc);
//
// Fills all elements of "array" where the mask is true with a sequence
// starting with 0 and incremented by one for each element
// where the mask is true.
// The first axis varies most rapidly.
template void indgen(MaskedArray &a);
//
// Fills all elements of "array" where the mask is true with a sequence
// starting with "start" and incremented by one for each element
// where the mask is true.
// The first axis varies most rapidly.
template void indgen(MaskedArray &a, T start);
//
// ArrayError
//
//
// Sum of every element of the MaskedArray where the Mask is true.
template T sum(const MaskedArray &a);
//
// Sum of the squares of every element of the MaskedArray where the Mask is true.
template T sumsquares(const MaskedArray &a);
//
// Product of every element of the MaskedArray where the Mask is true.
// This could of course easily overflow.
template T product(const MaskedArray &a);
//
// The mean of "a" is the sum of all elements of "a" divided by the number
// of elements of "a".
template T mean(const MaskedArray &a);
//
// The variance of "a" is the sum of (a(i) - mean(a))**2/(a.nelements() - ddof).
// Similar to numpy the argument ddof tells if the population variance (ddof=0)
// or the sample variance (ddof=1) is taken.
// The variance functions proper use ddof=1.
//
Note that for a complex valued T the absolute values are used; in that way
// the variance is equal to the sum of the variances of the real and imaginary parts.
// Hence the imaginary part in the return value is 0.
template T variance(const MaskedArray &a);
template T pvariance(const MaskedArray &a, size_t ddof=0);
// Rather than using a computed mean, use the supplied value.
template T variance(const MaskedArray &a, T mean);
template T pvariance(const MaskedArray &a, T mean, size_t ddof=0);
// The standard deviation of "a" is the square root of its variance.
template T stddev(const MaskedArray &a);
template T pstddev(const MaskedArray &a, size_t ddof=0);
template T stddev(const MaskedArray &a, T mean);
template T pstddev(const MaskedArray &a, T mean, size_t ddof=0);
//
// The average deviation of "a" is the sum of abs(a(i) - mean(a))/N. (N.B.
// N, not N-1 in the denominator).
template T avdev(const MaskedArray &a);
//
// The average deviation of "a" is the sum of abs(a(i) - mean(a))/N. (N.B.
// N, not N-1 in the denominator).
// Rather than using a computed mean, use the supplied value.
template T avdev(const MaskedArray &a,T mean);
//
// The root-mean-square of "a" is the sqrt of sum(a*a)/N.
template T rms(const MaskedArray &a);
//
// The median of "a" is a(n/2).
// When a has an even number of elements and the switch takeEvenMean is set,
// the median is 0.5*(a(n/2) + a((n+1)/2)).
// According to Numerical Recipes (2nd edition) it makes little sense to take
// the mean when the array is large enough (> 100 elements). Therefore
// the default for takeEvenMean is false when the array has > 100 elements,
// otherwise it is true.
//
If "sorted"==true we assume the data is already sorted and we
// compute the median directly. Otherwise the function GenSort::kthLargest
// is used to find the median (kthLargest is about 6 times faster
// than a full quicksort).
//
template inline T median(const MaskedArray &a, bool sorted=false)
{ return median (a, sorted, (a.nelements() <= 100)); }
template T median(const MaskedArray &a, bool sorted,
bool takeEvenMean);
//
// The median absolute deviation from the median. Interface is as for
// the median functions
//
template inline T madfm(const MaskedArray &a, bool sorted=false)
{ return madfm (a, sorted, (a.nelements() <= 100)); }
template T madfm(const MaskedArray &a, bool sorted,
bool takeEvenMean);
//
// Returns a MaskedArray where every element is squared.
template MaskedArray square(const MaskedArray &val);
// Returns a MaskedArray where every element is cubed.
template MaskedArray cube(const MaskedArray &val);
//
template class MaskedSumFunc {
public:
T operator() (const MaskedArray& arr) const { return sum(arr); }
};
template class MaskedProductFunc {
public:
T operator() (const MaskedArray& arr) const { return product(arr); }
};
template class MaskedMinFunc {
public:
T operator() (const MaskedArray& arr) const { return min(arr); }
};
template class MaskedMaxFunc {
public:
T operator() (const MaskedArray& arr) const { return max(arr); }
};
template class MaskedMeanFunc {
public:
T operator() (const MaskedArray& arr) const { return mean(arr); }
};
template class MaskedVarianceFunc {
public:
T operator() (const MaskedArray& arr) const { return variance(arr); }
};
template class MaskedStddevFunc {
public:
T operator() (const MaskedArray& arr) const { return stddev(arr); }
};
template class MaskedAvdevFunc {
public:
T operator() (const MaskedArray& arr) const { return avdev(arr); }
};
template class MaskedRmsFunc {
public:
T operator() (const MaskedArray& arr) const { return rms(arr); }
};
template class MaskedMedianFunc {
public:
explicit MaskedMedianFunc (bool sorted=false, bool takeEvenMean=true)
: itsSorted(sorted), itsTakeEvenMean(takeEvenMean) {}
T operator() (const MaskedArray& arr) const
{ return median(arr, itsSorted, itsTakeEvenMean); }
private:
bool itsSorted;
bool itsTakeEvenMean;
bool itsInPlace;
};
template class MaskedMadfmFunc {
public:
explicit MaskedMadfmFunc(bool sorted=false, bool takeEvenMean=true)
: itsSorted(sorted), itsTakeEvenMean(takeEvenMean) {}
float operator()(const MaskedArray& arr) const
{ return madfm(arr, itsSorted, itsTakeEvenMean); }
private:
bool itsSorted;
bool itsTakeEvenMean;
bool itsInPlace;
};
// Apply the given ArrayMath reduction function objects
// to each box in the array.
//
// Downsample an array by taking the mean of every [25,25] elements.
//
// Array downArr = boxedArrayMath(in, IPosition(2,25,25),
// MaskedMeanFunc());
//
//
// The dimensionality of the array can be larger than the box; in that
// case the missing axes of the box are assumed to have length 1.
// A box axis length <= 0 means the full array axis.
template
MaskedArray boxedArrayMath (const MaskedArray& array,
const IPosition& boxSize,
const FuncType& funcObj);
// Apply for each element in the array the given ArrayMath reduction function
// object to the box around that element. The full box is 2*halfBoxSize + 1.
// It can be used for arrays and boxes of any dimensionality; missing
// halfBoxSize values are set to 1.
//
// Determine for each element in the array the median of a box
// with size [51,51] around that element:
//
// Array medians = slidingArrayMath(in, IPosition(2,25,25),
// MaskedMedianFunc());
//
// This is a potentially expensive operation. On a high-end PC it took
// appr. 27 seconds to get the medians for an array of [1000,1000] using
// a halfBoxSize of [50,50].
//
//
The fillEdge argument determines how the edge is filled where
// no full boxes can be made. true means it is set to zero; false means
// that the edge is removed, thus the output array is smaller than the
// input array.
// This brute-force method of determining the medians outperforms
// all kinds of smart implementations. For a vector it is about as fast
// as the casacore class MedianSlider, for a 2D array
// it is much, much faster.
//
template
Array slidingArrayMath (const MaskedArray& array,
const IPosition& halfBoxSize,
const FuncType& funcObj,
bool fillEdge=true);
} //# NAMESPACE CASACORE - END
#include "MaskArrMath.tcc"
#endif