"""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)
syntax highlighted by Code2HTML, v. 0.9.1