path: root/numpy/ma/core.py

# pylint: disable-msg=E1002
"""MA: a facility for dealing with missing observations
MA is generally used as a numpy.array look-alike.
by Paul F. Dubois.

Copyright 1999, 2000, 2001 Regents of the University of California.
Released for unlimited redistribution.
Adapted for numpy_core 2005 by Travis Oliphant and
(mainly) Paul Dubois.

Subclassing of the base ndarray 2006 by Pierre Gerard-Marchant.
pgmdevlist_AT_gmail_DOT_com
Improvements suggested by Reggie Dugard (reggie_AT_merfinllc_DOT_com)

:author: Pierre Gerard-Marchant
:contact: pierregm_at_uga_dot_edu
"""
__author__ = "Pierre GF Gerard-Marchant"
__docformat__ = "restructuredtext en"

__all__ = ['MAError', 'MaskType', 'MaskedArray',
           'bool_', 'complex_', 'float_', 'int_', 'object_',
           'abs', 'absolute', 'add', 'all', 'allclose', 'allequal', 'alltrue',
           'amax', 'amin', 'anom', 'anomalies', 'any', 'arange',
           'arccos', 'arccosh', 'arcsin', 'arcsinh', 'arctan', 'arctan2',
           'arctanh', 'argmax', 'argmin', 'argsort', 'around',
           'array', 'asarray','asanyarray',
           'bitwise_and', 'bitwise_or', 'bitwise_xor',
           'ceil', 'choose', 'clip', 'common_fill_value', 'compress',
           'compressed', 'concatenate', 'conjugate', 'cos', 'cosh', 'count',
           'default_fill_value', 'diagonal', 'divide', 'dump', 'dumps',
           'empty', 'empty_like', 'equal', 'exp',
           'fabs', 'fmod', 'filled', 'floor', 'floor_divide','fix_invalid',
           'frombuffer', 'fromfunction',
           'getdata','getmask', 'getmaskarray', 'greater', 'greater_equal',
           'hypot',
           'identity', 'ids', 'indices', 'inner', 'innerproduct',
           'isMA', 'isMaskedArray', 'is_mask', 'is_masked', 'isarray',
           'left_shift', 'less', 'less_equal', 'load', 'loads', 'log', 'log10',
           'logical_and', 'logical_not', 'logical_or', 'logical_xor',
           'make_mask', 'make_mask_none', 'mask_or', 'masked',
           'masked_array', 'masked_equal', 'masked_greater',
           'masked_greater_equal', 'masked_inside', 'masked_invalid',
           'masked_less','masked_less_equal', 'masked_not_equal',
           'masked_object','masked_outside', 'masked_print_option',
           'masked_singleton','masked_values', 'masked_where', 'max', 'maximum',
           'mean', 'min', 'minimum', 'multiply',
           'negative', 'nomask', 'nonzero', 'not_equal',
           'ones', 'outer', 'outerproduct',
           'power', 'product', 'ptp', 'put', 'putmask',
           'rank', 'ravel', 'remainder', 'repeat', 'reshape', 'resize',
           'right_shift', 'round_',
           'set_fill_value', 'shape', 'sin', 'sinh', 'size', 'sometrue', 'sort',
           'sqrt', 'std', 'subtract', 'sum', 'swapaxes',
           'take', 'tan', 'tanh', 'transpose', 'true_divide',
           'var', 'where',
           'zeros']

import sys
import types
import cPickle
import operator

import numpy
from numpy.core import bool_, complex_, float_, int_, object_, str_

import numpy.core.umath as umath
import numpy.core.fromnumeric  as fromnumeric
import numpy.core.numeric as numeric
import numpy.core.numerictypes as ntypes
from numpy import bool_, dtype, typecodes, amax, amin, ndarray, iscomplexobj
from numpy import expand_dims as n_expand_dims
from numpy import array as narray
import warnings


MaskType = bool_
nomask = MaskType(0)

divide_tolerance = 1.e-35
numpy.seterr(all='ignore')

def doc_note(note):
    return "\nNotes\n-----\n%s" % note

#####--------------------------------------------------------------------------
#---- --- Exceptions ---
#####--------------------------------------------------------------------------
class MAError(Exception):
    "Class for MA related errors."
    def __init__ (self, args=None):
        "Creates an exception."
        Exception.__init__(self,args)
        self.args = args
    def __str__(self):
        "Calculates the string representation."
        return str(self.args)
    __repr__ = __str__

#####--------------------------------------------------------------------------
#---- --- Filling options ---
#####--------------------------------------------------------------------------
# b: boolean - c: complex - f: floats - i: integer - O: object - S: string
default_filler = {'b': True,
                  'c' : 1.e20 + 0.0j,
                  'f' : 1.e20,
                  'i' : 999999,
                  'O' : '?',
                  'S' : 'N/A',
                  'u' : 999999,
                  'V' : '???',
                  }
max_filler = ntypes._minvals
max_filler.update([(k,-numpy.inf) for k in [numpy.float32, numpy.float64]])
min_filler = ntypes._maxvals
min_filler.update([(k,numpy.inf) for k in [numpy.float32, numpy.float64]])
if 'float128' in ntypes.typeDict:
    max_filler.update([(numpy.float128,-numpy.inf)])
    min_filler.update([(numpy.float128, numpy.inf)])

def default_fill_value(obj):
    """Calculate the default fill value for the argument object.

    """
    if hasattr(obj,'dtype'):
        defval = default_filler[obj.dtype.kind]
    elif isinstance(obj, numeric.dtype):
        defval = default_filler[obj.kind]
    elif isinstance(obj, float):
        defval = default_filler['f']
    elif isinstance(obj, int) or isinstance(obj, long):
        defval = default_filler['i']
    elif isinstance(obj, str):
        defval = default_filler['S']
    elif isinstance(obj, complex):
        defval = default_filler['c']
    else:
        defval = default_filler['O']
    return defval

def minimum_fill_value(obj):
    """Calculate the default fill value suitable for taking the
    minimum of ``obj``.

    """
    if hasattr(obj, 'dtype'):
        objtype = obj.dtype
        filler = min_filler[objtype]
        if filler is None:
            raise TypeError, 'Unsuitable type for calculating minimum.'
        return filler
    elif isinstance(obj, float):
        return min_filler[ntypes.typeDict['float_']]
    elif isinstance(obj, int):
        return min_filler[ntypes.typeDict['int_']]
    elif isinstance(obj, long):
        return min_filler[ntypes.typeDict['uint']]
    elif isinstance(obj, numeric.dtype):
        return min_filler[obj]
    else:
        raise TypeError, 'Unsuitable type for calculating minimum.'

def maximum_fill_value(obj):
    """Calculate the default fill value suitable for taking the maximum
    of ``obj``.

    """
    if hasattr(obj, 'dtype'):
        objtype = obj.dtype
        filler = max_filler[objtype]
        if filler is None:
            raise TypeError, 'Unsuitable type for calculating minimum.'
        return filler
    elif isinstance(obj, float):
        return max_filler[ntypes.typeDict['float_']]
    elif isinstance(obj, int):
        return max_filler[ntypes.typeDict['int_']]
    elif isinstance(obj, long):
        return max_filler[ntypes.typeDict['uint']]
    elif isinstance(obj, numeric.dtype):
        return max_filler[obj]
    else:
        raise TypeError, 'Unsuitable type for calculating minimum.'


def _check_fill_value(fill_value, dtype):
    descr = numpy.dtype(dtype).descr
    if fill_value is None:
        if len(descr) > 1:
            fill_value = [default_fill_value(numeric.dtype(d[1]))
                          for d in descr]
        else:
            fill_value = default_fill_value(dtype)
    else:
        fill_value = narray(fill_value).tolist()
        fval = numpy.resize(fill_value, len(descr))
        if len(descr) > 1:
            fill_value = [numpy.asarray(f).astype(d[1]).item()
                          for (f,d) in zip(fval, descr)]
        else:
            fill_value = narray(fval, copy=False, dtype=dtype).item()
    return fill_value


def set_fill_value(a, fill_value):
    """Set the filling value of a, if a is a masked array.  Otherwise,
    do nothing.

    Returns
    -------
    None

    """
    if isinstance(a, MaskedArray):
        a._fill_value = _check_fill_value(fill_value, a.dtype)
    return

def get_fill_value(a):
    """Return the filling value of a, if any.  Otherwise, returns the
    default filling value for that type.

    """
    if isinstance(a, MaskedArray):
        result = a.fill_value
    else:
        result = default_fill_value(a)
    return result

def common_fill_value(a, b):
    """Return the common filling value of a and b, if any.
    If a and b have different filling values, returns None.

    """
    t1 = get_fill_value(a)
    t2 = get_fill_value(b)
    if t1 == t2:
        return t1
    return None


#####--------------------------------------------------------------------------
def filled(a, value = None):
    """Return a as an array with masked data replaced by value.  If
    value is None, get_fill_value(a) is used instead.  If a is already
    a ndarray, a itself is returned.

    Parameters
    ----------
    a : maskedarray or array_like
        An input object.
    value : {var}, optional
        Filling value. If not given, the output of get_fill_value(a)
        is used instead.

    Returns
    -------
    a : array_like

    """
    if hasattr(a, 'filled'):
        return a.filled(value)
    elif isinstance(a, ndarray):
        # Should we check for contiguity ? and a.flags['CONTIGUOUS']:
        return a
    elif isinstance(a, dict):
        return narray(a, 'O')
    else:
        return narray(a)

#####--------------------------------------------------------------------------
def get_masked_subclass(*arrays):
    """Return the youngest subclass of MaskedArray from a list of
    (masked) arrays.  In case of siblings, the first takes over.

    """
    if len(arrays) == 1:
        arr = arrays[0]
        if isinstance(arr, MaskedArray):
            rcls = type(arr)
        else:
            rcls = MaskedArray
    else:
        arrcls = [type(a) for a in arrays]
        rcls = arrcls[0]
        if not issubclass(rcls, MaskedArray):
            rcls = MaskedArray
        for cls in arrcls[1:]:
            if issubclass(cls, rcls):
                rcls = cls
    return rcls

#####--------------------------------------------------------------------------
def get_data(a, subok=True):
    """Return the _data part of a (if any), or a as a ndarray.

    Parameters
    ----------
    a : array_like
        A ndarray or a subclass of.
    subok : bool
        Whether to force the output to a 'pure' ndarray (False) or to
        return a subclass of ndarray if approriate (True).

    """
    data = getattr(a, '_data', numpy.array(a, subok=subok))
    if not subok:
        return data.view(ndarray)
    return data
getdata = get_data

def fix_invalid(a, copy=True, fill_value=None):
    """Return (a copy of) a where invalid data (nan/inf) are masked
    and replaced by fill_value.

    Note that a copy is performed by default (just in case...).

    Parameters
    ----------
    a : array_like
        A (subclass of) ndarray.
    copy : bool
        Whether to use a copy of a (True) or to fix a in place (False).
    fill_value : {var}, optional
        Value used for fixing invalid data.  If not given, the output
        of get_fill_value(a) is used instead.

    Returns
    -------
    b : MaskedArray

    """
    a = masked_array(a, copy=copy, subok=True)
    #invalid = (numpy.isnan(a._data) | numpy.isinf(a._data))
    invalid = numpy.logical_not(numpy.isfinite(a._data))
    if not invalid.any():
        return a
    a._mask |= invalid
    if fill_value is None:
        fill_value = a.fill_value
    a._data[invalid] = fill_value
    return a



#####--------------------------------------------------------------------------
#---- --- Ufuncs ---
#####--------------------------------------------------------------------------
ufunc_domain = {}
ufunc_fills = {}

class _DomainCheckInterval:
    """Define a valid interval, so that :

    ``domain_check_interval(a,b)(x) = true`` where
    ``x < a`` or ``x > b``.

    """
    def __init__(self, a, b):
        "domain_check_interval(a,b)(x) = true where x < a or y > b"
        if (a > b):
            (a, b) = (b, a)
        self.a = a
        self.b = b

    def __call__ (self, x):
        "Execute the call behavior."
        return umath.logical_or(umath.greater (x, self.b),
                                umath.less(x, self.a))
#............................
class _DomainTan:
    """Define a valid interval for the `tan` function, so that:

    ``domain_tan(eps) = True`` where ``abs(cos(x)) < eps``

    """
    def __init__(self, eps):
        "domain_tan(eps) = true where abs(cos(x)) < eps)"
        self.eps = eps
    def __call__ (self, x):
        "Executes the call behavior."
        return umath.less(umath.absolute(umath.cos(x)), self.eps)
#............................
class _DomainSafeDivide:
    """Define a domain for safe division."""
    def __init__ (self, tolerance=divide_tolerance):
        self.tolerance = tolerance
    def __call__ (self, a, b):
        return umath.absolute(a) * self.tolerance >= umath.absolute(b)
#............................
class _DomainGreater:
    "DomainGreater(v)(x) = true where x <= v"
    def __init__(self, critical_value):
        "DomainGreater(v)(x) = true where x <= v"
        self.critical_value = critical_value

    def __call__ (self, x):
        "Executes the call behavior."
        return umath.less_equal(x, self.critical_value)
#............................
class _DomainGreaterEqual:
    "DomainGreaterEqual(v)(x) = true where x < v"
    def __init__(self, critical_value):
        "DomainGreaterEqual(v)(x) = true where x < v"
        self.critical_value = critical_value

    def __call__ (self, x):
        "Executes the call behavior."
        return umath.less(x, self.critical_value)

#..............................................................................
class _MaskedUnaryOperation:
    """Defines masked version of unary operations, where invalid
    values are pre-masked.

    Parameters
    ----------
    f : callable
    fill :
        Default filling value (0).
    domain :
        Default domain (None).

    """
    def __init__ (self, mufunc, fill=0, domain=None):
        """ _MaskedUnaryOperation(aufunc, fill=0, domain=None)
            aufunc(fill) must be defined
            self(x) returns aufunc(x)
            with masked values where domain(x) is true or getmask(x) is true.
        """
        self.f = mufunc
        self.fill = fill
        self.domain = domain
        self.__doc__ = getattr(mufunc, "__doc__", str(mufunc))
        self.__name__ = getattr(mufunc, "__name__", str(mufunc))
        ufunc_domain[mufunc] = domain
        ufunc_fills[mufunc] = fill
    #
    def __call__ (self, a, *args, **kwargs):
        "Execute the call behavior."
        #
        m = getmask(a)
        d1 = get_data(a)
        #
        if self.domain is not None:
            dm = narray(self.domain(d1), copy=False)
            m = numpy.logical_or(m, dm)
            # The following two lines control the domain filling methods.
            d1 = d1.copy()
            # We could use smart indexing : d1[dm] = self.fill ...
            # ... but numpy.putmask looks more efficient, despite the copy.
            numpy.putmask(d1, dm, self.fill)
        # Take care of the masked singletong first ...
        if not m.ndim and m:
            return masked
        # Get the result class .......................
        if isinstance(a, MaskedArray):
            subtype = type(a)
        else:
            subtype = MaskedArray
        # Get the result  as a view of the subtype ...
        result = self.f(d1, *args, **kwargs).view(subtype)
        # Fix the mask if we don't have a scalar
        if result.ndim > 0:
            result._mask = m
            result._update_from(a)
        return result
    #
    def __str__ (self):
        return "Masked version of %s. [Invalid values are masked]" % str(self.f)

#..............................................................................
class _MaskedBinaryOperation:
    """Define masked version of binary operations, where invalid
    values are pre-masked.

    Parameters
    ----------
    f : callable
    fillx :
        Default filling value for the first argument (0).
    filly :
        Default filling value for the second argument (0).
    domain :
        Default domain (None).

    """
    def __init__ (self, mbfunc, fillx=0, filly=0):
        """abfunc(fillx, filly) must be defined.
           abfunc(x, filly) = x for all x to enable reduce.
        """
        self.f = mbfunc
        self.fillx = fillx
        self.filly = filly
        self.__doc__ = getattr(mbfunc, "__doc__", str(mbfunc))
        self.__name__ = getattr(mbfunc, "__name__", str(mbfunc))
        ufunc_domain[mbfunc] = None
        ufunc_fills[mbfunc] = (fillx, filly)
    #
    def __call__ (self, a, b, *args, **kwargs):
        "Execute the call behavior."
        m = mask_or(getmask(a), getmask(b))
        (d1, d2) = (get_data(a), get_data(b))
        result = self.f(d1, d2, *args, **kwargs).view(get_masked_subclass(a,b))
        if result.size > 1:
            if m is not nomask:
                result._mask = make_mask_none(result.shape)
                result._mask.flat = m
            if isinstance(a,MaskedArray):
                result._update_from(a)
            if isinstance(b,MaskedArray):
                result._update_from(b)
        elif m:
            return masked
        return result
    #
    def reduce (self, target, axis=0, dtype=None):
        """Reduce `target` along the given `axis`."""
        if isinstance(target, MaskedArray):
            tclass = type(target)
        else:
            tclass = MaskedArray
        m = getmask(target)
        t = filled(target, self.filly)
        if t.shape == ():
            t = t.reshape(1)
            if m is not nomask:
                m = make_mask(m, copy=1)
                m.shape = (1,)
        if m is nomask:
            return self.f.reduce(t, axis).view(tclass)
        t = t.view(tclass)
        t._mask = m
        tr = self.f.reduce(getdata(t), axis, dtype=dtype or t.dtype)
        mr = umath.logical_and.reduce(m, axis)
        tr = tr.view(tclass)
        if mr.ndim > 0:
            tr._mask = mr
            return tr
        elif mr:
            return masked
        return tr

    def outer (self, a, b):
        """Return the function applied to the outer product of a and b.

        """
        ma = getmask(a)
        mb = getmask(b)
        if ma is nomask and mb is nomask:
            m = nomask
        else:
            ma = getmaskarray(a)
            mb = getmaskarray(b)
            m = umath.logical_or.outer(ma, mb)
        if (not m.ndim) and m:
            return masked
        rcls = get_masked_subclass(a,b)
        # We could fill the arguments first, butis it useful ?
        # d = self.f.outer(filled(a, self.fillx), filled(b, self.filly)).view(rcls)
        d = self.f.outer(getdata(a), getdata(b)).view(rcls)
        if d.ndim > 0:
            d._mask = m
        return d

    def accumulate (self, target, axis=0):
        """Accumulate `target` along `axis` after filling with y fill
        value.

        """
        if isinstance(target, MaskedArray):
            tclass = type(target)
        else:
            tclass = masked_array
        t = filled(target, self.filly)
        return self.f.accumulate(t, axis).view(tclass)

    def __str__ (self):
        return "Masked version of " + str(self.f)

#..............................................................................
class _DomainedBinaryOperation:
    """Define binary operations that have a domain, like divide.

    They have no reduce, outer or accumulate.

    Parameters
    ----------
    f : function.
    domain : Default domain.
    fillx : Default filling value for the first argument (0).
    filly : Default filling value for the second argument (0).

    """
    def __init__ (self, dbfunc, domain, fillx=0, filly=0):
        """abfunc(fillx, filly) must be defined.
           abfunc(x, filly) = x for all x to enable reduce.
        """
        self.f = dbfunc
        self.domain = domain
        self.fillx = fillx
        self.filly = filly
        self.__doc__ = getattr(dbfunc, "__doc__", str(dbfunc))
        self.__name__ = getattr(dbfunc, "__name__", str(dbfunc))
        ufunc_domain[dbfunc] = domain
        ufunc_fills[dbfunc] = (fillx, filly)

    def __call__(self, a, b):
        "Execute the call behavior."
        ma = getmask(a)
        mb = getmask(b)
        d1 = getdata(a)
        d2 = get_data(b)
        t = narray(self.domain(d1, d2), copy=False)
        if t.any(None):
            mb = mask_or(mb, t)
            # The following line controls the domain filling
            d2 = numpy.where(t,self.filly,d2)
        m = mask_or(ma, mb)
        if (not m.ndim) and m:
            return masked
        result =  self.f(d1, d2).view(get_masked_subclass(a,b))
        if result.ndim > 0:
            result._mask = m
            if isinstance(a,MaskedArray):
                result._update_from(a)
            if isinstance(b,MaskedArray):
                result._update_from(b)
        return result

    def __str__ (self):
        return "Masked version of " + str(self.f)

#..............................................................................
# Unary ufuncs
exp = _MaskedUnaryOperation(umath.exp)
conjugate = _MaskedUnaryOperation(umath.conjugate)
sin = _MaskedUnaryOperation(umath.sin)
cos = _MaskedUnaryOperation(umath.cos)
tan = _MaskedUnaryOperation(umath.tan)
arctan = _MaskedUnaryOperation(umath.arctan)
arcsinh = _MaskedUnaryOperation(umath.arcsinh)
sinh = _MaskedUnaryOperation(umath.sinh)
cosh = _MaskedUnaryOperation(umath.cosh)
tanh = _MaskedUnaryOperation(umath.tanh)
abs = absolute = _MaskedUnaryOperation(umath.absolute)
fabs = _MaskedUnaryOperation(umath.fabs)
negative = _MaskedUnaryOperation(umath.negative)
floor = _MaskedUnaryOperation(umath.floor)
ceil = _MaskedUnaryOperation(umath.ceil)
around = _MaskedUnaryOperation(fromnumeric.round_)
logical_not = _MaskedUnaryOperation(umath.logical_not)
# Domained unary ufuncs .......................................................
sqrt = _MaskedUnaryOperation(umath.sqrt, 0.0,
                             _DomainGreaterEqual(0.0))
log = _MaskedUnaryOperation(umath.log, 1.0,
                            _DomainGreater(0.0))
log10 = _MaskedUnaryOperation(umath.log10, 1.0,
                              _DomainGreater(0.0))
tan = _MaskedUnaryOperation(umath.tan, 0.0,
                            _DomainTan(1.e-35))
arcsin = _MaskedUnaryOperation(umath.arcsin, 0.0,
                               _DomainCheckInterval(-1.0, 1.0))
arccos = _MaskedUnaryOperation(umath.arccos, 0.0,
                               _DomainCheckInterval(-1.0, 1.0))
arccosh = _MaskedUnaryOperation(umath.arccosh, 1.0,
                                _DomainGreaterEqual(1.0))
arctanh = _MaskedUnaryOperation(umath.arctanh, 0.0,
                                _DomainCheckInterval(-1.0+1e-15, 1.0-1e-15))
# Binary ufuncs ...............................................................
add = _MaskedBinaryOperation(umath.add)
subtract = _MaskedBinaryOperation(umath.subtract)
multiply = _MaskedBinaryOperation(umath.multiply, 1, 1)
arctan2 = _MaskedBinaryOperation(umath.arctan2, 0.0, 1.0)
equal = _MaskedBinaryOperation(umath.equal)
equal.reduce = None
not_equal = _MaskedBinaryOperation(umath.not_equal)
not_equal.reduce = None
less_equal = _MaskedBinaryOperation(umath.less_equal)
less_equal.reduce = None
greater_equal = _MaskedBinaryOperation(umath.greater_equal)
greater_equal.reduce = None
less = _MaskedBinaryOperation(umath.less)
less.reduce = None
greater = _MaskedBinaryOperation(umath.greater)
greater.reduce = None
logical_and = _MaskedBinaryOperation(umath.logical_and)
alltrue = _MaskedBinaryOperation(umath.logical_and, 1, 1).reduce
logical_or = _MaskedBinaryOperation(umath.logical_or)
sometrue = logical_or.reduce
logical_xor = _MaskedBinaryOperation(umath.logical_xor)
bitwise_and = _MaskedBinaryOperation(umath.bitwise_and)
bitwise_or = _MaskedBinaryOperation(umath.bitwise_or)
bitwise_xor = _MaskedBinaryOperation(umath.bitwise_xor)
hypot = _MaskedBinaryOperation(umath.hypot)
# Domained binary ufuncs ......................................................
divide = _DomainedBinaryOperation(umath.divide, _DomainSafeDivide(), 0, 1)
true_divide = _DomainedBinaryOperation(umath.true_divide,
                                        _DomainSafeDivide(), 0, 1)
floor_divide = _DomainedBinaryOperation(umath.floor_divide,
                                         _DomainSafeDivide(), 0, 1)
remainder = _DomainedBinaryOperation(umath.remainder,
                                      _DomainSafeDivide(), 0, 1)
fmod = _DomainedBinaryOperation(umath.fmod, _DomainSafeDivide(), 0, 1)


#####--------------------------------------------------------------------------
#---- --- Mask creation functions ---
#####--------------------------------------------------------------------------
def get_mask(a):
    """Return the mask of a, if any, or nomask.

    To get a full array of booleans of the same shape as a, use
    getmaskarray.

    """
    return getattr(a, '_mask', nomask)
getmask = get_mask

def getmaskarray(a):
    """Return the mask of a, if any, or a boolean array of the shape
    of a, full of False.

    """
    m = getmask(a)
    if m is nomask:
        m = make_mask_none(fromnumeric.shape(a))
    return m

def is_mask(m):
    """Return True if m is a legal mask.

    Does not check contents, only type.

    """
    try:
        return m.dtype.type is MaskType
    except AttributeError:
        return False
#
def make_mask(m, copy=False, shrink=True, flag=None):
    """Return m as a mask, creating a copy if necessary or requested.

    The function can accept any sequence of integers or nomask.  Does
    not check that contents must be 0s and 1s.

    Parameters
    ----------
    m : array_like
        Potential mask.
    copy : bool
        Whether to return a copy of m (True) or m itself (False).
    shrink : bool
        Whether to shrink m to nomask if all its values are False.

    """
    if flag is not None:
        warnings.warn("The flag 'flag' is now called 'shrink'!",
                      DeprecationWarning)
        shrink = flag
    if m is nomask:
        return nomask
    elif isinstance(m, ndarray):
        m = filled(m, True)
        if m.dtype.type is MaskType:
            if copy:
                result = narray(m, dtype=MaskType, copy=copy)
            else:
                result = m
        else:
            result = narray(m, dtype=MaskType)
    else:
        result = narray(filled(m, True), dtype=MaskType)
    # Bas les masques !
    if shrink and not result.any():
        return nomask
    else:
        return result

def make_mask_none(s):
    """Return a mask of shape s, filled with False.

    Parameters
    ----------
    s : tuple
        A tuple indicating the shape of the final mask.

    """
    result = numeric.zeros(s, dtype=MaskType)
    return result

def mask_or (m1, m2, copy=False, shrink=True):
    """Return the combination of two masks m1 and m2.

    The masks are combined with the *logical_or* operator, treating
    nomask as False.  The result may equal m1 or m2 if the other is
    nomask.

    Parameters
    ----------
    m1 : array_like
        First mask.
    m2 : array_like
        Second mask
    copy : bool
        Whether to return a copy.
    shrink : bool
        Whether to shrink m to nomask if all its values are False.

     """
    if m1 is nomask:
        return make_mask(m2, copy=copy, shrink=shrink)
    if m2 is nomask:
        return make_mask(m1, copy=copy, shrink=shrink)
    if m1 is m2 and is_mask(m1):
        return m1
    return make_mask(umath.logical_or(m1, m2), copy=copy, shrink=shrink)

#####--------------------------------------------------------------------------
#--- --- Masking functions ---
#####--------------------------------------------------------------------------
def masked_where(condition, a, copy=True):
    """Return a as an array masked where condition is true.

    Masked values of a or condition are kept.

    Parameters
    ----------
    condition : array_like
        Masking condition.
    a : array_like
        Array to mask.
    copy : bool
        Whether to return a copy of a (True) or modify a in place.

    """
    cond = make_mask(condition)
    a = narray(a, copy=copy, subok=True)
    if hasattr(a, '_mask'):
        cond = mask_or(cond, a._mask)
        cls = type(a)
    else:
        cls = MaskedArray
    result = a.view(cls)
    result._mask = cond
    return result

def masked_greater(x, value, copy=True):
    "Shortcut to masked_where, with condition = (x > value)."
    return masked_where(greater(x, value), x, copy=copy)

def masked_greater_equal(x, value, copy=True):
    "Shortcut to masked_where, with condition = (x >= value)."
    return masked_where(greater_equal(x, value), x, copy=copy)

def masked_less(x, value, copy=True):
    "Shortcut to masked_where, with condition = (x < value)."
    return masked_where(less(x, value), x, copy=copy)

def masked_less_equal(x, value, copy=True):
    "Shortcut to masked_where, with condition = (x <= value)."
    return masked_where(less_equal(x, value), x, copy=copy)

def masked_not_equal(x, value, copy=True):
    "Shortcut to masked_where, with condition = (x != value)."
    return masked_where((x != value), x, copy=copy)

#
def masked_equal(x, value, copy=True):
    """Shortcut to masked_where, with condition = (x == value).  For
    floating point, consider `masked_values(x, value)` instead.

    """
    # An alternative implementation relies on filling first: probably not needed.
    # d = filled(x, 0)
    # c = umath.equal(d, value)
    # m = mask_or(c, getmask(x))
    # return array(d, mask=m, copy=copy)
    return masked_where((x == value), x, copy=copy)

def masked_inside(x, v1, v2, copy=True):
    """Shortcut to masked_where, where condition is True for x inside
    the interval [v1,v2] (v1 <= x <= v2).  The boundaries v1 and v2
    can be given in either order.

    Notes
    -----
    The array x is prefilled with its filling value.

    """
    if v2 < v1:
        (v1, v2) = (v2, v1)
    xf = filled(x)
    condition = (xf >= v1) & (xf <= v2)
    return masked_where(condition, x, copy=copy)

def masked_outside(x, v1, v2, copy=True):
    """Shortcut to masked_where, where condition is True for x outside
    the interval [v1,v2] (x < v1)|(x > v2).  The boundaries v1 and v2
    can be given in either order.

    Notes
    -----
    The array x is prefilled with its filling value.

    """
    if v2 < v1:
        (v1, v2) = (v2, v1)
    xf = filled(x)
    condition = (xf < v1) | (xf > v2)
    return masked_where(condition, x, copy=copy)

#
def masked_object(x, value, copy=True):
    """Mask the array x where the data are exactly equal to value.

    This function is suitable only for object arrays: for floating
    point, please use ``masked_values`` instead.

    Notes
    -----
    The mask is set to `nomask` if posible.

    """
    if isMaskedArray(x):
        condition = umath.equal(x._data, value)
        mask = x._mask
    else:
        condition = umath.equal(fromnumeric.asarray(x), value)
        mask = nomask
    mask = mask_or(mask, make_mask(condition, shrink=True))
    return masked_array(x, mask=mask, copy=copy, fill_value=value)

def masked_values(x, value, rtol=1.e-5, atol=1.e-8, copy=True):
    """Mask the array x where the data are approximately equal in
    value, i.e.

    (abs(x - value) <= atol+rtol*abs(value))

    Suitable only for floating points. For integers, please use
    ``masked_equal``.  The mask is set to nomask if posible.

    Parameters
    ----------
    x : array_like
        Array to fill.
    value : float
        Masking value.
    rtol : float
        Tolerance parameter.
    atol : float
        Tolerance parameter (1e-8).
    copy : bool
        Whether to return a copy of x.

    """
    abs = umath.absolute
    xnew = filled(x, value)
    if issubclass(xnew.dtype.type, numeric.floating):
        condition = umath.less_equal(abs(xnew-value), atol+rtol*abs(value))
        mask = getattr(x, '_mask', nomask)
    else:
        condition = umath.equal(xnew, value)
        mask = nomask
    mask = mask_or(mask, make_mask(condition, shrink=True))
    return masked_array(xnew, mask=mask, copy=copy, fill_value=value)

def masked_invalid(a, copy=True):
    """Mask the array for invalid values (nans or infs).  Any
    preexisting mask is conserved.

    """
    a = narray(a, copy=copy, subok=True)
    condition = ~(numpy.isfinite(a))
    if hasattr(a, '_mask'):
        condition = mask_or(condition, a._mask)
        cls = type(a)
    else:
        cls = MaskedArray
    result = a.view(cls)
    result._mask = condition
    return result


#####--------------------------------------------------------------------------
#---- --- Printing options ---
#####--------------------------------------------------------------------------
class _MaskedPrintOption:
    """Handle the string used to represent missing data in a masked
    array.

    """
    def __init__ (self, display):
        "Create the masked_print_option object."
        self._display = display
        self._enabled = True

    def display(self):
        "Display the string to print for masked values."
        return self._display

    def set_display (self, s):
        "Set the string to print for masked values."
        self._display = s

    def enabled(self):
        "Is the use of the display value enabled?"
        return self._enabled

    def enable(self, shrink=1):
        "Set the enabling shrink to `shrink`."
        self._enabled = shrink

    def __str__ (self):
        return str(self._display)

    __repr__ = __str__

#if you single index into a masked location you get this object.
masked_print_option = _MaskedPrintOption('--')

#####--------------------------------------------------------------------------
#---- --- MaskedArray class ---
#####--------------------------------------------------------------------------

#...............................................................................
class _arraymethod(object):
    """Define a wrapper for basic array methods.

    Upon call, returns a masked array, where the new _data array is
    the output of the corresponding method called on the original
    _data.

    If onmask is True, the new mask is the output of the method called
    on the initial mask. Otherwise, the new mask is just a reference
    to the initial mask.

    Parameters
    ----------
    _name : String
        Name of the function to apply on data.
    _onmask : bool
        Whether the mask must be processed also (True) or left
        alone (False). Default: True.
    obj : Object
        The object calling the arraymethod.

    """
    def __init__(self, funcname, onmask=True):
        self._name = funcname
        self._onmask = onmask
        self.obj = None
        self.__doc__ = self.getdoc()
    #
    def getdoc(self):
        "Return the doc of the function (from the doc of the method)."
        methdoc = getattr(ndarray, self._name, None)
        methdoc = getattr(numpy, self._name, methdoc)
        if methdoc is not None:
            return methdoc.__doc__
    #
    def __get__(self, obj, objtype=None):
        self.obj = obj
        return self
    #
    def __call__(self, *args, **params):
        methodname = self._name
        data = self.obj._data
        mask = self.obj._mask
        cls = type(self.obj)
        result = getattr(data, methodname)(*args, **params).view(cls)
        result._update_from(self.obj)
        if result.ndim:
            if not self._onmask:
                result.__setmask__(mask)
            elif mask is not nomask:
                result.__setmask__(getattr(mask, methodname)(*args, **params))
        else:
            if mask.ndim and mask.all():
                return masked
        return result
#..........................................................

class FlatIter(object):
    "Define an interator."
    def __init__(self, ma):
        self.ma = ma
        self.ma_iter = numpy.asarray(ma).flat

        if ma._mask is nomask:
            self.maskiter = None
        else:
            self.maskiter = ma._mask.flat

    def __iter__(self):
        return self

    ### This won't work is ravel makes a copy
    def __setitem__(self, index, value):
        a = self.ma.ravel()
        a[index] = value

    def next(self):
        d = self.ma_iter.next()
        if self.maskiter is not None and self.maskiter.next():
            d = masked
        return d


class MaskedArray(numeric.ndarray):
    """Arrays with possibly masked values.  Masked values of True
    exclude the corresponding element from any computation.

    Construction:
        x = MaskedArray(data, mask=nomask, dtype=None, copy=True,
        fill_value=None, keep_mask=True, hard_mask=False, shrink=True)

    Parameters
    ----------
    data : {var}
        Input data.
    mask : {nomask, sequence}
        Mask.  Must be convertible to an array of booleans with
        the same shape as data: True indicates a masked (eg.,
        invalid) data.
    dtype : dtype
        Data type of the output. If None, the type of the data
        argument is used.  If dtype is not None and different from
        data.dtype, a copy is performed.
    copy : bool
        Whether to copy the input data (True), or to use a
        reference instead.  Note: data are NOT copied by default.
    subok : {True, boolean}
        Whether to return a subclass of MaskedArray (if possible)
        or a plain MaskedArray.
    ndmin : {0, int}
        Minimum number of dimensions
    fill_value : {var}
        Value used to fill in the masked values when necessary. If
        None, a default based on the datatype is used.
    keep_mask : {True, boolean}
        Whether to combine mask with the mask of the input data,
        if any (True), or to use only mask for the output (False).
    hard_mask : {False, boolean}
        Whether to use a hard mask or not. With a hard mask,
        masked values cannot be unmasked.
    shrink : {True, boolean}
        Whether to force compression of an empty mask.

    """

    __array_priority__ = 15
    _defaultmask = nomask
    _defaulthardmask = False
    _baseclass =  numeric.ndarray

    def __new__(cls, data=None, mask=nomask, dtype=None, copy=False,
                subok=True, ndmin=0, fill_value=None,
                keep_mask=True, hard_mask=False, flag=None,shrink=True,
                **options):
        """Create a new masked array from scratch.

        Note: you can also create an array with the .view(MaskedArray)
        method.

        """
        if flag is not None:
            warnings.warn("The flag 'flag' is now called 'shrink'!",
                          DeprecationWarning)
            shrink = flag
        # Process data............
        _data = narray(data, dtype=dtype, copy=copy, subok=True, ndmin=ndmin)
        _baseclass = getattr(data, '_baseclass', type(_data))
        _basedict = getattr(data, '_basedict', getattr(data, '__dict__', {}))
        if not isinstance(data, MaskedArray) or not subok:
            _data = _data.view(cls)
        else:
            _data = _data.view(type(data))
        # Backwards compatibility w/ numpy.core.ma .......
        if hasattr(data,'_mask') and not isinstance(data, ndarray):
            _data._mask = data._mask
            _sharedmask = True
        # Process mask ...........
        if mask is nomask:
            if not keep_mask:
                if shrink:
                    _data._mask = nomask
                else:
                    _data._mask = make_mask_none(_data)
            if copy:
                _data._mask = _data._mask.copy()
                _data._sharedmask = False
            else:
                _data._sharedmask = True
        else:
            mask = narray(mask, dtype=MaskType, copy=copy)
            if mask.shape != _data.shape:
                (nd, nm) = (_data.size, mask.size)
                if nm == 1:
                    mask = numeric.resize(mask, _data.shape)
                elif nm == nd:
                    mask = fromnumeric.reshape(mask, _data.shape)
                else:
                    msg = "Mask and data not compatible: data size is %i, "+\
                          "mask size is %i."
                    raise MAError, msg % (nd, nm)
                copy = True
            if _data._mask is nomask:
                _data._mask = mask
                _data._sharedmask = not copy
            else:
                if not keep_mask:
                    _data._mask = mask
                    _data._sharedmask = not copy
                else:
                    _data._mask = umath.logical_or(mask, _data._mask)
                    _data._sharedmask = False

        # Update fill_value.......
        _data._fill_value = _check_fill_value(fill_value, _data.dtype)
        # Process extra options ..
        _data._hardmask = hard_mask
        _data._baseclass = _baseclass
        _data._basedict = _basedict
        return _data
    #
    def _update_from(self, obj):
        """Copies some attributes of obj to self.
        """  
        if obj is not None and isinstance(obj,ndarray):
            _baseclass = type(obj)
        else:
            _baseclass = ndarray
        _basedict = getattr(obj,'_basedict',getattr(obj,'__dict__',{}))
        _dict = dict(_fill_value=getattr(obj, '_fill_value', None),
                     _hardmask=getattr(obj, '_hardmask', False),
                     _sharedmask=getattr(obj, '_sharedmask', False),
                     _baseclass=getattr(obj,'_baseclass',_baseclass),
                     _basedict=_basedict,)
        self.__dict__.update(_dict)
        self.__dict__.update(_basedict)        
        return
    #........................
    def __array_finalize__(self,obj):
        """Finalizes the masked array.
        """
        # Get main attributes .........
        self._update_from(obj)
        self._mask = getattr(obj, '_mask', nomask)
        # Finalize the mask ...........
        if self._mask is not nomask:
            self._mask.shape = self.shape
        return
    #..................................
    def __array_wrap__(self, obj, context=None):
        """Special hook for ufuncs.
        Wraps the numpy array and sets the mask according to context.
        """
        result = obj.view(type(self))
        result._update_from(self)
        #..........
        if context is not None:
            result._mask = result._mask.copy()
            (func, args, _) = context
            m = reduce(mask_or, [getmaskarray(arg) for arg in args])
            # Get the domain mask................
            domain = ufunc_domain.get(func, None)
            if domain is not None:
                if len(args) > 2:
                    d = reduce(domain, args)
                else:
                    d = domain(*args)
                # Fill the result where the domain is wrong
                try:
                    # Binary domain: take the last value
                    fill_value = ufunc_fills[func][-1]
                except TypeError:
                    # Unary domain: just use this one
                    fill_value = ufunc_fills[func]
                except KeyError:
                    # Domain not recognized, use fill_value instead
                    fill_value = self.fill_value
                result = result.copy()
                numpy.putmask(result, d, fill_value)
                # Update the mask
                if m is nomask:
                    if d is not nomask:
                        m = d
                else:
                    m |= d
            # Make sure the mask has the proper size
            if result.shape == () and m:
                return masked
            else:
                result._mask = m
                result._sharedmask = False
        #....
        return result
    #.............................................
    def __getitem__(self, indx):
        """x.__getitem__(y) <==> x[y]

        Return the item described by i, as a masked array.

        """
        # This test is useful, but we should keep things light...
#        if getmask(indx) is not nomask:
#            msg = "Masked arrays must be filled before they can be used as indices!"
#            raise IndexError, msg
        dout = ndarray.__getitem__(self.view(ndarray), indx)
        # We could directly use ndarray.__getitem__ on self...
        # But then we would have to modify __array_finalize__ to prevent the
        # mask of being reshaped if it hasn't been set up properly yet...
        # So it's easier to stick to the current version
        m = self._mask
        if not getattr(dout,'ndim', False):
            # Just a scalar............
            if m is not nomask and m[indx]:
                return masked
        else:
            # Force dout to MA ........
            dout = dout.view(type(self))
            # Inherit attributes from self
            dout._update_from(self)
            # Check the fill_value ....
            if isinstance(indx, basestring):
                fvindx = list(self.dtype.names).index(indx)
                dout._fill_value = self.fill_value[fvindx]
            # Update the mask if needed
            if m is not nomask:
                if isinstance(indx, basestring):
                    dout._mask = m.reshape(dout.shape)
                else:
                    dout._mask = ndarray.__getitem__(m, indx).reshape(dout.shape)
#               Note: Don't try to check for m.any(), that'll take too long...
#                mask = ndarray.__getitem__(m, indx).reshape(dout.shape)
#                if self._shrinkmask and not m.any():
#                    dout._mask = nomask
#                else:
#                    dout._mask = mask
        return dout
    #........................
    def __setitem__(self, indx, value):
        """x.__setitem__(i, y) <==> x[i]=y

        Set item described by index. If value is masked, masks those
        locations.

        """
        if self is masked:
            raise MAError, 'Cannot alter the masked element.'
        # This test is useful, but we should keep things light...
#        if getmask(indx) is not nomask:
#            msg = "Masked arrays must be filled before they can be used as indices!"
#            raise IndexError, msg
        if isinstance(indx, basestring):
            ndarray.__setitem__(self._data,indx, getdata(value))
            warnings.warn("MaskedArray.__setitem__ on fields: "\
                          "The mask is NOT affected!")
            return
        #....
        if value is masked:
            m = self._mask
            if m is nomask:
                m = numpy.zeros(self.shape, dtype=MaskType)
            m[indx] = True
            self._mask = m
            self._sharedmask = False
            return
        #....
        dval = narray(value, copy=False, dtype=self.dtype)
        valmask = getmask(value)        
        if self._mask is nomask:
            # Set the data, then the mask
            ndarray.__setitem__(self._data,indx,dval)
            if valmask is not nomask:
                self._mask = numpy.zeros(self.shape, dtype=MaskType)
                self._mask[indx] = valmask
        elif not self._hardmask:
            # Unshare the mask if necessary to avoid propagation
            self.unshare_mask()
            # Set the data, then the mask
            ndarray.__setitem__(self._data,indx,dval)
            self._mask[indx] = valmask
        elif hasattr(indx, 'dtype') and (indx.dtype==bool_):
            indx = indx * umath.logical_not(self._mask)
            ndarray.__setitem__(self._data,indx,dval)
        else:
            mindx = mask_or(self._mask[indx], valmask, copy=True)
            dindx = self._data[indx]
            if dindx.size > 1:
                dindx[~mindx] = dval
            elif mindx is nomask:
                dindx = dval
            ndarray.__setitem__(self._data,indx,dindx)
            self._mask[indx] = mindx
    #............................................
    def __getslice__(self, i, j):
        """x.__getslice__(i, j) <==> x[i:j]

        Return the slice described by (i, j).  The use of negative
        indices is not supported.

        """
        return self.__getitem__(slice(i,j))
    #........................
    def __setslice__(self, i, j, value):
        """x.__setslice__(i, j, value) <==> x[i:j]=value

        Set the slice (i,j) of a to value. If value is masked, mask
        those locations.

        """
        self.__setitem__(slice(i,j), value)
    #............................................
    def __setmask__(self, mask, copy=False):
        """Set the mask.

        """
        if mask is not nomask:
            mask = narray(mask, copy=copy, dtype=MaskType)
            # We could try to check whether shrinking is needed..
            # ... but we would waste some precious time
#            if self._shrinkmask and not mask.any():
#                mask = nomask
        if self._mask is nomask:
            self._mask = mask
        elif self._hardmask:
            if mask is not nomask:
                self._mask.__ior__(mask)
        else:
            # This one is tricky: if we set the mask that way, we may break the
            # propagation. But if we don't, we end up with a mask full of False
            # and a test on nomask fails...
            if mask is nomask:
                self._mask = nomask
            else:
                self.unshare_mask()
                self._mask.flat = mask
        if self._mask.shape:
            self._mask = numeric.reshape(self._mask, self.shape)
    _set_mask = __setmask__
    #....
    def _get_mask(self):
        """Return the current mask.

        """
        # We could try to force a reshape, but that wouldn't work in some cases.
#        return self._mask.reshape(self.shape)
        return self._mask
    mask = property(fget=_get_mask, fset=__setmask__, doc="Mask")
    #............................................
    def harden_mask(self):
        """Force the mask to hard.

        """
        self._hardmask = True

    def soften_mask(self):
        """Force the mask to soft.

        """
        self._hardmask = False

    def unshare_mask(self):
        """Copy the mask and set the sharedmask flag to False.

        """
        if self._sharedmask:
            self._mask = self._mask.copy()
            self._sharedmask = False

    def shrink_mask(self):
        """Reduce a mask to nomask when possible.

        """
        m = self._mask
        if m.ndim and not m.any():
            self._mask = nomask

    #............................................
    def _get_data(self):
        """Return the current data, as a view of the original
        underlying data.

        """
        return self.view(self._baseclass)
    _data = property(fget=_get_data)
    data = property(fget=_get_data)

    def raw_data(self):
        """Return the _data part of the MaskedArray.

        DEPRECATED: You should really use ``.data`` instead...

        """
        warnings.warn('Use .data instead.', DeprecationWarning)
        return self._data
    #............................................
    def _get_flat(self):
        """Return a flat iterator.

        """
        return FlatIter(self)
    #
    def _set_flat (self, value):
        """Set a flattened version of self to value.

        """
        y = self.ravel()
        y[:] = value
    #
    flat = property(fget=_get_flat, fset=_set_flat,
                    doc="Flat version of the array.")
    #............................................
    def get_fill_value(self):
        """Return the filling value.

        """
        if self._fill_value is None:
            self._fill_value = _check_fill_value(None, self.dtype)
        return self._fill_value

    def set_fill_value(self, value=None):
        """Set the filling value to value.

        If value is None, use a default based on the data type.

        """
        self._fill_value = _check_fill_value(value,self.dtype)

    fill_value = property(fget=get_fill_value, fset=set_fill_value,
                          doc="Filling value.")

    def filled(self, fill_value=None):
        """Return a copy of self._data, where masked values are filled
        with fill_value.

        If fill_value is None, self.fill_value is used instead.

        Notes
        -----
        + Subclassing is preserved
        + The result is NOT a MaskedArray !

        Examples
        --------
        >>> x = array([1,2,3,4,5], mask=[0,0,1,0,1], fill_value=-999)
        >>> x.filled()
        array([1,2,-999,4,-999])
        >>> type(x.filled())
        <type 'numpy.ndarray'>

        """
        m = self._mask
        if m is nomask or not m.any():
            return self._data
        #
        if fill_value is None:
            fill_value = self.fill_value
        #
        if self is masked_singleton:
            result = numeric.asanyarray(fill_value)
        else:
            result = self._data.copy()
            try:
                numpy.putmask(result, m, fill_value)
            except (TypeError, AttributeError):
                fill_value = narray(fill_value, dtype=object)
                d = result.astype(object)
                result = fromnumeric.choose(m, (d, fill_value))
            except IndexError:
                #ok, if scalar
                if self._data.shape:
                    raise
                elif m:
                    result = narray(fill_value, dtype=self.dtype)
                else:
                    result = self._data
        return result

    def compressed(self):
        """Return a 1-D array of all the non-masked data.

        """
        data = ndarray.ravel(self._data)
        if self._mask is not nomask:
            data = data.compress(numpy.logical_not(ndarray.ravel(self._mask)))
        return data


    def compress(self, condition, axis=None, out=None):
        """Return a where condition is True.
        If condition is a MaskedArray, missing values are considered as False.

        Returns
        -------
        A MaskedArray object.

        Notes
        -----
        Please note the difference with compressed() !
        The output of compress has a mask, the output of compressed does not.

        """
        # Get the basic components
        (_data, _mask) = (self._data, self._mask)
        # Force the condition to a regular ndarray (forget the missing values...)
        condition = narray(condition, copy=False, subok=False)
        #
        _new = _data.compress(condition, axis=axis, out=out).view(type(self))
        _new._update_from(self)
        if _mask is not nomask:
            _new._mask = _mask.compress(condition, axis=axis)
        return _new

    #............................................
    def __str__(self):
        """String representation.

        """
        if masked_print_option.enabled():
            f = masked_print_option
            if self is masked:
                return str(f)
            m = self._mask
            if m is nomask:
                res = self._data
            else:
                if m.shape == ():
                    if m:
                        return str(f)
                    else:
                        return str(self._data)
                # convert to object array to make filled work
#CHECK: the two lines below seem more robust than the self._data.astype
#                res = numeric.empty(self._data.shape, object_)
#                numeric.putmask(res,~m,self._data)
                res = self._data.astype("|O8")
                res[m] = f
        else:
            res = self.filled(self.fill_value)
        return str(res)

    def __repr__(self):
        """Literal string representation.

        """
        with_mask = """\
masked_%(name)s(data =
 %(data)s,
      mask =
 %(mask)s,
      fill_value=%(fill)s)
"""
        with_mask1 = """\
masked_%(name)s(data = %(data)s,
      mask = %(mask)s,
      fill_value=%(fill)s)
"""
        n = len(self.shape)
        name = repr(self._data).split('(')[0]
        if n <= 1:
            return with_mask1 % {
                'name': name,
                'data': str(self),
                'mask': str(self._mask),
                'fill': str(self.fill_value),
                }
        return with_mask % {
            'name': name,
            'data': str(self),
            'mask': str(self._mask),
            'fill': str(self.fill_value),
            }
    #............................................
    def __add__(self, other):
        "Add other to self, and return a new masked array."
        return add(self, other)
    #
    def __sub__(self, other):
        "Subtract other to self, and return a new masked array."
        return subtract(self, other)
    #
    def __mul__(self, other):
        "Multiply other by self, and return a new masked array."
        return multiply(self, other)
    #
    def __div__(self, other):
        "Divide other into self, and return a new masked array."
        return divide(self, other)
    #
    def __truediv__(self, other):
        "Divide other into self, and return a new masked array."
        return true_divide(self, other)
    #
    def __floordiv__(self, other):
        "Divide other into self, and return a new masked array."
        return floor_divide(self, other)
    #
    def __pow__(self, other):
        "Raise self to the power other, masking the potential NaNs/Infs"
        return power(self, other)
    #............................................
    def __iadd__(self, other):
        "Add other to self in-place."
        ndarray.__iadd__(self._data, getdata(other))
        m = getmask(other)
        if self._mask is nomask:
            self._mask = m
        elif m is not nomask:
            self._mask += m
        return self
    #....
    def __isub__(self, other):
        "Subtract other from self in-place."
        ndarray.__isub__(self._data, getdata(other))
        m = getmask(other)
        if self._mask is nomask:
            self._mask = m
        elif m is not nomask:
            self._mask += m
        return self
    #....
    def __imul__(self, other):
        "Multiply self by other in-place."
        ndarray.__imul__(self._data, getdata(other))
        m = getmask(other)
        if self._mask is nomask:
            self._mask = m
        elif m is not nomask:
            self._mask += m
        return self
    #....
    def __idiv__(self, other):
        "Divide self by other in-place."
        other_data = getdata(other)
        dom_mask = _DomainSafeDivide().__call__(self._data, other_data)
        other_mask = getmask(other)
        new_mask = mask_or(other_mask, dom_mask)
        # The following 3 lines control the domain filling
        if dom_mask.any():
            other_data = other_data.copy()
            numpy.putmask(other_data, dom_mask, 1)
        ndarray.__idiv__(self._data, other_data)
        self._mask = mask_or(self._mask, new_mask)
        return self
    #...
    def __ipow__(self, other):
        "Raise self to the power other, in place"
        _data = self._data
        other_data = getdata(other)
        other_mask = getmask(other)
        ndarray.__ipow__(_data, other_data)
        invalid = numpy.logical_not(numpy.isfinite(_data))
        new_mask = mask_or(other_mask,invalid)
        self._mask = mask_or(self._mask, new_mask)
        # The following line is potentially problematic, as we change _data...
        numpy.putmask(self._data,invalid,self.fill_value)
        return self
    #............................................
    def __float__(self):
        "Convert to float."
        if self.size > 1:
            raise TypeError,\
                   "Only length-1 arrays can be converted to Python scalars"
        elif self._mask:
            warnings.warn("Warning: converting a masked element to nan.")
            return numpy.nan
        return float(self.item())

    def __int__(self):
        "Convert to int."
        if self.size > 1:
            raise TypeError,\
                   "Only length-1 arrays can be converted to Python scalars"
        elif self._mask:
            raise MAError, 'Cannot convert masked element to a Python int.'
        return int(self.item())
    #............................................
    def get_imag(self):
        result = self._data.imag.view(type(self))
        result.__setmask__(self._mask)
        return result
    imag = property(fget=get_imag,doc="Imaginary part")

    def get_real(self):
        result = self._data.real.view(type(self))
        result.__setmask__(self._mask)
        return result
    real = property(fget=get_real,doc="Real part")


    #............................................
    def count(self, axis=None):
        """Count the non-masked elements of the array along the given
        axis.

        Parameters
        ----------
        axis : int, optional
            Axis along which to count the non-masked elements. If
            not given, all the non masked elements are counted.

        Returns
        -------
        A masked array where the mask is True where all data are
        masked.  If axis is None, returns either a scalar ot the
        masked singleton if all values are masked.

        """
        m = self._mask
        s = self.shape
        ls = len(s)
        if m is nomask:
            if ls == 0:
                return 1
            if ls == 1:
                return s[0]
            if axis is None:
                return self.size
            else:
                n = s[axis]
                t = list(s)
                del t[axis]
                return numeric.ones(t) * n
        n1 = numpy.size(m, axis)
        n2 = m.astype(int_).sum(axis)
        if axis is None:
            return (n1-n2)
        else:
            return narray(n1 - n2)
    #............................................
    flatten = _arraymethod('flatten')
    #
    def ravel(self):
        """Returns a 1D version of self, as a view."""
        r = ndarray.ravel(self._data).view(type(self))
        r._update_from(self)
        if self._mask is not nomask:
            r._mask = ndarray.ravel(self._mask).reshape(r.shape)
        else:
            r._mask = nomask
        return r
    #
    repeat = _arraymethod('repeat')
    #
    def reshape (self, *s):
        """Reshape the array to shape s.

        Returns
        -------
        A new masked array.

        Notes
        -----
        If you want to modify the shape in place, please use
        ``a.shape = s``

        """
        result = self._data.reshape(*s).view(type(self))
        result.__dict__.update(self.__dict__)
        if result._mask is not nomask:
            result._mask = self._mask.copy()
            result._mask.shape = result.shape
        return result
    #
    def resize(self, newshape, refcheck=True, order=False):
        """Attempt to modify the size and the shape of the array in place.

        The array must own its own memory and not be referenced by
        other arrays.

        Returns
        -------
        None.

        """
        try:
            self._data.resize(newshape, refcheck, order)
            if self.mask is not nomask:
                self._mask.resize(newshape, refcheck, order)
        except ValueError:
            raise ValueError("Cannot resize an array that has been referenced "
                             "or is referencing another array in this way.\n"
                             "Use the resize function.")
        return None
    #
    def put(self, indices, values, mode='raise'):
        """Set storage-indexed locations to corresponding values.

        a.put(values, indices, mode) sets a.flat[n] = values[n] for
        each n in indices.  If ``values`` is shorter than ``indices``
        then it will repeat.  If ``values`` has some masked values, the
        initial mask is updated in consequence, else the corresponding
        values are unmasked.

        """
        m = self._mask
        # Hard mask: Get rid of the values/indices that fall on masked data
        if self._hardmask and self._mask is not nomask:
            mask = self._mask[indices]
            indices = narray(indices, copy=False)
            values = narray(values, copy=False, subok=True)
            values.resize(indices.shape)
            indices = indices[~mask]
            values = values[~mask]
        #....
        self._data.put(indices, values, mode=mode)
        #....
        if m is nomask:
            m = getmask(values)
        else:
            m = m.copy()
            if getmask(values) is nomask:
                m.put(indices, False, mode=mode)
            else:
                m.put(indices, values._mask, mode=mode)
            m = make_mask(m, copy=False, shrink=True)
        self._mask = m
    #............................................
    def ids (self):
        """Return the addresses of the data and mask areas."""
        if self._mask is nomask:
            return (self.ctypes.data, id(nomask))
        return (self.ctypes.data, self._mask.ctypes.data)
    #............................................
    def all(self, axis=None, out=None):
        """Return True if all entries along the given axis are True,
        False otherwise.  Masked values are considered as True during
        computation.

        Parameter
        ----------
            axis : int, optional
                Axis along which the operation is performed.  If None,
                the operation is performed on a flatten array
            out : {MaskedArray}, optional
                Alternate optional output.  If not None, out should be
                a valid MaskedArray of the same shape as the output of
                self._data.all(axis).

        Returns            A masked array, where the mask is True if all data along
        -------
        the axis are masked.

        Notes
        -----
        An exception is raised if ``out`` is not None and not of the
        same type as self.

        """
        if out is None:
            d = self.filled(True).all(axis=axis).view(type(self))
            if d.ndim > 0:
                d.__setmask__(self._mask.all(axis))
            return d
        elif type(out) is not type(self):
            raise TypeError("The external array should have " \
                            "a type %s (got %s instead)" %\
                            (type(self), type(out)))
        self.filled(True).all(axis=axis, out=out)
        if out.ndim:
            out.__setmask__(self._mask.all(axis))
        return out


    def any(self, axis=None, out=None):
        """Returns True if at least one entry along the given axis is
        True.

        Returns False if all entries are False.
        Masked values are considered as True during computation.

        Parameter
        ----------
            axis : int, optional
                Axis along which the operation is performed.
                If None, the operation is performed on a flatten array
            out : {MaskedArray}, optional
                Alternate optional output.  If not None, out should be
                a valid MaskedArray of the same shape as the output of
                self._data.all(axis).

        Returns            A masked array, where the mask is True if all data along
        -------
        the axis are masked.

        Notes
        -----
        An exception is raised if ``out`` is not None and not of the
        same type as self.

        """
        if out is None:
            d = self.filled(False).any(axis=axis).view(type(self))
            if d.ndim > 0:
                d.__setmask__(self._mask.all(axis))
            return d
        elif type(out) is not type(self):
            raise TypeError("The external array should have a type %s "\
                            "(got %s instead)" %\
                            (type(self), type(out)))
        self.filled(False).any(axis=axis, out=out)
        if out.ndim:
            out.__setmask__(self._mask.all(axis))
        return out


    def nonzero(self):
        """Return the indices of the elements of a that are not zero
        nor masked, as a tuple of arrays.

        There are as many tuples as dimensions of a, each tuple
        contains the indices of the non-zero elements in that
        dimension.  The corresponding non-zero values can be obtained
        with ``a[a.nonzero()]``.

        To group the indices by element, rather than dimension, use
        instead: ``transpose(a.nonzero())``.

        The result of this is always a 2d array, with a row for each
        non-zero element.

        """
        return narray(self.filled(0), copy=False).nonzero()
    #............................................
    def trace(self, offset=0, axis1=0, axis2=1, dtype=None, out=None):
        """a.trace(offset=0, axis1=0, axis2=1, dtype=None, out=None)

        Return the sum along the offset diagonal of the array's
        indicated `axis1` and `axis2`.

        """
        # TODO: What are we doing with `out`?
        m = self._mask
        if m is nomask:
            result = super(MaskedArray, self).trace(offset=offset, axis1=axis1,
                                                    axis2=axis2, out=out)
            return result.astype(dtype)
        else:
            D = self.diagonal(offset=offset, axis1=axis1, axis2=axis2)
            return D.astype(dtype).filled(0).sum(axis=None)
    #............................................
    def sum(self, axis=None, dtype=None):
        """Sum the array over the given axis.

        Masked elements are set to 0 internally.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        dtype : dtype, optional
            Datatype for the intermediary computation. If not given,
            the current dtype is used instead.

        """
        if self._mask is nomask:
            mask = nomask
        else:
            mask = self._mask.all(axis)
            if (not mask.ndim) and mask:
                return masked
        result = self.filled(0).sum(axis, dtype=dtype).view(type(self))
        if result.ndim > 0:
            result.__setmask__(mask)
        return result

    def cumsum(self, axis=None, dtype=None):
        """Return the cumulative sum of the elements of the array
        along the given axis.

        Masked values are set to 0 internally.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        dtype : {dtype}, optional
            Datatype for the intermediary computation. If not
            given, the current dtype is used instead.

        """
        result = self.filled(0).cumsum(axis=axis, dtype=dtype).view(type(self))
        result.__setmask__(self.mask)
        return result

    def prod(self, axis=None, dtype=None):
        """Return the product of the elements of the array along the
        given axis.

        Masked elements are set to 1 internally.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        dtype : {dtype}, optional
            Datatype for the intermediary computation. If not
            given, the current dtype is used instead.

        """
        if self._mask is nomask:
            mask = nomask
        else:
            mask = self._mask.all(axis)
            if (not mask.ndim) and mask:
                return masked
        result = self.filled(1).prod(axis=axis, dtype=dtype).view(type(self))
        if result.ndim:
            result.__setmask__(mask)
        return result

    product = prod

    def cumprod(self, axis=None, dtype=None):
        """Return the cumulative product of the elements of the array
        along the given axis.

        Masked values are set to 1 internally.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        dtype : {dtype}, optional
            Datatype for the intermediary computation. If not
            given, the current dtype is used instead.

        """
        result = self.filled(1).cumprod(axis=axis, dtype=dtype).view(type(self))
        result.__setmask__(self.mask)
        return result

    def mean(self, axis=None, dtype=None, out=None):
        """Average the array over the given axis.  Equivalent to

        a.sum(axis, dtype) / a.size(axis).

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        dtype : {dtype}, optional
            Datatype for the intermediary computation. If not
            given, the current dtype is used instead.

        """
        if self._mask is nomask:
            result = super(MaskedArray, self).mean(axis=axis, dtype=dtype)
        else:
            dsum = self.sum(axis=axis, dtype=dtype)
            cnt = self.count(axis=axis)
            result = dsum*1./cnt
        if out is not None:
            out.flat = result.ravel()
        return result

    def anom(self, axis=None, dtype=None):
        """Return the anomalies (deviations from the average) along
        the given axis.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        dtype : {dtype}, optional
            Datatype for the intermediary computation. If not
            given, the current dtype is used instead.

        """
        m = self.mean(axis, dtype)
        if not axis:
            return (self - m)
        else:
            return (self - expand_dims(m,axis))

    def var(self, axis=None, dtype=None, ddof=0):
        """Return the variance, a measure of the spread of a distribution.

        The variance is the average of the squared deviations from the
        mean, i.e. var = mean(abs(x - x.mean())**2).

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        dtype : {dtype}, optional
            Datatype for the intermediary computation. If not
            given, the current dtype is used instead.

        Notes
        -----
        The value returned is by default a biased estimate of the
        true variance, since the mean is computed by dividing by N-ddof.
        For the (more standard) unbiased estimate, use ddof=1 or.
        Note that for complex numbers the absolute value is taken before
        squaring, so that the result is always real and nonnegative.

        """
        if self._mask is nomask:
            # TODO: Do we keep super, or var _data and take a view ?
            return super(MaskedArray, self).var(axis=axis, dtype=dtype,
                                                ddof=ddof)
        else:
            cnt = self.count(axis=axis)-ddof
            danom = self.anom(axis=axis, dtype=dtype)
            if iscomplexobj(self):
                danom = umath.absolute(danom)**2
            else:
                danom *= danom
            dvar = narray(danom.sum(axis) / cnt).view(type(self))
            if axis is not None:
                dvar._mask = mask_or(self._mask.all(axis), (cnt==1))
            dvar._update_from(self)
            return dvar

    def std(self, axis=None, dtype=None, ddof=0):
        """Return the standard deviation, a measure of the spread of a
        distribution.

        The standard deviation is the square root of the average of
        the squared deviations from the mean, i.e.

        std = sqrt(mean(abs(x - x.mean())**2)).

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        dtype : {dtype}, optional
            Datatype for the intermediary computation.
            If not given, the current dtype is used instead.

        Notes
        -----
        The value returned is by default a biased estimate of the
        true standard deviation, since the mean is computed by dividing
        by N-ddof.  For the more standard unbiased estimate, use ddof=1.
        Note that for complex numbers the absolute value is taken before
        squaring, so that the result is always real and nonnegative.
        """
        dvar = self.var(axis,dtype,ddof=ddof)
        if axis is not None or dvar is not masked:
            dvar = sqrt(dvar)
        return dvar

    #............................................
    def round(self, decimals=0, out=None):
        result = self._data.round(decimals).view(type(self))
        result._mask = self._mask
        result._update_from(self)
        if out is None:
            return result
        out[:] = result
        return
    round.__doc__ = ndarray.round.__doc__

    #............................................
    def argsort(self, axis=None, fill_value=None, kind='quicksort',
                order=None):
        """Return an ndarray of indices that sort the array along the
        specified axis.  Masked values are filled beforehand to
        fill_value.

        Parameters
        ----------
        axis : int, optional
            Axis to be indirectly sorted.
            If not given, uses a flatten version of the array.
        fill_value : {var}
            Value used to fill in the masked values.
            If not given, self.fill_value is used instead.
        kind : {string}
            Sorting algorithm (default 'quicksort')
            Possible values: 'quicksort', 'mergesort', or 'heapsort'

        Notes
        -----
        This method executes an indirect sort along the given axis
        using the algorithm specified by the kind keyword. It returns
        an array of indices of the same shape as 'a' that index data
        along the given axis in sorted order.

        The various sorts are characterized by average speed, worst
        case performance need for work space, and whether they are
        stable.  A stable sort keeps items with the same key in the
        same relative order. The three available algorithms have the
        following properties:

        |------------------------------------------------------|
        |    kind   | speed |  worst case | work space | stable|
        |------------------------------------------------------|
        |'quicksort'|   1   | O(n^2)      |     0      |   no  |
        |'mergesort'|   2   | O(n*log(n)) |    ~n/2    |   yes |
        |'heapsort' |   3   | O(n*log(n)) |     0      |   no  |
        |------------------------------------------------------|

        All the sort algorithms make temporary copies of the data when
        the sort is not along the last axis. Consequently, sorts along
        the last axis are faster and use less space than sorts along
        other axis.

        """
        if fill_value is None:
            fill_value = default_fill_value(self)
        d = self.filled(fill_value).view(ndarray)
        return d.argsort(axis=axis, kind=kind, order=order)
    #........................
    def argmin(self, axis=None, fill_value=None):
        """Return an ndarray of indices for the minimum values of a
        along the specified axis.

        Masked values are treated as if they had the value fill_value.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        fill_value : {var}, optional
            Value used to fill in the masked values.  If None, the
            output of minimum_fill_value(self._data) is used.

        """
        if fill_value is None:
            fill_value = minimum_fill_value(self)
        d = self.filled(fill_value).view(ndarray)
        return d.argmin(axis)
    #........................
    def argmax(self, axis=None, fill_value=None):
        """Returns the array of indices for the maximum values of `a`
        along the specified axis.

        Masked values are treated as if they had the value fill_value.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        fill_value : {var}, optional
            Value used to fill in the masked values.  If None, the
            output of maximum_fill_value(self._data) is used.

        """
        if fill_value is None:
            fill_value = maximum_fill_value(self._data)
        d = self.filled(fill_value).view(ndarray)
        return d.argmax(axis)

    def sort(self, axis=-1, kind='quicksort', order=None,
             endwith=True, fill_value=None):
        """Sort along the given axis.

        Parameters
        ----------
        axis : int
            Axis to be indirectly sorted.
        kind : {string}
            Sorting algorithm (default 'quicksort')
            Possible values: 'quicksort', 'mergesort', or 'heapsort'.
        order : {var}
            If a has fields defined, then the order keyword can be
            the field name to sort on or a list (or tuple) of
            field names to indicate the order that fields should
            be used to define the sort.
        fill_value : {var}
            Value used to fill in the masked values.  If None, use
            the the output of minimum_fill_value().
        endwith : bool
            Whether missing values (if any) should be forced in
            the upper indices (at the end of the array) (True) or
            lower indices (at the beginning).

        Returns
        -------
        When used as method, returns None.
            When used as a function, returns an array.

        Notes
        -----
        This method sorts 'a' in place along the given axis using
        the algorithm specified by the kind keyword.

        The various sorts may characterized by average speed,
        worst case performance need for work space, and whether
        they are stable.  A stable sort keeps items with the same
        key in the same relative order and is most useful when
        used w/ argsort where the key might differ from the items
        being sorted.  The three available algorithms have the
        following properties:

        |------------------------------------------------------|
        |    kind   | speed |  worst case | work space | stable|
        |------------------------------------------------------|
        |'quicksort'|   1   | O(n^2)      |     0      |   no  |
        |'mergesort'|   2   | O(n*log(n)) |    ~n/2    |   yes |
        |'heapsort' |   3   | O(n*log(n)) |     0      |   no  |
        |------------------------------------------------------|

        """
        if self._mask is nomask:
            ndarray.sort(self,axis=axis, kind=kind, order=order)
        else:
            if fill_value is None:
                if endwith:
                    filler = minimum_fill_value(self)
                else:
                    filler = maximum_fill_value(self)
            else:
                filler = fill_value
            idx = numpy.indices(self.shape)
            idx[axis] = self.filled(filler).argsort(axis=axis,kind=kind,order=order)
            idx_l = idx.tolist()
            tmp_mask = self._mask[idx_l].flat
            tmp_data = self._data[idx_l].flat
            self.flat = tmp_data
            self._mask.flat = tmp_mask
        return

    #............................................
    def min(self, axis=None, fill_value=None):
        """Return the minimum of a along the given axis.

        Masked values are filled with fill_value.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        fill_value : {var}, optional
            Value used to fill in the masked values.
            If None, use the the output of minimum_fill_value().

    """
        mask = self._mask
        # Check all/nothing case ......
        if mask is nomask:
            return super(MaskedArray, self).min(axis=axis)
        elif (not mask.ndim) and mask:
            return masked
        # Get the mask ................
        if axis is None:
            mask = umath.logical_and.reduce(mask.flat)
        else:
            mask = umath.logical_and.reduce(mask, axis=axis)
        # Skip if all masked ..........
        if not mask.ndim and mask:
            return masked
        # Get the fill value ...........
        if fill_value is None:
            fill_value = minimum_fill_value(self)
        # Get the data ................
        result = self.filled(fill_value).min(axis=axis).view(type(self))
        if result.ndim > 0:
            result._mask = mask
        return result

    def mini(self, axis=None):
        if axis is None:
            return minimum(self)
        else:
            return minimum.reduce(self, axis)

    #........................
    def max(self, axis=None, fill_value=None):
        """Return the maximum/a along the given axis.

        Masked values are filled with fill_value.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        fill_value : {var}, optional
            Value used to fill in the masked values.
            If None, use the the output of maximum_fill_value().
        """
        mask = self._mask
        # Check all/nothing case ......
        if mask is nomask:
            return super(MaskedArray, self).max(axis=axis)
        elif (not mask.ndim) and mask:
            return masked
        # Check the mask ..............
        if axis is None:
            mask = umath.logical_and.reduce(mask.flat)
        else:
            mask = umath.logical_and.reduce(mask, axis=axis)
        # Skip if all masked ..........
        if not mask.ndim and mask:
            return masked
        # Get the fill value ..........
        if fill_value is None:
            fill_value = maximum_fill_value(self)
        # Get the data ................
        result = self.filled(fill_value).max(axis=axis).view(type(self))
        if result.ndim > 0:
            result._mask = mask
        return result
    #........................
    def ptp(self, axis=None, fill_value=None):
        """Return the visible data range (max-min) along the given axis.

        Parameters
        ----------
        axis : int, optional
            Axis along which to perform the operation.
            If None, applies to a flattened version of the array.
        fill_value : {var}, optional
            Value used to fill in the masked values.  If None, the
            maximum uses the maximum default, the minimum uses the
            minimum default.

        """
        return self.max(axis, fill_value) - self.min(axis, fill_value)

    # Array methods ---------------------------------------
    copy = _arraymethod('copy')
    diagonal = _arraymethod('diagonal')
    take = _arraymethod('take')
    transpose = _arraymethod('transpose')
    T = property(fget=lambda self:self.transpose())
    swapaxes = _arraymethod('swapaxes')
    clip = _arraymethod('clip', onmask=False)
    copy = _arraymethod('copy')
    squeeze = _arraymethod('squeeze')
    #--------------------------------------------
    def tolist(self, fill_value=None):
        """Copy the data portion of the array to a hierarchical python
        list and returns that list.

        Data items are converted to the nearest compatible Python
        type.  Masked values are converted to fill_value. If
        fill_value is None, the corresponding entries in the output
        list will be ``None``.

        """
        if fill_value is not None:
            return self.filled(fill_value).tolist()
        result = self.filled().tolist()
        # Set temps to save time when dealing w/ mrecarrays...
        _mask = self._mask
        if _mask is nomask:
            return result
        nbdims = self.ndim
        dtypesize = len(self.dtype)
        if nbdims == 0:
            return tuple([None]*dtypesize)
        elif nbdims == 1:
            maskedidx = _mask.nonzero()[0].tolist()
            if dtypesize:
                nodata = tuple([None]*dtypesize)
            else:
                nodata = None
            [operator.setitem(result,i,nodata) for i in maskedidx]
        else:
            for idx in zip(*[i.tolist() for i in _mask.nonzero()]):
                tmp = result
                for i in idx[:-1]:
                    tmp = tmp[i]
                tmp[idx[-1]] = None
        return result
    #........................
    def tostring(self, fill_value=None, order='C'):
        """Return a copy of array data as a Python string containing the raw
        bytes in the array.

        Parameters
        ----------
        fill_value : {var}, optional
            Value used to fill in the masked values.
            If None, uses self.fill_value instead.
        order : {string}
            Order of the data item in the copy {"C","F","A"}.
            "C"       -- C order (row major)
            "Fortran" -- Fortran order (column major)
            "Any"     -- Current order of array.
            None      -- Same as "Any"

        """
        return self.filled(fill_value).tostring(order=order)
    #........................
    def tofile(self, fid, sep="", format="%s"):
        raise NotImplementedError("Not implemented yet, sorry...")

    #--------------------------------------------
    # Pickling
    def __getstate__(self):
        """Return the internal state of the masked array, for pickling
        purposes.

        """
        state = (1,
                 self.shape,
                 self.dtype,
                 self.flags.fnc,
                 self._data.tostring(),
                 getmaskarray(self).tostring(),
                 self._fill_value,
                 )
        return state
    #
    def __setstate__(self, state):
        """Restore the internal state of the masked array, for
        pickling purposes.  ``state`` is typically the output of the
        ``__getstate__`` output, and is a 5-tuple:

        - class name
        - a tuple giving the shape of the data
        - a typecode for the data
        - a binary string for the data
        - a binary string for the mask.

        """
        (ver, shp, typ, isf, raw, msk, flv) = state
        ndarray.__setstate__(self, (shp, typ, isf, raw))
        self._mask.__setstate__((shp, dtype(bool), isf, msk))
        self.fill_value = flv
    #
    def __reduce__(self):
        """Return a 3-tuple for pickling a MaskedArray.

        """
        return (_mareconstruct,
                (self.__class__, self._baseclass, (0,), 'b', ),
                self.__getstate__())


def _mareconstruct(subtype, baseclass, baseshape, basetype,):
    """Internal function that builds a new MaskedArray from the
    information stored in a pickle.

    """
    _data = ndarray.__new__(baseclass, baseshape, basetype)
    _mask = ndarray.__new__(ndarray, baseshape, 'b1')
    return subtype.__new__(subtype, _data, mask=_mask, dtype=basetype,)


#####--------------------------------------------------------------------------
#---- --- Shortcuts ---
#####---------------------------------------------------------------------------
def isMaskedArray(x):
    "Is x a masked array, that is, an instance of MaskedArray?"
    return isinstance(x, MaskedArray)
isarray = isMaskedArray
isMA = isMaskedArray  #backward compatibility
# We define the masked singleton as a float for higher precedence...
# Note that it can be tricky sometimes w/ type comparison
masked_singleton = MaskedArray(0, dtype=float_, mask=True)
masked = masked_singleton

masked_array = MaskedArray

def array(data, dtype=None, copy=False, order=False,
          mask=nomask, fill_value=None,
          keep_mask=True, hard_mask=False, shrink=True, subok=True, ndmin=0,
          ):
    """array(data, dtype=None, copy=False, order=False, mask=nomask,
             fill_value=None, keep_mask=True, hard_mask=False, shrink=True,
             subok=True, ndmin=0)

    Acts as shortcut to MaskedArray, with options in a different order
    for convenience.  And backwards compatibility...

    """
    #TODO: we should try to put 'order' somwehere
    return MaskedArray(data, mask=mask, dtype=dtype, copy=copy, subok=subok,
                       keep_mask=keep_mask, hard_mask=hard_mask,
                       fill_value=fill_value, ndmin=ndmin, shrink=shrink)
array.__doc__ = masked_array.__doc__

def is_masked(x):
    """Does x have masked values?"""
    m = getmask(x)
    if m is nomask:
        return False
    elif m.any():
        return True
    return False


#####---------------------------------------------------------------------------
#---- --- Extrema functions ---
#####---------------------------------------------------------------------------
class _extrema_operation(object):
    "Generic class for maximum/minimum functions."
    def __call__(self, a, b=None):
        "Executes the call behavior."
        if b is None:
            return self.reduce(a)
        return where(self.compare(a, b), a, b)
    #.........
    def reduce(self, target, axis=None):
        "Reduce target along the given axis."
        target = narray(target, copy=False, subok=True)
        m = getmask(target)
        if axis is not None:
            kargs = { 'axis' : axis }
        else:
            kargs = {}
            target = target.ravel()
            if not (m is nomask):
                m = m.ravel()
        if m is nomask:
            t = self.ufunc.reduce(target, **kargs)
        else:
            target = target.filled(self.fill_value_func(target)).view(type(target))
            t = self.ufunc.reduce(target, **kargs)
            m = umath.logical_and.reduce(m, **kargs)
            if hasattr(t, '_mask'):
                t._mask = m
            elif m:
                t = masked
        return t
    #.........
    def outer (self, a, b):
        "Return the function applied to the outer product of a and b."
        ma = getmask(a)
        mb = getmask(b)
        if ma is nomask and mb is nomask:
            m = nomask
        else:
            ma = getmaskarray(a)
            mb = getmaskarray(b)
            m = logical_or.outer(ma, mb)
        result = self.ufunc.outer(filled(a), filled(b))
        result._mask = m
        return result

#............................
class _minimum_operation(_extrema_operation):
    "Object to calculate minima"
    def __init__ (self):
        """minimum(a, b) or minimum(a)
In one argument case, returns the scalar minimum.
        """
        self.ufunc = umath.minimum
        self.afunc = amin
        self.compare = less
        self.fill_value_func = minimum_fill_value

#............................
class _maximum_operation(_extrema_operation):
    "Object to calculate maxima"
    def __init__ (self):
        """maximum(a, b) or maximum(a)
           In one argument case returns the scalar maximum.
        """
        self.ufunc = umath.maximum
        self.afunc = amax
        self.compare = greater
        self.fill_value_func = maximum_fill_value

#..........................................................
def min(array, axis=None, out=None):
    """Return the minima along the given axis.

    If `axis` is None, applies to the flattened array.

    """
    if out is not None:
        raise TypeError("Output arrays Unsupported for masked arrays")
    if axis is None:
        return minimum(array)
    else:
        return minimum.reduce(array, axis)
min.__doc__ = MaskedArray.min.__doc__
#............................
def max(obj, axis=None, out=None):
    if out is not None:
        raise TypeError("Output arrays Unsupported for masked arrays")
    if axis is None:
        return maximum(obj)
    else:
        return maximum.reduce(obj, axis)
max.__doc__ = MaskedArray.max.__doc__
#.............................
def ptp(obj, axis=None):
    """a.ptp(axis=None) =  a.max(axis)-a.min(axis)"""
    try:
        return obj.max(axis)-obj.min(axis)
    except AttributeError:
        return max(obj, axis=axis) - min(obj, axis=axis)
ptp.__doc__ = MaskedArray.ptp.__doc__


#####---------------------------------------------------------------------------
#---- --- Definition of functions from the corresponding methods ---
#####---------------------------------------------------------------------------
class _frommethod:
    """Define functions from existing MaskedArray methods.

    Parameters
    ----------
        _methodname : string
            Name of the method to transform.

    """
    def __init__(self, methodname):
        self._methodname = methodname
        self.__doc__ = self.getdoc()
    def getdoc(self):
        "Return the doc of the function (from the doc of the method)."
        try:
            return getattr(MaskedArray, self._methodname).__doc__
        except:
            return getattr(numpy, self._methodname).__doc__
    def __call__(self, a, *args, **params):
        if isinstance(a, MaskedArray):
            return getattr(a, self._methodname).__call__(*args, **params)
        #FIXME ----
        #As x is not a MaskedArray, we transform it to a ndarray with asarray
        #... and call the corresponding method.
        #Except that sometimes it doesn't work (try reshape([1,2,3,4],(2,2)))
        #we end up with a "SystemError: NULL result without error in PyObject_Call"
        #A dirty trick is then to call the initial numpy function...
        method = getattr(narray(a, copy=False), self._methodname)
        try:
            return method(*args, **params)
        except SystemError:
            return getattr(numpy,self._methodname).__call__(a, *args, **params)

all = _frommethod('all')
anomalies = anom = _frommethod('anom')
any = _frommethod('any')
conjugate = _frommethod('conjugate')
ids = _frommethod('ids')
nonzero = _frommethod('nonzero')
diagonal = _frommethod('diagonal')
maximum = _maximum_operation()
mean = _frommethod('mean')
minimum = _minimum_operation ()
product = _frommethod('prod')
ptp = _frommethod('ptp')
ravel = _frommethod('ravel')
repeat = _frommethod('repeat')
round = _frommethod('round')
std = _frommethod('std')
sum = _frommethod('sum')
swapaxes = _frommethod('swapaxes')
take = _frommethod('take')
trace = _frommethod('trace')
var = _frommethod('var')
compress = _frommethod('compress')

#..............................................................................
def power(a, b, third=None):
    """Computes a**b elementwise.

    """
    if third is not None:
        raise MAError, "3-argument power not supported."
    # Get the masks
    ma = getmask(a)
    mb = getmask(b)
    m = mask_or(ma, mb)
    # Get the rawdata
    fa = getdata(a)
    fb = getdata(b)
    # Get the type of the result (so that we preserve subclasses)
    if isinstance(a,MaskedArray):
        basetype = type(a)
    else:
        basetype = MaskedArray
    # Get the result and view it as a (subclass of) MaskedArray
    result = umath.power(fa,fb).view(basetype)
    # Find where we're in trouble w/ NaNs and Infs
    invalid = numpy.logical_not(numpy.isfinite(result.view(ndarray)))
    # Retrieve some extra attributes if needed
    if isinstance(result,MaskedArray):
        result._update_from(a)
    # Add the initial mask
    if m is not nomask:
        if numpy.isscalar(result):
            return masked
        result._mask = m
    # Fix the invalid parts
    if invalid.any():
        if not result.ndim:
            return masked
        result[invalid] = masked
        result._data[invalid] = result.fill_value
    return result
    
#    if fb.dtype.char in typecodes["Integer"]:
#        return masked_array(umath.power(fa, fb), m)
#    m = mask_or(m, (fa < 0) & (fb != fb.astype(int))) 
#    if m is nomask:
#        return masked_array(umath.power(fa, fb))
#    else:
#        fa = fa.copy()
#        if m.all():
#            fa.flat = 1
#        else: 
#            numpy.putmask(fa,m,1)
#        return masked_array(umath.power(fa, fb), m)

#..............................................................................
def argsort(a, axis=None, kind='quicksort', order=None, fill_value=None):
    "Function version of the eponymous method."
    if fill_value is None:
        fill_value = default_fill_value(a)
    d = filled(a, fill_value)
    if axis is None:
        return d.argsort(kind=kind, order=order)
    return d.argsort(axis, kind=kind, order=order)
argsort.__doc__ = MaskedArray.argsort.__doc__

def argmin(a, axis=None, fill_value=None):
    "Function version of the eponymous method."
    if fill_value is None:
        fill_value = default_fill_value(a)
    d = filled(a, fill_value)
    return d.argmin(axis=axis)
argmin.__doc__ = MaskedArray.argmin.__doc__

def argmax(a, axis=None, fill_value=None):
    "Function version of the eponymous method."
    if fill_value is None:
        fill_value = default_fill_value(a)
        try:
            fill_value = - fill_value
        except:
            pass
    d = filled(a, fill_value)
    return d.argmax(axis=axis)
argmin.__doc__ = MaskedArray.argmax.__doc__

def sort(a, axis=-1, kind='quicksort', order=None, endwith=True, fill_value=None):
    "Function version of the eponymous method."
    a = narray(a, copy=True, subok=True)
    if axis is None:
        a = a.flatten()
        axis = 0
    if fill_value is None:
        if endwith:
            filler = minimum_fill_value(a)
        else:
            filler = maximum_fill_value(a)
    else:
        filler = fill_value
#    return
    indx = numpy.indices(a.shape).tolist()
    indx[axis] = filled(a,filler).argsort(axis=axis,kind=kind,order=order)
    return a[indx]
sort.__doc__ = MaskedArray.sort.__doc__


def compressed(x):
    """Return a 1-D array of all the non-masked data."""
    if getmask(x) is nomask:
        return numpy.asanyarray(x)
    else:
        return x.compressed()

def concatenate(arrays, axis=0):
    "Concatenate the arrays along the given axis."
    d = numpy.concatenate([getdata(a) for a in arrays], axis)
    rcls = get_masked_subclass(*arrays)
    data = d.view(rcls)
    # Check whether one of the arrays has a non-empty mask...
    for x in arrays:
        if getmask(x) is not nomask:
            break
    else:
        return data
    # OK, so we have to concatenate the masks
    dm = numpy.concatenate([getmaskarray(a) for a in arrays], axis)
    # If we decide to keep a '_shrinkmask' option, we want to check that ...
    # ... all of them are True, and then check for dm.any()
#    shrink = numpy.logical_or.reduce([getattr(a,'_shrinkmask',True) for a in arrays])
#    if shrink and not dm.any():
    if not dm.any():
        data._mask = nomask
    else:
        data._mask = dm.reshape(d.shape)
    return data

def count(a, axis = None):
    return masked_array(a, copy=False).count(axis)
count.__doc__ = MaskedArray.count.__doc__


def expand_dims(x,axis):
    """Expand the shape of the array by including a new axis before
    the given one.

    """
    result = n_expand_dims(x,axis)
    if isinstance(x, MaskedArray):
        new_shape = result.shape
        result = x.view()
        result.shape = new_shape
        if result._mask is not nomask:
            result._mask.shape = new_shape
    return result

#......................................
def left_shift (a, n):
    "Left shift n bits."
    m = getmask(a)
    if m is nomask:
        d = umath.left_shift(filled(a), n)
        return masked_array(d)
    else:
        d = umath.left_shift(filled(a, 0), n)
        return masked_array(d, mask=m)

def right_shift (a, n):
    "Right shift n bits."
    m = getmask(a)
    if m is nomask:
        d = umath.right_shift(filled(a), n)
        return masked_array(d)
    else:
        d = umath.right_shift(filled(a, 0), n)
        return masked_array(d, mask=m)

#......................................
def put(a, indices, values, mode='raise'):
    """Set storage-indexed locations to corresponding values.

    Values and indices are filled if necessary.

    """
    # We can't use 'frommethod', the order of arguments is different
    try:
        return a.put(indices, values, mode=mode)
    except AttributeError:
        return narray(a, copy=False).put(indices, values, mode=mode)

def putmask(a, mask, values): #, mode='raise'):
    """Set a.flat[n] = values[n] for each n where mask.flat[n] is true.

    If values is not the same size of a and mask then it will repeat
    as necessary.  This gives different behavior than
    a[mask] = values.

    Note: Using a masked array as values will NOT transform a ndarray in
          a maskedarray.

    """
    # We can't use 'frommethod', the order of arguments is different
    if not isinstance(a, MaskedArray):
        a = a.view(MaskedArray)
    (valdata, valmask) = (getdata(values), getmask(values))
    if getmask(a) is nomask:
        if valmask is not nomask:
            a._sharedmask = True
            a.mask = numpy.zeros(a.shape, dtype=bool_)
            numpy.putmask(a._mask, mask, valmask)
    elif a._hardmask:
        if valmask is not nomask:
            m = a._mask.copy()
            numpy.putmask(m, mask, valmask)
            a.mask |= m
    else:
        if valmask is nomask:
            valmask = getmaskarray(values)
        numpy.putmask(a._mask, mask, valmask)
    numpy.putmask(a._data, mask, valdata)
    return

def transpose(a,axes=None):
    """Return a view of the array with dimensions permuted according to axes,
    as a masked array.

    If ``axes`` is None (default), the output view has reversed
    dimensions compared to the original.

    """
    #We can't use 'frommethod', as 'transpose' doesn't take keywords
    try:
        return a.transpose(axes)
    except AttributeError:
        return narray(a, copy=False).transpose(axes).view(MaskedArray)

def reshape(a, new_shape):
    """Change the shape of the array a to new_shape."""
    #We can't use 'frommethod', it whine about some parameters. Dmmit.
    try:
        return a.reshape(new_shape)
    except AttributeError:
        return narray(a, copy=False).reshape(new_shape).view(MaskedArray)

def resize(x, new_shape):
    """Return a new array with the specified shape.

    The total size of the original array can be any size.  The new
    array is filled with repeated copies of a. If a was masked, the
    new array will be masked, and the new mask will be a repetition of
    the old one.

    """
    # We can't use _frommethods here, as N.resize is notoriously whiny.
    m = getmask(x)
    if m is not nomask:
        m = numpy.resize(m, new_shape)
    result = numpy.resize(x, new_shape).view(get_masked_subclass(x))
    if result.ndim:
        result._mask = m
    return result


#................................................
def rank(obj):
    "maskedarray version of the numpy function."
    return fromnumeric.rank(getdata(obj))
rank.__doc__ = numpy.rank.__doc__
#
def shape(obj):
    "maskedarray version of the numpy function."
    return fromnumeric.shape(getdata(obj))
shape.__doc__ = numpy.shape.__doc__
#
def size(obj, axis=None):
    "maskedarray version of the numpy function."
    return fromnumeric.size(getdata(obj), axis)
size.__doc__ = numpy.size.__doc__
#................................................

#####--------------------------------------------------------------------------
#---- --- Extra functions ---
#####--------------------------------------------------------------------------
def where (condition, x=None, y=None):
    """where(condition | x, y)

    Returns a (subclass of) masked array, shaped like condition, where
    the elements are x when condition is True, and y otherwise.  If
    neither x nor y are given, returns a tuple of indices where
    condition is True (a la condition.nonzero()).

    Parameters
    ----------
    condition : {var}
        The condition to meet. Must be convertible to an integer
        array.
    x : {var}, optional
        Values of the output when the condition is met
    y : {var}, optional
        Values of the output when the condition is not met.

    """
    if x is None and y is None:
        return filled(condition, 0).nonzero()
    elif x is None or y is None:
        raise ValueError, "Either both or neither x and y should be given."
    # Get the condition ...............
    fc = filled(condition, 0).astype(bool_)
    notfc = numpy.logical_not(fc)
    # Get the data ......................................
    xv = getdata(x)
    yv = getdata(y)
    if x is masked:
        ndtype = yv.dtype
        xm = numpy.ones(fc.shape, dtype=MaskType)
    elif y is masked:
        ndtype = xv.dtype
        ym = numpy.ones(fc.shape, dtype=MaskType)
    else:
        ndtype = numpy.max([xv.dtype, yv.dtype])
        xm = getmask(x)
    d = numpy.empty(fc.shape, dtype=ndtype).view(MaskedArray)
    numpy.putmask(d._data, fc, xv.astype(ndtype))
    numpy.putmask(d._data, notfc, yv.astype(ndtype))
    d._mask = numpy.zeros(fc.shape, dtype=MaskType)
    numpy.putmask(d._mask, fc, getmask(x))
    numpy.putmask(d._mask, notfc, getmask(y))
    d._mask |= getmaskarray(condition)
    if not d._mask.any():
        d._mask = nomask
    return d

def choose (indices, t, out=None, mode='raise'):
    "Return array shaped like indices with elements chosen from t"
    #TODO: implement options `out` and `mode`, if possible.
    def fmask (x):
        "Returns the filled array, or True if masked."
        if x is masked:
            return 1
        return filled(x)
    def nmask (x):
        "Returns the mask, True if ``masked``, False if ``nomask``."
        if x is masked:
            return 1
        m = getmask(x)
        if m is nomask:
            return 0
        return m
    c = filled(indices, 0)
    masks = [nmask(x) for x in t]
    a = [fmask(x) for x in t]
    d = numpy.choose(c, a)
    m = numpy.choose(c, masks)
    m = make_mask(mask_or(m, getmask(indices)), copy=0, shrink=True)
    return masked_array(d, mask=m)

def round_(a, decimals=0, out=None):
    """Return a copy of a, rounded to 'decimals' places.

    When 'decimals' is negative, it specifies the number of positions
    to the left of the decimal point.  The real and imaginary parts of
    complex numbers are rounded separately. Nothing is done if the
    array is not of float type and 'decimals' is greater than or equal
    to 0.

    Parameters
    ----------
    decimals : int
        Number of decimals to round to. May be negative.
    out : array_like
        Existing array to use for output.
        If not given, returns a default copy of a.

    Notes
    -----
    If out is given and does not have a mask attribute, the mask of a
    is lost!

    """
    if out is None:
        return numpy.round_(a, decimals, out)
    else:
        numpy.round_(getdata(a), decimals, out)
        if hasattr(out, '_mask'):
            out._mask = getmask(a)
        return out

def arange(stop, start=None, step=1, dtype=None):
    "maskedarray version of the numpy function."
    return numpy.arange(stop, start, step, dtype).view(MaskedArray)
arange.__doc__ = numpy.arange.__doc__

def inner(a, b):
    "maskedarray version of the numpy function."
    fa = filled(a, 0)
    fb = filled(b, 0)
    if len(fa.shape) == 0:
        fa.shape = (1,)
    if len(fb.shape) == 0:
        fb.shape = (1,)
    return numpy.inner(fa, fb).view(MaskedArray)
inner.__doc__ = numpy.inner.__doc__
inner.__doc__ += doc_note("Masked values are replaced by 0.")
innerproduct = inner

def outer(a, b):
    "maskedarray version of the numpy function."
    fa = filled(a, 0).ravel()
    fb = filled(b, 0).ravel()
    d = numeric.outer(fa, fb)
    ma = getmask(a)
    mb = getmask(b)
    if ma is nomask and mb is nomask:
        return masked_array(d)
    ma = getmaskarray(a)
    mb = getmaskarray(b)
    m = make_mask(1-numeric.outer(1-ma, 1-mb), copy=0)
    return masked_array(d, mask=m)
outer.__doc__ = numpy.outer.__doc__
outer.__doc__ += doc_note("Masked values are replaced by 0.")
outerproduct = outer

def allequal (a, b, fill_value=True):
    """Return True if all entries of a and b are equal, using
    fill_value as a truth value where either or both are masked.

    """
    m = mask_or(getmask(a), getmask(b))
    if m is nomask:
        x = getdata(a)
        y = getdata(b)
        d = umath.equal(x, y)
        return d.all()
    elif fill_value:
        x = getdata(a)
        y = getdata(b)
        d = umath.equal(x, y)
        dm = array(d, mask=m, copy=False)
        return dm.filled(True).all(None)
    else:
        return False

def allclose (a, b, fill_value=True, rtol=1.e-5, atol=1.e-8):
    """ Return True if all elements of a and b are equal subject to
    given tolerances.

    If fill_value is True, masked values are considered equal.
    If fill_value is False, masked values considered unequal.
    The relative error rtol should be positive and << 1.0
    The absolute error atol comes into play for those elements of b
    that are very small or zero; it says how small `a` must be also.

    """
    m = mask_or(getmask(a), getmask(b))
    d1 = getdata(a)
    d2 = getdata(b)
    x = filled(array(d1, copy=0, mask=m), fill_value).astype(float)
    y = filled(array(d2, copy=0, mask=m), 1).astype(float)
    d = umath.less_equal(umath.absolute(x-y), atol + rtol * umath.absolute(y))
    return fromnumeric.alltrue(fromnumeric.ravel(d))

#..............................................................................
def asarray(a, dtype=None):
    """asarray(data, dtype) = array(data, dtype, copy=0, subok=0)

    Return a as a MaskedArray object of the given dtype.
    If dtype is not given or None, is is set to the dtype of a.
    No copy is performed if a is already an array.
    Subclasses are converted to the base class MaskedArray.

    """
    return masked_array(a, dtype=dtype, copy=False, keep_mask=True, subok=False)

def asanyarray(a, dtype=None):
    """asanyarray(data, dtype) = array(data, dtype, copy=0, subok=1)

    Return a as an masked array.
    If dtype is not given or None, is is set to the dtype of a.
    No copy is performed if a is already an array.
    Subclasses are conserved.

    """
    return masked_array(a, dtype=dtype, copy=False, keep_mask=True, subok=True)


def empty(new_shape, dtype=float):
    "maskedarray version of the numpy function."
    return numpy.empty(new_shape, dtype).view(MaskedArray)
empty.__doc__ = numpy.empty.__doc__

def empty_like(a):
    "maskedarray version of the numpy function."
    return numpy.empty_like(a).view(MaskedArray)
empty_like.__doc__ = numpy.empty_like.__doc__

def ones(new_shape, dtype=float):
    "maskedarray version of the numpy function."
    return numpy.ones(new_shape, dtype).view(MaskedArray)
ones.__doc__ = numpy.ones.__doc__

def zeros(new_shape, dtype=float):
    "maskedarray version of the numpy function."
    return numpy.zeros(new_shape, dtype).view(MaskedArray)
zeros.__doc__ = numpy.zeros.__doc__

#####--------------------------------------------------------------------------
#---- --- Pickling ---
#####--------------------------------------------------------------------------
def dump(a,F):
    """Pickle the MaskedArray `a` to the file `F`.  `F` can either be
    the handle of an exiting file, or a string representing a file
    name.

    """
    if not hasattr(F,'readline'):
        F = open(F,'w')
    return cPickle.dump(a,F)

def dumps(a):
    """Return a string corresponding to the pickling of the
    MaskedArray.

    """
    return cPickle.dumps(a)

def load(F):
    """Wrapper around ``cPickle.load`` which accepts either a
    file-like object or a filename.

    """
    if not hasattr(F, 'readline'):
        F = open(F,'r')
    return cPickle.load(F)

def loads(strg):
    "Load a pickle from the current string."""
    return cPickle.loads(strg)

################################################################################
def fromfile(file, dtype=float, count=-1, sep=''):
    raise NotImplementedError("Not yet implemented. Sorry")


class _convert2ma:
    """Convert functions from numpy to numpy.ma.

    Parameters
    ----------
        _methodname : string
            Name of the method to transform.

    """
    __doc__ = None
    def __init__(self, funcname):
        self._func = getattr(numpy, funcname)
        self.__doc__ = self.getdoc()
    def getdoc(self):
        "Return the doc of the function (from the doc of the method)."
        return self._func.__doc__
    def __call__(self, a, *args, **params):
        return self._func.__call__(a, *args, **params).view(MaskedArray)

frombuffer = _convert2ma('frombuffer')
fromfunction = _convert2ma('fromfunction')
identity = _convert2ma('identity')
indices = numpy.indices
clip = numpy.clip

###############################################################################