# pylint: disable-msg=E1002
"""MA: a facility for dealing with missing observations
MA is generally used as a numpy.array look-alike.
by Paul F. Dubois.
Copyright 1999, 2000, 2001 Regents of the University of California.
Released for unlimited redistribution.
Adapted for numpy_core 2005 by Travis Oliphant and
(mainly) Paul Dubois.
Subclassing of the base ndarray 2006 by Pierre Gerard-Marchant.
pgmdevlist_AT_gmail_DOT_com
Improvements suggested by Reggie Dugard (reggie_AT_merfinllc_DOT_com)
:author: Pierre Gerard-Marchant
:contact: pierregm_at_uga_dot_edu
"""
__author__ = "Pierre GF Gerard-Marchant"
__docformat__ = "restructuredtext en"
__all__ = ['MAError', 'MaskType', 'MaskedArray',
'bool_', 'complex_', 'float_', 'int_', 'object_',
'abs', 'absolute', 'add', 'all', 'allclose', 'allequal', 'alltrue',
'amax', 'amin', 'anom', 'anomalies', 'any', 'arange',
'arccos', 'arccosh', 'arcsin', 'arcsinh', 'arctan', 'arctan2',
'arctanh', 'argmax', 'argmin', 'argsort', 'around',
'array', 'asarray','asanyarray',
'bitwise_and', 'bitwise_or', 'bitwise_xor',
'ceil', 'choose', 'clip', 'common_fill_value', 'compress',
'compressed', 'concatenate', 'conjugate', 'cos', 'cosh', 'count',
'default_fill_value', 'diagonal', 'divide', 'dump', 'dumps',
'empty', 'empty_like', 'equal', 'exp',
'fabs', 'fmod', 'filled', 'floor', 'floor_divide','fix_invalid',
'frombuffer', 'fromfunction',
'getdata','getmask', 'getmaskarray', 'greater', 'greater_equal',
'hypot',
'identity', 'ids', 'indices', 'inner', 'innerproduct',
'isMA', 'isMaskedArray', 'is_mask', 'is_masked', 'isarray',
'left_shift', 'less', 'less_equal', 'load', 'loads', 'log', 'log10',
'logical_and', 'logical_not', 'logical_or', 'logical_xor',
'make_mask', 'make_mask_none', 'mask_or', 'masked',
'masked_array', 'masked_equal', 'masked_greater',
'masked_greater_equal', 'masked_inside', 'masked_invalid',
'masked_less','masked_less_equal', 'masked_not_equal',
'masked_object','masked_outside', 'masked_print_option',
'masked_singleton','masked_values', 'masked_where', 'max', 'maximum',
'mean', 'min', 'minimum', 'multiply',
'negative', 'nomask', 'nonzero', 'not_equal',
'ones', 'outer', 'outerproduct',
'power', 'product', 'ptp', 'put', 'putmask',
'rank', 'ravel', 'remainder', 'repeat', 'reshape', 'resize',
'right_shift', 'round_',
'set_fill_value', 'shape', 'sin', 'sinh', 'size', 'sometrue', 'sort',
'sqrt', 'std', 'subtract', 'sum', 'swapaxes',
'take', 'tan', 'tanh', 'transpose', 'true_divide',
'var', 'where',
'zeros']
import sys
import types
import cPickle
import operator
import numpy
from numpy.core import bool_,complex_,float_,int_,object_,str_
import numpy.core.umath as umath
import numpy.core.fromnumeric as fromnumeric
import numpy.core.numeric as numeric
import numpy.core.numerictypes as ntypes
from numpy import bool_,dtype,typecodes,amax,amin,ndarray,iscomplexobj
from numpy import expand_dims
from numpy import array
import warnings
MaskType = bool_
nomask = MaskType(0)
divide_tolerance = 1.e-35
numpy.seterr(all='ignore')
def doc_note(note):
return "\nNotes\n-----\n%s" % note
#####--------------------------------------------------------------------------
#---- --- Exceptions ---
#####--------------------------------------------------------------------------
class MAError(Exception):
"Class for MA related errors."
def __init__ (self, args=None):
"Creates an exception."
Exception.__init__(self,args)
self.args = args
def __str__(self):
"Calculates the string representation."
return str(self.args)
__repr__ = __str__
#####--------------------------------------------------------------------------
#---- --- Filling options ---
#####--------------------------------------------------------------------------
# b: boolean - c: complex - f: floats - i: integer - O: object - S: string
default_filler = {'b': True,
'c' : 1.e20 + 0.0j,
'f' : 1.e20,
'i' : 999999,
'O' : '?',
'S' : 'N/A',
'u' : 999999,
'V' : '???',
}
max_filler = ntypes._minvals
max_filler.update([(k,-numpy.inf) for k in [numpy.float32, numpy.float64]])
min_filler = ntypes._maxvals
min_filler.update([(k,numpy.inf) for k in [numpy.float32, numpy.float64]])
if 'float128' in ntypes.typeDict:
max_filler.update([(numpy.float128,-numpy.inf)])
min_filler.update([(numpy.float128, numpy.inf)])
def default_fill_value(obj):
"""Calculate the default fill value for the argument object.
"""
if hasattr(obj,'dtype'):
defval = default_filler[obj.dtype.kind]
elif isinstance(obj, numeric.dtype):
defval = default_filler[obj.kind]
elif isinstance(obj, float):
defval = default_filler['f']
elif isinstance(obj, int) or isinstance(obj, long):
defval = default_filler['i']
elif isinstance(obj, str):
defval = default_filler['S']
elif isinstance(obj, complex):
defval = default_filler['c']
else:
defval = default_filler['O']
return defval
def minimum_fill_value(obj):
"""Calculate the default fill value suitable for taking the
minimum of ``obj``.
"""
if hasattr(obj, 'dtype'):
objtype = obj.dtype
filler = min_filler[objtype]
if filler is None:
raise TypeError, 'Unsuitable type for calculating minimum.'
return filler
elif isinstance(obj, float):
return min_filler[ntypes.typeDict['float_']]
elif isinstance(obj, int):
return min_filler[ntypes.typeDict['int_']]
elif isinstance(obj, long):
return min_filler[ntypes.typeDict['uint']]
elif isinstance(obj, numeric.dtype):
return min_filler[obj]
else:
raise TypeError, 'Unsuitable type for calculating minimum.'
def maximum_fill_value(obj):
"""Calculate the default fill value suitable for taking the maximum
of ``obj``.
"""
if hasattr(obj, 'dtype'):
objtype = obj.dtype
filler = max_filler[objtype]
if filler is None:
raise TypeError, 'Unsuitable type for calculating minimum.'
return filler
elif isinstance(obj, float):
return max_filler[ntypes.typeDict['float_']]
elif isinstance(obj, int):
return max_filler[ntypes.typeDict['int_']]
elif isinstance(obj, long):
return max_filler[ntypes.typeDict['uint']]
elif isinstance(obj, numeric.dtype):
return max_filler[obj]
else:
raise TypeError, 'Unsuitable type for calculating minimum.'
def _check_fill_value(fill_value, dtype):
descr = numpy.dtype(dtype).descr
if fill_value is None:
if len(descr) > 1:
fill_value = [default_fill_value(numeric.dtype(d[1]))
for d in descr]
else:
fill_value = default_fill_value(dtype)
else:
fill_value = narray(fill_value).tolist()
fval = numpy.resize(fill_value, len(descr))
if len(descr) > 1:
fill_value = [numpy.asarray(f).astype(d[1]).item()
for (f,d) in zip(fval, descr)]
else:
fill_value = narray(fval, copy=False, dtype=dtype).item()
return fill_value
def set_fill_value(a, fill_value):
"""Set the filling value of a, if a is a masked array. Otherwise,
do nothing.
Returns
-------
None
"""
if isinstance(a, MaskedArray):
a._fill_value = _check_fill_value(fill_value, a.dtype)
return
def get_fill_value(a):
"""Return the filling value of a, if any. Otherwise, returns the
default filling value for that type.
"""
if isinstance(a, MaskedArray):
result = a.fill_value
else:
result = default_fill_value(a)
return result
def common_fill_value(a, b):
"""Return the common filling value of a and b, if any.
If a and b have different filling values, returns None.
"""
t1 = get_fill_value(a)
t2 = get_fill_value(b)
if t1 == t2:
return t1
return None
#####--------------------------------------------------------------------------
def filled(a, value = None):
"""Return a as an array with masked data replaced by value. If
value is None, get_fill_value(a) is used instead. If a is already
a ndarray, a itself is returned.
Parameters
----------
a : maskedarray or array_like
An input object.
value : {var}, optional
Filling value. If not given, the output of get_fill_value(a)
is used instead.
Returns
-------
a : array_like
"""
if hasattr(a, 'filled'):
return a.filled(value)
elif isinstance(a, ndarray):
# Should we check for contiguity ? and a.flags['CONTIGUOUS']:
return a
elif isinstance(a, dict):
return narray(a, 'O')
else:
return narray(a)
#####--------------------------------------------------------------------------
def get_masked_subclass(*arrays):
"""Return the youngest subclass of MaskedArray from a list of
(masked) arrays. In case of siblings, the first takes over.
"""
if len(arrays) == 1:
arr = arrays[0]
if isinstance(arr, MaskedArray):
rcls = type(arr)
else:
rcls = MaskedArray
else:
arrcls = [type(a) for a in arrays]
rcls = arrcls[0]
if not issubclass(rcls, MaskedArray):
rcls = MaskedArray
for cls in arrcls[1:]:
if issubclass(cls, rcls):
rcls = cls
return rcls
#####--------------------------------------------------------------------------
def get_data(a, subok=True):
"""Return the _data part of a (if any), or a as a ndarray.
Parameters
----------
a : array_like
A ndarray or a subclass of.
subok : bool
Whether to force the output to a 'pure' ndarray (False) or to
return a subclass of ndarray if approriate (True).
"""
data = getattr(a, '_data', numpy.array(a, subok=subok))
if not subok:
return data.view(ndarray)
return data
getdata = get_data
def fix_invalid(a, copy=True, fill_value=None):
"""Return (a copy of) a where invalid data (nan/inf) are masked
and replaced by fill_value.
Note that a copy is performed by default (just in case...).
Parameters
----------
a : array_like
A (subclass of) ndarray.
copy : bool
Whether to use a copy of a (True) or to fix a in place (False).
fill_value : {var}, optional
Value used for fixing invalid data. If not given, the output
of get_fill_value(a) is used instead.
Returns
-------
b : MaskedArray
"""
a = masked_array(a, copy=copy, subok=True)
#invalid = (numpy.isnan(a._data) | numpy.isinf(a._data))
invalid = numpy.logical_not(numpy.isfinite(a._data))
if not invalid.any():
return a
a._mask |= invalid
if fill_value is None:
fill_value = a.fill_value
a._data[invalid] = fill_value
return a
#####--------------------------------------------------------------------------
#---- --- Ufuncs ---
#####--------------------------------------------------------------------------
ufunc_domain = {}
ufunc_fills = {}
class _DomainCheckInterval:
"""Define a valid interval, so that :
``domain_check_interval(a,b)(x) = true`` where
``x < a`` or ``x > b``.
"""
def __init__(self, a, b):
"domain_check_interval(a,b)(x) = true where x < a or y > b"
if (a > b):
(a, b) = (b, a)
self.a = a
self.b = b
def __call__ (self, x):
"Execute the call behavior."
return umath.logical_or(umath.greater (x, self.b),
umath.less(x, self.a))
#............................
class _DomainTan:
"""Define a valid interval for the `tan` function, so that:
``domain_tan(eps) = True`` where ``abs(cos(x)) < eps``
"""
def __init__(self, eps):
"domain_tan(eps) = true where abs(cos(x)) < eps)"
self.eps = eps
def __call__ (self, x):
"Executes the call behavior."
return umath.less(umath.absolute(umath.cos(x)), self.eps)
#............................
class _DomainSafeDivide:
"""Define a domain for safe division."""
def __init__ (self, tolerance=divide_tolerance):
self.tolerance = tolerance
def __call__ (self, a, b):
return umath.absolute(a) * self.tolerance >= umath.absolute(b)
#............................
class _DomainGreater:
"DomainGreater(v)(x) = true where x <= v"
def __init__(self, critical_value):
"DomainGreater(v)(x) = true where x <= v"
self.critical_value = critical_value
def __call__ (self, x):
"Executes the call behavior."
return umath.less_equal(x, self.critical_value)
#............................
class _DomainGreaterEqual:
"DomainGreaterEqual(v)(x) = true where x < v"
def __init__(self, critical_value):
"DomainGreaterEqual(v)(x) = true where x < v"
self.critical_value = critical_value
def __call__ (self, x):
"Executes the call behavior."
return umath.less(x, self.critical_value)
#..............................................................................
class _MaskedUnaryOperation:
"""Defines masked version of unary operations, where invalid
values are pre-masked.
Parameters
----------
f : callable
fill :
Default filling value (0).
domain :
Default domain (None).
"""
def __init__ (self, mufunc, fill=0, domain=None):
""" _MaskedUnaryOperation(aufunc, fill=0, domain=None)
aufunc(fill) must be defined
self(x) returns aufunc(x)
with masked values where domain(x) is true or getmask(x) is true.
"""
self.f = mufunc
self.fill = fill
self.domain = domain
self.__doc__ = getattr(mufunc, "__doc__", str(mufunc))
self.__name__ = getattr(mufunc, "__name__", str(mufunc))
ufunc_domain[mufunc] = domain
ufunc_fills[mufunc] = fill
#
def __call__ (self, a, *args, **kwargs):
"Execute the call behavior."
#
m = getmask(a)
d1 = get_data(a)
#
if self.domain is not None:
dm = narray(self.domain(d1), copy=False)
m = numpy.logical_or(m, dm)
# The following two lines control the domain filling methods.
d1 = d1.copy()
# We could use smart indexing : d1[dm] = self.fill ...
# ... but numpy.putmask looks more efficient, despite the copy.
numpy.putmask(d1, dm, self.fill)
# Take care of the masked singletong first ...
if not m.ndim and m:
return masked
# Get the result class .......................
if isinstance(a, MaskedArray):
subtype = type(a)
else:
subtype = MaskedArray
# Get the result as a view of the subtype ...
result = self.f(d1, *args, **kwargs).view(subtype)
# Fix the mask if we don't have a scalar
if result.ndim > 0:
result._mask = m
result._update_from(a)
return result
#
def __str__ (self):
return "Masked version of %s. [Invalid values are masked]" % str(self.f)
#..............................................................................
class _MaskedBinaryOperation:
"""Define masked version of binary operations, where invalid
values are pre-masked.
Parameters
----------
f : callable
fillx :
Default filling value for the first argument (0).
filly :
Default filling value for the second argument (0).
domain :
Default domain (None).
"""
def __init__ (self, mbfunc, fillx=0, filly=0):
"""abfunc(fillx, filly) must be defined.
abfunc(x, filly) = x for all x to enable reduce.
"""
self.f = mbfunc
self.fillx = fillx
self.filly = filly
self.__doc__ = getattr(mbfunc, "__doc__", str(mbfunc))
self.__name__ = getattr(mbfunc, "__name__", str(mbfunc))
ufunc_domain[mbfunc] = None
ufunc_fills[mbfunc] = (fillx, filly)
#
def __call__ (self, a, b, *args, **kwargs):
"Execute the call behavior."
m = mask_or(getmask(a), getmask(b))
(d1, d2) = (get_data(a), get_data(b))
result = self.f(d1, d2, *args, **kwargs).view(get_masked_subclass(a,b))
if result.size > 1:
if m is not nomask:
result._mask = make_mask_none(result.shape)
result._mask.flat = m
if isinstance(a,MaskedArray):
result._update_from(a)
if isinstance(b,MaskedArray):
result._update_from(b)
elif m:
return masked
return result
#
def reduce (self, target, axis=0, dtype=None):
"""Reduce `target` along the given `axis`."""
if isinstance(target, MaskedArray):
tclass = type(target)
else:
tclass = MaskedArray
m = getmask(target)
t = filled(target, self.filly)
if t.shape == ():
t = t.reshape(1)
if m is not nomask:
m = make_mask(m, copy=1)
m.shape = (1,)
if m is nomask:
return self.f.reduce(t, axis).view(tclass)
t = t.view(tclass)
t._mask = m
tr = self.f.reduce(getdata(t), axis, dtype=dtype or t.dtype)
mr = umath.logical_and.reduce(m, axis)
tr = tr.view(tclass)
if mr.ndim > 0:
tr._mask = mr
return tr
elif mr:
return masked
return tr
def outer (self, a, b):
"""Return the function applied to the outer product of a and b.
"""
ma = getmask(a)
mb = getmask(b)
if ma is nomask and mb is nomask:
m = nomask
else:
ma = getmaskarray(a)
mb = getmaskarray(b)
m = umath.logical_or.outer(ma, mb)
if (not m.ndim) and m:
return masked
rcls = get_masked_subclass(a,b)
# We could fill the arguments first, butis it useful ?
# d = self.f.outer(filled(a, self.fillx), filled(b, self.filly)).view(rcls)
d = self.f.outer(getdata(a), getdata(b)).view(rcls)
if d.ndim > 0:
d._mask = m
return d
def accumulate (self, target, axis=0):
"""Accumulate `target` along `axis` after filling with y fill
value.
"""
if isinstance(target, MaskedArray):
tclass = type(target)
else:
tclass = masked_array
t = filled(target, self.filly)
return self.f.accumulate(t, axis).view(tclass)
def __str__ (self):
return "Masked version of " + str(self.f)
#..............................................................................
class _DomainedBinaryOperation:
"""Define binary operations that have a domain, like divide.
They have no reduce, outer or accumulate.
Parameters
----------
f : function.
domain : Default domain.
fillx : Default filling value for the first argument (0).
filly : Default filling value for the second argument (0).
"""
def __init__ (self, dbfunc, domain, fillx=0, filly=0):
"""abfunc(fillx, filly) must be defined.
abfunc(x, filly) = x for all x to enable reduce.
"""
self.f = dbfunc
self.domain = domain
self.fillx = fillx
self.filly = filly
self.__doc__ = getattr(dbfunc, "__doc__", str(dbfunc))
self.__name__ = getattr(dbfunc, "__name__", str(dbfunc))
ufunc_domain[dbfunc] = domain
ufunc_fills[dbfunc] = (fillx, filly)
def __call__(self, a, b):
"Execute the call behavior."
ma = getmask(a)
mb = getmask(b)
d1 = getdata(a)
d2 = get_data(b)
t = narray(self.domain(d1, d2), copy=False)
if t.any(None):
mb = mask_or(mb, t)
# The following line controls the domain filling
d2 = numpy.where(t,self.filly,d2)
m = mask_or(ma, mb)
if (not m.ndim) and m:
return masked
result = self.f(d1, d2).view(get_masked_subclass(a,b))
if result.ndim > 0:
result._mask = m
if isinstance(a,MaskedArray):
result._update_from(a)
if isinstance(b,MaskedArray):
result._update_from(b)
return result
def __str__ (self):
return "Masked version of " + str(self.f)
#..............................................................................
# Unary ufuncs
exp = _MaskedUnaryOperation(umath.exp)
conjugate = _MaskedUnaryOperation(umath.conjugate)
sin = _MaskedUnaryOperation(umath.sin)
cos = _MaskedUnaryOperation(umath.cos)
tan = _MaskedUnaryOperation(umath.tan)
arctan = _MaskedUnaryOperation(umath.arctan)
arcsinh = _MaskedUnaryOperation(umath.arcsinh)
sinh = _MaskedUnaryOperation(umath.sinh)
cosh = _MaskedUnaryOperation(umath.cosh)
tanh = _MaskedUnaryOperation(umath.tanh)
abs = absolute = _MaskedUnaryOperation(umath.absolute)
fabs = _MaskedUnaryOperation(umath.fabs)
negative = _MaskedUnaryOperation(umath.negative)
floor = _MaskedUnaryOperation(umath.floor)
ceil = _MaskedUnaryOperation(umath.ceil)
around = _MaskedUnaryOperation(fromnumeric.round_)
logical_not = _MaskedUnaryOperation(umath.logical_not)
# Domained unary ufuncs .......................................................
sqrt = _MaskedUnaryOperation(umath.sqrt, 0.0,
_DomainGreaterEqual(0.0))
log = _MaskedUnaryOperation(umath.log, 1.0,
_DomainGreater(0.0))
log10 = _MaskedUnaryOperation(umath.log10, 1.0,
_DomainGreater(0.0))
tan = _MaskedUnaryOperation(umath.tan, 0.0,
_DomainTan(1.e-35))
arcsin = _MaskedUnaryOperation(umath.arcsin, 0.0,
_DomainCheckInterval(-1.0, 1.0))
arccos = _MaskedUnaryOperation(umath.arccos, 0.0,
_DomainCheckInterval(-1.0, 1.0))
arccosh = _MaskedUnaryOperation(umath.arccosh, 1.0,
_DomainGreaterEqual(1.0))
arctanh = _MaskedUnaryOperation(umath.arctanh, 0.0,
_DomainCheckInterval(-1.0+1e-15, 1.0-1e-15))
# Binary ufuncs ...............................................................
add = _MaskedBinaryOperation(umath.add)
subtract = _MaskedBinaryOperation(umath.subtract)
multiply = _MaskedBinaryOperation(umath.multiply, 1, 1)
arctan2 = _MaskedBinaryOperation(umath.arctan2, 0.0, 1.0)
equal = _MaskedBinaryOperation(umath.equal)
equal.reduce = None
not_equal = _MaskedBinaryOperation(umath.not_equal)
not_equal.reduce = None
less_equal = _MaskedBinaryOperation(umath.less_equal)
less_equal.reduce = None
greater_equal = _MaskedBinaryOperation(umath.greater_equal)
greater_equal.reduce = None
less = _MaskedBinaryOperation(umath.less)
less.reduce = None
greater = _MaskedBinaryOperation(umath.greater)
greater.reduce = None
logical_and = _MaskedBinaryOperation(umath.logical_and)
alltrue = _MaskedBinaryOperation(umath.logical_and, 1, 1).reduce
logical_or = _MaskedBinaryOperation(umath.logical_or)
sometrue = logical_or.reduce
logical_xor = _MaskedBinaryOperation(umath.logical_xor)
bitwise_and = _MaskedBinaryOperation(umath.bitwise_and)
bitwise_or = _MaskedBinaryOperation(umath.bitwise_or)
bitwise_xor = _MaskedBinaryOperation(umath.bitwise_xor)
hypot = _MaskedBinaryOperation(umath.hypot)
# Domained binary ufuncs ......................................................
divide = _DomainedBinaryOperation(umath.divide, _DomainSafeDivide(), 0, 1)
true_divide = _DomainedBinaryOperation(umath.true_divide,
_DomainSafeDivide(), 0, 1)
floor_divide = _DomainedBinaryOperation(umath.floor_divide,
_DomainSafeDivide(), 0, 1)
remainder = _DomainedBinaryOperation(umath.remainder,
_DomainSafeDivide(), 0, 1)
fmod = _DomainedBinaryOperation(umath.fmod, _DomainSafeDivide(), 0, 1)
#####--------------------------------------------------------------------------
#---- --- Mask creation functions ---
#####--------------------------------------------------------------------------
def get_mask(a):
"""Return the mask of a, if any, or nomask.
To get a full array of booleans of the same shape as a, use
getmaskarray.
"""
return getattr(a, '_mask', nomask)
getmask = get_mask
def getmaskarray(a):
"""Return the mask of a, if any, or a boolean array of the shape
of a, full of False.
"""
m = getmask(a)
if m is nomask:
m = make_mask_none(fromnumeric.shape(a))
return m
def is_mask(m):
"""Return True if m is a legal mask.
Does not check contents, only type.
"""
try:
return m.dtype.type is MaskType
except AttributeError:
return False
#
def make_mask(m, copy=False, shrink=True, flag=None):
"""Return m as a mask, creating a copy if necessary or requested.
The function can accept any sequence of integers or nomask. Does
not check that contents must be 0s and 1s.
Parameters
----------
m : array_like
Potential mask.
copy : bool
Whether to return a copy of m (True) or m itself (False).
shrink : bool
Whether to shrink m to nomask if all its values are False.
"""
if flag is not None:
warnings.warn("The flag 'flag' is now called 'shrink'!",
DeprecationWarning)
shrink = flag
if m is nomask:
return nomask
elif isinstance(m, ndarray):
m = filled(m, True)
if m.dtype.type is MaskType:
if copy:
result = narray(m, dtype=MaskType, copy=copy)
else:
result = m
else:
result = narray(m, dtype=MaskType)
else:
result = narray(filled(m, True), dtype=MaskType)
# Bas les masques !
if shrink and not result.any():
return nomask
else:
return result
def make_mask_none(s):
"""Return a mask of shape s, filled with False.
Parameters
----------
s : tuple
A tuple indicating the shape of the final mask.
"""
result = numeric.zeros(s, dtype=MaskType)
return result
def mask_or (m1, m2, copy=False, shrink=True):
"""Return the combination of two masks m1 and m2.
The masks are combined with the *logical_or* operator, treating
nomask as False. The result may equal m1 or m2 if the other is
nomask.
Parameters
----------
m1 : array_like
First mask.
m2 : array_like
Second mask
copy : bool
Whether to return a copy.
shrink : bool
Whether to shrink m to nomask if all its values are False.
"""
if m1 is nomask:
return make_mask(m2, copy=copy, shrink=shrink)
if m2 is nomask:
return make_mask(m1, copy=copy, shrink=shrink)
if m1 is m2 and is_mask(m1):
return m1
return make_mask(umath.logical_or(m1, m2), copy=copy, shrink=shrink)
#####--------------------------------------------------------------------------
#--- --- Masking functions ---
#####--------------------------------------------------------------------------
def masked_where(condition, a, copy=True):
"""Return a as an array masked where condition is true.
Masked values of a or condition are kept.
Parameters
----------
condition : array_like
Masking condition.
a : array_like
Array to mask.
copy : bool
Whether to return a copy of a (True) or modify a in place.
"""
cond = make_mask(condition)
a = narray(a, copy=copy, subok=True)
if hasattr(a, '_mask'):
cond = mask_or(cond, a._mask)
cls = type(a)
else:
cls = MaskedArray
result = a.view(cls)
result._mask = cond
return result
def masked_greater(x, value, copy=True):
"Shortcut to masked_where, with condition = (x > value)."
return masked_where(greater(x, value), x, copy=copy)
def masked_greater_equal(x, value, copy=True):
"Shortcut to masked_where, with condition = (x >= value)."
return masked_where(greater_equal(x, value), x, copy=copy)
def masked_less(x, value, copy=True):
"Shortcut to masked_where, with condition = (x < value)."
return masked_where(less(x, value), x, copy=copy)
def masked_less_equal(x, value, copy=True):
"Shortcut to masked_where, with condition = (x <= value)."
return masked_where(less_equal(x, value), x, copy=copy)
def masked_not_equal(x, value, copy=True):
"Shortcut to masked_where, with condition = (x != value)."
return masked_where((x != value), x, copy=copy)
#
def masked_equal(x, value, copy=True):
"""Shortcut to masked_where, with condition = (x == value). For
floating point, consider `masked_values(x, value)` instead.
"""
# An alternative implementation relies on filling first: probably not needed.
# d = filled(x, 0)
# c = umath.equal(d, value)
# m = mask_or(c, getmask(x))
# return array(d, mask=m, copy=copy)
return masked_where((x == value), x, copy=copy)
def masked_inside(x, v1, v2, copy=True):
"""Shortcut to masked_where, where condition is True for x inside
the interval [v1,v2] (v1 <= x <= v2). The boundaries v1 and v2
can be given in either order.
Notes
-----
The array x is prefilled with its filling value.
"""
if v2 < v1:
(v1, v2) = (v2, v1)
xf = filled(x)
condition = (xf >= v1) & (xf <= v2)
return masked_where(condition, x, copy=copy)
def masked_outside(x, v1, v2, copy=True):
"""Shortcut to masked_where, where condition is True for x outside
the interval [v1,v2] (x < v1)|(x > v2). The boundaries v1 and v2
can be given in either order.
Notes
-----
The array x is prefilled with its filling value.
"""
if v2 < v1:
(v1, v2) = (v2, v1)
xf = filled(x)
condition = (xf < v1) | (xf > v2)
return masked_where(condition, x, copy=copy)
#
def masked_object(x, value, copy=True):
"""Mask the array x where the data are exactly equal to value.
This function is suitable only for object arrays: for floating
point, please use ``masked_values`` instead.
Notes
-----
The mask is set to `nomask` if posible.
"""
if isMaskedArray(x):
condition = umath.equal(x._data, value)
mask = x._mask
else:
condition = umath.equal(fromnumeric.asarray(x), value)
mask = nomask
mask = mask_or(mask, make_mask(condition, shrink=True))
return masked_array(x, mask=mask, copy=copy, fill_value=value)
def masked_values(x, value, rtol=1.e-5, atol=1.e-8, copy=True):
"""Mask the array x where the data are approximately equal in
value, i.e.
(abs(x - value) <= atol+rtol*abs(value))
Suitable only for floating points. For integers, please use
``masked_equal``. The mask is set to nomask if posible.
Parameters
----------
x : array_like
Array to fill.
value : float
Masking value.
rtol : float
Tolerance parameter.
atol : float
Tolerance parameter (1e-8).
copy : bool
Whether to return a copy of x.
"""
abs = umath.absolute
xnew = filled(x, value)
if issubclass(xnew.dtype.type, numeric.floating):
condition = umath.less_equal(abs(xnew-value), atol+rtol*abs(value))
mask = getattr(x, '_mask', nomask)
else:
condition = umath.equal(xnew, value)
mask = nomask
mask = mask_or(mask, make_mask(condition, shrink=True))
return masked_array(xnew, mask=mask, copy=copy, fill_value=value)
def masked_invalid(a, copy=True):
"""Mask the array for invalid values (nans or infs). Any
preexisting mask is conserved.
"""
a = narray(a, copy=copy, subok=True)
condition = ~(numpy.isfinite(a))
if hasattr(a, '_mask'):
condition = mask_or(condition, a._mask)
cls = type(a)
else:
cls = MaskedArray
result = a.view(cls)
result._mask = condition
return result
#####--------------------------------------------------------------------------
#---- --- Printing options ---
#####--------------------------------------------------------------------------
class _MaskedPrintOption:
"""Handle the string used to represent missing data in a masked
array.
"""
def __init__ (self, display):
"Create the masked_print_option object."
self._display = display
self._enabled = True
def display(self):
"Display the string to print for masked values."
return self._display
def set_display (self, s):
"Set the string to print for masked values."
self._display = s
def enabled(self):
"Is the use of the display value enabled?"
return self._enabled
def enable(self, shrink=1):
"Set the enabling shrink to `shrink`."
self._enabled = shrink
def __str__ (self):
return str(self._display)
__repr__ = __str__
#if you single index into a masked location you get this object.
masked_print_option = _MaskedPrintOption('--')
#####--------------------------------------------------------------------------
#---- --- MaskedArray class ---
#####--------------------------------------------------------------------------
#...............................................................................
class _arraymethod(object):
"""Define a wrapper for basic array methods.
Upon call, returns a masked array, where the new _data array is
the output of the corresponding method called on the original
_data.
If onmask is True, the new mask is the output of the method called
on the initial mask. Otherwise, the new mask is just a reference
to the initial mask.
Parameters
----------
_name : String
Name of the function to apply on data.
_onmask : bool
Whether the mask must be processed also (True) or left
alone (False). Default: True.
obj : Object
The object calling the arraymethod.
"""
def __init__(self, funcname, onmask=True):
self._name = funcname
self._onmask = onmask
self.obj = None
self.__doc__ = self.getdoc()
#
def getdoc(self):
"Return the doc of the function (from the doc of the method)."
methdoc = getattr(ndarray, self._name, None)
methdoc = getattr(numpy, self._name, methdoc)
if methdoc is not None:
return methdoc.__doc__
#
def __get__(self, obj, objtype=None):
self.obj = obj
return self
#
def __call__(self, *args, **params):
methodname = self._name
data = self.obj._data
mask = self.obj._mask
cls = type(self.obj)
result = getattr(data, methodname)(*args, **params).view(cls)
result._update_from(self.obj)
if result.ndim:
if not self._onmask:
result.__setmask__(mask)
elif mask is not nomask:
result.__setmask__(getattr(mask, methodname)(*args, **params))
else:
if mask.ndim and mask.all():
return masked
return result
#..........................................................
class FlatIter(object):
"Define an interator."
def __init__(self, ma):
self.ma = ma
self.ma_iter = numpy.asarray(ma).flat
if ma._mask is nomask:
self.maskiter = None
else:
self.maskiter = ma._mask.flat
def __iter__(self):
return self
### This won't work is ravel makes a copy
def __setitem__(self, index, value):
a = self.ma.ravel()
a[index] = value
def next(self):
d = self.ma_iter.next()
if self.maskiter is not None and self.maskiter.next():
d = masked
return d
class MaskedArray(numeric.ndarray):
"""Arrays with possibly masked values. Masked values of True
exclude the corresponding element from any computation.
Construction:
x = MaskedArray(data, mask=nomask, dtype=None, copy=True,
fill_value=None, keep_mask=True, hard_mask=False, shrink=True)
Parameters
----------
data : {var}
Input data.
mask : {nomask, sequence}
Mask. Must be convertible to an array of booleans with
the same shape as data: True indicates a masked (eg.,
invalid) data.
dtype : dtype
Data type of the output. If None, the type of the data
argument is used. If dtype is not None and different from
data.dtype, a copy is performed.
copy : bool
Whether to copy the input data (True), or to use a
reference instead. Note: data are NOT copied by default.
subok : {True, boolean}
Whether to return a subclass of MaskedArray (if possible)
or a plain MaskedArray.
ndmin : {0, int}
Minimum number of dimensions
fill_value : {var}
Value used to fill in the masked values when necessary. If
None, a default based on the datatype is used.
keep_mask : {True, boolean}
Whether to combine mask with the mask of the input data,
if any (True), or to use only mask for the output (False).
hard_mask : {False, boolean}
Whether to use a hard mask or not. With a hard mask,
masked values cannot be unmasked.
shrink : {True, boolean}
Whether to force compression of an empty mask.
"""
__array_priority__ = 15
_defaultmask = nomask
_defaulthardmask = False
_baseclass = numeric.ndarray
def __new__(cls, data=None, mask=nomask, dtype=None, copy=False,
subok=True, ndmin=0, fill_value=None,
keep_mask=True, hard_mask=False, flag=None,shrink=True,
**options):
"""Create a new masked array from scratch.
Note: you can also create an array with the .view(MaskedArray)
method.
"""
if flag is not None:
warnings.warn("The flag 'flag' is now called 'shrink'!",
DeprecationWarning)
shrink = flag
# Process data............
_data = narray(data, dtype=dtype, copy=copy, subok=True, ndmin=ndmin)
_baseclass = getattr(data, '_baseclass', type(_data))
_basedict = getattr(data, '_basedict', getattr(data, '__dict__', {}))
if not isinstance(data, MaskedArray) or not subok:
_data = _data.view(cls)
else:
_data = _data.view(type(data))
# Backwards compatibility w/ numpy.core.ma .......
if hasattr(data,'_mask') and not isinstance(data, ndarray):
_data._mask = data._mask
_sharedmask = True
# Process mask ...........
if mask is nomask:
if not keep_mask:
if shrink:
_data._mask = nomask
else:
_data._mask = make_mask_none(_data)
if copy:
_data._mask = _data._mask.copy()
_data._sharedmask = False
else:
_data._sharedmask = True
else:
mask = narray(mask, dtype=MaskType, copy=copy)
if mask.shape != _data.shape:
(nd, nm) = (_data.size, mask.size)
if nm == 1:
mask = numeric.resize(mask, _data.shape)
elif nm == nd:
mask = fromnumeric.reshape(mask, _data.shape)
else:
msg = "Mask and data not compatible: data size is %i, "+\
"mask size is %i."
raise MAError, msg % (nd, nm)
copy = True
if _data._mask is nomask:
_data._mask = mask
_data._sharedmask = not copy
else:
if not keep_mask:
_data._mask = mask
_data._sharedmask = not copy
else:
_data._mask = umath.logical_or(mask, _data._mask)
_data._sharedmask = False
# Update fill_value.......
_data._fill_value = _check_fill_value(fill_value, _data.dtype)
# Process extra options ..
_data._hardmask = hard_mask
_data._baseclass = _baseclass
_data._basedict = _basedict
return _data
#
def _update_from(self, obj):
"""Copies some attributes of obj to self.
"""
if obj is not None and isinstance(obj,ndarray):
_baseclass = type(obj)
else:
_baseclass = ndarray
_basedict = getattr(obj,'_basedict',getattr(obj,'__dict__',{}))
_dict = dict(_fill_value=getattr(obj, '_fill_value', None),
_hardmask=getattr(obj, '_hardmask', False),
_sharedmask=getattr(obj, '_sharedmask', False),
_baseclass=getattr(obj,'_baseclass',_baseclass),
_basedict=_basedict,)
self.__dict__.update(_dict)
self.__dict__.update(_basedict)
return
#........................
def __array_finalize__(self,obj):
"""Finalizes the masked array.
"""
# Get main attributes .........
self._update_from(obj)
self._mask = getattr(obj, '_mask', nomask)
# Finalize the mask ...........
if self._mask is not nomask:
self._mask.shape = self.shape
return
#..................................
def __array_wrap__(self, obj, context=None):
"""Special hook for ufuncs.
Wraps the numpy array and sets the mask according to context.
"""
result = obj.view(type(self))
result._update_from(self)
#..........
if context is not None:
result._mask = result._mask.copy()
(func, args, _) = context
m = reduce(mask_or, [getmaskarray(arg) for arg in args])
# Get the domain mask................
domain = ufunc_domain.get(func, None)
if domain is not None:
if len(args) > 2:
d = reduce(domain, args)
else:
d = domain(*args)
# Fill the result where the domain is wrong
try:
# Binary domain: take the last value
fill_value = ufunc_fills[func][-1]
except TypeError:
# Unary domain: just use this one
fill_value = ufunc_fills[func]
except KeyError:
# Domain not recognized, use fill_value instead
fill_value = self.fill_value
result = result.copy()
numpy.putmask(result, d, fill_value)
# Update the mask
if m is nomask:
if d is not nomask:
m = d
else:
m |= d
# Make sure the mask has the proper size
if result.shape == () and m:
return masked
else:
result._mask = m
result._sharedmask = False
#....
return result
#.............................................
def __getitem__(self, indx):
"""x.__getitem__(y) <==> x[y]
Return the item described by i, as a masked array.
"""
# This test is useful, but we should keep things light...
# if getmask(indx) is not nomask:
# msg = "Masked arrays must be filled before they can be used as indices!"
# raise IndexError, msg
dout = ndarray.__getitem__(self.view(ndarray), indx)
# We could directly use ndarray.__getitem__ on self...
# But then we would have to modify __array_finalize__ to prevent the
# mask of being reshaped if it hasn't been set up properly yet...
# So it's easier to stick to the current version
m = self._mask
if not getattr(dout,'ndim', False):
# Just a scalar............
if m is not nomask and m[indx]:
return masked
else:
# Force dout to MA ........
dout = dout.view(type(self))
# Inherit attributes from self
dout._update_from(self)
# Check the fill_value ....
if isinstance(indx, basestring):
fvindx = list(self.dtype.names).index(indx)
dout._fill_value = self.fill_value[fvindx]
# Update the mask if needed
if m is not nomask:
if isinstance(indx, basestring):
dout._mask = m.reshape(dout.shape)
else:
dout._mask = ndarray.__getitem__(m, indx).reshape(dout.shape)
# Note: Don't try to check for m.any(), that'll take too long...
# mask = ndarray.__getitem__(m, indx).reshape(dout.shape)
# if self._shrinkmask and not m.any():
# dout._mask = nomask
# else:
# dout._mask = mask
return dout
#........................
def __setitem__(self, indx, value):
"""x.__setitem__(i, y) <==> x[i]=y
Set item described by index. If value is masked, masks those
locations.
"""
if self is masked:
raise MAError, 'Cannot alter the masked element.'
# This test is useful, but we should keep things light...
# if getmask(indx) is not nomask:
# msg = "Masked arrays must be filled before they can be used as indices!"
# raise IndexError, msg
if isinstance(indx, basestring):
ndarray.__setitem__(self._data,indx, getdata(value))
warnings.warn("MaskedArray.__setitem__ on fields: "\
"The mask is NOT affected!")
return
#....
if value is masked:
m = self._mask
if m is nomask:
m = numpy.zeros(self.shape, dtype=MaskType)
m[indx] = True
self._mask = m
self._sharedmask = False
return
#....
dval = narray(value, copy=False, dtype=self.dtype)
valmask = getmask(value)
if self._mask is nomask:
# Set the data, then the mask
ndarray.__setitem__(self._data,indx,dval)
if valmask is not nomask:
self._mask = numpy.zeros(self.shape, dtype=MaskType)
self._mask[indx] = valmask
elif not self._hardmask:
# Unshare the mask if necessary to avoid propagation
self.unshare_mask()
# Set the data, then the mask
ndarray.__setitem__(self._data,indx,dval)
self._mask[indx] = valmask
elif hasattr(indx, 'dtype') and (indx.dtype==bool_):
indx = indx * umath.logical_not(self._mask)
ndarray.__setitem__(self._data,indx,dval)
else:
mindx = mask_or(self._mask[indx], valmask, copy=True)
dindx = self._data[indx]
if dindx.size > 1:
dindx[~mindx] = dval
elif mindx is nomask:
dindx = dval
ndarray.__setitem__(self._data,indx,dindx)
self._mask[indx] = mindx
#............................................
def __getslice__(self, i, j):
"""x.__getslice__(i, j) <==> x[i:j]
Return the slice described by (i, j). The use of negative
indices is not supported.
"""
return self.__getitem__(slice(i,j))
#........................
def __setslice__(self, i, j, value):
"""x.__setslice__(i, j, value) <==> x[i:j]=value
Set the slice (i,j) of a to value. If value is masked, mask
those locations.
"""
self.__setitem__(slice(i,j), value)
#............................................
def __setmask__(self, mask, copy=False):
"""Set the mask.
"""
if mask is not nomask:
mask = narray(mask, copy=copy, dtype=MaskType)
# We could try to check whether shrinking is needed..
# ... but we would waste some precious time
# if self._shrinkmask and not mask.any():
# mask = nomask
if self._mask is nomask:
self._mask = mask
elif self._hardmask:
if mask is not nomask:
self._mask.__ior__(mask)
else:
# This one is tricky: if we set the mask that way, we may break the
# propagation. But if we don't, we end up with a mask full of False
# and a test on nomask fails...
if mask is nomask:
self._mask = nomask
else:
self.unshare_mask()
self._mask.flat = mask
if self._mask.shape:
self._mask = numeric.reshape(self._mask, self.shape)
_set_mask = __setmask__
#....
def _get_mask(self):
"""Return the current mask.
"""
# We could try to force a reshape, but that wouldn't work in some cases.
# return self._mask.reshape(self.shape)
return self._mask
mask = property(fget=_get_mask, fset=__setmask__, doc="Mask")
#............................................
def harden_mask(self):
"""Force the mask to hard.
"""
self._hardmask = True
def soften_mask(self):
"""Force the mask to soft.
"""
self._hardmask = False
def unshare_mask(self):
"""Copy the mask and set the sharedmask flag to False.
"""
if self._sharedmask:
self._mask = self._mask.copy()
self._sharedmask = False
def shrink_mask(self):
"""Reduce a mask to nomask when possible.
"""
m = self._mask
if m.ndim and not m.any():
self._mask = nomask
#............................................
def _get_data(self):
"""Return the current data, as a view of the original
underlying data.
"""
return self.view(self._baseclass)
_data = property(fget=_get_data)
data = property(fget=_get_data)
def raw_data(self):
"""Return the _data part of the MaskedArray.
DEPRECATED: You should really use ``.data`` instead...
"""
warnings.warn('Use .data instead.', DeprecationWarning)
return self._data
#............................................
def _get_flat(self):
"""Return a flat iterator.
"""
return FlatIter(self)
#
def _set_flat (self, value):
"""Set a flattened version of self to value.
"""
y = self.ravel()
y[:] = value
#
flat = property(fget=_get_flat, fset=_set_flat,
doc="Flat version of the array.")
#............................................
def get_fill_value(self):
"""Return the filling value.
"""
if self._fill_value is None:
self._fill_value = _check_fill_value(None, self.dtype)
return self._fill_value
def set_fill_value(self, value=None):
"""Set the filling value to value.
If value is None, use a default based on the data type.
"""
self._fill_value = _check_fill_value(value,self.dtype)
fill_value = property(fget=get_fill_value, fset=set_fill_value,
doc="Filling value.")
def filled(self, fill_value=None):
"""Return a copy of self._data, where masked values are filled
with fill_value.
If fill_value is None, self.fill_value is used instead.
Notes
-----
+ Subclassing is preserved
+ The result is NOT a MaskedArray !
Examples
--------
>>> x = array([1,2,3,4,5], mask=[0,0,1,0,1], fill_value=-999)
>>> x.filled()
array([1,2,-999,4,-999])
>>> type(x.filled())
<type 'numpy.ndarray'>
"""
m = self._mask
if m is nomask or not m.any():
return self._data
#
if fill_value is None:
fill_value = self.fill_value
#
if self is masked_singleton:
result = numeric.asanyarray(fill_value)
else:
result = self._data.copy()
try:
numpy.putmask(result, m, fill_value)
except (TypeError, AttributeError):
fill_value = narray(fill_value, dtype=object)
d = result.astype(object)
result = fromnumeric.choose(m, (d, fill_value))
except IndexError:
#ok, if scalar
if self._data.shape:
raise
elif m:
result = narray(fill_value, dtype=self.dtype)
else:
result = self._data
return result
def compressed(self):
"""Return a 1-D array of all the non-masked data.
"""
data = ndarray.ravel(self._data)
if self._mask is not nomask:
data = data.compress(numpy.logical_not(ndarray.ravel(self._mask)))
return data
def compress(self, condition, axis=None, out=None):
"""Return a where condition is True.
If condition is a MaskedArray, missing values are considered as False.
Returns
-------
A MaskedArray object.
Notes
-----
Please note the difference with compressed() !
The output of compress has a mask, the output of compressed does not.
"""
# Get the basic components
(_data, _mask) = (self._data, self._mask)
# Force the condition to a regular ndarray (forget the missing values...)
condition = narray(condition, copy=False, subok=False)
#
_new = _data.compress(condition, axis=axis, out=out).view(type(self))
_new._update_from(self)
if _mask is not nomask:
_new._mask = _mask.compress(condition, axis=axis)
return _new
#............................................
def __str__(self):
"""String representation.
"""
if masked_print_option.enabled():
f = masked_print_option
if self is masked:
return str(f)
m = self._mask
if m is nomask:
res = self._data
else:
if m.shape == ():
if m:
return str(f)
else:
return str(self._data)
# convert to object array to make filled work
#CHECK: the two lines below seem more robust than the self._data.astype
# res = numeric.empty(self._data.shape, object_)
# numeric.putmask(res,~m,self._data)
res = self._data.astype("|O8")
res[m] = f
else:
res = self.filled(self.fill_value)
return str(res)
def __repr__(self):
"""Literal string representation.
"""
with_mask = """\
masked_%(name)s(data =
%(data)s,
mask =
%(mask)s,
fill_value=%(fill)s)
"""
with_mask1 = """\
masked_%(name)s(data = %(data)s,
mask = %(mask)s,
fill_value=%(fill)s)
"""
n = len(self.shape)
name = repr(self._data).split('(')[0]
if n <= 1:
return with_mask1 % {
'name': name,
'data': str(self),
'mask': str(self._mask),
'fill': str(self.fill_value),
}
return with_mask % {
'name': name,
'data': str(self),
'mask': str(self._mask),
'fill': str(self.fill_value),
}
#............................................
def __add__(self, other):
"Add other to self, and return a new masked array."
return add(self, other)
#
def __sub__(self, other):
"Subtract other to self, and return a new masked array."
return subtract(self, other)
#
def __mul__(self, other):
"Multiply other by self, and return a new masked array."
return multiply(self, other)
#
def __div__(self, other):
"Divide other into self, and return a new masked array."
return divide(self, other)
#
def __truediv__(self, other):
"Divide other into self, and return a new masked array."
return true_divide(self, other)
#
def __floordiv__(self, other):
"Divide other into self, and return a new masked array."
return floor_divide(self, other)
#
def __pow__(self, other):
"Raise self to the power other, masking the potential NaNs/Infs"
return power(self, other)
#............................................
def __iadd__(self, other):
"Add other to self in-place."
ndarray.__iadd__(self._data, getdata(other))
m = getmask(other)
if self._mask is nomask:
self._mask = m
elif m is not nomask:
self._mask += m
return self
#....
def __isub__(self, other):
"Subtract other from self in-place."
ndarray.__isub__(self._data, getdata(other))
m = getmask(other)
if self._mask is nomask:
self._mask = m
elif m is not nomask:
self._mask += m
return self
#....
def __imul__(self, other):
"Multiply self by other in-place."
ndarray.__imul__(self._data, getdata(other))
m = getmask(other)
if self._mask is nomask:
self._mask = m
elif m is not nomask:
self._mask += m
return self
#....
def __idiv__(self, other):
"Divide self by other in-place."
other_data = getdata(other)
dom_mask = _DomainSafeDivide().__call__(self._data, other_data)
other_mask = getmask(other)
new_mask = mask_or(other_mask, dom_mask)
# The following 3 lines control the domain filling
if dom_mask.any():
other_data = other_data.copy()
numpy.putmask(other_data, dom_mask, 1)
ndarray.__idiv__(self._data, other_data)
self._mask = mask_or(self._mask, new_mask)
return self
#...
def __ipow__(self, other):
"Raise self to the power other, in place"
_data = self._data
other_data = getdata(other)
other_mask = getmask(other)
ndarray.__ipow__(_data, other_data)
invalid = numpy.logical_not(numpy.isfinite(_data))
new_mask = mask_or(other_mask,invalid)
self._mask = mask_or(self._mask, new_mask)
# The following line is potentially problematic, as we change _data...
numpy.putmask(self._data,invalid,self.fill_value)
return self
#............................................
def __float__(self):
"Convert to float."
if self.size > 1:
raise TypeError,\
"Only length-1 arrays can be converted to Python scalars"
elif self._mask:
warnings.warn("Warning: converting a masked element to nan.")
return numpy.nan
return float(self.item())
def __int__(self):
"Convert to int."
if self.size > 1:
raise TypeError,\
"Only length-1 arrays can be converted to Python scalars"
elif self._mask:
raise MAError, 'Cannot convert masked element to a Python int.'
return int(self.item())
#............................................
def get_imag(self):
result = self._data.imag.view(type(self))
result.__setmask__(self._mask)
return result
imag = property(fget=get_imag,doc="Imaginary part")
def get_real(self):
result = self._data.real.view(type(self))
result.__setmask__(self._mask)
return result
real = property(fget=get_real,doc="Real part")
#............................................
def count(self, axis=None):
"""Count the non-masked elements of the array along the given
axis.
Parameters
----------
axis : int, optional
Axis along which to count the non-masked elements. If
not given, all the non masked elements are counted.
Returns
-------
A masked array where the mask is True where all data are
masked. If axis is None, returns either a scalar ot the
masked singleton if all values are masked.
"""
m = self._mask
s = self.shape
ls = len(s)
if m is nomask:
if ls == 0:
return 1
if ls == 1:
return s[0]
if axis is None:
return self.size
else:
n = s[axis]
t = list(s)
del t[axis]
return numeric.ones(t) * n
n1 = numpy.size(m, axis)
n2 = m.astype(int_).sum(axis)
if axis is None:
return (n1-n2)
else:
return narray(n1 - n2)
#............................................
flatten = _arraymethod('flatten')
#
def ravel(self):
"""Returns a 1D version of self, as a view."""
r = ndarray.ravel(self._data).view(type(self))
r._update_from(self)
if self._mask is not nomask:
r._mask = ndarray.ravel(self._mask).reshape(r.shape)
else:
r._mask = nomask
return r
#
repeat = _arraymethod('repeat')
#
def reshape (self, *s):
"""Reshape the array to shape s.
Returns
-------
A new masked array.
Notes
-----
If you want to modify the shape in place, please use
``a.shape = s``
"""
result = self._data.reshape(*s).view(type(self))
result.__dict__.update(self.__dict__)
if result._mask is not nomask:
result._mask = self._mask.copy()
result._mask.shape = result.shape
return result
#
def resize(self, newshape, refcheck=True, order=False):
"""Attempt to modify the size and the shape of the array in place.
The array must own its own memory and not be referenced by
other arrays.
Returns
-------
None.
"""
try:
self._data.resize(newshape, refcheck, order)
if self.mask is not nomask:
self._mask.resize(newshape, refcheck, order)
except ValueError:
raise ValueError("Cannot resize an array that has been referenced "
"or is referencing another array in this way.\n"
"Use the resize function.")
return None
#
def put(self, indices, values, mode='raise'):
"""Set storage-indexed locations to corresponding values.
a.put(values, indices, mode) sets a.flat[n] = values[n] for
each n in indices. If ``values`` is shorter than ``indices``
then it will repeat. If ``values`` has some masked values, the
initial mask is updated in consequence, else the corresponding
values are unmasked.
"""
m = self._mask
# Hard mask: Get rid of the values/indices that fall on masked data
if self._hardmask and self._mask is not nomask:
mask = self._mask[indices]
indices = narray(indices, copy=False)
values = narray(values, copy=False, subok=True)
values.resize(indices.shape)
indices = indices[~mask]
values = values[~mask]
#....
self._data.put(indices, values, mode=mode)
#....
if m is nomask:
m = getmask(values)
else:
m = m.copy()
if getmask(values) is nomask:
m.put(indices, False, mode=mode)
else:
m.put(indices, values._mask, mode=mode)
m = make_mask(m, copy=False, shrink=True)
self._mask = m
#............................................
def ids (self):
"""Return the addresses of the data and mask areas."""
if self._mask is nomask:
return (self.ctypes.data, id(nomask))
return (self.ctypes.data, self._mask.ctypes.data)
#............................................
def all(self, axis=None, out=None):
"""Return True if all entries along the given axis are True,
False otherwise. Masked values are considered as True during
computation.
Parameter
----------
axis : int, optional
Axis along which the operation is performed. If None,
the operation is performed on a flatten array
out : {MaskedArray}, optional
Alternate optional output. If not None, out should be
a valid MaskedArray of the same shape as the output of
self._data.all(axis).
Returns A masked array, where the mask is True if all data along
-------
the axis are masked.
Notes
-----
An exception is raised if ``out`` is not None and not of the
same type as self.
"""
if out is None:
d = self.filled(True).all(axis=axis).view(type(self))
if d.ndim > 0:
d.__setmask__(self._mask.all(axis))
return d
elif type(out) is not type(self):
raise TypeError("The external array should have " \
"a type %s (got %s instead)" %\
(type(self), type(out)))
self.filled(True).all(axis=axis, out=out)
if out.ndim:
out.__setmask__(self._mask.all(axis))
return out
def any(self, axis=None, out=None):
"""Returns True if at least one entry along the given axis is
True.
Returns False if all entries are False.
Masked values are considered as True during computation.
Parameter
----------
axis : int, optional
Axis along which the operation is performed.
If None, the operation is performed on a flatten array
out : {MaskedArray}, optional
Alternate optional output. If not None, out should be
a valid MaskedArray of the same shape as the output of
self._data.all(axis).
Returns A masked array, where the mask is True if all data along
-------
the axis are masked.
Notes
-----
An exception is raised if ``out`` is not None and not of the
same type as self.
"""
if out is None:
d = self.filled(False).any(axis=axis).view(type(self))
if d.ndim > 0:
d.__setmask__(self._mask.all(axis))
return d
elif type(out) is not type(self):
raise TypeError("The external array should have a type %s "\
"(got %s instead)" %\
(type(self), type(out)))
self.filled(False).any(axis=axis, out=out)
if out.ndim:
out.__setmask__(self._mask.all(axis))
return out
def nonzero(self):
"""Return the indices of the elements of a that are not zero
nor masked, as a tuple of arrays.
There are as many tuples as dimensions of a, each tuple
contains the indices of the non-zero elements in that
dimension. The corresponding non-zero values can be obtained
with ``a[a.nonzero()]``.
To group the indices by element, rather than dimension, use
instead: ``transpose(a.nonzero())``.
The result of this is always a 2d array, with a row for each
non-zero element.
"""
return narray(self.filled(0), copy=False).nonzero()
#............................................
def trace(self, offset=0, axis1=0, axis2=1, dtype=None, out=None):
"""a.trace(offset=0, axis1=0, axis2=1, dtype=None, out=None)
Return the sum along the offset diagonal of the array's
indicated `axis1` and `axis2`.
"""
# TODO: What are we doing with `out`?
m = self._mask
if m is nomask:
result = super(MaskedArray, self).trace(offset=offset, axis1=axis1,
axis2=axis2, out=out)
return result.astype(dtype)
else:
D = self.diagonal(offset=offset, axis1=axis1, axis2=axis2)
return D.astype(dtype).filled(0).sum(axis=None)
#............................................
def sum(self, axis=None, dtype=None):
"""Sum the array over the given axis.
Masked elements are set to 0 internally.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
dtype : dtype, optional
Datatype for the intermediary computation. If not given,
the current dtype is used instead.
"""
if self._mask is nomask:
mask = nomask
else:
mask = self._mask.all(axis)
if (not mask.ndim) and mask:
return masked
result = self.filled(0).sum(axis, dtype=dtype).view(type(self))
if result.ndim > 0:
result.__setmask__(mask)
return result
def cumsum(self, axis=None, dtype=None):
"""Return the cumulative sum of the elements of the array
along the given axis.
Masked values are set to 0 internally.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
dtype : {dtype}, optional
Datatype for the intermediary computation. If not
given, the current dtype is used instead.
"""
result = self.filled(0).cumsum(axis=axis, dtype=dtype).view(type(self))
result.__setmask__(self.mask)
return result
def prod(self, axis=None, dtype=None):
"""Return the product of the elements of the array along the
given axis.
Masked elements are set to 1 internally.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
dtype : {dtype}, optional
Datatype for the intermediary computation. If not
given, the current dtype is used instead.
"""
if self._mask is nomask:
mask = nomask
else:
mask = self._mask.all(axis)
if (not mask.ndim) and mask:
return masked
result = self.filled(1).prod(axis=axis, dtype=dtype).view(type(self))
if result.ndim:
result.__setmask__(mask)
return result
product = prod
def cumprod(self, axis=None, dtype=None):
"""Return the cumulative product of the elements of the array
along the given axis.
Masked values are set to 1 internally.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
dtype : {dtype}, optional
Datatype for the intermediary computation. If not
given, the current dtype is used instead.
"""
result = self.filled(1).cumprod(axis=axis, dtype=dtype).view(type(self))
result.__setmask__(self.mask)
return result
def mean(self, axis=None, dtype=None, out=None):
"""Average the array over the given axis. Equivalent to
a.sum(axis, dtype) / a.size(axis).
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
dtype : {dtype}, optional
Datatype for the intermediary computation. If not
given, the current dtype is used instead.
"""
if self._mask is nomask:
result = super(MaskedArray, self).mean(axis=axis, dtype=dtype)
else:
dsum = self.sum(axis=axis, dtype=dtype)
cnt = self.count(axis=axis)
result = dsum*1./cnt
if out is not None:
out.flat = result.ravel()
return result
def anom(self, axis=None, dtype=None):
"""Return the anomalies (deviations from the average) along
the given axis.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
dtype : {dtype}, optional
Datatype for the intermediary computation. If not
given, the current dtype is used instead.
"""
m = self.mean(axis, dtype)
if not axis:
return (self - m)
else:
return (self - expand_dims(m,axis))
def var(self, axis=None, dtype=None, ddof=0):
"""Return the variance, a measure of the spread of a distribution.
The variance is the average of the squared deviations from the
mean, i.e. var = mean(abs(x - x.mean())**2).
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
dtype : {dtype}, optional
Datatype for the intermediary computation. If not
given, the current dtype is used instead.
Notes
-----
The value returned is by default a biased estimate of the
true variance, since the mean is computed by dividing by N-ddof.
For the (more standard) unbiased estimate, use ddof=1 or.
Note that for complex numbers the absolute value is taken before
squaring, so that the result is always real and nonnegative.
"""
if self._mask is nomask:
# TODO: Do we keep super, or var _data and take a view ?
return super(MaskedArray, self).var(axis=axis, dtype=dtype,
ddof=ddof)
else:
cnt = self.count(axis=axis)-ddof
danom = self.anom(axis=axis, dtype=dtype)
if iscomplexobj(self):
danom = umath.absolute(danom)**2
else:
danom *= danom
dvar = narray(danom.sum(axis) / cnt).view(type(self))
if axis is not None:
dvar._mask = mask_or(self._mask.all(axis), (cnt==1))
dvar._update_from(self)
return dvar
def std(self, axis=None, dtype=None, ddof=0):
"""Return the standard deviation, a measure of the spread of a
distribution.
The standard deviation is the square root of the average of
the squared deviations from the mean, i.e.
std = sqrt(mean(abs(x - x.mean())**2)).
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
dtype : {dtype}, optional
Datatype for the intermediary computation.
If not given, the current dtype is used instead.
Notes
-----
The value returned is by default a biased estimate of the
true standard deviation, since the mean is computed by dividing
by N-ddof. For the more standard unbiased estimate, use ddof=1.
Note that for complex numbers the absolute value is taken before
squaring, so that the result is always real and nonnegative.
"""
dvar = self.var(axis,dtype,ddof=ddof)
if axis is not None or dvar is not masked:
dvar = sqrt(dvar)
return dvar
#............................................
def round(self, decimals=0, out=None):
result = self._data.round(decimals).view(type(self))
result._mask = self._mask
result._update_from(self)
if out is None:
return result
out[:] = result
return
round.__doc__ = ndarray.round.__doc__
#............................................
def argsort(self, axis=None, fill_value=None, kind='quicksort',
order=None):
"""Return an ndarray of indices that sort the array along the
specified axis. Masked values are filled beforehand to
fill_value.
Parameters
----------
axis : int, optional
Axis to be indirectly sorted.
If not given, uses a flatten version of the array.
fill_value : {var}
Value used to fill in the masked values.
If not given, self.fill_value is used instead.
kind : {string}
Sorting algorithm (default 'quicksort')
Possible values: 'quicksort', 'mergesort', or 'heapsort'
Notes
-----
This method executes an indirect sort along the given axis
using the algorithm specified by the kind keyword. It returns
an array of indices of the same shape as 'a' that index data
along the given axis in sorted order.
The various sorts are characterized by average speed, worst
case performance need for work space, and whether they are
stable. A stable sort keeps items with the same key in the
same relative order. The three available algorithms have the
following properties:
|------------------------------------------------------|
| kind | speed | worst case | work space | stable|
|------------------------------------------------------|
|'quicksort'| 1 | O(n^2) | 0 | no |
|'mergesort'| 2 | O(n*log(n)) | ~n/2 | yes |
|'heapsort' | 3 | O(n*log(n)) | 0 | no |
|------------------------------------------------------|
All the sort algorithms make temporary copies of the data when
the sort is not along the last axis. Consequently, sorts along
the last axis are faster and use less space than sorts along
other axis.
"""
if fill_value is None:
fill_value = default_fill_value(self)
d = self.filled(fill_value).view(ndarray)
return d.argsort(axis=axis, kind=kind, order=order)
#........................
def argmin(self, axis=None, fill_value=None):
"""Return an ndarray of indices for the minimum values of a
along the specified axis.
Masked values are treated as if they had the value fill_value.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
fill_value : {var}, optional
Value used to fill in the masked values. If None, the
output of minimum_fill_value(self._data) is used.
"""
if fill_value is None:
fill_value = minimum_fill_value(self)
d = self.filled(fill_value).view(ndarray)
return d.argmin(axis)
#........................
def argmax(self, axis=None, fill_value=None):
"""Returns the array of indices for the maximum values of `a`
along the specified axis.
Masked values are treated as if they had the value fill_value.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
fill_value : {var}, optional
Value used to fill in the masked values. If None, the
output of maximum_fill_value(self._data) is used.
"""
if fill_value is None:
fill_value = maximum_fill_value(self._data)
d = self.filled(fill_value).view(ndarray)
return d.argmax(axis)
def sort(self, axis=-1, kind='quicksort', order=None,
endwith=True, fill_value=None):
"""Sort along the given axis.
Parameters
----------
axis : int
Axis to be indirectly sorted.
kind : {string}
Sorting algorithm (default 'quicksort')
Possible values: 'quicksort', 'mergesort', or 'heapsort'.
order : {var}
If a has fields defined, then the order keyword can be
the field name to sort on or a list (or tuple) of
field names to indicate the order that fields should
be used to define the sort.
fill_value : {var}
Value used to fill in the masked values. If None, use
the the output of minimum_fill_value().
endwith : bool
Whether missing values (if any) should be forced in
the upper indices (at the end of the array) (True) or
lower indices (at the beginning).
Returns
-------
When used as method, returns None.
When used as a function, returns an array.
Notes
-----
This method sorts 'a' in place along the given axis using
the algorithm specified by the kind keyword.
The various sorts may characterized by average speed,
worst case performance need for work space, and whether
they are stable. A stable sort keeps items with the same
key in the same relative order and is most useful when
used w/ argsort where the key might differ from the items
being sorted. The three available algorithms have the
following properties:
|------------------------------------------------------|
| kind | speed | worst case | work space | stable|
|------------------------------------------------------|
|'quicksort'| 1 | O(n^2) | 0 | no |
|'mergesort'| 2 | O(n*log(n)) | ~n/2 | yes |
|'heapsort' | 3 | O(n*log(n)) | 0 | no |
|------------------------------------------------------|
"""
if self._mask is nomask:
ndarray.sort(self,axis=axis, kind=kind, order=order)
else:
if fill_value is None:
if endwith:
filler = minimum_fill_value(self)
else:
filler = maximum_fill_value(self)
else:
filler = fill_value
idx = numpy.indices(self.shape)
idx[axis] = self.filled(filler).argsort(axis=axis,kind=kind,order=order)
idx_l = idx.tolist()
tmp_mask = self._mask[idx_l].flat
tmp_data = self._data[idx_l].flat
self.flat = tmp_data
self._mask.flat = tmp_mask
return
#............................................
def min(self, axis=None, fill_value=None):
"""Return the minimum of a along the given axis.
Masked values are filled with fill_value.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
fill_value : {var}, optional
Value used to fill in the masked values.
If None, use the the output of minimum_fill_value().
"""
mask = self._mask
# Check all/nothing case ......
if mask is nomask:
return super(MaskedArray, self).min(axis=axis)
elif (not mask.ndim) and mask:
return masked
# Get the mask ................
if axis is None:
mask = umath.logical_and.reduce(mask.flat)
else:
mask = umath.logical_and.reduce(mask, axis=axis)
# Skip if all masked ..........
if not mask.ndim and mask:
return masked
# Get the fill value ...........
if fill_value is None:
fill_value = minimum_fill_value(self)
# Get the data ................
result = self.filled(fill_value).min(axis=axis).view(type(self))
if result.ndim > 0:
result._mask = mask
return result
def mini(self, axis=None):
if axis is None:
return minimum(self)
else:
return minimum.reduce(self, axis)
#........................
def max(self, axis=None, fill_value=None):
"""Return the maximum/a along the given axis.
Masked values are filled with fill_value.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
fill_value : {var}, optional
Value used to fill in the masked values.
If None, use the the output of maximum_fill_value().
"""
mask = self._mask
# Check all/nothing case ......
if mask is nomask:
return super(MaskedArray, self).max(axis=axis)
elif (not mask.ndim) and mask:
return masked
# Check the mask ..............
if axis is None:
mask = umath.logical_and.reduce(mask.flat)
else:
mask = umath.logical_and.reduce(mask, axis=axis)
# Skip if all masked ..........
if not mask.ndim and mask:
return masked
# Get the fill value ..........
if fill_value is None:
fill_value = maximum_fill_value(self)
# Get the data ................
result = self.filled(fill_value).max(axis=axis).view(type(self))
if result.ndim > 0:
result._mask = mask
return result
#........................
def ptp(self, axis=None, fill_value=None):
"""Return the visible data range (max-min) along the given axis.
Parameters
----------
axis : int, optional
Axis along which to perform the operation.
If None, applies to a flattened version of the array.
fill_value : {var}, optional
Value used to fill in the masked values. If None, the
maximum uses the maximum default, the minimum uses the
minimum default.
"""
return self.max(axis, fill_value) - self.min(axis, fill_value)
# Array methods ---------------------------------------
copy = _arraymethod('copy')
diagonal = _arraymethod('diagonal')
take = _arraymethod('take')
transpose = _arraymethod('transpose')
T = property(fget=lambda self:self.transpose())
swapaxes = _arraymethod('swapaxes')
clip = _arraymethod('clip', onmask=False)
copy = _arraymethod('copy')
squeeze = _arraymethod('squeeze')
#--------------------------------------------
def tolist(self, fill_value=None):
"""Copy the data portion of the array to a hierarchical python
list and returns that list.
Data items are converted to the nearest compatible Python
type. Masked values are converted to fill_value. If
fill_value is None, the corresponding entries in the output
list will be ``None``.
"""
if fill_value is not None:
return self.filled(fill_value).tolist()
result = self.filled().tolist()
# Set temps to save time when dealing w/ mrecarrays...
_mask = self._mask
if _mask is nomask:
return result
nbdims = self.ndim
dtypesize = len(self.dtype)
if nbdims == 0:
return tuple([None]*dtypesize)
elif nbdims == 1:
maskedidx = _mask.nonzero()[0].tolist()
if dtypesize:
nodata = tuple([None]*dtypesize)
else:
nodata = None
[operator.setitem(result,i,nodata) for i in maskedidx]
else:
for idx in zip(*[i.tolist() for i in _mask.nonzero()]):
tmp = result
for i in idx[:-1]:
tmp = tmp[i]
tmp[idx[-1]] = None
return result
#........................
def tostring(self, fill_value=None, order='C'):
"""Return a copy of array data as a Python string containing the raw
bytes in the array.
Parameters
----------
fill_value : {var}, optional
Value used to fill in the masked values.
If None, uses self.fill_value instead.
order : {string}
Order of the data item in the copy {"C","F","A"}.
"C" -- C order (row major)
"Fortran" -- Fortran order (column major)
"Any" -- Current order of array.
None -- Same as "Any"
"""
return self.filled(fill_value).tostring(order=order)
#........................
def tofile(self, fid, sep="", format="%s"):
raise NotImplementedError("Not implemented yet, sorry...")
#--------------------------------------------
# Pickling
def __getstate__(self):
"""Return the internal state of the masked array, for pickling
purposes.
"""
state = (1,
self.shape,
self.dtype,
self.flags.fnc,
self._data.tostring(),
getmaskarray(self).tostring(),
self._fill_value,
)
return state
#
def __setstate__(self, state):
"""Restore the internal state of the masked array, for
pickling purposes. ``state`` is typically the output of the
``__getstate__`` output, and is a 5-tuple:
- class name
- a tuple giving the shape of the data
- a typecode for the data
- a binary string for the data
- a binary string for the mask.
"""
(ver, shp, typ, isf, raw, msk, flv) = state
ndarray.__setstate__(self, (shp, typ, isf, raw))
self._mask.__setstate__((shp, dtype(bool), isf, msk))
self.fill_value = flv
#
def __reduce__(self):
"""Return a 3-tuple for pickling a MaskedArray.
"""
return (_mareconstruct,
(self.__class__, self._baseclass, (0,), 'b', ),
self.__getstate__())
def _mareconstruct(subtype, baseclass, baseshape, basetype,):
"""Internal function that builds a new MaskedArray from the
information stored in a pickle.
"""
_data = ndarray.__new__(baseclass, baseshape, basetype)
_mask = ndarray.__new__(ndarray, baseshape, 'b1')
return subtype.__new__(subtype, _data, mask=_mask, dtype=basetype,)
#####--------------------------------------------------------------------------
#---- --- Shortcuts ---
#####---------------------------------------------------------------------------
def isMaskedArray(x):
"Is x a masked array, that is, an instance of MaskedArray?"
return isinstance(x, MaskedArray)
isarray = isMaskedArray
isMA = isMaskedArray #backward compatibility
# We define the masked singleton as a float for higher precedence...
# Note that it can be tricky sometimes w/ type comparison
masked_singleton = MaskedArray(0, dtype=float_, mask=True)
masked = masked_singleton
masked_array = MaskedArray
def array(data, dtype=None, copy=False, order=False,
mask=nomask, fill_value=None,
keep_mask=True, hard_mask=False, shrink=True, subok=True, ndmin=0,
):
"""array(data, dtype=None, copy=False, order=False, mask=nomask,
fill_value=None, keep_mask=True, hard_mask=False, shrink=True,
subok=True, ndmin=0)
Acts as shortcut to MaskedArray, with options in a different order
for convenience. And backwards compatibility...
"""
#TODO: we should try to put 'order' somwehere
return MaskedArray(data, mask=mask, dtype=dtype, copy=copy, subok=subok,
keep_mask=keep_mask, hard_mask=hard_mask,
fill_value=fill_value, ndmin=ndmin, shrink=shrink)
array.__doc__ = masked_array.__doc__
def is_masked(x):
"""Does x have masked values?"""
m = getmask(x)
if m is nomask:
return False
elif m.any():
return True
return False
#####---------------------------------------------------------------------------
#---- --- Extrema functions ---
#####---------------------------------------------------------------------------
class _extrema_operation(object):
"Generic class for maximum/minimum functions."
def __call__(self, a, b=None):
"Executes the call behavior."
if b is None:
return self.reduce(a)
return where(self.compare(a, b), a, b)
#.........
def reduce(self, target, axis=None):
"Reduce target along the given axis."
target = narray(target, copy=False, subok=True)
m = getmask(target)
if axis is not None:
kargs = { 'axis' : axis }
else:
kargs = {}
target = target.ravel()
if not (m is nomask):
m = m.ravel()
if m is nomask:
t = self.ufunc.reduce(target, **kargs)
else:
target = target.filled(self.fill_value_func(target)).view(type(target))
t = self.ufunc.reduce(target, **kargs)
m = umath.logical_and.reduce(m, **kargs)
if hasattr(t, '_mask'):
t._mask = m
elif m:
t = masked
return t
#.........
def outer (self, a, b):
"Return the function applied to the outer product of a and b."
ma = getmask(a)
mb = getmask(b)
if ma is nomask and mb is nomask:
m = nomask
else:
ma = getmaskarray(a)
mb = getmaskarray(b)
m = logical_or.outer(ma, mb)
result = self.ufunc.outer(filled(a), filled(b))
result._mask = m
return result
#............................
class _minimum_operation(_extrema_operation):
"Object to calculate minima"
def __init__ (self):
"""minimum(a, b) or minimum(a)
In one argument case, returns the scalar minimum.
"""
self.ufunc = umath.minimum
self.afunc = amin
self.compare = less
self.fill_value_func = minimum_fill_value
#............................
class _maximum_operation(_extrema_operation):
"Object to calculate maxima"
def __init__ (self):
"""maximum(a, b) or maximum(a)
In one argument case returns the scalar maximum.
"""
self.ufunc = umath.maximum
self.afunc = amax
self.compare = greater
self.fill_value_func = maximum_fill_value
#..........................................................
def min(array, axis=None, out=None):
"""Return the minima along the given axis.
If `axis` is None, applies to the flattened array.
"""
if out is not None:
raise TypeError("Output arrays Unsupported for masked arrays")
if axis is None:
return minimum(array)
else:
return minimum.reduce(array, axis)
min.__doc__ = MaskedArray.min.__doc__
#............................
def max(obj, axis=None, out=None):
if out is not None:
raise TypeError("Output arrays Unsupported for masked arrays")
if axis is None:
return maximum(obj)
else:
return maximum.reduce(obj, axis)
max.__doc__ = MaskedArray.max.__doc__
#.............................
def ptp(obj, axis=None):
"""a.ptp(axis=None) = a.max(axis)-a.min(axis)"""
try:
return obj.max(axis)-obj.min(axis)
except AttributeError:
return max(obj, axis=axis) - min(obj, axis=axis)
ptp.__doc__ = MaskedArray.ptp.__doc__
#####---------------------------------------------------------------------------
#---- --- Definition of functions from the corresponding methods ---
#####---------------------------------------------------------------------------
class _frommethod:
"""Define functions from existing MaskedArray methods.
Parameters
----------
_methodname : string
Name of the method to transform.
"""
def __init__(self, methodname):
self._methodname = methodname
self.__doc__ = self.getdoc()
def getdoc(self):
"Return the doc of the function (from the doc of the method)."
try:
return getattr(MaskedArray, self._methodname).__doc__
except:
return getattr(numpy, self._methodname).__doc__
def __call__(self, a, *args, **params):
if isinstance(a, MaskedArray):
return getattr(a, self._methodname).__call__(*args, **params)
#FIXME ----
#As x is not a MaskedArray, we transform it to a ndarray with asarray
#... and call the corresponding method.
#Except that sometimes it doesn't work (try reshape([1,2,3,4],(2,2)))
#we end up with a "SystemError: NULL result without error in PyObject_Call"
#A dirty trick is then to call the initial numpy function...
method = getattr(narray(a, copy=False), self._methodname)
try:
return method(*args, **params)
except SystemError:
return getattr(numpy,self._methodname).__call__(a, *args, **params)
all = _frommethod('all')
anomalies = anom = _frommethod('anom')
any = _frommethod('any')
conjugate = _frommethod('conjugate')
ids = _frommethod('ids')
nonzero = _frommethod('nonzero')
diagonal = _frommethod('diagonal')
maximum = _maximum_operation()
mean = _frommethod('mean')
minimum = _minimum_operation ()
product = _frommethod('prod')
ptp = _frommethod('ptp')
ravel = _frommethod('ravel')
repeat = _frommethod('repeat')
round = _frommethod('round')
std = _frommethod('std')
sum = _frommethod('sum')
swapaxes = _frommethod('swapaxes')
take = _frommethod('take')
trace = _frommethod('trace')
var = _frommethod('var')
compress = _frommethod('compress')
#..............................................................................
def power(a, b, third=None):
"""Computes a**b elementwise.
"""
if third is not None:
raise MAError, "3-argument power not supported."
# Get the masks
ma = getmask(a)
mb = getmask(b)
m = mask_or(ma, mb)
# Get the rawdata
fa = getdata(a)
fb = getdata(b)
# Get the type of the result (so that we preserve subclasses)
if isinstance(a,MaskedArray):
basetype = type(a)
else:
basetype = MaskedArray
# Get the result and view it as a (subclass of) MaskedArray
result = umath.power(fa,fb).view(basetype)
# Find where we're in trouble w/ NaNs and Infs
invalid = numpy.logical_not(numpy.isfinite(result.view(ndarray)))
# Retrieve some extra attributes if needed
if isinstance(result,MaskedArray):
result._update_from(a)
# Add the initial mask
if m is not nomask:
if numpy.isscalar(result):
return masked
result._mask = m
# Fix the invalid parts
if invalid.any():
if not result.ndim:
return masked
result[invalid] = masked
result._data[invalid] = result.fill_value
return result
# if fb.dtype.char in typecodes["Integer"]:
# return masked_array(umath.power(fa, fb), m)
# m = mask_or(m, (fa < 0) & (fb != fb.astype(int)))
# if m is nomask:
# return masked_array(umath.power(fa, fb))
# else:
# fa = fa.copy()
# if m.all():
# fa.flat = 1
# else:
# numpy.putmask(fa,m,1)
# return masked_array(umath.power(fa, fb), m)
#..............................................................................
def argsort(a, axis=None, kind='quicksort', order=None, fill_value=None):
"Function version of the eponymous method."
if fill_value is None:
fill_value = default_fill_value(a)
d = filled(a, fill_value)
if axis is None:
return d.argsort(kind=kind, order=order)
return d.argsort(axis, kind=kind, order=order)
argsort.__doc__ = MaskedArray.argsort.__doc__
def argmin(a, axis=None, fill_value=None):
"Function version of the eponymous method."
if fill_value is None:
fill_value = default_fill_value(a)
d = filled(a, fill_value)
return d.argmin(axis=axis)
argmin.__doc__ = MaskedArray.argmin.__doc__
def argmax(a, axis=None, fill_value=None):
"Function version of the eponymous method."
if fill_value is None:
fill_value = default_fill_value(a)
try:
fill_value = - fill_value
except:
pass
d = filled(a, fill_value)
return d.argmax(axis=axis)
argmin.__doc__ = MaskedArray.argmax.__doc__
def sort(a, axis=-1, kind='quicksort', order=None, endwith=True, fill_value=None):
"Function version of the eponymous method."
a = narray(a, copy=True, subok=True)
if axis is None:
a = a.flatten()
axis = 0
if fill_value is None:
if endwith:
filler = minimum_fill_value(a)
else:
filler = maximum_fill_value(a)
else:
filler = fill_value
# return
indx = numpy.indices(a.shape).tolist()
indx[axis] = filled(a,filler).argsort(axis=axis,kind=kind,order=order)
return a[indx]
sort.__doc__ = MaskedArray.sort.__doc__
def compressed(x):
"""Return a 1-D array of all the non-masked data."""
if getmask(x) is nomask:
return numpy.asanyarray(x)
else:
return x.compressed()
def concatenate(arrays, axis=0):
"Concatenate the arrays along the given axis."
d = numpy.concatenate([getdata(a) for a in arrays], axis)
rcls = get_masked_subclass(*arrays)
data = d.view(rcls)
# Check whether one of the arrays has a non-empty mask...
for x in arrays:
if getmask(x) is not nomask:
break
else:
return data
# OK, so we have to concatenate the masks
dm = numpy.concatenate([getmaskarray(a) for a in arrays], axis)
# If we decide to keep a '_shrinkmask' option, we want to check that ...
# ... all of them are True, and then check for dm.any()
# shrink = numpy.logical_or.reduce([getattr(a,'_shrinkmask',True) for a in arrays])
# if shrink and not dm.any():
if not dm.any():
data._mask = nomask
else:
data._mask = dm.reshape(d.shape)
return data
def count(a, axis = None):
return masked_array(a, copy=False).count(axis)
count.__doc__ = MaskedArray.count.__doc__
def expand_dims(x,axis):
"""Expand the shape of the array by including a new axis before
the given one.
"""
result = n_expand_dims(x,axis)
if isinstance(x, MaskedArray):
new_shape = result.shape
result = x.view()
result.shape = new_shape
if result._mask is not nomask:
result._mask.shape = new_shape
return result
#......................................
def left_shift (a, n):
"Left shift n bits."
m = getmask(a)
if m is nomask:
d = umath.left_shift(filled(a), n)
return masked_array(d)
else:
d = umath.left_shift(filled(a, 0), n)
return masked_array(d, mask=m)
def right_shift (a, n):
"Right shift n bits."
m = getmask(a)
if m is nomask:
d = umath.right_shift(filled(a), n)
return masked_array(d)
else:
d = umath.right_shift(filled(a, 0), n)
return masked_array(d, mask=m)
#......................................
def put(a, indices, values, mode='raise'):
"""Set storage-indexed locations to corresponding values.
Values and indices are filled if necessary.
"""
# We can't use 'frommethod', the order of arguments is different
try:
return a.put(indices, values, mode=mode)
except AttributeError:
return narray(a, copy=False).put(indices, values, mode=mode)
def putmask(a, mask, values): #, mode='raise'):
"""Set a.flat[n] = values[n] for each n where mask.flat[n] is true.
If values is not the same size of a and mask then it will repeat
as necessary. This gives different behavior than
a[mask] = values.
Note: Using a masked array as values will NOT transform a ndarray in
a maskedarray.
"""
# We can't use 'frommethod', the order of arguments is different
if not isinstance(a, MaskedArray):
a = a.view(MaskedArray)
(valdata, valmask) = (getdata(values), getmask(values))
if getmask(a) is nomask:
if valmask is not nomask:
a._sharedmask = True
a.mask = numpy.zeros(a.shape, dtype=bool_)
numpy.putmask(a._mask, mask, valmask)
elif a._hardmask:
if valmask is not nomask:
m = a._mask.copy()
numpy.putmask(m, mask, valmask)
a.mask |= m
else:
if valmask is nomask:
valmask = getmaskarray(values)
numpy.putmask(a._mask, mask, valmask)
numpy.putmask(a._data, mask, valdata)
return
def transpose(a,axes=None):
"""Return a view of the array with dimensions permuted according to axes,
as a masked array.
If ``axes`` is None (default), the output view has reversed
dimensions compared to the original.
"""
#We can't use 'frommethod', as 'transpose' doesn't take keywords
try:
return a.transpose(axes)
except AttributeError:
return narray(a, copy=False).transpose(axes).view(MaskedArray)
def reshape(a, new_shape):
"""Change the shape of the array a to new_shape."""
#We can't use 'frommethod', it whine about some parameters. Dmmit.
try:
return a.reshape(new_shape)
except AttributeError:
return narray(a, copy=False).reshape(new_shape).view(MaskedArray)
def resize(x, new_shape):
"""Return a new array with the specified shape.
The total size of the original array can be any size. The new
array is filled with repeated copies of a. If a was masked, the
new array will be masked, and the new mask will be a repetition of
the old one.
"""
# We can't use _frommethods here, as N.resize is notoriously whiny.
m = getmask(x)
if m is not nomask:
m = numpy.resize(m, new_shape)
result = numpy.resize(x, new_shape).view(get_masked_subclass(x))
if result.ndim:
result._mask = m
return result
#................................................
def rank(obj):
"maskedarray version of the numpy function."
return fromnumeric.rank(getdata(obj))
rank.__doc__ = numpy.rank.__doc__
#
def shape(obj):
"maskedarray version of the numpy function."
return fromnumeric.shape(getdata(obj))
shape.__doc__ = numpy.shape.__doc__
#
def size(obj, axis=None):
"maskedarray version of the numpy function."
return fromnumeric.size(getdata(obj), axis)
size.__doc__ = numpy.size.__doc__
#................................................
#####--------------------------------------------------------------------------
#---- --- Extra functions ---
#####--------------------------------------------------------------------------
def where (condition, x=None, y=None):
"""where(condition | x, y)
Returns a (subclass of) masked array, shaped like condition, where
the elements are x when condition is True, and y otherwise. If
neither x nor y are given, returns a tuple of indices where
condition is True (a la condition.nonzero()).
Parameters
----------
condition : {var}
The condition to meet. Must be convertible to an integer
array.
x : {var}, optional
Values of the output when the condition is met
y : {var}, optional
Values of the output when the condition is not met.
"""
if x is None and y is None:
return filled(condition, 0).nonzero()
elif x is None or y is None:
raise ValueError, "Either both or neither x and y should be given."
# Get the condition ...............
fc = filled(condition, 0).astype(bool_)
notfc = numpy.logical_not(fc)
# Get the data ......................................
xv = getdata(x)
yv = getdata(y)
if x is masked:
ndtype = yv.dtype
xm = numpy.ones(fc.shape, dtype=MaskType)
elif y is masked:
ndtype = xv.dtype
ym = numpy.ones(fc.shape, dtype=MaskType)
else:
ndtype = numpy.max([xv.dtype, yv.dtype])
xm = getmask(x)
d = numpy.empty(fc.shape, dtype=ndtype).view(MaskedArray)
numpy.putmask(d._data, fc, xv.astype(ndtype))
numpy.putmask(d._data, notfc, yv.astype(ndtype))
d._mask = numpy.zeros(fc.shape, dtype=MaskType)
numpy.putmask(d._mask, fc, getmask(x))
numpy.putmask(d._mask, notfc, getmask(y))
d._mask |= getmaskarray(condition)
if not d._mask.any():
d._mask = nomask
return d
def choose (indices, t, out=None, mode='raise'):
"Return array shaped like indices with elements chosen from t"
#TODO: implement options `out` and `mode`, if possible.
def fmask (x):
"Returns the filled array, or True if masked."
if x is masked:
return 1
return filled(x)
def nmask (x):
"Returns the mask, True if ``masked``, False if ``nomask``."
if x is masked:
return 1
m = getmask(x)
if m is nomask:
return 0
return m
c = filled(indices, 0)
masks = [nmask(x) for x in t]
a = [fmask(x) for x in t]
d = numpy.choose(c, a)
m = numpy.choose(c, masks)
m = make_mask(mask_or(m, getmask(indices)), copy=0, shrink=True)
return masked_array(d, mask=m)
def round_(a, decimals=0, out=None):
"""Return a copy of a, rounded to 'decimals' places.
When 'decimals' is negative, it specifies the number of positions
to the left of the decimal point. The real and imaginary parts of
complex numbers are rounded separately. Nothing is done if the
array is not of float type and 'decimals' is greater than or equal
to 0.
Parameters
----------
decimals : int
Number of decimals to round to. May be negative.
out : array_like
Existing array to use for output.
If not given, returns a default copy of a.
Notes
-----
If out is given and does not have a mask attribute, the mask of a
is lost!
"""
if out is None:
return numpy.round_(a, decimals, out)
else:
numpy.round_(getdata(a), decimals, out)
if hasattr(out, '_mask'):
out._mask = getmask(a)
return out
def arange(stop, start=None, step=1, dtype=None):
"maskedarray version of the numpy function."
return numpy.arange(stop, start, step, dtype).view(MaskedArray)
arange.__doc__ = numpy.arange.__doc__
def inner(a, b):
"maskedarray version of the numpy function."
fa = filled(a, 0)
fb = filled(b, 0)
if len(fa.shape) == 0:
fa.shape = (1,)
if len(fb.shape) == 0:
fb.shape = (1,)
return numpy.inner(fa, fb).view(MaskedArray)
inner.__doc__ = numpy.inner.__doc__
inner.__doc__ += doc_note("Masked values are replaced by 0.")
innerproduct = inner
def outer(a, b):
"maskedarray version of the numpy function."
fa = filled(a, 0).ravel()
fb = filled(b, 0).ravel()
d = numeric.outer(fa, fb)
ma = getmask(a)
mb = getmask(b)
if ma is nomask and mb is nomask:
return masked_array(d)
ma = getmaskarray(a)
mb = getmaskarray(b)
m = make_mask(1-numeric.outer(1-ma, 1-mb), copy=0)
return masked_array(d, mask=m)
outer.__doc__ = numpy.outer.__doc__
outer.__doc__ += doc_note("Masked values are replaced by 0.")
outerproduct = outer
def allequal (a, b, fill_value=True):
"""Return True if all entries of a and b are equal, using
fill_value as a truth value where either or both are masked.
"""
m = mask_or(getmask(a), getmask(b))
if m is nomask:
x = getdata(a)
y = getdata(b)
d = umath.equal(x, y)
return d.all()
elif fill_value:
x = getdata(a)
y = getdata(b)
d = umath.equal(x, y)
dm = array(d, mask=m, copy=False)
return dm.filled(True).all(None)
else:
return False
def allclose (a, b, fill_value=True, rtol=1.e-5, atol=1.e-8):
""" Return True if all elements of a and b are equal subject to
given tolerances.
If fill_value is True, masked values are considered equal.
If fill_value is False, masked values considered unequal.
The relative error rtol should be positive and << 1.0
The absolute error atol comes into play for those elements of b
that are very small or zero; it says how small `a` must be also.
"""
m = mask_or(getmask(a), getmask(b))
d1 = getdata(a)
d2 = getdata(b)
x = filled(array(d1, copy=0, mask=m), fill_value).astype(float)
y = filled(array(d2, copy=0, mask=m), 1).astype(float)
d = umath.less_equal(umath.absolute(x-y), atol + rtol * umath.absolute(y))
return fromnumeric.alltrue(fromnumeric.ravel(d))
#..............................................................................
def asarray(a, dtype=None):
"""asarray(data, dtype) = array(data, dtype, copy=0, subok=0)
Return a as a MaskedArray object of the given dtype.
If dtype is not given or None, is is set to the dtype of a.
No copy is performed if a is already an array.
Subclasses are converted to the base class MaskedArray.
"""
return masked_array(a, dtype=dtype, copy=False, keep_mask=True, subok=False)
def asanyarray(a, dtype=None):
"""asanyarray(data, dtype) = array(data, dtype, copy=0, subok=1)
Return a as an masked array.
If dtype is not given or None, is is set to the dtype of a.
No copy is performed if a is already an array.
Subclasses are conserved.
"""
return masked_array(a, dtype=dtype, copy=False, keep_mask=True, subok=True)
def empty(new_shape, dtype=float):
"maskedarray version of the numpy function."
return numpy.empty(new_shape, dtype).view(MaskedArray)
empty.__doc__ = numpy.empty.__doc__
def empty_like(a):
"maskedarray version of the numpy function."
return numpy.empty_like(a).view(MaskedArray)
empty_like.__doc__ = numpy.empty_like.__doc__
def ones(new_shape, dtype=float):
"maskedarray version of the numpy function."
return numpy.ones(new_shape, dtype).view(MaskedArray)
ones.__doc__ = numpy.ones.__doc__
def zeros(new_shape, dtype=float):
"maskedarray version of the numpy function."
return numpy.zeros(new_shape, dtype).view(MaskedArray)
zeros.__doc__ = numpy.zeros.__doc__
#####--------------------------------------------------------------------------
#---- --- Pickling ---
#####--------------------------------------------------------------------------
def dump(a,F):
"""Pickle the MaskedArray `a` to the file `F`. `F` can either be
the handle of an exiting file, or a string representing a file
name.
"""
if not hasattr(F,'readline'):
F = open(F,'w')
return cPickle.dump(a,F)
def dumps(a):
"""Return a string corresponding to the pickling of the
MaskedArray.
"""
return cPickle.dumps(a)
def load(F):
"""Wrapper around ``cPickle.load`` which accepts either a
file-like object or a filename.
"""
if not hasattr(F, 'readline'):
F = open(F,'r')
return cPickle.load(F)
def loads(strg):
"Load a pickle from thecurrentstring.Convertfunctionsnumpytonumpy.ma. import
Parameters
----------
_methodname : string
Name of the method to transform.
"""
__doc__ = None
def __init__(self, funcname):
self._func = getattr(numpy, funcname)
self.__doc__ = self.getdoc()
def getdoc(self):
"Return the doc of the function (from the doc of the method)."
return self._func.__doc__
def __call__(self, a, *args, **params):
return self._func.__call__(a, *args, **params).view(MaskedArray)
frombuffer = _convert2ma('frombuffer')
fromfunction = _convert2ma('fromfunction')
identity = _convert2ma('identity')
indices = numpy.indices
clip = numpy.clip
###############################################################################
|