diff options
Diffstat (limited to 'numpy/base/ma.py')
| -rw-r--r-- | numpy/base/ma.py | 2062 |
1 files changed, 2062 insertions, 0 deletions
diff --git a/numpy/base/ma.py b/numpy/base/ma.py new file mode 100644 index 000000000..245351349 --- /dev/null +++ b/numpy/base/ma.py @@ -0,0 +1,2062 @@ +"""MA: a facility for dealing with missing observations +MA is generally used as a scipy.array look-alike. +by Paul F. Dubois. + +Copyright 1999, 2000, 2001 Regents of the University of California. +Released for unlimited redistribution. +Adapted for scipy_core 2005 by Travis Oliphant and +(mainly) Paul Dubois. +""" +import string, types, sys + +import umath +import oldnumeric +import function_base +from numeric import e, pi, newaxis, ndarray, inf +from oldnumeric import typecodes, amax, amin +from numerictypes import * +import numeric + + +MaskType=bool_ +divide_tolerance = 1.e-35 + +class MAError (Exception): + def __init__ (self, args=None): + "Create an exception" + self.args = args + def __str__(self): + "Calculate the string representation" + return str(self.args) + __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) + +#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.dtypechar + 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.dtypechar + 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.dtypechar + 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 None. + Returns None 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 None + +def getmaskarray (a): + """Mask of values in a; an array of zeros if mask is None + 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 None: + return make_mask_none(shape(a)) + else: + return m + +def is_mask (m): + """Is m a legal mask? Does not check contents, only type. + """ + if m is None or (isinstance(m, ndarray) and \ + m.dtype is MaskType): + return 1 + else: + return 0 + +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 None. Does not check + that contents must be 0s and 1s. + if flag, return None if m contains no true elements. + """ + if m is None: + return None + elif isinstance(m, ndarray): + if m.dtype 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 oldnumeric.sometrue(oldnumeric.ravel(result)): + return None + 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 None as false. + Result may equal m1 or m2 if the other is None. + """ + if m1 is None: return make_mask(m2) + if m2 is None: 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)) + + 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)) + if m is None: + result = self.f(d1, *args, **kwargs) + if type(result) is ndarray: + return masked_array (result) + else: + return result + else: + dx = masked_array(d1, m) + result = self.f(filled(dx, self.fill), *args, **kwargs) + if type(result) is ndarray: + return masked_array(result, m) + elif m[...]: + return masked + else: + return result + + 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)) + + 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 oldnumeric.sometrue(t, None): + d2 = where(t, self.filly, d2) + mb = mask_or(mb, t) + m = mask_or(ma, mb) + if m is None: + result = self.f(d1, d2) + if type(result) is ndarray: + return masked_array(result) + else: + return result + result = self.f(d1, d2) + if type(result) is ndarray: + if m.shape != result.shape: + m = mask_or(getmaskarray(a), getmaskarray(b)) + return masked_array(result, m) + elif m[...]: + return masked + else: + return result + 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)) + + def __call__ (self, a, b, *args, **kwargs): + "Execute the call behavior." + m = mask_or(getmask(a), getmask(b)) + if m is None: + d1 = filled(a, self.fillx) + d2 = filled(b, self.filly) + result = self.f(d1, d2, *args, **kwargs) + if type(result) is ndarray: + return masked_array(result) + else: + return result + d1 = filled(a, self.fillx) + d2 = filled(b, self.filly) + result = self.f(d1, d2, *args, **kwargs) + if type(result) is ndarray: + if m.shape != result.shape: + m = mask_or(getmaskarray(a), getmaskarray(b)) + return masked_array(result, m) + elif m[...]: + return masked + else: + return result + + def reduce (self, target, axis=0): + """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 None: + m = make_mask(m, copy=1) + m.shape = (1,) + if m is None: + return masked_array (self.f.reduce (t, axis)) + else: + t = masked_array (t, m) + t = self.f.reduce(filled(t, self.filly), axis) + 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 None and mb is None: + m = None + 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) +arctanh = masked_unary_operation(umath.arctanh) +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) +nonzero = masked_unary_operation(oldnumeric.nonzero) +around = masked_unary_operation(function_base.round_) +floor = masked_unary_operation(umath.floor) +ceil = masked_unary_operation(umath.ceil) +sometrue = masked_unary_operation(oldnumeric.sometrue) +alltrue = masked_unary_operation(oldnumeric.alltrue, 1) +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) +logical_or = masked_binary_operation(umath.logical_or) +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 oldnumeric.rank(filled(object)) + +def shape (object): + return oldnumeric.shape(filled(object)) + +def size (object, axis=None): + return oldnumeric.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, fortran=False, + mask = None, fill_value=None) + + If copy=False, every effort is made not to copy the data: + If data is a MaskedArray, and argument mask=None, + 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 dtypechar is not None and + is != data.dtypechar 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=dtypechar, copy=copy) + + If mask is None 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. + """ + def __init__(self, data, dtype=None, copy=True, fortran=False, + mask=None, fill_value=None): + """array(data, dtype=None, copy=True, fortran=False, mask=None, fill_value=None) + If data already a numeric array, its dtype becomes the default value of dtype. + """ + tc = dtype + need_data_copied = copy + if isinstance(data, MaskedArray): + c = data.data + ctc = c.dtypechar + if tc is None: + tc = ctc + elif dtype2char(tc) != ctc: + need_data_copied = True + if mask is None: + mask = data.mask + elif mask is not None: #attempting to change the mask + need_data_copied = True + + elif isinstance(data, ndarray): + c = data + ctc = c.dtypechar + if tc is None: + tc = ctc + elif dtype2char(tc) != ctc: + need_data_copied = True + else: + need_data_copied = False #because I'll do it now + c = numeric.array(data, dtype=tc, copy=True, fortran=fortran) + + if need_data_copied: + if tc == ctc: + self._data = numeric.array(c, dtype=tc, copy=True, fortran=fortran) + else: + self._data = c.astype(tc) + else: + self._data = c + + if mask is None: + self._mask = None + self._shared_mask = 0 + else: + self._mask = make_mask (mask) + if self._mask is None: + 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 = oldnumeric.resize(self._mask, self._data.shape) + self._shared_mask = 0 + elif nd == 1: + self._data = oldnumeric.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): + "Special hook for numeric. Converts to numeric if possible." + if self._mask is not None: + if oldnumeric.ravel(self._mask).any(): + raise MAError, \ + """Cannot automatically convert masked array to numeric because data + is masked in one or more locations. + """ + else: # Mask is all false + # Optimize to avoid future invocations of this section. + self._mask = None + self._shared_mask = 0 + if t: + return self._data.astype(t) + else: + return self._data + + 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 None: + self._mask = self._mask.copy() + self._mask.shape = newshape + + def _get_flat(self): + """Calculate the flat value. + """ + if self._mask is None: + return masked_array(self._data.ravel(), mask=None, + 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 None: + return masked_array(self._data.real, mask=None, + fill_value = self.fill_value()) + else: + return masked_array(self._data.real, mask=self._mask.ravel(), + 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 None: + return masked_array(self._data.imag, mask=None, + fill_value = self.fill_value()) + else: + return masked_array(self._data.imag, mask=self._mask.ravel(), + 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 + 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 None: + 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 None: + 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 None: + 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 None: + 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 None: + 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) + + def __getslice__(self, i, j): + "Get slice described by i, j" + self.unshare_mask() + m = self._mask + dout = self._data[i:j] + if m is None: + return masked_array(dout, fill_value=self._fill_value) + else: + return masked_array(dout, mask = m[i:j], 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 the masked element.' + if value is masked: + if self._mask is None: + 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 None: + if self._mask is not None: + self.unshare_mask() + self._mask[index] = False + else: + if self._mask is None: + self._mask = make_mask_none(d.shape) + self._shared_mask = True + else: + self.unshare_mask() + self._mask[index] = m + + def __setslice__(self, i, j, value): + "Set slice i:j; if value is masked, mask those locations." + d = self._data + if self is masked: + raise MAError, "Cannot alter the 'masked' object." + if value is masked: + if self._mask is None: + self._mask = make_mask_none(d.shape) + self._shared_mask = False + self._mask[i:j] = True + return + m = getmask(value) + value = filled(value).astype(d.dtype) + d[i:j] = value + if m is None: + if self._mask is not None: + self.unshare_mask() + self._mask[i:j] = False + else: + if self._mask is None: + self._mask = make_mask_none(self._data.shape) + self._shared_mask = False + self._mask[i:j] = m + + 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.dtypechar + f = filled(other,0) + t1 = f.dtypechar + 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 None: + self._data += f + m = getmask(other) + self._mask = m + self._shared_mask = m is not None + 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.dtypechar + f = filled(other,0) + t1 = f.dtypechar + 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 None: + self._data *= f + m = getmask(other) + self._mask = m + self._shared_mask = m is not None + 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.dtypechar + f = filled(other,0) + t1 = f.dtypechar + 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 None: + self._data -= f + m = getmask(other) + self._mask = m + self._shared_mask = m is not None + 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.dtypechar + f = filled(other,0) + t1 = f.dtypechar + 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 = oldnumeric.ravel(self._data) + if self._mask is None: + return array(d) + else: + m = 1 - oldnumeric.ravel(self._mask) + c = oldnumeric.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 None: + 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 = oldnumeric.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 None, 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 None: + 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).reshape(*d.shape) + else: + try: + result = numeric.array(d, dtype=d.dtype, copy=1) + result[m] = value + except: + #ok, can't put that value in here + value = numeric.array(value, dtype=object) + d = d.astype(object) + result = oldnumeric.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 None: + ind = iota + else: + ind = oldnumeric.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 None. + """ + d = self._data + if self._mask is not None: + d[self._mask] = filled(values).astype(d.dtype) + self._shared_mask = 0 + self._mask = None + + def ravel (self): + """Return a 1-D view of self.""" + if self._mask is None: + 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 None. 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 None: + 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_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_dtypechar(self): + return self._data.dtypechar + dtypechar = property(fget=_get_dtypechar, doc="type character of the array.") + + def _get_dtype(self): + return self._data.dtype + dtype = property(fget=_get_dtype, doc="type of the array elements.") + + def item(self): + "Return Python scalar if possible." + if self._mask is not None: + m = oldnumeric.ravel(self._mask) + try: + if m[0]: + return masked + except IndexError: + return masked + return self._data.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 None if possible." + if self._mask is None: return + m = make_mask(self._mask, flag=1) + if m is None: + self._mask = None + 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 + + 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 + +#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 oldnumeric.alltrue(oldnumeric.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 None: + x = filled(a) + y = filled(b) + d = umath.equal(x, y) + return oldnumeric.alltrue(oldnumeric.ravel(d)) + elif fill_value: + x = filled(a) + y = filled(b) + d = umath.equal(x, y) + dm = array(d, mask=m, copy=0) + return oldnumeric.alltrue(oldnumeric.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 None 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, 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 arrayrange(start, stop=None, step=1, dtype=None): + """Just like range() except it returns a array whose type can be specified + by the keyword argument dtypechar. + """ + return array(numeric.arrayrange(start, stop, step, dtype)) + +arange = arrayrange + +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 None: + 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 None: + 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 None: + m = oldnumeric.resize(m, new_shape) + result = array(oldnumeric.resize(filled(a), new_shape), mask=m) + result.set_fill_value(get_fill_value(a)) + return result + +def repeat(a, repeats, axis=0): + """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): + repeats = tuple([repeats]*(shape(af)[axis])) + + m = getmask(a) + if m is not None: + m = oldnumeric.repeat(m, repeats, axis) + d = oldnumeric.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=int): + """zeros(n, dtype=int) = + an array of all zeros of the given length or shape.""" + return array(numeric.zeros(shape, dtype)) + +def ones (shape, dtype=int): + """ones(n, dtype=int) = + 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.dtypechar in typecodes["Integer"]: + return masked_array(umath.power(fa, fb), m) + md = make_mask(umath.less_equal (fa, 0), flag=1) + m = mask_or(m, md) + if m is None: + 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=None, fill_value=None): + """masked_array(a, mask=None) = + array(a, mask=mask, copy=0, fill_value=fill_value) + """ + return array(a, mask=mask, copy=0, fill_value=fill_value) + +sum = add.reduce +product = multiply.reduce + +def average (a, axis=0, weights=None, returned = 0): + """average(a, axis=0, 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)/(sum(weights)*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 None: + 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 = oldnumeric.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 None: + 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],): + ni = 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],): + ni = 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) + if x is masked: + xv = 0 + xm = 1 + else: + xv = filled(x) + xm = getmask(x) + if xm is None: xm = 0 + if y is masked: + yv = 0 + ym = 1 + else: + yv = filled(y) + ym = getmask(y) + if ym is None: ym = 0 + 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): + "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 None: 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 None: + 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 = oldnumeric.ravel(filled(a)) + if m is None: + 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 None: 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 take (a, indices, axis=0): + "take(a, indices, axis=0) returns selection of items from a." + m = getmask(a) + d = masked_array(a).raw_data() + if m is None: + 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): + "transpose(a, axes=None) reorder dimensions per tuple axes" + m = getmask(a) + d = filled(a) + if m is None: + 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): + """put(a, indices, values) 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 None: + 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 None: + return + numeric.putmask(a.raw_data(), mask, values) + m = getmask(a) + if m is None: return + a.unshare_mask() + numeric.putmask(a.raw_mask(), mask, 0) + +def innerproduct(a,b): + """innerproduct(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.innerproduct(fa, fb)) + +def outerproduct(a, b): + """outerproduct(a,b) = {a[i]*b[j]}, has shape (len(a),len(b))""" + fa = filled(a,0).ravel() + fb = filled(b,0).ravel() + d = numeric.outerproduct(fa, fb) + ma = getmask(a) + mb = getmask(b) + if ma is None and mb is None: + return masked_array(d) + ma = getmaskarray(a) + mb = getmaskarray(b) + m = make_mask(1-numeric.outerproduct(1-ma,1-mb), copy=0) + return masked_array(d, m) + +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): + """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 None: + 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 None: + 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 None: + 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 None and mb is None: + m = None + 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 None: + 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 None: + 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 None and mb is None: + m = None + 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 = oldnumeric.sort(d, axis) + if getmask(x) is None: + 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 = oldnumeric.diagonal(filled(a), k, axis1, axis2) + m = getmask(a) + if m is None: + return masked_array(d, m) + else: + return masked_array(d, oldnumeric.diagonal(m, k, axis1, axis2)) + +def argsort (x, axis = -1, 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 scipy array. + """ + d = filled(x, fill_value) + return oldnumeric.argsort(d, axis) + +def argmin (x, axis = -1, 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 scipy array if x has more than one dimension. + Otherwise, returns a scalar index. + """ + d = filled(x, fill_value) + return oldnumeric.argmin(d, axis) + +def argmax (x, axis = -1, 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 scipy 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 oldnumeric.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) + +masked = MaskedArray([0], int, mask=[1])[0:0] +masked = masked[0:0] |