"""MA: a facility for dealing with missing observations MA is generally used as a Numeric.array look-alike. There are some differences in semantics, see manual. In particular note that slices are copies, not references. by Paul F. Dubois L-264 Lawrence Livermore National Laboratory dubois@users.sourceforge.net Copyright 1999, 2000 Regents of the University of California. Released for unlimited redistribution; see file Legal.htm """ import Numeric import string, types from MA_version import version, version_info __version__ = version __version_info__ = version_info del version, version_info from Precision import * from Numeric import e, pi, NewAxis from activeattr import ActiveAttributes MaskType=Int0 class MAError (Exception): def __init__ (self, args=None): self.args = args def __str__(self): return str(self.args) __repr__ = __str__ class _MaskedValue: "One instance of this class, masked, is created." def __init__ (self, display): self.set_display(display) self.__dict__['_enabled'] = 1 def display (self): "Show what prints for masked values." return self.__dict__['_display'] def set_display (self, s): "set_display(s) sets what prints for masked values." if type(s) != types.StringType: raise MAError, "masked.set_display(s), s must be a string." self.__dict__['_display'] = s def enabled (self): "Is the use of the display value enabled?" return self.__dict__['_enabled'] def enable(flag=1): "Set the enabling flag to flag." self.__dict__['_enabled'] = flag def __str__(self): return self.__dict__['_display'] def __repr__(self): d = self.__dict__['_display'] return "MaskedValue('" + self.__dict__['_display'] +"')" def nope(self, *other): raise MAError, 'Cannot do requested operation with a masked value.' __add__ = nope __radd__ = nope __sub__ = nope __rsub__ = nope __mult__ = nope __rmult__ = nope __div__ = nope __rdiv__ = nope __pow__ = nope __sqrt__ = nope __int__ = nope __neg__ = nope __float__ = nope __abs__ = nope __setitem__ = nope __setslice__ = nope __getitem__ = nope __getslice__ = nope __len__ = nope __getattr__ = nope __setattr__ = nope #if you single index into a masked location you get this object. masked = _MaskedValue('--') # Use single element arrays or scalars. default_real_fill_value = Numeric.array([1.0e20]).astype(Float32) default_complex_fill_value = Numeric.array([1.0e20 + 0.0j]).astype(Complex32) default_character_fill_value = masked default_integer_fill_value = Numeric.array([0]).astype(UnsignedInt8) default_object_fill_value = masked def default_fill_value (obj): "Function to calculate default fill value for an object." try: x = obj.typecode() except AttributeError: t = type(obj) if t is types.FloatType: return default_real_fill_value elif t is types.IntType or t is types.LongType: return default_integer_fill_value elif t is types.StringType: return default_character_fill_value elif t is types.ComplexType: return default_complex_fill_value else: return default_object_fill_value 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 Numeric.absolute(default_integer_fill_value) else: return default_object_fill_value def set_fill_value (a, fill_value): "Set fill value of a if it is a masked array." if isMA(a): a.set_fill_value (fill_value) def getmask (a): """Mask of values in a; could be None. To get an array for sure use getmaskarray.""" if isMA(a): return a.mask() else: return None def getmaskarray (a): "Mask of values in a; an array of zeros if mask is None." 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 (type(m) == Numeric.ArrayType and m.typecode() == 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 is_mask(m): if m is None: return None result = Numeric.array(m, MaskType, copy, 1) else: result = Numeric.array(filled(m,1), MaskType, copy, 1) if flag and not Numeric.sometrue(Numeric.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, MaskType) result.savespace(1) 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(Numeric.logical_or(m1, m2)) def filled (a, value = None, copy=0): """a as a contiguous Numeric array with any masked areas replaced by value No data is copied if a has no masked area and data is already a contiguous array unless copy is 1. if value is None or the special element "masked", fill_value(a) is used instead. If a is already a Numeric array, and copy == 0, then 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 isMA(a): m = a.mask() d = a.raw_data() if m is None: if copy or not d.iscontiguous(): result = Numeric.array(d, typecode=d.typecode(), copy=1, savespace = d.spacesaver()) else: result = d elif value is masked: result = Numeric.choose(m, (d, masked)) else: if value is None: v = a.fill_value() else: v = value tc = d.typecode() result = Numeric.array(d, typecode=tc, copy=1, savespace = d.spacesaver()) Numeric.putmask(result, m, v) result.savespace(d.spacesaver()) elif type(a) == Numeric.ArrayType: if copy or not a.iscontiguous(): result = Numeric.array(a, typecode=a.typecode(), copy=1, savespace=a.spacesaver()) else: result = a else: result = Numeric.array(a) assert type(result) is Numeric.ArrayType assert result.iscontiguous() return result def fill_value (a): """ The fill value of a, if it has one; otherwise, the default fill value for that type. """ if isMA(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 = fill_value(a) t2 = 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): self.y1 = y1 self.y2 = y2 def __call__ (self, x): return Numeric.logical_or(Numeric.greater (x, self.y2), Numeric.less(x, self.y1) ) class domain_tan: "domain_tan(eps) = true where abs(cos(x)) < eps)" def __init__(self, eps): self.eps = eps def __call__ (self, x): return Numeric.less(Numeric.absolute(Numeric.cos(x)), self.eps) class domain_greater: "domain_greater(v)(x) = true where x <= v" def __init__(self, critical_value): self.critical_value = critical_value def __call__ (self, x): return Numeric.less_equal (x, self.critical_value) class domain_greater_equal: "domain_greater_equal(v)(x) = true where x < v" def __init__(self, critical_value): self.critical_value = critical_value def __call__ (self, x): return Numeric.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): # 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) if type(result) is Numeric.ArrayType: return masked_array (result) else: return result else: dx = masked_array(d1, m) result = self.f(filled(dx, self.fill)) if type(result) is Numeric.ArrayType: return masked_array(result, m) elif m[...]: return masked else: return result class domain_safe_divide: """Tester for a domain for x/y Testing for exact type of result seems too expensive so I use a test that is right for single precision float/complex or integer For double precision may mask out some things that would have been ok Better ideas welcome. -- Paul Dubois """ def __init__(self): self.eps = 1.e-35 def __call__ (self, x, y): return Numeric.greater_equal(Numeric.absolute(x)*self.eps, Numeric.absolute(y)) class masked_binary_operation: def __init__ (self, abfunc, fillx=0, filly=0, domain=None): """abfunc(fillx, filly) must be defined. domain(x, y) returns a mask where abfunc(x,y) undefined. abfunc(x, filly) = x for all x to enable reduce. """ self.f = abfunc self.fillx = fillx self.filly = filly self.domain = domain self.__doc__ = getattr(abfunc, "__doc__", str(abfunc)) def __call__ (self, a, b): m = mask_or(getmask(a), getmask(b)) if m is None and self.domain is None: d1 = filled(a, self.fillx) d2 = filled(b, self.filly) result = self.f(d1, d2) if type(result) is Numeric.ArrayType: return masked_array(result) else: return result d1 = filled(a, self.fillx) d2 = filled(b, self.filly) if self.domain is not None: m = mask_or(m, self.domain(d1, d2)) if Numeric.sometrue(Numeric.ravel(m)): a1 = masked_array(d1, m) b1 = masked_array(d2, m) d1 = filled(a1, self.fillx) d2 = filled(b1, self.filly) else: m = None result = self.f(d1, d2) if type(result) is Numeric.ArrayType: return masked_array(result, m) elif m is None: return result elif m[...]: return masked else: return result def reduce (self, target, axis=0): m = getmask(target) t = filled(target, self.filly) if self.domain is not None: m = mask_or(m, self.domain(self.fillx, t)) 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 = Numeric.logical_and.reduce(m, axis) return masked_array(t, m, fill_value(target)) def accumulate (self, target, axis=0): m = getmask(target) t = filled(target, self.filly) if self.domain is not None: m = mask_or(m, self.domain(self.fillx, t)) if m is None: return masked_array (self.f.accumulate (t, axis)) else: t = array(t, typecode=t.typecode(), mask = m, copy=0) t = self.f.accumulate(filled(t, self.filly), axis) return masked_array(t, m, fill_value(target)) sqrt = masked_unary_operation(Numeric.sqrt, 0.0, domain_greater_equal(0.0)) log = masked_unary_operation(Numeric.log, 1.0, domain_greater(0.0)) log10 = masked_unary_operation(Numeric.log10, 1.0, domain_greater(0.0)) exp = masked_unary_operation(Numeric.exp) conjugate = masked_unary_operation(Numeric.conjugate) sin = masked_unary_operation(Numeric.sin) cos = masked_unary_operation(Numeric.cos) tan = masked_unary_operation(Numeric.tan, 0.0, domain_tan(1.e-35)) arcsin = masked_unary_operation(Numeric.arcsin, 0.0, domain_check_interval(-1.0, 1.0)) arccos = masked_unary_operation(Numeric.arccos, 0.0, domain_check_interval(-1.0, 1.0)) arctan = masked_unary_operation(Numeric.arctan) # Missing from Numeric # arcsinh = masked_unary_operation(Numeric.arcsinh) # arccosh = masked_unary_operation(Numeric.arccosh) # arctanh = masked_unary_operation(Numeric.arctanh) sinh = masked_unary_operation(Numeric.sinh) cosh = masked_unary_operation(Numeric.cosh) tanh = masked_unary_operation(Numeric.tanh) absolute = masked_unary_operation(Numeric.absolute) fabs = masked_unary_operation(Numeric.fabs) negative = masked_unary_operation(Numeric.negative) nonzero = masked_unary_operation(Numeric.nonzero) around = masked_unary_operation(Numeric.around) floor = masked_unary_operation(Numeric.floor) sometrue = masked_unary_operation(Numeric.sometrue) alltrue = masked_unary_operation(Numeric.alltrue, 1) logical_not = masked_unary_operation(Numeric.logical_not) add = masked_binary_operation(Numeric.add) subtract = masked_binary_operation(Numeric.subtract) subtract.reduce = None multiply = masked_binary_operation(Numeric.multiply, 1, 1) divide = masked_binary_operation(Numeric.divide, 0, 1, domain_safe_divide()) divide.reduce = None remainder = masked_binary_operation(Numeric.remainder, 0, 1, domain_safe_divide()) remainder.reduce = None fmod = masked_binary_operation(Numeric.fmod, 0, 1, domain_safe_divide()) fmod.reduce = None hypot = masked_binary_operation(Numeric.hypot) hypot.reduce = None arctan2 = masked_binary_operation(Numeric.arctan2, 0.0, 1.0) arctan2.reduce = None equal = masked_binary_operation(Numeric.equal) equal.reduce = None not_equal = masked_binary_operation(Numeric.not_equal) not_equal.reduce = None less_equal = masked_binary_operation(Numeric.less_equal) less_equal.reduce = None greater_equal = masked_binary_operation(Numeric.greater_equal) greater_equal.reduce = None less = masked_binary_operation(Numeric.less) less.reduce = None greater = masked_binary_operation(Numeric.greater) greater.reduce = None logical_and = masked_binary_operation(Numeric.logical_and) logical_or = masked_binary_operation(Numeric.logical_or) logical_xor = masked_binary_operation(Numeric.logical_xor) bitwise_and = masked_binary_operation(Numeric.bitwise_and) bitwise_or = masked_binary_operation(Numeric.bitwise_or) bitwise_xor = masked_binary_operation(Numeric.bitwise_xor) class MA (ActiveAttributes): """Arrays with possibly masked values. Masked values of 1 exclude element from the computation. Construction: x = array(data, typecode=None, copy=1, savespace=0, mask = None, fill_value=None) raw data stored is the result of: Numeric.array(data, typecode=typecode, copy=copy, savespace=savespace) If mask is None there are no masked values. Otherwise mask must be convertible to an array of integers of typecode MaskType, with values 1 or 0, and of 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. If copy == 0, an attempt is made to not copy the data. If savespace is 1, the data is given the spacesaver property, and the mask is replaced by None if all its elements are true. """ def __init__(self, data, typecode=None, copy=1, savespace=None, mask=None, fill_value=None, ): """array(data, typecode=None,copy=1, savespace=None, mask=None, fill_value=None) If data already a Numeric array, its typecode and spacesaver() become the default values for typecode and savespace. """ ActiveAttributes.__init__ (self) tc = typecode ss = savespace if tc is None and type(data) is Numeric.ArrayType: tc = data.typecode() if ss is None: if type(data) is Numeric.ArrayType: ss = data.spacesaver() else: ss = 0 self.__data = Numeric.array(filled(data, copy=copy), typecode=tc, copy=0, savespace=ss) self.__mask = make_mask (mask, flag=ss) self.set_fill_value(fill_value) 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 = Numeric.resize(self.__mask, self.__data.shape) self.__shared_mask = 0 elif nd == 1: self.__data = Numeric.resize(self.__data, self.__mask.shape) else: raise MAError, "Mask and data not compatible." self.set_active_attribute_handler ('shape', MA._get_shape, MA._set_shape) self.set_active_attribute_handler ('flat', MA._get_flat, None) self.set_active_attribute_handler ('real', MA._get_real, None) self.set_active_attribute_handler ('imaginary', MA._get_imaginary, None) def __array__ (self): return self.filled () def _get_shape(self): return self.__data.shape def _set_shape (self, newshape): self.__data.shape = newshape if self.__mask is not None: self.unshare_mask() self.__mask.shape = newshape def _get_flat(self): if self.__mask is None: return masked_array(self.__data.flat, mask=None, fill_value = self.fill_value()) else: return masked_array(self.__data.flat, mask=self.__mask.flat, fill_value = self.fill_value()) def _get_real(self): 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.flat, fill_value = self.fill_value()) def _get_imaginary(self): if self.__mask is None: return masked_array(self.__data.imaginary, mask=None, fill_value = self.fill_value()) else: return masked_array(self.__data.imaginary, mask=self.__mask.flat, fill_value = self.fill_value()) def __str__(self): if masked.enabled(): f = masked else: f = self.fill_value() return str(filled(self, f)) def __repr__(self): with_mask = """*** Masked array, mask present *** Data: %(data)s Mask (fill value %(fill)s) %(mask)s """ without_mask = """*** Masked array, no mask *** Data: %(data)s """ if self.__mask: return with_mask % { 'data': str(self), 'mask': str(self.mask()), 'fill': str(self.fill_value()) } else: return without_mask % { 'data': str(self), } def __float__(self): x = filled(self, 0) return masked_array(float(x), self.__mask)[...] def __int__(self): x = filled(self, 0) return masked_array(int(x), self.__mask)[...] # Note copy semantics here differ from Numeric def __getitem__(self, i): m = self.__mask dout = self.__data[i] ss = self.__data.spacesaver() tc =self.__data.typecode() if type(dout) is Numeric.ArrayType: if m is None: result = array(dout, typecode=tc, copy = 1, savespace=ss) else: result = array(dout, typecode=tc, copy = 1, savespace=ss, mask = m[i], fill_value=self.fill_value()) return result elif m is None or not m[i]: return dout #scalar else: #scalar but masked return masked def __getslice__(self, i, j): m = self.__mask dout = self.__data[i:j] ss = self.__data.spacesaver() tc =self.__data.typecode() if m is None: return array(dout, typecode=tc, copy = 1, savespace=ss) else: return array(dout, typecode=tc, copy = 1, savespace=ss, 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): self.unshare_mask() if value is masked: if self.__mask is None: self.__mask = make_mask_none(self.__data.shape) assert is_mask(self.__mask) self.__mask[index] = 1 return self.__data[index] = filled(value, self.fill_value()) m = getmask(value) if m is None: if self.__mask is not None: self.__mask[index] = 0 else: if self.__mask is None: self.__mask = make_mask_none(self.__data.shape) self.__mask[index] = m def __setslice__(self, i, j, value): self.__setitem__(slice(i,j), value) def __len__ (self): """Return length of first dimension. This is weird but Python's slicing behavior depends on it.""" return len(self.__data) def __abs__(self): return absolute(self) def __neg__(self): return negative(self) def __add__(self, other): return add(self, other) __radd__ = __add__ def __sub__(self, other): return subtract(self, other) def __rsub__(self, other): return subtract(other, self) def __mul__(self, other): return multiply(self, other) __rmul__ = __mul__ def __div__(self, other): return divide(self, other) def __rdiv__(self, other): return divide(other, self) def __pow__(self,other, third=None): return power(self, other, third) def __sqrt__(self): return sqrt(self) def compressed (self): "A 1-D array of all the non-masked data." d = Numeric.ravel(self.__data) if self.__mask is None: return array(d) else: m = 1 - Numeric.ravel(self.__mask) c = Numeric.compress(m, d) return array(c, copy=0) def put (self, values): """Set the non-masked entries of self to filled(values). No change to mask """ iota = Numeric.arange(self.size()) if self.__mask is None: ind = iota else: ind = Numeric.compress(1 - self.__mask, iota) if len(ind) != size(values): raise MAError, "x.put(values) incorrect count of values." Numeric.put (self.__data, ind, filled(values)) def putmask (self, values): """Set the masked entries of self to filled(values). Mask changed to None. """ if self.__mask is not None: iota = Numeric.arange(self.size()) ind = Numeric.compress(self.__mask, iota) if len(ind) != size(values): raise MAError, "x.put(values) incorrect count of values." Numeric.put (self.__data, ind, filled(values)) self.__mask = None self.__shared_mask = 0 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() """ d = self.__data if self.__mask is None and d.iscontiguous(): return d else: if fill_value is None: fill_value = self.fill_value() return filled(self, fill_value) def ids (self): """Return the ids of the data and mask areas""" return (id(self.__data), id(self.__mask)) def mask(self): "Return the data mask, or None." return self.__mask def raw_data (self): """ The non-filled data; portions may be meaningless. Expert use only.""" return self.__data def spacesaver (self): "Get the spacesaver attribute." return self.__data.spacesaver() def savespace (self, value): "Set the spacesaver attribute to value" self.__data.savespace(value) 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.__data) self.__fill_value = v def size (self, axis = None): "Number of elements in array, or in a particular axis." s = self.__data.shape rank = len(s) if rank ==0: s=(1,) if axis is None: return Numeric.multiply.reduce(s) else: return s[axis] def spacesaver (self): "spacesaver() queries the spacesaver attribute." return self.__data.spacesaver() def typecode(self): return self.__data.typecode() 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 astype (self, tc): "return self as array of given type." d = self.__data.astype(tc) d.savespace(self.__data.spacesaver()) return array(d, mask=self.__mask) def iscontiguous (self): return self.__data.iscontiguous() def byte_swapped(self): return self.__data.byte_swapped() array = MA def isMA (x): "Is x a masked array, that is, an instance of MA?" return isinstance(x, MA) isarray = isMA 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 must be positive and << 1.0 The absolute error atol comes into play for those elements of y that are very small or zero; it says how small x must be also. """ m = mask_or(getmask(a), getmask(b)) x = filled(array(a, mask=m, copy=0), fill_value) y = filled(array(b, mask=m, copy=0), 1) d = Numeric.less(Numeric.absolute(x-y), atol + rtol * Numeric.absolute(y)) return Numeric.alltrue(Numeric.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 = Numeric.equal(x, y) return Numeric.alltrue(Numeric.ravel(d)) elif fill_value: x = filled(a) y = filled(b) d = Numeric.equal(x, y) dm = array(d, mask=m, copy=0) return Numeric.alltrue(Numeric.ravel(filled(dm, 1))) else: return 0 def masked_values (data, value, rtol=1.e-5, atol=1.e-8, copy=1, savespace=0): """ masked_values(data, value, rtol=1.e-5, atol=1.e-8) Create a masked array; mask is None if possible. May share data values with original array, but not recommended. """ d = Numeric.array(filled(data, value, copy=copy)) dm = Numeric.less(Numeric.absolute(d-value), atol+rtol*Numeric.absolute(d)) result = array(d, mask = dm, savespace=savespace, copy=0, fill_value=value) return result def masked_object (data, value, copy=1, savespace=0): "Create array masked where exactly data equal to value" d = filled(data, value) dm = Numeric.equal(d, value) result = array(d, mask=dm, copy=copy, savespace=savespace, fill_value=value) return result def rank (a): "Get the rank of a" return len(shape(a)) def shape (a): "Get the shape of a" return array(a,copy=0).shape def arrayrange(start, stop=None, step=1, typecode=None): """Just like range() except it returns a array whose type can be specfied by the keyword argument typecode. """ return array(Numeric.arrayrange(start, stop, step, typecode)) arange = arrayrange 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.""" result = array(Numeric.resize(a, new_shape)) return result def identity(n): """identity(n) returns the identity matrix of shape n x n. """ return array(Numeric.identity(n)) def indices (dimensions, typecode=None): """indices(dimensions,typecode=None) returns an array representing a grid of indices with row-only, and column-only variation. """ return array(Numeric.indices(dimensions, typecode)) def zeros (shape, typecode=Int, savespace=0): """zeros(n, typecode=Int, savespace=0) = an array of all zeros of the given length or shape.""" return array(Numeric.zeros(shape, typecode, savespace)) def ones (shape, typecode=Int, savespace=0): """ones(n, typecode=Int, savespace=0) = an array of all ones of the given length or shape.""" return array(Numeric.ones(shape, typecode, savespace)) def size (a, axis=None): "Get the number of elements in a, or along a certain axis." return masked_array(a).size(axis) def count (a, axis = None): "Count of the non-masked elements." m = getmask(a) if m is None: return size(a, axis) else: n1 = size(a, axis) if axis is None: axis = 0 n2 = sum (m, axis) return n1 - n2 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.typecode() in typecodes["Integer"]: return masked_array(Numeric.power(fa, fb), m) md = make_mask(Numeric.less_equal (fa, 0), flag=1) m = mask_or(m, md) if m is None: return masked_array(Numeric.power(fa, fb)) else: fa = Numeric.where(m, 1, fa) return masked_array(Numeric.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) Use fill_value(a) if None. """ # # This is an unfortunate copy of what is in fill_value # but I want the name fill_value as a parameter. # if fill_value is None: if isMA(a): fill_value = a.fill_value() else: fill_value = default_fill_value(a) return array(a, mask=mask, copy=0, fill_value=fill_value) def sum (a, axis = 0, fill_value=0): "Sum of elements along a certain axis using fill_value for missing." s = Numeric.add.reduce(filled(a, fill_value), axis) if type(s) is Numeric.ArrayType: return masked_array(s) else: return s def product (a, axis = 0, fill_value=1): "Product of elements along axis using fill_value for missing elements." s = Numeric.multiply.reduce(filled(a, fill_value), axis) if type(s) is Numeric.ArrayType: return array(s, copy=0) else: return s def average (a, axis=0, weights=None): """average(a, axis=0, weights=None) Computes average along indicated axis. Masked elements are ignored. Result may equal masked if average cannot be computed. If weights are given, result is sum(a*weights)/sum(weights), with all elements masked in a or in weights ignored. weights if given must have a's shape. Denominator is multiplied by 1.0 to prevent truncation for integers. """ if weights is None: d = 1.0 * count(a, axis) n = sum(a, axis) else: wm = mask_or(getmask(weights), getmask(a)) w = masked_array(weights, wm) d = 1.0 * sum(w, axis) n = sum(a*w, axis) return n / d def choose (condition, t): """Shaped like condition, values t[0] where condition not 0, t[1] otherwise. If t[0] or t[1] is masked, special treatment to preserve type. """ c = not_equal(condition, 0) if t[0] is masked: cm = make_mask(filled(c, 1)) m = mask_or(cm, getmask(t[1])) return masked_array(filled(t[1]), m) elif t[1] is masked: cm = make_mask(filled(c, 1)) m = mask_or(cm, getmask(t[0])) return masked_array(filled(t[0]), m) else: d = Numeric.choose(filled(c, 0), (filled(t[0]), filled(t[1]))) m = mask_or(getmask(c), mask_or(getmask(t[0]), getmask(t[1]))) return masked_array(d, m) def where (condition, x, y): "where(condition, x, y) is x where condition is true, y otherwise" return choose(condition, (y, x)) def masked_where(condition, x, copy=1): """Return x as an array masked where condition is true. Also masked where x or condition masked. """ c = not_equal(condition, 0) cm = make_mask(filled(c,1)) m = mask_or(getmask(x), cm) return array(filled(x), copy=copy, mask=m) def masked_greater(x, value): "masked_greater(x, value) = x masked where x > value" return masked_where(greater(x, value), x) def masked_greater_equal(x, value): "masked_greater_equal(x, value) = x masked where x >= value" return masked_where(greater_equal(x, value), x) def masked_less(x, value): "masked_less(x, value) = x masked where x < value" return masked_where(less(x, value), x) def masked_less_equal(x, value): "masked_less_equal(x, value) = x masked where x <= value" return masked_where(less_equal(x, value), x) def masked_not_equal(x, value): "masked_not_equal(x, value) = x masked where x != value" return masked_where(not_equal(x, value), x) def masked_equal(x, value): """masked_equal(x, value) = x masked where x == value For floating point consider masked_values(x, value) instead. """ return masked_where(equal(x, value), x) def masked_outside(x, v1, v2): "x with mask of all values of x that are outside [v1,v2]" if v2 < v1: t = v2 v2 = v1 v1 = t d = filled(x) c = logical_or(less(d, v1), greater(d, v2)) m = mask_or(c, getmask(x)) return array(d, mask = c) def reshape (a, newshape): "Copy of a with a new shape." m = getmask(a) d = Numeric.reshape(filled(a), 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 = Numeric.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 = filled(a) 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.mask(), ind, 0) def putmask (a, mask, values): "put (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.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. """ 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).flat fb = filled(b,0).flat 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. """ return innerproduct(filled(a,0), swapaxes(filled(b,0), -1, -2)) def compress (condition, x, dimension=-1): cm = array(condition, mask=getmask(x), copy=0) c = cm.filled(0) return array(Numeric.compress(c, filled(x), dimension), copy=0) def maximum (a, b = None): "returns maximum element of a single array, or elementwise" if b is None: xm = getmask(a) if xm is None: return max(filled(a)) else: xc = a.compressed() if(len(xc) == 0): raise MAError, "maximum element does not exist" return max(filled(xc)) else: return choose(less(b,a), (b, a)) def minimum (a, b = None): "returns minimum element of a single array, or elementwise" if b is None: xm = getmask(a) if xm is None: return min(filled(a)) else: xc = a.compressed() if(len(xc) == 0): raise MAError, "minimum element does not exist" return min(filled(xc)) else: return choose(less(b,a), (a, b)) def sort (x, axis = -1, fill_value=None): """Sort x along the given axis, using fill_value if given as the sort value for any masked elements. Shape of returned value is the same as that of x. """ d = filled(x, fill_value) s = Numeric.argsort(d, axis) ds = Numeric.take (d, s, axis) m = getmask(x) if m is None: return array(ds, copy=0) else: return array(ds, mask=Numeric.take(m, s, axis), copy=0, fill_value=fill_value) 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 fill_value(x) """ d = filled(x, fill_value) return Numeric.argsort(d, axis) def asarray(data, typecode=None): """asarray(data, typecode=None) = array(data, typecode=None, copy=0)""" return array(data, typecode=typecode, copy=0) # This section is stolen from a post about how to limit array printing. __MaxElements = 300 #Maximum size for printing def limitedArrayRepr(a, max_line_width = None, precision = None, suppress_small = None): global __MaxElements s = a.shape elems = Numeric.multiply.reduce(s) if elems > __MaxElements: if len(s) > 1: return 'array (%s) , type = %s, has %d elements' % \ (string.join(map(str, s), ","), a.typecode(), elems) else: return Numeric.array2string (a[:__MaxElements], max_line_width, precision, suppress_small,',',0) + \ ('\n + %d more elements' % (elems - __MaxElements)) else: return Numeric.array2string (a, max_line_width, precision, suppress_small,',',0) __original_str = Numeric.array_str __original_repr = Numeric.array_repr def set_print_limit (m=0): "Set the maximum # of elements for printing arrays. <=0 = no limit" import multiarray global __MaxElements n = int(m) __MaxElements = n if n <= 0: Numeric.array_str = __original_str Numeric.array_repr = __original_repr multiarray.set_string_function(__original_str, 0) multiarray.set_string_function(__original_repr, 1) else: Numeric.array_str = limitedArrayRepr Numeric.array_repr = limitedArrayRepr multiarray.set_string_function(limitedArrayRepr, 0) multiarray.set_string_function(limitedArrayRepr, 1) def get_print_limit (): "Get the maximum # of elements for printing arrays. " return __MaxElements set_print_limit(__MaxElements)