diff options
Diffstat (limited to 'numpy/oldnumeric/ma.py')
| -rw-r--r-- | numpy/oldnumeric/ma.py | 2269 |
1 files changed, 2262 insertions, 7 deletions
diff --git a/numpy/oldnumeric/ma.py b/numpy/oldnumeric/ma.py index 81555570a..af83cb0a4 100644 --- a/numpy/oldnumeric/ma.py +++ b/numpy/oldnumeric/ma.py @@ -1,14 +1,2269 @@ -# Incompatibility in that getmask and a.mask returns nomask -# instead of None +"""MA: a facility for dealing with missing observations +MA is generally used as a numpy.array look-alike. +by Paul F. Dubois. -from numpy.ma import * -import numpy.ma as nca +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. + +""" +import types, sys + +import umath +import fromnumeric +from numeric import newaxis, ndarray, inf +from fromnumeric import amax, amin +from numerictypes import bool_, typecodes +import numeric +import warnings + +# Ufunc domain lookup for __array_wrap__ +ufunc_domain = {} +# Ufunc fills lookup for __array__ +ufunc_fills = {} + +MaskType = bool_ +nomask = MaskType(0) +divide_tolerance = 1.e-35 + +class MAError (Exception): + def __init__ (self, args=None): + "Create an exception" + + # The .args attribute must be a tuple. + if not isinstance(args, tuple): + args = (args,) + self.args = args + def __str__(self): + "Calculate the string representation" + return str(self.args[0]) + __repr__ = __str__ + +class _MaskedPrintOption: + "One instance of this class, masked_print_option, is created." + def __init__ (self, display): + "Create the masked print option object." + self.set_display(display) + self._enabled = 1 + + def display (self): + "Show what prints for masked values." + return self._display + + def set_display (self, s): + "set_display(s) sets what prints for masked values." + self._display = s + + def enabled (self): + "Is the use of the display value enabled?" + return self._enabled + + def enable(self, flag=1): + "Set the enabling flag to flag." + self._enabled = flag + + 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('--') + +# Use single element arrays or scalars. +default_real_fill_value = 1.e20 +default_complex_fill_value = 1.e20 + 0.0j +default_character_fill_value = '-' +default_integer_fill_value = 999999 +default_object_fill_value = '?' + +def default_fill_value (obj): + "Function to calculate default fill value for an object." + if isinstance(obj, types.FloatType): + return default_real_fill_value + elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType): + return default_integer_fill_value + elif isinstance(obj, types.StringType): + return default_character_fill_value + elif isinstance(obj, types.ComplexType): + return default_complex_fill_value + elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray): + x = obj.dtype.char + if x in typecodes['Float']: + return default_real_fill_value + if x in typecodes['Integer']: + return default_integer_fill_value + if x in typecodes['Complex']: + return default_complex_fill_value + if x in typecodes['Character']: + return default_character_fill_value + if x in typecodes['UnsignedInteger']: + return umath.absolute(default_integer_fill_value) + return default_object_fill_value + else: + return default_object_fill_value + +def minimum_fill_value (obj): + "Function to calculate default fill value suitable for taking minima." + if isinstance(obj, types.FloatType): + return numeric.inf + elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType): + return sys.maxint + elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray): + x = obj.dtype.char + if x in typecodes['Float']: + return numeric.inf + if x in typecodes['Integer']: + return sys.maxint + if x in typecodes['UnsignedInteger']: + return sys.maxint + else: + raise TypeError, 'Unsuitable type for calculating minimum.' + +def maximum_fill_value (obj): + "Function to calculate default fill value suitable for taking maxima." + if isinstance(obj, types.FloatType): + return -inf + elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType): + return -sys.maxint + elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray): + x = obj.dtype.char + if x in typecodes['Float']: + return -inf + if x in typecodes['Integer']: + return -sys.maxint + if x in typecodes['UnsignedInteger']: + return 0 + else: + raise TypeError, 'Unsuitable type for calculating maximum.' + +def set_fill_value (a, fill_value): + "Set fill value of a if it is a masked array." + if isMaskedArray(a): + a.set_fill_value (fill_value) + +def getmask (a): + """Mask of values in a; could be nomask. + Returns nomask if a is not a masked array. + To get an array for sure use getmaskarray.""" + if isinstance(a, MaskedArray): + return a.raw_mask() + else: + return nomask + +def getmaskarray (a): + """Mask of values in a; an array of zeros if mask is nomask + or not a masked array, and is a byte-sized integer. + Do not try to add up entries, for example. + """ + m = getmask(a) + if m is nomask: + return make_mask_none(shape(a)) + else: + return m + +def is_mask (m): + """Is m 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=0, flag=0): + """make_mask(m, copy=0, flag=0) + return m as a mask, creating a copy if necessary or requested. + Can accept any sequence of integers or nomask. Does not check + that contents must be 0s and 1s. + if flag, return nomask if m contains no true elements. + """ + if m is nomask: + return nomask + elif isinstance(m, ndarray): + if m.dtype.type is MaskType: + if copy: + result = numeric.array(m, dtype=MaskType, copy=copy) + else: + result = m + else: + result = m.astype(MaskType) + else: + result = filled(m, True).astype(MaskType) + + if flag and not fromnumeric.sometrue(fromnumeric.ravel(result)): + return nomask + else: + return result + +def make_mask_none (s): + "Return a mask of all zeros of shape s." + result = numeric.zeros(s, dtype=MaskType) + result.shape = s + return result + +def mask_or (m1, m2): + """Logical or of the mask candidates m1 and m2, treating nomask as false. + Result may equal m1 or m2 if the other is nomask. + """ + if m1 is nomask: return make_mask(m2) + if m2 is nomask: return make_mask(m1) + if m1 is m2 and is_mask(m1): return m1 + return make_mask(umath.logical_or(m1, m2)) + +def filled (a, value = None): + """a as a contiguous numeric array with any masked areas replaced by value + if value is None or the special element "masked", get_fill_value(a) + is used instead. + + If a is already a contiguous numeric array, a itself is returned. + + filled(a) can be used to be sure that the result is numeric when + passing an object a to other software ignorant of MA, in particular to + numeric itself. + """ + if isinstance(a, MaskedArray): + return a.filled(value) + elif isinstance(a, ndarray) and a.flags['CONTIGUOUS']: + return a + elif isinstance(a, types.DictType): + return numeric.array(a, 'O') + else: + return numeric.array(a) + +def get_fill_value (a): + """ + The fill value of a, if it has one; otherwise, the default fill value + for that type. + """ + if isMaskedArray(a): + result = a.fill_value() + else: + result = default_fill_value(a) + return result + +def common_fill_value (a, b): + "The common fill_value of a and b, if there is one, or None" + t1 = get_fill_value(a) + t2 = get_fill_value(b) + if t1 == t2: return t1 + return None + +# Domain functions return 1 where the argument(s) are not in the domain. +class domain_check_interval: + "domain_check_interval(a,b)(x) = true where x < a or y > b" + def __init__(self, y1, y2): + "domain_check_interval(a,b)(x) = true where x < a or y > b" + self.y1 = y1 + self.y2 = y2 + + def __call__ (self, x): + "Execute the call behavior." + return umath.logical_or(umath.greater (x, self.y2), + umath.less(x, self.y1) + ) + +class domain_tan: + "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): + "Execute the call behavior." + return umath.less(umath.absolute(umath.cos(x)), self.eps) + +class domain_greater: + "domain_greater(v)(x) = true where x <= v" + def __init__(self, critical_value): + "domain_greater(v)(x) = true where x <= v" + self.critical_value = critical_value + + def __call__ (self, x): + "Execute the call behavior." + return umath.less_equal (x, self.critical_value) + +class domain_greater_equal: + "domain_greater_equal(v)(x) = true where x < v" + def __init__(self, critical_value): + "domain_greater_equal(v)(x) = true where x < v" + self.critical_value = critical_value + + def __call__ (self, x): + "Execute the call behavior." + return umath.less (x, self.critical_value) + +class masked_unary_operation: + def __init__ (self, aufunc, fill=0, domain=None): + """ masked_unary_operation(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 = aufunc + self.fill = fill + self.domain = domain + self.__doc__ = getattr(aufunc, "__doc__", str(aufunc)) + self.__name__ = getattr(aufunc, "__name__", str(aufunc)) + ufunc_domain[aufunc] = domain + ufunc_fills[aufunc] = fill, + + def __call__ (self, a, *args, **kwargs): + "Execute the call behavior." +# numeric tries to return scalars rather than arrays when given scalars. + m = getmask(a) + d1 = filled(a, self.fill) + if self.domain is not None: + m = mask_or(m, self.domain(d1)) + result = self.f(d1, *args, **kwargs) + return masked_array(result, m) + + def __str__ (self): + return "Masked version of " + str(self.f) + + +class domain_safe_divide: + def __init__ (self, tolerance=divide_tolerance): + self.tolerance = tolerance + def __call__ (self, a, b): + return umath.absolute(a) * self.tolerance >= umath.absolute(b) + +class domained_binary_operation: + """Binary operations that have a domain, like divide. These are complicated + so they are a separate class. They have no reduce, outer or accumulate. + """ + def __init__ (self, abfunc, domain, fillx=0, filly=0): + """abfunc(fillx, filly) must be defined. + abfunc(x, filly) = x for all x to enable reduce. + """ + self.f = abfunc + self.domain = domain + self.fillx = fillx + self.filly = filly + self.__doc__ = getattr(abfunc, "__doc__", str(abfunc)) + self.__name__ = getattr(abfunc, "__name__", str(abfunc)) + ufunc_domain[abfunc] = domain + ufunc_fills[abfunc] = fillx, filly + + def __call__(self, a, b): + "Execute the call behavior." + ma = getmask(a) + mb = getmask(b) + d1 = filled(a, self.fillx) + d2 = filled(b, self.filly) + t = self.domain(d1, d2) + + if fromnumeric.sometrue(t, None): + d2 = where(t, self.filly, d2) + mb = mask_or(mb, t) + m = mask_or(ma, mb) + result = self.f(d1, d2) + return masked_array(result, m) + + def __str__ (self): + return "Masked version of " + str(self.f) + +class masked_binary_operation: + def __init__ (self, abfunc, fillx=0, filly=0): + """abfunc(fillx, filly) must be defined. + abfunc(x, filly) = x for all x to enable reduce. + """ + self.f = abfunc + self.fillx = fillx + self.filly = filly + self.__doc__ = getattr(abfunc, "__doc__", str(abfunc)) + ufunc_domain[abfunc] = None + ufunc_fills[abfunc] = fillx, filly + + def __call__ (self, a, b, *args, **kwargs): + "Execute the call behavior." + m = mask_or(getmask(a), getmask(b)) + d1 = filled(a, self.fillx) + d2 = filled(b, self.filly) + result = self.f(d1, d2, *args, **kwargs) + if isinstance(result, ndarray) \ + and m.ndim != 0 \ + and m.shape != result.shape: + m = mask_or(getmaskarray(a), getmaskarray(b)) + return masked_array(result, m) + + def reduce (self, target, axis=0, dtype=None): + """Reduce target along the given axis with this function.""" + 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: + t = self.f.reduce(t, axis) + else: + t = masked_array (t, m) + # XXX: "or t.dtype" below is a workaround for what appears + # XXX: to be a bug in reduce. + t = self.f.reduce(filled(t, self.filly), axis, + dtype=dtype or t.dtype) + m = umath.logical_and.reduce(m, axis) + if isinstance(t, ndarray): + return masked_array(t, m, get_fill_value(target)) + elif m: + return masked + else: + 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) + d = self.f.outer(filled(a, self.fillx), filled(b, self.filly)) + return masked_array(d, m) + + def accumulate (self, target, axis=0): + """Accumulate target along axis after filling with y fill value.""" + t = filled(target, self.filly) + return masked_array (self.f.accumulate (t, axis)) + def __str__ (self): + return "Masked version of " + str(self.f) + +sqrt = masked_unary_operation(umath.sqrt, 0.0, domain_greater_equal(0.0)) +log = masked_unary_operation(umath.log, 1.0, domain_greater(0.0)) +log10 = masked_unary_operation(umath.log10, 1.0, domain_greater(0.0)) +exp = masked_unary_operation(umath.exp) +conjugate = masked_unary_operation(umath.conjugate) +sin = masked_unary_operation(umath.sin) +cos = masked_unary_operation(umath.cos) +tan = masked_unary_operation(umath.tan, 0.0, domain_tan(1.e-35)) +arcsin = masked_unary_operation(umath.arcsin, 0.0, domain_check_interval(-1.0, 1.0)) +arccos = masked_unary_operation(umath.arccos, 0.0, domain_check_interval(-1.0, 1.0)) +arctan = masked_unary_operation(umath.arctan) +# Missing from numeric +arcsinh = masked_unary_operation(umath.arcsinh) +arccosh = masked_unary_operation(umath.arccosh, 1.0, domain_greater_equal(1.0)) +arctanh = masked_unary_operation(umath.arctanh, 0.0, domain_check_interval(-1.0+1e-15, 1.0-1e-15)) +sinh = masked_unary_operation(umath.sinh) +cosh = masked_unary_operation(umath.cosh) +tanh = masked_unary_operation(umath.tanh) +absolute = masked_unary_operation(umath.absolute) +fabs = masked_unary_operation(umath.fabs) +negative = masked_unary_operation(umath.negative) + +def nonzero(a): + """returns the indices of the elements of a which are not zero + and not masked + """ + return numeric.asarray(filled(a, 0).nonzero()) + +around = masked_unary_operation(fromnumeric.round_) +floor = masked_unary_operation(umath.floor) +ceil = masked_unary_operation(umath.ceil) +logical_not = masked_unary_operation(umath.logical_not) + +add = masked_binary_operation(umath.add) +subtract = masked_binary_operation(umath.subtract) +subtract.reduce = None +multiply = masked_binary_operation(umath.multiply, 1, 1) +divide = domained_binary_operation(umath.divide, domain_safe_divide(), 0, 1) +true_divide = domained_binary_operation(umath.true_divide, domain_safe_divide(), 0, 1) +floor_divide = domained_binary_operation(umath.floor_divide, domain_safe_divide(), 0, 1) +remainder = domained_binary_operation(umath.remainder, domain_safe_divide(), 0, 1) +fmod = domained_binary_operation(umath.fmod, domain_safe_divide(), 0, 1) +hypot = masked_binary_operation(umath.hypot) +arctan2 = masked_binary_operation(umath.arctan2, 0.0, 1.0) +arctan2.reduce = None +equal = masked_binary_operation(umath.equal) +equal.reduce = None +not_equal = masked_binary_operation(umath.not_equal) +not_equal.reduce = None +less_equal = masked_binary_operation(umath.less_equal) +less_equal.reduce = None +greater_equal = masked_binary_operation(umath.greater_equal) +greater_equal.reduce = None +less = masked_binary_operation(umath.less) +less.reduce = None +greater = masked_binary_operation(umath.greater) +greater.reduce = None +logical_and = masked_binary_operation(umath.logical_and) +alltrue = masked_binary_operation(umath.logical_and, 1, 1).reduce +logical_or = masked_binary_operation(umath.logical_or) +sometrue = logical_or.reduce +logical_xor = masked_binary_operation(umath.logical_xor) +bitwise_and = masked_binary_operation(umath.bitwise_and) +bitwise_or = masked_binary_operation(umath.bitwise_or) +bitwise_xor = masked_binary_operation(umath.bitwise_xor) + +def rank (object): + return fromnumeric.rank(filled(object)) + +def shape (object): + return fromnumeric.shape(filled(object)) + +def size (object, axis=None): + return fromnumeric.size(filled(object), axis) + +class MaskedArray (object): + """Arrays with possibly masked values. + Masked values of 1 exclude the corresponding element from + any computation. + + Construction: + x = array(data, dtype=None, copy=True, order=False, + mask = nomask, fill_value=None) + + If copy=False, every effort is made not to copy the data: + If data is a MaskedArray, and argument mask=nomask, + then the candidate data is data.data and the + mask used is data.mask. If data is a numeric array, + it is used as the candidate raw data. + If dtype is not None and + is != data.dtype.char then a data copy is required. + Otherwise, the candidate is used. + + If a data copy is required, raw data stored is the result of: + numeric.array(data, dtype=dtype.char, copy=copy) + + If mask is nomask there are no masked values. Otherwise mask must + be convertible to an array of booleans with the same shape as x. + + fill_value is used to fill in masked values when necessary, + such as when printing and in method/function filled(). + The fill_value is not used for computation within this module. + """ + __array_priority__ = 10.1 + def __init__(self, data, dtype=None, copy=True, order=False, + mask=nomask, fill_value=None): + """array(data, dtype=None, copy=True, order=False, mask=nomask, fill_value=None) + If data already a numeric array, its dtype becomes the default value of dtype. + """ + if dtype is None: + tc = None + else: + tc = numeric.dtype(dtype) + need_data_copied = copy + if isinstance(data, MaskedArray): + c = data.data + if tc is None: + tc = c.dtype + elif tc != c.dtype: + need_data_copied = True + if mask is nomask: + mask = data.mask + elif mask is not nomask: #attempting to change the mask + need_data_copied = True + + elif isinstance(data, ndarray): + c = data + if tc is None: + tc = c.dtype + elif tc != c.dtype: + need_data_copied = True + else: + need_data_copied = False #because I'll do it now + c = numeric.array(data, dtype=tc, copy=True, order=order) + tc = c.dtype + + if need_data_copied: + if tc == c.dtype: + self._data = numeric.array(c, dtype=tc, copy=True, order=order) + else: + self._data = c.astype(tc) + else: + self._data = c + + if mask is nomask: + self._mask = nomask + self._shared_mask = 0 + else: + self._mask = make_mask (mask) + if self._mask is nomask: + self._shared_mask = 0 + else: + self._shared_mask = (self._mask is mask) + nm = size(self._mask) + nd = size(self._data) + if nm != nd: + if nm == 1: + self._mask = fromnumeric.resize(self._mask, self._data.shape) + self._shared_mask = 0 + elif nd == 1: + self._data = fromnumeric.resize(self._data, self._mask.shape) + self._data.shape = self._mask.shape + else: + raise MAError, "Mask and data not compatible." + elif nm == 1 and shape(self._mask) != shape(self._data): + self.unshare_mask() + self._mask.shape = self._data.shape + + self.set_fill_value(fill_value) + + def __array__ (self, t=None, context=None): + "Special hook for numeric. Converts to numeric if possible." + if self._mask is not nomask: + if fromnumeric.ravel(self._mask).any(): + if context is None: + warnings.warn("Cannot automatically convert masked array to "\ + "numeric because data\n is masked in one or "\ + "more locations."); + return self._data + #raise MAError, \ + # """Cannot automatically convert masked array to numeric because data + # is masked in one or more locations. + # """ + else: + func, args, i = context + fills = ufunc_fills.get(func) + if fills is None: + raise MAError, "%s not known to ma" % func + return self.filled(fills[i]) + else: # Mask is all false + # Optimize to avoid future invocations of this section. + self._mask = nomask + self._shared_mask = 0 + if t: + return self._data.astype(t) + else: + return self._data + + def __array_wrap__ (self, array, context=None): + """Special hook for ufuncs. + + Wraps the numpy array and sets the mask according to + context. + """ + if context is None: + return MaskedArray(array, copy=False, mask=nomask) + func, args = context[:2] + domain = ufunc_domain[func] + m = reduce(mask_or, [getmask(a) for a in args]) + if domain is not None: + m = mask_or(m, domain(*[getattr(a, '_data', a) + for a in args])) + if m is not nomask: + try: + shape = array.shape + except AttributeError: + pass + else: + if m.shape != shape: + m = reduce(mask_or, [getmaskarray(a) for a in args]) + + return MaskedArray(array, copy=False, mask=m) + + def _get_shape(self): + "Return the current shape." + return self._data.shape + + def _set_shape (self, newshape): + "Set the array's shape." + self._data.shape = newshape + if self._mask is not nomask: + self._mask = self._mask.copy() + self._mask.shape = newshape + + def _get_flat(self): + """Calculate the flat value. + """ + if self._mask is nomask: + return masked_array(self._data.ravel(), mask=nomask, + fill_value = self.fill_value()) + else: + return masked_array(self._data.ravel(), + mask=self._mask.ravel(), + fill_value = self.fill_value()) + + def _set_flat (self, value): + "x.flat = value" + y = self.ravel() + y[:] = value + + def _get_real(self): + "Get the real part of a complex array." + if self._mask is nomask: + return masked_array(self._data.real, mask=nomask, + fill_value = self.fill_value()) + else: + return masked_array(self._data.real, mask=self._mask, + fill_value = self.fill_value()) + + def _set_real (self, value): + "x.real = value" + y = self.real + y[...] = value + + def _get_imaginary(self): + "Get the imaginary part of a complex array." + if self._mask is nomask: + return masked_array(self._data.imag, mask=nomask, + fill_value = self.fill_value()) + else: + return masked_array(self._data.imag, mask=self._mask, + fill_value = self.fill_value()) + + def _set_imaginary (self, value): + "x.imaginary = value" + y = self.imaginary + y[...] = value + + def __str__(self): + """Calculate the str representation, using masked for fill if + it is enabled. Otherwise fill with fill value. + """ + if masked_print_option.enabled(): + f = masked_print_option + # XXX: Without the following special case masked + # XXX: would print as "[--]", not "--". Can we avoid + # XXX: checks for masked by choosing a different value + # XXX: for the masked singleton? 2005-01-05 -- sasha + if self is masked: + return str(f) + m = self._mask + if m is not nomask and m.shape == () and m: + return str(f) + # convert to object array to make filled work + self = self.astype(object) + else: + f = self.fill_value() + res = self.filled(f) + return str(res) + + def __repr__(self): + """Calculate the repr representation, using masked for fill if + it is enabled. Otherwise fill with fill value. + """ + with_mask = """\ +array(data = + %(data)s, + mask = + %(mask)s, + fill_value=%(fill)s) +""" + with_mask1 = """\ +array(data = %(data)s, + mask = %(mask)s, + fill_value=%(fill)s) +""" + without_mask = """array( + %(data)s)""" + without_mask1 = """array(%(data)s)""" + + n = len(self.shape) + if self._mask is nomask: + if n <= 1: + return without_mask1 % {'data':str(self.filled())} + return without_mask % {'data':str(self.filled())} + else: + if n <= 1: + return with_mask % { + 'data': str(self.filled()), + 'mask': str(self._mask), + 'fill': str(self.fill_value()) + } + return with_mask % { + 'data': str(self.filled()), + 'mask': str(self._mask), + 'fill': str(self.fill_value()) + } + without_mask1 = """array(%(data)s)""" + if self._mask is nomask: + return without_mask % {'data':str(self.filled())} + else: + return with_mask % { + 'data': str(self.filled()), + 'mask': str(self._mask), + 'fill': str(self.fill_value()) + } + + def __float__(self): + "Convert self to float." + self.unmask() + if self._mask is not nomask: + raise MAError, 'Cannot convert masked element to a Python float.' + return float(self.data.item()) + + def __int__(self): + "Convert self to int." + self.unmask() + if self._mask is not nomask: + raise MAError, 'Cannot convert masked element to a Python int.' + return int(self.data.item()) + + def __getitem__(self, i): + "Get item described by i. Not a copy as in previous versions." + self.unshare_mask() + m = self._mask + dout = self._data[i] + if m is nomask: + try: + if dout.size == 1: + return dout + else: + return masked_array(dout, fill_value=self._fill_value) + except AttributeError: + return dout + mi = m[i] + if mi.size == 1: + if mi: + return masked + else: + return dout + else: + return masked_array(dout, mi, fill_value=self._fill_value) + +# -------- +# setitem and setslice notes +# note that if value is masked, it means to mask those locations. +# setting a value changes the mask to match the value in those locations. + + def __setitem__(self, index, value): + "Set item described by index. If value is masked, mask those locations." + d = self._data + if self is masked: + raise MAError, 'Cannot alter masked elements.' + if value is masked: + if self._mask is nomask: + self._mask = make_mask_none(d.shape) + self._shared_mask = False + else: + self.unshare_mask() + self._mask[index] = True + return + m = getmask(value) + value = filled(value).astype(d.dtype) + d[index] = value + if m is nomask: + if self._mask is not nomask: + self.unshare_mask() + self._mask[index] = False + else: + if self._mask is nomask: + self._mask = make_mask_none(d.shape) + self._shared_mask = True + else: + self.unshare_mask() + self._mask[index] = m + + def __nonzero__(self): + """returns true if any element is non-zero or masked + + """ + # XXX: This changes bool conversion logic from MA. + # XXX: In MA bool(a) == len(a) != 0, but in numpy + # XXX: scalars do not have len + m = self._mask + d = self._data + return bool(m is not nomask and m.any() + or d is not nomask and d.any()) + + def __len__ (self): + """Return length of first dimension. This is weird but Python's + slicing behavior depends on it.""" + return len(self._data) + + def __and__(self, other): + "Return bitwise_and" + return bitwise_and(self, other) + + def __or__(self, other): + "Return bitwise_or" + return bitwise_or(self, other) + + def __xor__(self, other): + "Return bitwise_xor" + return bitwise_xor(self, other) + + __rand__ = __and__ + __ror__ = __or__ + __rxor__ = __xor__ + + def __abs__(self): + "Return absolute(self)" + return absolute(self) + + def __neg__(self): + "Return negative(self)" + return negative(self) + + def __pos__(self): + "Return array(self)" + return array(self) + + def __add__(self, other): + "Return add(self, other)" + return add(self, other) + + __radd__ = __add__ + + def __mod__ (self, other): + "Return remainder(self, other)" + return remainder(self, other) + + def __rmod__ (self, other): + "Return remainder(other, self)" + return remainder(other, self) + + def __lshift__ (self, n): + return left_shift(self, n) + + def __rshift__ (self, n): + return right_shift(self, n) + + def __sub__(self, other): + "Return subtract(self, other)" + return subtract(self, other) + + def __rsub__(self, other): + "Return subtract(other, self)" + return subtract(other, self) + + def __mul__(self, other): + "Return multiply(self, other)" + return multiply(self, other) + + __rmul__ = __mul__ + + def __div__(self, other): + "Return divide(self, other)" + return divide(self, other) + + def __rdiv__(self, other): + "Return divide(other, self)" + return divide(other, self) + + def __truediv__(self, other): + "Return divide(self, other)" + return true_divide(self, other) + + def __rtruediv__(self, other): + "Return divide(other, self)" + return true_divide(other, self) + + def __floordiv__(self, other): + "Return divide(self, other)" + return floor_divide(self, other) + + def __rfloordiv__(self, other): + "Return divide(other, self)" + return floor_divide(other, self) + + def __pow__(self, other, third=None): + "Return power(self, other, third)" + return power(self, other, third) + + def __sqrt__(self): + "Return sqrt(self)" + return sqrt(self) + + def __iadd__(self, other): + "Add other to self in place." + t = self._data.dtype.char + f = filled(other, 0) + t1 = f.dtype.char + if t == t1: + pass + elif t in typecodes['Integer']: + if t1 in typecodes['Integer']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Float']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Complex']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + elif t1 in typecodes['Complex']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + else: + raise TypeError, 'Incorrect type for in-place operation.' + + if self._mask is nomask: + self._data += f + m = getmask(other) + self._mask = m + self._shared_mask = m is not nomask + else: + result = add(self, masked_array(f, mask=getmask(other))) + self._data = result.data + self._mask = result.mask + self._shared_mask = 1 + return self + + def __imul__(self, other): + "Add other to self in place." + t = self._data.dtype.char + f = filled(other, 0) + t1 = f.dtype.char + if t == t1: + pass + elif t in typecodes['Integer']: + if t1 in typecodes['Integer']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Float']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Complex']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + elif t1 in typecodes['Complex']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + else: + raise TypeError, 'Incorrect type for in-place operation.' + + if self._mask is nomask: + self._data *= f + m = getmask(other) + self._mask = m + self._shared_mask = m is not nomask + else: + result = multiply(self, masked_array(f, mask=getmask(other))) + self._data = result.data + self._mask = result.mask + self._shared_mask = 1 + return self + + def __isub__(self, other): + "Subtract other from self in place." + t = self._data.dtype.char + f = filled(other, 0) + t1 = f.dtype.char + if t == t1: + pass + elif t in typecodes['Integer']: + if t1 in typecodes['Integer']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Float']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Complex']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + elif t1 in typecodes['Complex']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + else: + raise TypeError, 'Incorrect type for in-place operation.' + + if self._mask is nomask: + self._data -= f + m = getmask(other) + self._mask = m + self._shared_mask = m is not nomask + else: + result = subtract(self, masked_array(f, mask=getmask(other))) + self._data = result.data + self._mask = result.mask + self._shared_mask = 1 + return self + + + + def __idiv__(self, other): + "Divide self by other in place." + t = self._data.dtype.char + f = filled(other, 0) + t1 = f.dtype.char + if t == t1: + pass + elif t in typecodes['Integer']: + if t1 in typecodes['Integer']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Float']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + elif t in typecodes['Complex']: + if t1 in typecodes['Integer']: + f = f.astype(t) + elif t1 in typecodes['Float']: + f = f.astype(t) + elif t1 in typecodes['Complex']: + f = f.astype(t) + else: + raise TypeError, 'Incorrect type for in-place operation.' + else: + raise TypeError, 'Incorrect type for in-place operation.' + mo = getmask(other) + result = divide(self, masked_array(f, mask=mo)) + self._data = result.data + dm = result.raw_mask() + if dm is not self._mask: + self._mask = dm + self._shared_mask = 1 + return self + + def __eq__(self, other): + return equal(self,other) + + def __ne__(self, other): + return not_equal(self,other) + + def __lt__(self, other): + return less(self,other) + + def __le__(self, other): + return less_equal(self,other) + + def __gt__(self, other): + return greater(self,other) + + def __ge__(self, other): + return greater_equal(self,other) + + def astype (self, tc): + "return self as array of given type." + d = self._data.astype(tc) + return array(d, mask=self._mask) + + def byte_swapped(self): + """Returns the raw data field, byte_swapped. Included for consistency + with numeric but doesn't make sense in this context. + """ + return self._data.byte_swapped() + + def compressed (self): + "A 1-D array of all the non-masked data." + d = fromnumeric.ravel(self._data) + if self._mask is nomask: + return array(d) + else: + m = 1 - fromnumeric.ravel(self._mask) + c = fromnumeric.compress(m, d) + return array(c, copy=0) + + def count (self, axis = None): + "Count of the non-masked elements in a, or along a certain axis." + m = self._mask + s = self._data.shape + ls = len(s) + if m is nomask: + if ls == 0: + return 1 + if ls == 1: + return s[0] + if axis is None: + return reduce(lambda x, y:x*y, s) + else: + n = s[axis] + t = list(s) + del t[axis] + return ones(t) * n + if axis is None: + w = fromnumeric.ravel(m).astype(int) + n1 = size(w) + if n1 == 1: + n2 = w[0] + else: + n2 = umath.add.reduce(w) + return n1 - n2 + else: + n1 = size(m, axis) + n2 = sum(m.astype(int), axis) + return n1 - n2 + + def dot (self, other): + "s.dot(other) = innerproduct(s, other)" + return innerproduct(self, other) + + def fill_value(self): + "Get the current fill value." + return self._fill_value + + def filled (self, fill_value=None): + """A numeric array with masked values filled. If fill_value is None, + use self.fill_value(). + + If mask is nomask, copy data only if not contiguous. + Result is always a contiguous, numeric array. +# Is contiguous really necessary now? + """ + d = self._data + m = self._mask + if m is nomask: + if d.flags['CONTIGUOUS']: + return d + else: + return d.copy() + else: + if fill_value is None: + value = self._fill_value + else: + value = fill_value + + if self is masked: + result = numeric.array(value) + else: + try: + result = numeric.array(d, dtype=d.dtype, copy=1) + result[m] = value + except (TypeError, AttributeError): + #ok, can't put that value in here + value = numeric.array(value, dtype=object) + d = d.astype(object) + result = fromnumeric.choose(m, (d, value)) + return result + + def ids (self): + """Return the ids of the data and mask areas""" + return (id(self._data), id(self._mask)) + + def iscontiguous (self): + "Is the data contiguous?" + return self._data.flags['CONTIGUOUS'] + + def itemsize(self): + "Item size of each data item." + return self._data.itemsize + + + def outer(self, other): + "s.outer(other) = outerproduct(s, other)" + return outerproduct(self, other) + + def put (self, values): + """Set the non-masked entries of self to filled(values). + No change to mask + """ + iota = numeric.arange(self.size) + d = self._data + if self._mask is nomask: + ind = iota + else: + ind = fromnumeric.compress(1 - self._mask, iota) + d[ind] = filled(values).astype(d.dtype) + + def putmask (self, values): + """Set the masked entries of self to filled(values). + Mask changed to nomask. + """ + d = self._data + if self._mask is not nomask: + d[self._mask] = filled(values).astype(d.dtype) + self._shared_mask = 0 + self._mask = nomask + + def ravel (self): + """Return a 1-D view of self.""" + if self._mask is nomask: + return masked_array(self._data.ravel()) + else: + return masked_array(self._data.ravel(), self._mask.ravel()) + + def raw_data (self): + """ Obsolete; use data property instead. + The raw data; portions may be meaningless. + May be noncontiguous. Expert use only.""" + return self._data + data = property(fget=raw_data, + doc="The data, but values at masked locations are meaningless.") + + def raw_mask (self): + """ Obsolete; use mask property instead. + May be noncontiguous. Expert use only. + """ + return self._mask + mask = property(fget=raw_mask, + doc="The mask, may be nomask. Values where mask true are meaningless.") + + def reshape (self, *s): + """This array reshaped to shape s""" + d = self._data.reshape(*s) + if self._mask is nomask: + return masked_array(d) + else: + m = self._mask.reshape(*s) + return masked_array(d, m) + + def set_fill_value (self, v=None): + "Set the fill value to v. Omit v to restore default." + if v is None: + v = default_fill_value (self.raw_data()) + self._fill_value = v + + def _get_ndim(self): + return self._data.ndim + ndim = property(_get_ndim, doc=numeric.ndarray.ndim.__doc__) + + def _get_size (self): + return self._data.size + size = property(fget=_get_size, doc="Number of elements in the array.") +## CHECK THIS: signature of numeric.array.size? + + def _get_dtype(self): + return self._data.dtype + dtype = property(fget=_get_dtype, doc="type of the array elements.") + + def item(self, *args): + "Return Python scalar if possible" + if self._mask is not nomask: + m = self._mask.item(*args) + try: + if m[0]: + return masked + except IndexError: + return masked + return self._data.item(*args) + + def itemset(self, *args): + "Set Python scalar into array" + item = args[-1] + args = args[:-1] + self[args] = item + + def tolist(self, fill_value=None): + "Convert to list" + return self.filled(fill_value).tolist() + + def tostring(self, fill_value=None): + "Convert to string" + return self.filled(fill_value).tostring() + + def unmask (self): + "Replace the mask by nomask if possible." + if self._mask is nomask: return + m = make_mask(self._mask, flag=1) + if m is nomask: + self._mask = nomask + self._shared_mask = 0 + + def unshare_mask (self): + "If currently sharing mask, make a copy." + if self._shared_mask: + self._mask = make_mask (self._mask, copy=1, flag=0) + self._shared_mask = 0 + + def _get_ctypes(self): + return self._data.ctypes + + def _get_T(self): + if (self.ndim < 2): + return self + return self.transpose() + + shape = property(_get_shape, _set_shape, + doc = 'tuple giving the shape of the array') + + flat = property(_get_flat, _set_flat, + doc = 'Access array in flat form.') + + real = property(_get_real, _set_real, + doc = 'Access the real part of the array') + + imaginary = property(_get_imaginary, _set_imaginary, + doc = 'Access the imaginary part of the array') + + imag = imaginary + + ctypes = property(_get_ctypes, None, doc="ctypes") + + T = property(_get_T, None, doc="get transpose") + +#end class MaskedArray + +array = MaskedArray + +def isMaskedArray (x): + "Is x a masked array, that is, an instance of MaskedArray?" + return isinstance(x, MaskedArray) + +isarray = isMaskedArray +isMA = isMaskedArray #backward compatibility + +def allclose (a, b, fill_value=1, rtol=1.e-5, atol=1.e-8): + """ Returns true if all components of a and b are equal + subject to given tolerances. + If fill_value is 1, masked values considered equal. + If fill_value is 0, 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 = filled(a) + d2 = filled(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 allequal (a, b, fill_value=1): + """ + 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 = filled(a) + y = filled(b) + d = umath.equal(x, y) + return fromnumeric.alltrue(fromnumeric.ravel(d)) + elif fill_value: + x = filled(a) + y = filled(b) + d = umath.equal(x, y) + dm = array(d, mask=m, copy=0) + return fromnumeric.alltrue(fromnumeric.ravel(filled(dm, 1))) + else: + return 0 + +def masked_values (data, value, rtol=1.e-5, atol=1.e-8, copy=1): + """ + masked_values(data, value, rtol=1.e-5, atol=1.e-8) + Create a masked array; mask is nomask if possible. + If copy==0, and otherwise possible, result + may share data values with original array. + Let d = filled(data, value). Returns d + masked where abs(data-value)<= atol + rtol * abs(value) + if d is of a floating point type. Otherwise returns + masked_object(d, value, copy) + """ + abs = umath.absolute + d = filled(data, value) + if issubclass(d.dtype.type, numeric.floating): + m = umath.less_equal(abs(d-value), atol+rtol*abs(value)) + m = make_mask(m, flag=1) + return array(d, mask = m, copy=copy, + fill_value=value) + else: + return masked_object(d, value, copy=copy) + +def masked_object (data, value, copy=1): + "Create array masked where exactly data equal to value" + d = filled(data, value) + dm = make_mask(umath.equal(d, value), flag=1) + return array(d, mask=dm, copy=copy, fill_value=value) + +def arange(start, stop=None, step=1, dtype=None): + """Just like range() except it returns a array whose type can be specified + by the keyword argument dtype. + """ + return array(numeric.arange(start, stop, step, dtype)) + +arrayrange = arange + +def fromstring (s, t): + "Construct a masked array from a string. Result will have no mask." + return masked_array(numeric.fromstring(s, t)) + +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, 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, m) + +def resize (a, new_shape): + """resize(a, new_shape) returns a new array with the specified shape. + The original array's total size can be any size.""" + m = getmask(a) + if m is not nomask: + m = fromnumeric.resize(m, new_shape) + result = array(fromnumeric.resize(filled(a), new_shape), mask=m) + result.set_fill_value(get_fill_value(a)) + return result + +def new_repeat(a, repeats, axis=None): + """repeat elements of a repeats times along axis + repeats is a sequence of length a.shape[axis] + telling how many times to repeat each element. + """ + af = filled(a) + if isinstance(repeats, types.IntType): + if axis is None: + num = af.size + else: + num = af.shape[axis] + repeats = tuple([repeats]*num) + + m = getmask(a) + if m is not nomask: + m = fromnumeric.repeat(m, repeats, axis) + d = fromnumeric.repeat(af, repeats, axis) + result = masked_array(d, m) + result.set_fill_value(get_fill_value(a)) + return result + + + +def identity(n): + """identity(n) returns the identity matrix of shape n x n. + """ + return array(numeric.identity(n)) + +def indices (dimensions, dtype=None): + """indices(dimensions,dtype=None) returns an array representing a grid + of indices with row-only, and column-only variation. + """ + return array(numeric.indices(dimensions, dtype)) + +def zeros (shape, dtype=float): + """zeros(n, dtype=float) = + an array of all zeros of the given length or shape.""" + return array(numeric.zeros(shape, dtype)) + +def ones (shape, dtype=float): + """ones(n, dtype=float) = + an array of all ones of the given length or shape.""" + return array(numeric.ones(shape, dtype)) + +def count (a, axis = None): + "Count of the non-masked elements in a, or along a certain axis." + a = masked_array(a) + return a.count(axis) + +def power (a, b, third=None): + "a**b" + if third is not None: + raise MAError, "3-argument power not supported." + ma = getmask(a) + mb = getmask(b) + m = mask_or(ma, mb) + fa = filled(a, 1) + fb = filled(b, 1) + if fb.dtype.char in typecodes["Integer"]: + return masked_array(umath.power(fa, fb), m) + md = make_mask(umath.less(fa, 0), flag=1) + m = mask_or(m, md) + if m is nomask: + return masked_array(umath.power(fa, fb)) + else: + fa = numeric.where(m, 1, fa) + return masked_array(umath.power(fa, fb), m) + +def masked_array (a, mask=nomask, fill_value=None): + """masked_array(a, mask=nomask) = + array(a, mask=mask, copy=0, fill_value=fill_value) + """ + return array(a, mask=mask, copy=0, fill_value=fill_value) + +def sum (target, axis=None, dtype=None): + if axis is None: + target = ravel(target) + axis = 0 + return add.reduce(target, axis, dtype) + +def product (target, axis=None, dtype=None): + if axis is None: + target = ravel(target) + axis = 0 + return multiply.reduce(target, axis, dtype) + +def new_average (a, axis=None, weights=None, returned = 0): + """average(a, axis=None, weights=None) + Computes average along indicated axis. + If axis is None, average over the entire array + Inputs can be integer or floating types; result is of type float. + + If weights are given, result is sum(a*weights,axis=0)/(sum(weights,axis=0)*1.0) + weights must have a's shape or be the 1-d with length the size + of a in the given axis. + + If returned, return a tuple: the result and the sum of the weights + or count of values. Results will have the same shape. + + masked values in the weights will be set to 0.0 + """ + a = masked_array(a) + mask = a.mask + ash = a.shape + if ash == (): + ash = (1,) + if axis is None: + if mask is nomask: + if weights is None: + n = add.reduce(a.raw_data().ravel()) + d = reduce(lambda x, y: x * y, ash, 1.0) + else: + w = filled(weights, 0.0).ravel() + n = umath.add.reduce(a.raw_data().ravel() * w) + d = umath.add.reduce(w) + del w + else: + if weights is None: + n = add.reduce(a.ravel()) + w = fromnumeric.choose(mask, (1.0, 0.0)).ravel() + d = umath.add.reduce(w) + del w + else: + w = array(filled(weights, 0.0), float, mask=mask).ravel() + n = add.reduce(a.ravel() * w) + d = add.reduce(w) + del w + else: + if mask is nomask: + if weights is None: + d = ash[axis] * 1.0 + n = umath.add.reduce(a.raw_data(), axis) + else: + w = filled(weights, 0.0) + wsh = w.shape + if wsh == (): + wsh = (1,) + if wsh == ash: + w = numeric.array(w, float, copy=0) + n = add.reduce(a*w, axis) + d = add.reduce(w, axis) + del w + elif wsh == (ash[axis],): + r = [newaxis]*len(ash) + r[axis] = slice(None, None, 1) + w = eval ("w["+ repr(tuple(r)) + "] * ones(ash, float)") + n = add.reduce(a*w, axis) + d = add.reduce(w, axis) + del w, r + else: + raise ValueError, 'average: weights wrong shape.' + else: + if weights is None: + n = add.reduce(a, axis) + w = numeric.choose(mask, (1.0, 0.0)) + d = umath.add.reduce(w, axis) + del w + else: + w = filled(weights, 0.0) + wsh = w.shape + if wsh == (): + wsh = (1,) + if wsh == ash: + w = array(w, float, mask=mask, copy=0) + n = add.reduce(a*w, axis) + d = add.reduce(w, axis) + elif wsh == (ash[axis],): + r = [newaxis]*len(ash) + r[axis] = slice(None, None, 1) + w = eval ("w["+ repr(tuple(r)) + "] * masked_array(ones(ash, float), mask)") + n = add.reduce(a*w, axis) + d = add.reduce(w, axis) + else: + raise ValueError, 'average: weights wrong shape.' + del w + #print n, d, repr(mask), repr(weights) + if n is masked or d is masked: return masked + result = divide (n, d) + del n + + if isinstance(result, MaskedArray): + result.unmask() + if returned: + if not isinstance(d, MaskedArray): + d = masked_array(d) + if not d.shape == result.shape: + d = ones(result.shape, float) * d + d.unmask() + if returned: + return result, d + else: + return result + +def where (condition, x, y): + """where(condition, x, y) is x where condition is nonzero, y otherwise. + condition must be convertible to an integer array. + Answer is always the shape of condition. + The type depends on x and y. It is integer if both x and y are + the value masked. + """ + fc = filled(not_equal(condition, 0), 0) + xv = filled(x) + xm = getmask(x) + yv = filled(y) + ym = getmask(y) + d = numeric.choose(fc, (yv, xv)) + md = numeric.choose(fc, (ym, xm)) + m = getmask(condition) + m = make_mask(mask_or(m, md), copy=0, flag=1) + return masked_array(d, m) + +def choose (indices, t, out=None, mode='raise'): + "Returns array shaped like indices with elements chosen from t" + def fmask (x): + if x is masked: return 1 + return filled(x) + def nmask (x): + 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 = numeric.choose(c, a) + m = numeric.choose(c, masks) + m = make_mask(mask_or(m, getmask(indices)), copy=0, flag=1) + return masked_array(d, m) + +def masked_where(condition, x, copy=1): + """Return x as an array masked where condition is true. + Also masked where x or condition masked. + """ + cm = filled(condition,1) + m = mask_or(getmask(x), cm) + return array(filled(x), copy=copy, mask=m) + +def masked_greater(x, value, copy=1): + "masked_greater(x, value) = x masked where x > value" + return masked_where(greater(x, value), x, copy) + +def masked_greater_equal(x, value, copy=1): + "masked_greater_equal(x, value) = x masked where x >= value" + return masked_where(greater_equal(x, value), x, copy) + +def masked_less(x, value, copy=1): + "masked_less(x, value) = x masked where x < value" + return masked_where(less(x, value), x, copy) + +def masked_less_equal(x, value, copy=1): + "masked_less_equal(x, value) = x masked where x <= value" + return masked_where(less_equal(x, value), x, copy) + +def masked_not_equal(x, value, copy=1): + "masked_not_equal(x, value) = x masked where x != value" + d = filled(x, 0) + c = umath.not_equal(d, value) + m = mask_or(c, getmask(x)) + return array(d, mask=m, copy=copy) + +def masked_equal(x, value, copy=1): + """masked_equal(x, value) = x masked where x == value + For floating point consider masked_values(x, value) instead. + """ + d = filled(x, 0) + c = umath.equal(d, value) + m = mask_or(c, getmask(x)) + return array(d, mask=m, copy=copy) + +def masked_inside(x, v1, v2, copy=1): + """x with mask of all values of x that are inside [v1,v2] + v1 and v2 can be given in either order. + """ + if v2 < v1: + t = v2 + v2 = v1 + v1 = t + d = filled(x, 0) + c = umath.logical_and(umath.less_equal(d, v2), umath.greater_equal(d, v1)) + m = mask_or(c, getmask(x)) + return array(d, mask = m, copy=copy) + +def masked_outside(x, v1, v2, copy=1): + """x with mask of all values of x that are outside [v1,v2] + v1 and v2 can be given in either order. + """ + if v2 < v1: + t = v2 + v2 = v1 + v1 = t + d = filled(x, 0) + c = umath.logical_or(umath.less(d, v1), umath.greater(d, v2)) + m = mask_or(c, getmask(x)) + return array(d, mask = m, copy=copy) + +def reshape (a, *newshape): + "Copy of a with a new shape." + m = getmask(a) + d = filled(a).reshape(*newshape) + if m is nomask: + return masked_array(d) + else: + return masked_array(d, mask=numeric.reshape(m, *newshape)) + +def ravel (a): + "a as one-dimensional, may share data and mask" + m = getmask(a) + d = fromnumeric.ravel(filled(a)) + if m is nomask: + return masked_array(d) + else: + return masked_array(d, mask=numeric.ravel(m)) + +def concatenate (arrays, axis=0): + "Concatenate the arrays along the given axis" + d = [] + for x in arrays: + d.append(filled(x)) + d = numeric.concatenate(d, axis) + for x in arrays: + if getmask(x) is not nomask: break + else: + return masked_array(d) + dm = [] + for x in arrays: + dm.append(getmaskarray(x)) + dm = numeric.concatenate(dm, axis) + return masked_array(d, mask=dm) + +def swapaxes (a, axis1, axis2): + m = getmask(a) + d = masked_array(a).data + if m is nomask: + return masked_array(data=numeric.swapaxes(d, axis1, axis2)) + else: + return masked_array(data=numeric.swapaxes(d, axis1, axis2), + mask=numeric.swapaxes(m, axis1, axis2),) + + +def new_take (a, indices, axis=None, out=None, mode='raise'): + "returns selection of items from a." + m = getmask(a) + # d = masked_array(a).raw_data() + d = masked_array(a).data + if m is nomask: + return masked_array(numeric.take(d, indices, axis)) + else: + return masked_array(numeric.take(d, indices, axis), + mask = numeric.take(m, indices, axis)) + +def transpose(a, axes=None): + "reorder dimensions per tuple axes" + m = getmask(a) + d = filled(a) + if m is nomask: + return masked_array(numeric.transpose(d, axes)) + else: + return masked_array(numeric.transpose(d, axes), + mask = numeric.transpose(m, axes)) + + +def put(a, indices, values, mode='raise'): + """sets storage-indexed locations to corresponding values. + + Values and indices are filled if necessary. + + """ + d = a.raw_data() + ind = filled(indices) + v = filled(values) + numeric.put (d, ind, v) + m = getmask(a) + if m is not nomask: + a.unshare_mask() + numeric.put(a.raw_mask(), ind, 0) + +def putmask(a, mask, values): + "putmask(a, mask, values) sets a where mask is true." + if mask is nomask: + return + numeric.putmask(a.raw_data(), mask, values) + m = getmask(a) + if m is nomask: return + a.unshare_mask() + numeric.putmask(a.raw_mask(), mask, 0) + +def inner(a, b): + """inner(a,b) returns the dot product of two arrays, which has + shape a.shape[:-1] + b.shape[:-1] with elements computed by summing the + product of the elements from the last dimensions of a and b. + Masked elements are replace by zeros. + """ + 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 masked_array(numeric.inner(fa, fb)) + +innerproduct = inner + +def outer(a, b): + """outer(a,b) = {a[i]*b[j]}, has shape (len(a),len(b))""" + 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, m) + +outerproduct = outer + +def dot(a, b): + """dot(a,b) returns matrix-multiplication between a and b. The product-sum + is over the last dimension of a and the second-to-last dimension of b. + Masked values are replaced by zeros. See also innerproduct. + """ + return innerproduct(filled(a, 0), numeric.swapaxes(filled(b, 0), -1, -2)) + +def compress(condition, x, dimension=-1, out=None): + """Select those parts of x for which condition is true. + Masked values in condition are considered false. + """ + c = filled(condition, 0) + m = getmask(x) + if m is not nomask: + m = numeric.compress(c, m, dimension) + d = numeric.compress(c, filled(x), dimension) + return masked_array(d, m) + +class _minimum_operation: + "Object to calculate minima" + def __init__ (self): + """minimum(a, b) or minimum(a) + In one argument case returns the scalar minimum. + """ + pass + + def __call__ (self, a, b=None): + "Execute the call behavior." + if b is None: + m = getmask(a) + if m is nomask: + d = amin(filled(a).ravel()) + return d + ac = a.compressed() + if len(ac) == 0: + return masked + else: + return amin(ac.raw_data()) + else: + return where(less(a, b), a, b) + + def reduce (self, target, axis=0): + """Reduce target along the given axis.""" + m = getmask(target) + if m is nomask: + t = filled(target) + return masked_array (umath.minimum.reduce (t, axis)) + else: + t = umath.minimum.reduce(filled(target, minimum_fill_value(target)), axis) + m = umath.logical_and.reduce(m, axis) + return masked_array(t, m, get_fill_value(target)) + + 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) + d = umath.minimum.outer(filled(a), filled(b)) + return masked_array(d, m) + +minimum = _minimum_operation () + +class _maximum_operation: + "Object to calculate maxima" + def __init__ (self): + """maximum(a, b) or maximum(a) + In one argument case returns the scalar maximum. + """ + pass + + def __call__ (self, a, b=None): + "Execute the call behavior." + if b is None: + m = getmask(a) + if m is nomask: + d = amax(filled(a).ravel()) + return d + ac = a.compressed() + if len(ac) == 0: + return masked + else: + return amax(ac.raw_data()) + else: + return where(greater(a, b), a, b) + + def reduce (self, target, axis=0): + """Reduce target along the given axis.""" + m = getmask(target) + if m is nomask: + t = filled(target) + return masked_array (umath.maximum.reduce (t, axis)) + else: + t = umath.maximum.reduce(filled(target, maximum_fill_value(target)), axis) + m = umath.logical_and.reduce(m, axis) + return masked_array(t, m, get_fill_value(target)) + + 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) + d = umath.maximum.outer(filled(a), filled(b)) + return masked_array(d, m) + +maximum = _maximum_operation () + +def sort (x, axis = -1, fill_value=None): + """If x does not have a mask, return a masked array formed from the + result of numeric.sort(x, axis). + Otherwise, fill x with fill_value. Sort it. + Set a mask where the result is equal to fill_value. + Note that this may have unintended consequences if the data contains the + fill value at a non-masked site. + + If fill_value is not given the default fill value for x's type will be + used. + """ + if fill_value is None: + fill_value = default_fill_value (x) + d = filled(x, fill_value) + s = fromnumeric.sort(d, axis) + if getmask(x) is nomask: + return masked_array(s) + return masked_values(s, fill_value, copy=0) + +def diagonal(a, k = 0, axis1=0, axis2=1): + """diagonal(a,k=0,axis1=0, axis2=1) = the k'th diagonal of a""" + d = fromnumeric.diagonal(filled(a), k, axis1, axis2) + m = getmask(a) + if m is nomask: + return masked_array(d, m) + else: + return masked_array(d, fromnumeric.diagonal(m, k, axis1, axis2)) + +def trace (a, offset=0, axis1=0, axis2=1, dtype=None, out=None): + """trace(a,offset=0, axis1=0, axis2=1) returns the sum along diagonals + (defined by the last two dimenions) of the array. + """ + return diagonal(a, offset, axis1, axis2).sum(dtype=dtype) + +def argsort (x, axis = -1, out=None, fill_value=None): + """Treating masked values as if they have the value fill_value, + return sort indices for sorting along given axis. + if fill_value is None, use get_fill_value(x) + Returns a numpy array. + """ + d = filled(x, fill_value) + return fromnumeric.argsort(d, axis) + +def argmin (x, axis = -1, out=None, fill_value=None): + """Treating masked values as if they have the value fill_value, + return indices for minimum values along given axis. + if fill_value is None, use get_fill_value(x). + Returns a numpy array if x has more than one dimension. + Otherwise, returns a scalar index. + """ + d = filled(x, fill_value) + return fromnumeric.argmin(d, axis) + +def argmax (x, axis = -1, out=None, fill_value=None): + """Treating masked values as if they have the value fill_value, + return sort indices for maximum along given axis. + if fill_value is None, use -get_fill_value(x) if it exists. + Returns a numpy array if x has more than one dimension. + Otherwise, returns a scalar index. + """ + if fill_value is None: + fill_value = default_fill_value (x) + try: + fill_value = - fill_value + except: + pass + d = filled(x, fill_value) + return fromnumeric.argmax(d, axis) + +def fromfunction (f, s): + """apply f to s to create array as in umath.""" + return masked_array(numeric.fromfunction(f, s)) + +def asarray(data, dtype=None): + """asarray(data, dtype) = array(data, dtype, copy=0) + """ + if isinstance(data, MaskedArray) and \ + (dtype is None or dtype == data.dtype): + return data + return array(data, dtype=dtype, copy=0) + +# Add methods to support ndarray interface +# XXX: I is better to to change the masked_*_operation adaptors +# XXX: to wrap ndarray methods directly to create ma.array methods. +from types import MethodType +def _m(f): + return MethodType(f, None, array) +def not_implemented(*args, **kwds): + raise NotImplementedError, "not yet implemented for numpy.ma arrays" +array.all = _m(alltrue) +array.any = _m(sometrue) +array.argmax = _m(argmax) +array.argmin = _m(argmin) +array.argsort = _m(argsort) +array.base = property(_m(not_implemented)) +array.byteswap = _m(not_implemented) + +def _choose(self, *args, **kwds): + return choose(self, args) +array.choose = _m(_choose) +del _choose + +def _clip(self,a_min,a_max,out=None): + return MaskedArray(data = self.data.clip(asarray(a_min).data, + asarray(a_max).data), + mask = mask_or(self.mask, + mask_or(getmask(a_min),getmask(a_max)))) +array.clip = _m(_clip) + +def _compress(self, cond, axis=None, out=None): + return compress(cond, self, axis) +array.compress = _m(_compress) +del _compress + +array.conj = array.conjugate = _m(conjugate) +array.copy = _m(not_implemented) + +def _cumprod(self, axis=None, dtype=None, out=None): + m = self.mask + if m is not nomask: + m = umath.logical_or.accumulate(self.mask, axis) + return MaskedArray(data = self.filled(1).cumprod(axis, dtype), mask=m) +array.cumprod = _m(_cumprod) + +def _cumsum(self, axis=None, dtype=None, out=None): + m = self.mask + if m is not nomask: + m = umath.logical_or.accumulate(self.mask, axis) + return MaskedArray(data=self.filled(0).cumsum(axis, dtype), mask=m) +array.cumsum = _m(_cumsum) + +array.diagonal = _m(diagonal) +array.dump = _m(not_implemented) +array.dumps = _m(not_implemented) +array.fill = _m(not_implemented) +array.flags = property(_m(not_implemented)) +array.flatten = _m(ravel) +array.getfield = _m(not_implemented) + +def _max(a, axis=None, out=None): + if out is not None: + raise TypeError("Output arrays Unsupported for masked arrays") + if axis is None: + return maximum(a) + else: + return maximum.reduce(a, axis) +array.max = _m(_max) +del _max +def _min(a, axis=None, out=None): + if out is not None: + raise TypeError("Output arrays Unsupported for masked arrays") + if axis is None: + return minimum(a) + else: + return minimum.reduce(a, axis) +array.min = _m(_min) +del _min +array.mean = _m(average) +array.nbytes = property(_m(not_implemented)) +array.newbyteorder = _m(not_implemented) +array.nonzero = _m(nonzero) +array.prod = _m(product) + +def _ptp(a,axis=None,out=None): + return a.max(axis,out)-a.min(axis) +array.ptp = _m(_ptp) +array.repeat = _m(repeat) +array.resize = _m(resize) +array.searchsorted = _m(not_implemented) +array.setfield = _m(not_implemented) +array.setflags = _m(not_implemented) +array.sort = _m(not_implemented) # NB: ndarray.sort is inplace + +def _squeeze(self): + try: + result = MaskedArray(data = self.data.squeeze(), + mask = self.mask.squeeze()) + except AttributeError: + result = _wrapit(self, 'squeeze') + return result +array.squeeze = _m(_squeeze) + +array.strides = property(_m(not_implemented)) +array.sum = _m(sum) +def _swapaxes(self,axis1,axis2): + return MaskedArray(data = self.data.swapaxes(axis1, axis2), + mask = self.mask.swapaxes(axis1, axis2)) +array.swapaxes = _m(_swapaxes) +array.take = _m(take) +array.tofile = _m(not_implemented) +array.trace = _m(trace) +array.transpose = _m(transpose) + +def _var(self,axis=None,dtype=None, out=None): + if axis is None: + return numeric.asarray(self.compressed()).var() + a = self.swapaxes(axis,0) + a = a - a.mean(axis=0) + a *= a + a /= a.count(axis=0) + return a.swapaxes(0,axis).sum(axis) +def _std(self,axis=None, dtype=None, out=None): + return (self.var(axis,dtype))**0.5 +array.var = _m(_var) +array.std = _m(_std) + +array.view = _m(not_implemented) +array.round = _m(around) +del _m, MethodType, not_implemented + + +masked = MaskedArray(0, int, mask=1) def repeat(a, repeats, axis=0): - return nca.repeat(a, repeats, axis) + return new_repeat(a, repeats, axis) def average(a, axis=0, weights=None, returned=0): - return nca.average(a, axis, weights, returned) + return new_average(a, axis, weights, returned) def take(a, indices, axis=0): - return nca.average(a, indices, axis=0) + return new_take(a, indices, axis=0) + + |