[Scipy-svn] r3371 - in trunk/scipy/sandbox/maskedarray: . alternative_versions tests

scipy-svn@scip... scipy-svn@scip...
Wed Sep 26 22:35:08 CDT 2007


Author: pierregm
Date: 2007-09-26 22:34:44 -0500 (Wed, 26 Sep 2007)
New Revision: 3371

Added:
   trunk/scipy/sandbox/maskedarray/alternative_versions/core_initial.py
   trunk/scipy/sandbox/maskedarray/bench.py
Modified:
   trunk/scipy/sandbox/maskedarray/__init__.py
   trunk/scipy/sandbox/maskedarray/core.py
   trunk/scipy/sandbox/maskedarray/extras.py
   trunk/scipy/sandbox/maskedarray/morestats.py
   trunk/scipy/sandbox/maskedarray/mrecords.py
   trunk/scipy/sandbox/maskedarray/mstats.py
   trunk/scipy/sandbox/maskedarray/tests/test_core.py
   trunk/scipy/sandbox/maskedarray/tests/test_mrecords.py
   trunk/scipy/sandbox/maskedarray/tests/test_subclassing.py
Log:
core : * simplified __getitem__ and __setitem__
	   * no automatic shrinking of the mask	
       * force prefilling on function domain
       * prevent unnecessary filling on function arguments
       * introduced fix_invalid and masked_invalid
       * reintroduced _sharedmask to manage mask propagation
       * updated doc
       
core._arraymethods: make sure the result has the same generic attributes as the initial data

Modified: trunk/scipy/sandbox/maskedarray/__init__.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/__init__.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/__init__.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -17,6 +17,7 @@
 import extras
 from extras import *
 
+import _nfcore
 
 __all__ = ['core', 'extras']
 __all__ += core.__all__

Added: trunk/scipy/sandbox/maskedarray/alternative_versions/core_initial.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/alternative_versions/core_initial.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/alternative_versions/core_initial.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -0,0 +1,2708 @@
+# 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
+:version: $Id: core.py 254 2007-08-15 04:11:52Z backtopop $
+"""
+__author__ = "Pierre GF Gerard-Marchant ($Author: backtopop $)"
+__version__ = '1.0'
+__revision__ = "$Revision: 254 $"
+__date__     = '$Date: 2007-08-15 00:11:52 -0400 (Wed, 15 Aug 2007) $'
+
+__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',
+           'bitwise_and', 'bitwise_or', 'bitwise_xor',
+           'ceil', 'choose', 'compressed', 'concatenate', 'conjugate',
+               'cos', 'cosh', 'count',
+           'diagonal', 'divide', 'dump', 'dumps',
+           'empty', 'empty_like', 'equal', 'exp',
+           'fabs', 'fmod', 'filled', 'floor', 'floor_divide',
+           'getmask', 'getmaskarray', 'greater', 'greater_equal', 'hypot',
+           'ids', '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_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_',
+           '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 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
+from numpy import expand_dims as n_expand_dims
+import warnings
+
+
+MaskType = bool_
+nomask = MaskType(0)
+
+divide_tolerance = 1.e-35
+numpy.seterr(all='ignore')
+
+# TODO: There's still a problem with N.add.reduce not working...
+# TODO: ...neither does N.add.accumulate
+
+#####--------------------------------------------------------------------------
+#---- --- 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,-numeric.inf) for k in [numpy.float32, numpy.float64]])
+min_filler = ntypes._maxvals
+min_filler.update([(k,numeric.inf) for k in [numpy.float32, numpy.float64]])
+if 'float128' in ntypes.typeDict:
+    max_filler.update([(numpy.float128,-numeric.inf)])
+    min_filler.update([(numpy.float128, numeric.inf)])
+
+
+def default_fill_value(obj):
+    "Calculates the default fill value for an object `obj`."
+    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):
+    "Calculates 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):
+    "Calculates 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 set_fill_value(a, fill_value):
+    "Sets the fill value of `a` if it is a masked array."
+    if isinstance(a, MaskedArray):
+        a.set_fill_value(fill_value)
+
+def get_fill_value(a):
+    """Returns the fill value of `a`, if any.
+    Otherwise, returns the default fill 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):
+    "Returns the common fill_value of `a` and `b`, if any, or `None`."
+    t1 = get_fill_value(a)
+    t2 = get_fill_value(b)
+    if t1 == t2:
+        return t1
+    return None
+
+#................................................
+def filled(a, value = None):
+    """Returns `a` as an array with masked data replaced by `value`.
+If `value` is `None` or the special element `masked`, `get_fill_value(a)`
+is used instead.
+
+If `a` is already a contiguous numeric array, `a` itself is returned.
+
+`filled(a)` can be used to be sure that the result is numeric when passing
+an object a to other software ignorant of MA, in particular to numpy itself.
+    """
+    if hasattr(a, 'filled'):
+        return a.filled(value)
+    elif isinstance(a, ndarray): # and a.flags['CONTIGUOUS']:
+        return a
+    elif isinstance(a, dict):
+        return numeric.array(a, 'O')
+    else:
+        return numeric.array(a)
+
+def get_masked_subclass(*arrays):
+    """Returns the youngest subclass of MaskedArray from a list of arrays,
+ or MaskedArray. 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
+
+#####--------------------------------------------------------------------------
+#---- --- Ufuncs ---
+#####--------------------------------------------------------------------------
+ufunc_domain = {}
+ufunc_fills = {}
+
+class domain_check_interval:
+    """Defines 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 domain_tan:
+    """Defines 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):
+        "Execute the call behavior."
+        return umath.less(umath.absolute(umath.cos(x)), self.eps)
+#............................
+class domain_safe_divide:
+    """defines 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 domain_greater:
+    "domain_greater(v)(x) = true where x <= v"
+    def __init__(self, critical_value):
+        "domain_greater(v)(x) = true where x <= v"
+        self.critical_value = critical_value
+
+    def __call__ (self, x):
+        "Execute the call behavior."
+        return umath.less_equal(x, self.critical_value)
+#............................
+class domain_greater_equal:
+    "domain_greater_equal(v)(x) = true where x < v"
+    def __init__(self, critical_value):
+        "domain_greater_equal(v)(x) = true where x < v"
+        self.critical_value = critical_value
+
+    def __call__ (self, x):
+        "Execute the call behavior."
+        return umath.less(x, self.critical_value)
+#..............................................................................
+class masked_unary_operation:
+    """Defines masked version of unary operations,
+where invalid values are pre-masked.
+
+:IVariables:
+    - `f` : function.
+    - `fill` : Default filling value *[0]*.
+    - `domain` : Default domain *[None]*.
+    """
+    def __init__ (self, mufunc, fill=0, domain=None):
+        """ masked_unary_operation(aufunc, fill=0, domain=None)
+            aufunc(fill) must be defined
+            self(x) returns aufunc(x)
+            with masked values where domain(x) is true or getmask(x) is true.
+        """
+        self.f = 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."
+# numeric tries to return scalars rather than arrays when given scalars.
+        m = getmask(a)
+        d1 = filled(a, self.fill)
+        if self.domain is not None:
+            m = mask_or(m, numeric.asarray(self.domain(d1)))
+        # Take care of the masked singletong first ...
+        if m.ndim == 0 and m:
+            return masked
+        # Get the result....
+        if isinstance(a, MaskedArray):
+            result = self.f(d1, *args, **kwargs).view(type(a))
+        else:
+            result = self.f(d1, *args, **kwargs).view(MaskedArray)
+        # Fix the mask if we don't have a scalar
+        if result.ndim > 0:
+            result._mask = m
+        return result
+    #
+    def __str__ (self):
+        return "Masked version of %s. [Invalid values are masked]" % str(self.f)
+#..............................................................................
+class masked_binary_operation:
+    """Defines masked version of binary operations,
+where invalid values are pre-masked.
+
+:IVariables:
+    - `f` : function.
+    - `fillx` : Default filling value for first array*[0]*.
+    - `filly` : Default filling value for second array*[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))
+        if (not m.ndim) and m:
+            return masked
+        d1 = filled(a, self.fillx)
+        d2 = filled(b, self.filly)
+# CHECK : Do we really need to fill the arguments ? Pro'ly not        
+#        result = self.f(a, b, *args, **kwargs).view(get_masked_subclass(a,b))
+        result = self.f(d1, d2, *args, **kwargs).view(get_masked_subclass(a,b))
+        if result.ndim > 0:
+            result._mask = m
+        return result
+    #
+    def reduce (self, target, axis=0, dtype=None):
+        """Reduces `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
+        # XXX: "or t.dtype" below is a workaround for what appears
+        # XXX: to be a bug in reduce.
+        tr = self.f.reduce(filled(t, self.filly), 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):
+        "Returns 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)
+        d = self.f.outer(filled(a, self.fillx), filled(b, self.filly)).view(rcls)
+        if d.ndim > 0:
+            d._mask = m
+        return d
+
+    def accumulate (self, target, axis=0):
+        """Accumulates `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 domained_binary_operation:
+    """Defines binary operations that have a domain, like divide.
+
+These are complicated so they are a separate class.
+They have no reduce, outer or accumulate.
+
+:IVariables:
+    - `f` : function.
+    - `fillx` : Default filling value for first array*[0]*.
+    - `filly` : Default filling value for second array*[0]*.
+    - `domain` : Default domain *[None]*.
+    """
+    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 = filled(a, self.fillx)
+        d2 = filled(b, self.filly)
+        t = numeric.asarray(self.domain(d1, d2))
+
+        if fromnumeric.sometrue(t, None):
+            d2 = numeric.where(t, self.filly, d2)
+            mb = mask_or(mb, t)
+        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
+        return result
+
+    def __str__ (self):
+        return "Masked version of " + str(self.f)
+
+#..............................................................................
+# Unary ufuncs
+exp = masked_unary_operation(umath.exp)
+conjugate = masked_unary_operation(umath.conjugate)
+sin = masked_unary_operation(umath.sin)
+cos = masked_unary_operation(umath.cos)
+tan = masked_unary_operation(umath.tan)
+arctan = masked_unary_operation(umath.arctan)
+arcsinh = masked_unary_operation(umath.arcsinh)
+sinh = masked_unary_operation(umath.sinh)
+cosh = masked_unary_operation(umath.cosh)
+tanh = masked_unary_operation(umath.tanh)
+abs = absolute = masked_unary_operation(umath.absolute)
+fabs = masked_unary_operation(umath.fabs)
+negative = masked_unary_operation(umath.negative)
+floor = masked_unary_operation(umath.floor)
+ceil = masked_unary_operation(umath.ceil)
+around = masked_unary_operation(fromnumeric.round_)
+logical_not = masked_unary_operation(umath.logical_not)
+# Domained unary ufuncs
+sqrt = masked_unary_operation(umath.sqrt, 0.0, domain_greater_equal(0.0))
+log = masked_unary_operation(umath.log, 1.0, domain_greater(0.0))
+log10 = masked_unary_operation(umath.log10, 1.0, domain_greater(0.0))
+tan = masked_unary_operation(umath.tan, 0.0, domain_tan(1.e-35))
+arcsin = masked_unary_operation(umath.arcsin, 0.0,
+                                domain_check_interval(-1.0, 1.0))
+arccos = masked_unary_operation(umath.arccos, 0.0,
+                                domain_check_interval(-1.0, 1.0))
+arccosh = masked_unary_operation(umath.arccosh, 1.0, domain_greater_equal(1.0))
+arctanh = masked_unary_operation(umath.arctanh, 0.0,
+                                 domain_check_interval(-1.0+1e-15, 1.0-1e-15))
+# Binary ufuncs
+add = masked_binary_operation(umath.add)
+subtract = masked_binary_operation(umath.subtract)
+multiply = masked_binary_operation(umath.multiply, 1, 1)
+arctan2 = masked_binary_operation(umath.arctan2, 0.0, 1.0)
+equal = masked_binary_operation(umath.equal)
+equal.reduce = None
+not_equal = masked_binary_operation(umath.not_equal)
+not_equal.reduce = None
+less_equal = masked_binary_operation(umath.less_equal)
+less_equal.reduce = None
+greater_equal = masked_binary_operation(umath.greater_equal)
+greater_equal.reduce = None
+less = masked_binary_operation(umath.less)
+less.reduce = None
+greater = masked_binary_operation(umath.greater)
+greater.reduce = None
+logical_and = masked_binary_operation(umath.logical_and)
+alltrue = masked_binary_operation(umath.logical_and, 1, 1).reduce
+logical_or = masked_binary_operation(umath.logical_or)
+sometrue = logical_or.reduce
+logical_xor = masked_binary_operation(umath.logical_xor)
+bitwise_and = masked_binary_operation(umath.bitwise_and)
+bitwise_or = masked_binary_operation(umath.bitwise_or)
+bitwise_xor = masked_binary_operation(umath.bitwise_xor)
+hypot = masked_binary_operation(umath.hypot)
+# Domained binary ufuncs
+divide = domained_binary_operation(umath.divide, domain_safe_divide(), 0, 1)
+true_divide = domained_binary_operation(umath.true_divide,
+                                        domain_safe_divide(), 0, 1)
+floor_divide = domained_binary_operation(umath.floor_divide,
+                                         domain_safe_divide(), 0, 1)
+remainder = domained_binary_operation(umath.remainder,
+                                      domain_safe_divide(), 0, 1)
+fmod = domained_binary_operation(umath.fmod, domain_safe_divide(), 0, 1)
+
+
+#####--------------------------------------------------------------------------
+#---- --- Mask creation functions ---
+#####--------------------------------------------------------------------------
+def getmask(a):
+    """Returns the mask of `a`, if any, or `nomask`.
+Returns `nomask` if `a` is not a masked array.
+To get an array for sure use getmaskarray."""
+    if hasattr(a, "_mask"):
+        return a._mask
+    else:
+        return nomask
+
+def getmaskarray(a):
+    """Returns the mask of `a`, if any.
+Otherwise, returns an array of `False`, with the same shape as `a`.
+    """
+    m = getmask(a)
+    if m is nomask:
+        return make_mask_none(fromnumeric.shape(a))
+    else:
+        return m
+
+def is_mask(m):
+    """Returns `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, small_mask=True, flag=None):
+    """make_mask(m, copy=0, small_mask=0)
+Returns `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.
+If `small_mask=True`, returns `nomask` if `m` contains no true elements.
+
+:Parameters:
+    - `m` (ndarray) : Mask.
+    - `copy` (boolean, *[False]*) : Returns a copy of `m` if true.
+    - `small_mask` (boolean, *[False]*): Flattens mask to `nomask` if `m` is all false.
+    """
+    if flag is not None:
+        warnings.warn("The flag 'flag' is now called 'small_mask'!",
+                      DeprecationWarning)
+        small_mask = flag
+    if m is nomask:
+        return nomask
+    elif isinstance(m, ndarray):
+        m = filled(m, True)
+        if m.dtype.type is MaskType:
+            if copy:
+                result = numeric.array(m, dtype=MaskType, copy=copy)
+            else:
+                result = m
+        else:
+            result = numeric.array(m, dtype=MaskType)
+    else:
+        result = numeric.array(filled(m, True), dtype=MaskType)
+    # Bas les masques !
+    if small_mask and not result.any():
+        return nomask
+    else:
+        return result
+
+def make_mask_none(s):
+    "Returns a mask of shape `s`, filled with `False`."
+    result = numeric.zeros(s, dtype=MaskType)
+    return result
+
+def mask_or (m1, m2, copy=False, small_mask=True):
+    """Returns 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:
+    - `m` (ndarray) : Mask.
+    - `copy` (boolean, *[False]*) : Returns a copy of `m` if true.
+    - `small_mask` (boolean, *[False]*): Flattens mask to `nomask` if `m` is all false.
+     """
+    if m1 is nomask:
+        return make_mask(m2, copy=copy, small_mask=small_mask)
+    if m2 is nomask:
+        return make_mask(m1, copy=copy, small_mask=small_mask)
+    if m1 is m2 and is_mask(m1):
+        return m1
+    return make_mask(umath.logical_or(m1, m2), copy=copy, small_mask=small_mask)
+
+#####--------------------------------------------------------------------------
+#--- --- Masking functions ---
+#####--------------------------------------------------------------------------
+def masked_where(condition, a, copy=True):
+    """Returns `x` as an array masked where `condition` is true.
+Masked values of `x` or `condition` are kept.
+
+:Parameters:
+    - `condition` (ndarray) : Masking condition.
+    - `x` (ndarray) : Array to mask.
+    - `copy` (boolean, *[False]*) : Returns a copy of `m` if true.
+    """
+    cond = filled(condition,1)
+    a = numeric.array(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=1):
+    "Shortcut to `masked_where`, with ``condition = (x > value)``."
+    return masked_where(greater(x, value), x, copy=copy)
+
+def masked_greater_equal(x, value, copy=1):
+    "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.
+    """
+    return masked_where((x == value), x, copy=copy)
+#    d = filled(x, 0)
+#    c = umath.equal(d, value)
+#    m = mask_or(c, getmask(x))
+#    return array(d, mask=m, copy=copy)
+
+def masked_inside(x, v1, v2, copy=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.
+    """
+    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.
+    """
+    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):
+    """Masks 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.
+The mask is set to `nomask` if posible.
+
+:parameter copy (Boolean, *[True]*):  Returns a copy of `x` if true. """
+    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, small_mask=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):
+    """Masks the array `x` where the data are approximately equal to `value`
+(that is, ``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:
+    - `rtol` (Float, *[1e-5]*): Tolerance parameter.
+    - `atol` (Float, *[1e-8]*): Tolerance parameter.
+    - `copy` (boolean, *[False]*) : Returns a copy of `x` if True.
+    """
+    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))
+        try:
+            mask = x._mask
+        except AttributeError:
+            mask = nomask
+    else:
+        condition = umath.equal(xnew, value)
+        mask = nomask
+    mask = mask_or(mask, make_mask(condition, small_mask=True))
+    return masked_array(xnew, mask=mask, copy=copy, fill_value=value)
+
+#####--------------------------------------------------------------------------
+#---- --- Printing options ---
+#####--------------------------------------------------------------------------
+class _MaskedPrintOption:
+    """Handles the string used to represent missing data in a masked array."""
+    def __init__ (self, display):
+        "Creates the masked_print_option object."
+        self._display = display
+        self._enabled = True
+
+    def display(self):
+        "Displays the string to print for masked values."
+        return self._display
+
+    def set_display (self, s):
+        "Sets 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, small_mask=1):
+        "Set the enabling small_mask to `small_mask`."
+        self._enabled = small_mask
+
+    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 ---
+#####--------------------------------------------------------------------------
+##def _getoptions(a_out, a_in):
+##    "Copies standards options of a_in to a_out."
+##    for att in [']
+#class _mathmethod(object):
+#    """Defines a wrapper for arithmetic methods.
+#Instead of directly calling a ufunc, the corresponding method of  the `array._data`
+#object is called instead.
+#    """
+#    def __init__ (self, methodname, fill_self=0, fill_other=0, domain=None):
+#        """
+#:Parameters:
+#    - `methodname` (String) : Method name.
+#    - `fill_self` (Float *[0]*) : Fill value for the instance.
+#    - `fill_other` (Float *[0]*) : Fill value for the target.
+#    - `domain` (Domain object *[None]*) : Domain of non-validity.
+#        """
+#        self.methodname = methodname
+#        self.fill_self = fill_self
+#        self.fill_other = fill_other
+#        self.domain = domain
+#        self.obj = None
+#        self.__doc__ = self.getdoc()
+#    #
+#    def getdoc(self):
+#        "Returns the doc of the function (from the doc of the method)."
+#        try:
+#            return getattr(MaskedArray, self.methodname).__doc__
+#        except:
+#            return getattr(ndarray, self.methodname).__doc__
+#    #
+#    def __get__(self, obj, objtype=None):
+#        self.obj = obj
+#        return self
+#    #
+#    def __call__ (self, other, *args):
+#        "Execute the call behavior."
+#        instance = self.obj
+#        m_self = instance._mask
+#        m_other = getmask(other)
+#        base = instance.filled(self.fill_self)
+#        target = filled(other, self.fill_other)
+#        if self.domain is not None:
+#            # We need to force the domain to a ndarray only.
+#            if self.fill_other > self.fill_self:
+#                domain = self.domain(base, target)
+#            else:
+#                domain = self.domain(target, base)
+#            if domain.any():
+#                #If `other` is a subclass of ndarray, `filled` must have the
+#                # same subclass, else we'll lose some info.
+#                #The easiest then is to fill `target` instead of creating
+#                # a pure ndarray.
+#                #Oh, and we better make a copy!
+#                if isinstance(other, ndarray):
+#                    # We don't want to modify other: let's copy target, then
+#                    target = target.copy()
+#                    target[fromnumeric.asarray(domain)] = self.fill_other
+#                else:
+#                    target = numeric.where(fromnumeric.asarray(domain),
+#                                           self.fill_other, target)
+#                m_other = mask_or(m_other, domain)
+#        m = mask_or(m_self, m_other)
+#        method = getattr(base, self.methodname)
+#        result = method(target, *args).view(type(instance))
+#        try:
+#            result._mask = m
+#        except AttributeError:
+#            if m:
+#                result = masked
+#        return result
+#...............................................................................
+class _arraymethod(object):
+    """Defines 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 calld on the initial mask.
+If `onmask` is False, the new mask is just a reference to the initial mask.
+
+:Parameters:
+    `funcname` : String
+        Name of the function to apply on data.
+    `onmask` : Boolean *[True]*
+        Whether the mask must be processed also (True) or left alone (False).
+    """
+    def __init__(self, funcname, onmask=True):
+        self._name = funcname
+        self._onmask = onmask
+        self.obj = None
+        self.__doc__ = self.getdoc()
+    #
+    def getdoc(self):
+        "Returns the doc of the function (from the doc of the method)."
+        methdoc = getattr(ndarray, self._name, None)
+        methdoc = getattr(numpy, self._name, methdoc)
+#        methdoc = getattr(MaskedArray, self._name, methdoc)
+        if methdoc is not None:
+            return methdoc.__doc__
+#        try:
+#            return getattr(MaskedArray, self._name).__doc__
+#        except:
+#            try:
+#                return getattr(numpy, self._name).__doc__
+#            except:
+#                return getattr(ndarray, self._name).__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._smallmask = self.obj._smallmask
+        if result.ndim:
+            if not self._onmask:
+                result._mask = mask
+            elif mask is not nomask:
+                result.__setmask__(getattr(mask, methodname)(*args, **params))
+        return result
+#..........................................................
+
+class flatiter(object):
+    "Defines 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 = array(data, dtype=None, copy=True, order=False,
+              mask = nomask, fill_value=None, small_mask=True)
+
+If copy=False, every effort is made not to copy the data:
+If `data` is a MaskedArray, and argument mask=nomask, then the candidate data
+is `data._data` and the mask used is `data._mask`.
+If `data` is a numeric array, it is used as the candidate raw data.
+If `dtype` is not None and is different from data.dtype.char then a data copy is required.
+Otherwise, the candidate is used.
+
+If a data copy is required, the raw (unmasked) data stored is the result of:
+numeric.array(data, dtype=dtype.char, copy=copy)
+
+If `mask` is `nomask` there are no masked values.
+Otherwise mask must be convertible to an array of booleans with the same shape as x.
+If `small_mask` is True, a mask consisting of zeros (False) only is compressed to `nomask`.
+Otherwise, the mask is not compressed.
+
+fill_value is used to fill in masked values when necessary, such as when
+printing and in method/function filled().
+The fill_value is not used for computation within this module.
+    """
+    __array_priority__ = 10.1
+    _defaultmask = nomask
+    _defaulthardmask = False
+    _baseclass =  numeric.ndarray
+    def __new__(cls, data=None, mask=nomask, dtype=None, copy=False, fill_value=None,
+                keep_mask=True, small_mask=True, hard_mask=False, flag=None,
+                subok=True, **options):
+        """array(data, dtype=None, copy=True, mask=nomask, fill_value=None)
+
+If `data` is already a ndarray, its dtype becomes the default value of dtype.
+        """
+        if flag is not None:
+            warnings.warn("The flag 'flag' is now called 'small_mask'!",
+                          DeprecationWarning)
+            small_mask = flag
+        # Process data............
+        _data = numeric.array(data, dtype=dtype, copy=copy, subok=subok)
+        _baseclass = getattr(data, '_baseclass', type(_data))
+        _basedict = getattr(data, '_basedict', getattr(data, '__dict__', None))
+        if not isinstance(data, MaskedArray): 
+            _data = _data.view(cls)
+        elif not subok:
+            _data = data.view(cls)
+        else:
+            _data = _data.view(type(data))
+        # Backwards compat .......
+        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:
+                _data._mask = nomask
+            if copy:
+                _data._mask = _data._mask.copy()
+        else:
+            mask = numeric.array(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)
+            if _data._mask is nomask:
+                _data._mask = mask
+                _data._sharedmask = True
+            else:
+                # Make a copy of the mask to avoid propagation
+                _data._sharedmask = False
+                if not keep_mask:
+                    _data._mask = mask
+                else:
+                    _data._mask = umath.logical_or(mask, _data._mask) 
+                    
+                    
+        # Update fill_value.......
+        _data._fill_value = getattr(data, '_fill_value', fill_value)
+        if _data._fill_value is None:
+            _data._fill_value = default_fill_value(_data)
+        # Process extra options ..
+        _data._hardmask = hard_mask
+        _data._smallmask = small_mask
+        _data._baseclass = _baseclass
+        _data._basedict = _basedict
+        return _data
+    #........................
+    def __array_finalize__(self,obj):
+        """Finalizes the masked array.
+        """
+        # Finalize mask ...............
+        self._mask = getattr(obj, '_mask', nomask)
+        if self._mask is not nomask:
+            self._mask.shape = self.shape
+        # Get the remaining options ...
+        self._hardmask = getattr(obj, '_hardmask', self._defaulthardmask)
+        self._smallmask = getattr(obj, '_smallmask', True)
+        self._sharedmask = True
+        self._baseclass = getattr(obj, '_baseclass', type(obj))
+        self._fill_value = getattr(obj, '_fill_value', None)
+        # Update special attributes ...
+        self._basedict = getattr(obj, '_basedict', getattr(obj, '__dict__', None))
+        if self._basedict is not None:
+            self.__dict__.update(self._basedict)
+        return
+    #..................................
+    def __array_wrap__(self, obj, context=None):
+        """Special hook for ufuncs.
+Wraps the numpy array and sets the mask according to context.
+        """
+        #TODO : Should we check for type result 
+        result = obj.view(type(self))
+        #..........
+        if context is not None:
+            result._mask = result._mask.copy()
+            (func, args, _) = context
+            m = reduce(mask_or, [getmask(arg) for arg in args])
+            # Get 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)
+                if m is nomask:
+                    if d is not nomask:
+                        m = d
+                else:
+                    m |= d
+            if not m.ndim and m:
+                if m:
+                    if result.shape == ():
+                        return masked
+                    result._mask = numeric.ones(result.shape, bool_)
+            else:
+                result._mask = m
+        #....
+#        result._mask = m
+        result._fill_value = self._fill_value
+        result._hardmask = self._hardmask
+        result._smallmask = self._smallmask
+        result._baseclass = self._baseclass
+        return result
+    #.............................................
+    def __getitem__(self, indx):
+        """x.__getitem__(y) <==> x[y]
+Returns the item described by i. Not a copy as in previous versions.
+        """
+        # 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
+        # super() can't work here if the underlying data is a matrix...
+        dout = (self._data).__getitem__(indx)
+        m = self._mask
+        if hasattr(dout, 'shape') and len(dout.shape) > 0:
+            # Not a scalar: make sure that dout is a MA
+            dout = dout.view(type(self))
+            dout._smallmask = self._smallmask
+            if m is not nomask:
+                # use _set_mask to take care of the shape
+                dout.__setmask__(m[indx])
+        elif m is not nomask and m[indx]:
+            return masked
+        return dout
+    #........................
+    def __setitem__(self, indx, value):
+        """x.__setitem__(i, y) <==> x[i]=y
+Sets item described by index. If value is masked, masks those locations.
+        """
+        if self is masked:
+            raise MAError, 'Cannot alter the masked element.'
+#        if getmask(indx) is not nomask:
+#            msg = "Masked arrays must be filled before they can be used as indices!"
+#            raise IndexError, msg
+        #....
+        if value is masked:
+            m = self._mask
+            if m is nomask:
+                m = make_mask_none(self.shape)
+#            else:
+#                m = m.copy()
+            m[indx] = True
+            self.__setmask__(m)
+            return
+        #....
+        dval = numeric.asarray(value).astype(self.dtype)
+        valmask = getmask(value)
+        if self._mask is nomask:
+            if valmask is not nomask:
+                self._mask = make_mask_none(self.shape)
+                self._mask[indx] = valmask
+        elif not self._hardmask:
+            _mask = self._mask.copy()
+            if valmask is nomask:
+                _mask[indx] = False
+            else:
+                _mask[indx] = valmask
+            self._set_mask(_mask)
+        elif hasattr(indx, 'dtype') and (indx.dtype==bool_):
+            indx = indx * umath.logical_not(self._mask)
+        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
+            dval = dindx
+            self._mask[indx] = mindx
+        # Set data ..........
+        #dval = filled(value).astype(self.dtype)
+        ndarray.__setitem__(self._data,indx,dval)
+    #............................................
+    def __getslice__(self, i, j):
+        """x.__getslice__(i, j) <==> x[i:j]
+Returns 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
+Sets a slice i:j to `value`.
+If `value` is masked, masks those locations."""
+        self.__setitem__(slice(i,j), value)
+    #............................................
+    def __setmask__(self, mask, copy=False):
+        newmask = make_mask(mask, copy=copy, small_mask=self._smallmask)
+#        self.unshare_mask()
+        if self._mask is nomask:
+            self._mask = newmask
+        elif self._hardmask:
+            if newmask is not nomask:
+                self._mask.__ior__(newmask)
+        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 newmask is nomask:
+                self._mask = nomask
+            else:
+                self._mask.flat = newmask
+        if self._mask.shape:
+            self._mask = numeric.reshape(self._mask, self.shape)
+    _set_mask = __setmask__
+    
+    def _get_mask(self):
+        """Returns the current mask."""
+        return self._mask
+
+    mask = property(fget=_get_mask, fset=__setmask__, doc="Mask")
+    #............................................
+    def harden_mask(self):
+        "Forces the mask to hard."
+        self._hardmask = True
+        
+    def soften_mask(self):
+        "Forces the mask to soft."
+        self._hardmask = False     
+        
+    def unshare_mask(self):
+        "Copies the mask and set the sharedmask flag to False."
+        if self._sharedmask:
+            self._mask = self._mask.copy()
+            self._sharedmask = False
+        
+    #............................................
+    def _get_data(self):
+        "Returns the current data (as a view of the original underlying data)>"
+        return self.view(self._baseclass)
+    _data = property(fget=_get_data)        
+    #............................................
+    def _get_flat(self):
+        """Calculates the flat value.
+        """
+        return flatiter(self)
+    #
+    def _set_flat (self, value):
+        "x.flat = value"
+        y = self.ravel()
+        y[:] = value
+    #
+    flat = property(fget=_get_flat, fset=_set_flat, doc="Flat version")
+    #............................................
+    def get_fill_value(self):
+        "Returns the filling value."
+        if self._fill_value is None:
+            self._fill_value = default_fill_value(self)
+        return self._fill_value
+
+    def set_fill_value(self, value=None):
+        """Sets the filling value to `value`.
+If None, uses the default, based on the data type."""
+        if value is None:
+            value = default_fill_value(self)
+        self._fill_value = value
+
+    fill_value = property(fget=get_fill_value, fset=set_fill_value,
+                          doc="Filling value")
+
+    def filled(self, fill_value=None):
+        """Returns an array of the same class as `_data`,
+ with masked values filled with `fill_value`.
+Subclassing is preserved.
+
+If `fill_value` is None, uses self.fill_value.
+        """
+        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)
+                #result[m] = fill_value
+            except (TypeError, AttributeError):
+                fill_value = numeric.array(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 = numeric.array(fill_value, dtype=self.dtype)
+                else:
+                    result = self._data
+        return result
+
+    def compressed(self):
+        "A 1-D array of all the non-masked data."
+        d = self.ravel()
+        if self._mask is nomask:
+            return d
+        elif not self._smallmask and not self._mask.any():
+            return d
+        else:
+            return d[numeric.logical_not(d._mask)]
+    #............................................
+    def __str__(self):
+        """x.__str__() <==> str(x)
+Calculates the string representation, using masked for fill if it is enabled.
+Otherwise, fills with fill value.
+        """
+        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):
+        """x.__repr__() <==> repr(x)
+Calculates the repr representation, using masked for fill if it is enabled.
+Otherwise fill with fill value.
+        """
+        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 __iadd__(self, other):
+        "Adds other to self in place."
+        ndarray.__iadd__(self._data,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):
+        "Subtracts other from self in place."
+        ndarray.__isub__(self._data,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):
+        "Multiplies self by other in place."
+        ndarray.__imul__(self._data,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):
+        "Divides self by other in place."
+        dom_mask = domain_safe_divide().__call__(self, filled(other,1))
+        other_mask = getmask(other)
+        new_mask = mask_or(other_mask, dom_mask)
+        ndarray.__idiv__(self._data, other)
+        self._mask = mask_or(self._mask, new_mask)
+        return self
+    #............................................
+    def __float__(self):
+        "Converts self to float."
+        if self._mask is not nomask:
+            warnings.warn("Warning: converting a masked element to nan.")
+            return numpy.nan
+            #raise MAError, 'Cannot convert masked element to a Python float.'
+        return float(self.item())
+
+    def __int__(self):
+        "Converts self to int."
+        if self._mask is not nomask:
+            raise MAError, 'Cannot convert masked element to a Python int.'
+        return int(self.item())
+    #............................................
+    def count(self, axis=None):
+        """Counts the non-masked elements of the array along a given axis,
+and returns a masked array where the mask is True where all data are masked.
+If `axis` is None, counts all the non-masked elements, and returns either a
+scalar or the masked singleton."""
+        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 = fromnumeric.size(m, axis)
+        n2 = m.astype(int_).sum(axis)
+        if axis is None:
+            return (n1-n2)
+        else:
+            return masked_array(n1 - n2)
+    #............................................
+    def reshape (self, *s):
+        """Reshapes the array to shape s.
+Returns a new masked array.
+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
+    #
+    repeat = _arraymethod('repeat')
+    #
+    def resize(self, newshape, refcheck=True, order=False):
+        """Attempts to modify size and shape of self inplace.
+        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
+    #
+    flatten = _arraymethod('flatten')
+    #
+    def put(self, indices, values, mode='raise'):
+        """Sets storage-indexed locations to corresponding values.
+a.put(values, indices, mode) sets a.flat[n] = values[n] for each n in indices.
+`values` can be scalar or an array shorter than indices, and it will be repeated,
+if necessary.
+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 = numeric.asarray(indices)
+            values = numeric.asanyarray(values)
+            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, small_mask=True)
+        self._mask = m
+    #............................................
+    def ids (self):
+        """Return the address of the data and mask areas."""
+        return (self.ctypes.data, self._mask.ctypes.data)    
+    #............................................
+    def all(self, axis=None, out=None):
+        """a.all(axis) returns True if all entries along the axis are True.
+    Returns False otherwise. If axis is None, uses the flatten array.
+    Masked data are considered as True during computation.
+    Outputs a masked array, where the mask is True if all data are masked along the axis.
+    Note: the out argument is not really operational...
+        """
+        d = self.filled(True).all(axis=axis, out=out).view(type(self))
+        if d.ndim > 0:
+            d.__setmask__(self._mask.all(axis))
+        return d
+
+    def any(self, axis=None, out=None):
+        """a.any(axis) returns True if some or all entries along the axis are True.
+    Returns False otherwise. If axis is None, uses the flatten array.
+    Masked data are considered as False during computation.
+    Outputs a masked array, where the mask is True if all data are masked along the axis.
+    Note: the out argument is not really operational...
+        """
+        d = self.filled(False).any(axis=axis, out=out).view(type(self))
+        if d.ndim > 0:
+            d.__setmask__(self._mask.all(axis))
+        return d
+    
+    def nonzero(self):
+        """a.nonzero() returns a tuple of arrays
+
+    Returns a tuple of arrays, one for each dimension of a,
+    containing 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
+        transpose(a.nonzero())
+    instead. The result of this is always a 2d array, with a row for
+    each non-zero element."""
+        return numeric.asarray(self.filled(0)).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)
+Returns 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).sum(axis=None)
+    #............................................
+    def sum(self, axis=None, dtype=None):
+        """a.sum(axis=None, dtype=None)
+Sums the array `a` over the given axis `axis`.
+Masked values are set to 0.
+If `axis` is None, applies to a flattened version of the array.
+    """
+        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):
+        """a.cumprod(axis=None, dtype=None)
+Returns the cumulative sum of the elements of array `a` along the given axis `axis`.
+Masked values are set to 0.
+If `axis` is None, applies to a flattened version of the array.
+        """
+        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):
+        """a.prod(axis=None, dtype=None)
+Returns the product of the elements of array `a` along the given axis `axis`.
+Masked elements are set to 1.
+If `axis` is None, applies to a flattened version of the array.
+        """
+        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):
+        """a.cumprod(axis=None, dtype=None)
+Returns the cumulative product of ethe lements of array `a` along the given axis `axis`.
+Masked values are set to 1.
+If `axis` is None, applies to a flattened version of the array.
+        """
+        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):
+        """a.mean(axis=None, dtype=None)
+
+    Averages the array over the given axis.  If the axis is None,
+    averages over all dimensions of the array.  Equivalent to
+
+      a.sum(axis, dtype) / size(a, axis).
+
+    The optional dtype argument is the data type for intermediate
+    calculations in the sum.
+
+    Returns a masked array, of the same class as a.
+        """
+        if self._mask is nomask:
+            return super(MaskedArray, self).mean(axis=axis, dtype=dtype)
+        else:
+            dsum = self.sum(axis=axis, dtype=dtype)
+            cnt = self.count(axis=axis)
+            return dsum*1./cnt
+
+    def anom(self, axis=None, dtype=None):
+        """a.anom(axis=None, dtype=None)
+    Returns the anomalies, or deviation from the average.
+            """
+        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):
+        """a.var(axis=None, dtype=None)
+Returns 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((x - x.mean())**2).
+        """
+        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)
+        else:
+            cnt = self.count(axis=axis)
+            danom = self.anom(axis=axis, dtype=dtype)
+            danom *= danom
+            dvar = numeric.array(danom.sum(axis) / cnt).view(type(self))
+            if axis is not None:
+                dvar._mask = mask_or(self._mask.all(axis), (cnt==1))
+            return dvar
+
+    def std(self, axis=None, dtype=None):
+        """a.std(axis=None, dtype=None)
+Returns 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((x - x.mean())**2)).
+        """
+        dvar = self.var(axis,dtype)
+        if axis is not None or dvar is not masked:
+            dvar = sqrt(dvar)
+        return dvar
+    #............................................
+    def argsort(self, axis=None, fill_value=None, kind='quicksort',
+                order=None):
+        """Returns an array of indices that sort 'a' along the specified axis.
+    Masked values are filled beforehand to `fill_value`.
+    If `fill_value` is None, uses the default for the data type.
+    Returns a numpy array.
+
+:Keywords:
+    `axis` : Integer *[None]*
+        Axis to be indirectly sorted (default -1)
+    `kind` : String *['quicksort']*
+        Sorting algorithm (default 'quicksort')
+        Possible values: 'quicksort', 'mergesort', or 'heapsort'
+
+    Returns: array of indices that sort 'a' along the specified axis.
+
+    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):
+        """Returns a 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`.
+    If `fill_value` is None, the default for the data type is used.
+    Returns a numpy array.
+
+:Keywords:
+    `axis` : Integer *[None]*
+        Axis to be indirectly sorted (default -1)
+    `fill_value` : var *[None]*
+        Default filling value. If None, uses the minimum default for the data type.
+        """
+        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`.
+    If `fill_value` is None, the maximum default for the data type is used.
+    Returns a numpy array.
+
+:Keywords:
+    `axis` : Integer *[None]*
+        Axis to be indirectly sorted (default -1)
+    `fill_value` : var *[None]*
+        Default filling value. If None, uses the data type default.
+        """
+        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 a along the given axis.
+
+    Keyword arguments:
+
+    axis  -- axis to be sorted (default -1)
+    kind  -- sorting algorithm (default 'quicksort')
+             Possible values: 'quicksort', 'mergesort', or 'heapsort'.
+    order -- 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.
+    endwith--Boolean flag indicating whether missing values (if any) should
+             be forced in the upper indices (at the end of the array) or
+             lower indices (at the beginning).
+
+    Returns: None.
+
+    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 with 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):
+        """Returns the minimum/a along the given axis.
+If `axis` is None, applies to the flattened array. Masked values are filled 
+with `fill_value` during processing. If `fill_value is None, it is set to the
+maximum_fill_value corresponding to the data type."""
+        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)
+        # Get the fil 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 max(self, axis=None, fill_value=None):
+        """Returns the maximum/a along the given axis.
+If `axis` is None, applies to the flattened array. Masked values are filled 
+with `fill_value` during processing. If `fill_value is None, it is set to the
+maximum_fill_value corresponding to the data type."""
+        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)
+        # 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):
+        """Returns the visible data range (max-min) along the given axis.
+If the axis is `None`, applies on a flattened array. Masked values are filled
+with `fill_value` for processing. If `fill_value` is None, the maximum is uses
+the maximum default, the minimum uses the minimum default."""
+        return self.max(axis, fill_value) - self.min(axis, fill_value)
+
+    # Array methods ---------------------------------------
+    conj = conjugate = _arraymethod('conjugate')
+    copy = _arraymethod('copy')
+    diagonal = _arraymethod('diagonal')
+    take = _arraymethod('take')
+    ravel = _arraymethod('ravel')
+    transpose = _arraymethod('transpose')
+    T = property(fget=lambda self:self.transpose())
+    swapaxes = _arraymethod('swapaxes')
+    clip = _arraymethod('clip', onmask=False)
+    compress = _arraymethod('compress')
+    copy = _arraymethod('copy')
+    squeeze = _arraymethod('squeeze')
+    #--------------------------------------------
+    def tolist(self, fill_value=None):
+        """Copies 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()
+        if self._mask is nomask:
+            return result
+        if self.ndim == 0:
+            return [None]
+        elif self.ndim == 1:
+            maskedidx = self._mask.nonzero()[0].tolist()
+            [operator.setitem(result,i,None) for i in maskedidx]
+        else:
+            for idx in zip(*[i.tolist() for i in self._mask.nonzero()]):
+                tmp = result
+                for i in idx[:-1]:
+                    tmp = tmp[i]
+                tmp[idx[-1]] = None
+        return result
+            
+            
+    #........................
+    def tostring(self, fill_value=None):
+        """a.tostring(order='C', fill_value=None) -> raw copy of array data as a Python string.
+
+    Keyword arguments:
+        order      : order of the data item in the copy {"C","F","A"} (default "C")
+        fill_value : value used in lieu of missing data 
+
+    Construct a Python string containing the raw bytes in the array. The order
+    of the data in arrays with ndim > 1 is specified by the 'order' keyword and
+    this keyword overrides the order of the array. The
+    choices are:
+
+        "C"       -- C order (row major)
+        "Fortran" -- Fortran order (column major)
+        "Any"     -- Current order of array.
+        None      -- Same as "Any"
+    
+    Masked data are filled with fill_value. If fill_value is None, the data-type-
+    dependent default is used."""
+        return self.filled(fill_value).tostring()   
+    #--------------------------------------------
+    # Backwards Compatibility. Heck...
+    @property
+    def data(self):
+        """Returns the `_data` part of the MaskedArray."""
+        return self._data
+    def raw_data(self):
+        """Returns the `_data` part of the MaskedArray.
+You should really use `data` instead..."""
+        return self._data
+    #--------------------------------------------
+    # Pickling
+    def __getstate__(self):
+        "Returns 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):
+        """Restores 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):
+        """Returns 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, small_mask=False)
+#MaskedArray.__dump__ = dump
+#MaskedArray.__dumps__ = dumps
+    
+    
+
+#####--------------------------------------------------------------------------
+#---- --- 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
+#masked = MaskedArray(0, int, mask=1)
+masked_singleton = MaskedArray(0, dtype=int_, mask=True)
+masked = masked_singleton
+
+masked_array = MaskedArray
+def array(data, dtype=None, copy=False, order=False, mask=nomask, subok=True,
+          keep_mask=True, small_mask=True, hard_mask=None, fill_value=None):
+    """array(data, dtype=None, copy=True, order=False, mask=nomask,
+             keep_mask=True, small_mask=True, fill_value=None)
+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, small_mask=small_mask,
+                       hard_mask=hard_mask, fill_value=fill_value)
+
+def is_masked(x):
+    """Returns whether x has some 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):
+        """Reduces target along the given axis."""
+        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):
+        "Returns 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):
+    """Returns 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)
+#............................
+def max(obj, axis=None, out=None):
+    """Returns the maxima 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 maximum(obj)
+    else:
+        return maximum.reduce(obj, axis)
+#.............................
+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)
+
+
+#####---------------------------------------------------------------------------
+#---- --- Definition of functions from the corresponding methods ---
+#####---------------------------------------------------------------------------
+class _frommethod:
+    """Defines functions from existing MaskedArray methods.
+:ivar _methodname (String): Name of the method to transform.
+    """
+    def __init__(self, methodname):
+        self._methodname = methodname
+        self.__doc__ = self.getdoc()
+    def getdoc(self):
+        "Returns 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(fromnumeric.asarray(a), 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')
+std = _frommethod('std')
+sum = _frommethod('sum')
+swapaxes = _frommethod('swapaxes')
+take = _frommethod('take')
+var = _frommethod('var')
+
+#..............................................................................
+def power(a, b, third=None):
+    """Computes a**b elementwise.
+    Masked values are set to 1."""
+    if third is not None:
+        raise MAError, "3-argument power not supported."
+    ma = getmask(a)
+    mb = getmask(b)
+    m = mask_or(ma, mb)
+    fa = filled(a, 1)
+    fb = filled(b, 1)
+    if fb.dtype.char in typecodes["Integer"]:
+        return masked_array(umath.power(fa, fb), m)
+    md = make_mask((fa < 0), small_mask=1)
+    m = mask_or(m, md)
+    if m is nomask:
+        return masked_array(umath.power(fa, fb))
+    else:
+        fa[m] = 1
+        return masked_array(umath.power(fa, fb), m)
+
+#..............................................................................
+def argsort(a, axis=None, kind='quicksort', order=None, fill_value=None):
+    """Returns an array of indices that sort 'a' along the specified axis.
+    Masked values are filled beforehand to `fill_value`.
+    If `fill_value` is None, uses the default for the data type.
+    Returns a numpy array.
+
+:Keywords:
+    `axis` : Integer *[None]*
+        Axis to be indirectly sorted (default -1)
+    `kind` : String *['quicksort']*
+        Sorting algorithm (default 'quicksort')
+        Possible values: 'quicksort', 'mergesort', or 'heapsort'
+
+    Returns: array of indices that sort 'a' along the specified axis.
+
+    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(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)
+
+def argmin(a, axis=None, fill_value=None):
+    """Returns the array of indices for the minimum values of `a` along the
+    specified axis.
+    Masked values are treated as if they had the value `fill_value`.
+    If `fill_value` is None, the default for the data type is used.
+    Returns a numpy array.
+
+:Keywords:
+    `axis` : Integer *[None]*
+        Axis to be indirectly sorted (default -1)
+    `fill_value` : var *[None]*
+        Default filling value. If None, uses the data type default.
+    """
+    if fill_value is None:
+        fill_value = default_fill_value(a)
+    d = filled(a, fill_value)
+    return d.argmin(axis=axis)
+
+def argmax(a, 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`.
+    If `fill_value` is None, the default for the data type is used.
+    Returns a numpy array.
+
+:Keywords:
+    `axis` : Integer *[None]*
+        Axis to be indirectly sorted (default -1)
+    `fill_value` : var *[None]*
+        Default filling value. If None, uses the data type default.
+    """
+    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)
+
+def sort(a, axis=-1, kind='quicksort', order=None, endwith=True, fill_value=None):
+    """
+    Sort a along the given axis.
+
+Keyword arguments:
+
+axis  -- axis to be sorted (default -1)
+kind  -- sorting algorithm (default 'quicksort')
+         Possible values: 'quicksort', 'mergesort', or 'heapsort'.
+order -- 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.
+endwith--Boolean flag indicating whether missing values (if any) should
+         be forced in the upper indices (at the end of the array) or
+         lower indices (at the beginning).
+
+Returns: None.
+
+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 with 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  |
+|------------------------------------------------------|
+
+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.
+
+"""
+    a = numeric.asanyarray(a)
+    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]
+
+def compressed(x):
+    """Returns a compressed version of a masked array (or just the array if it
+    wasn't masked first)."""
+    if getmask(x) is None:
+        return x
+    else:
+        return x.compressed()
+
+def count(a, axis = None):
+    "Count of the non-masked elements in a, or along a certain axis."
+    a = masked_array(a)
+    return a.count(axis)
+
+def concatenate(arrays, axis=0):
+    "Concatenates the arrays along the given axis"
+    d = numeric.concatenate([filled(a) for a in arrays], axis)
+    rcls = get_masked_subclass(*arrays)
+    data = d.view(rcls)
+    for x in arrays:
+        if getmask(x) is not nomask:
+            break
+    else:
+        return data
+    dm = numeric.concatenate([getmaskarray(a) for a in arrays], axis)
+    dm = make_mask(dm, copy=False, small_mask=True)
+    data._mask = dm
+    return data
+
+def expand_dims(x,axis):
+    """Expand the shape of a by including newaxis before given axis."""
+    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'):
+    """Sets 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 fromnumeric.asarray(a).put(indices, values, mode=mode)
+
+def putmask(a, mask, values): #, mode='raise'):
+    """`putmask(a, mask, v)` results in `a = v` for all places where `mask` is true.
+If `v` is shorter than `mask`, it will be repeated as necessary.
+In particular `v` can be a scalar or length 1 array."""
+    # We can't use 'frommethod', the order of arguments is different
+    try:
+        return a.putmask(values, mask)
+    except AttributeError:
+        return fromnumeric.asarray(a).putmask(values, mask)
+
+def transpose(a,axes=None):
+    """Returns a view of the array with dimensions permuted according to axes.
+If `axes` is None (default), returns array with dimensions reversed.
+    """
+    #We can't use 'frommethod', as 'transpose' doesn't take keywords
+    try:
+        return a.transpose(axes)
+    except AttributeError:
+        return fromnumeric.asarray(a).transpose(axes)
+
+def reshape(a, new_shape):
+    """Changes 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 fromnumeric.asarray(a).reshape(new_shape)
+
+def resize(x, new_shape):
+    """resize(a,new_shape) returns 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 = fromnumeric.resize(m, new_shape)
+    result = fromnumeric.resize(x, new_shape).view(get_masked_subclass(x))
+    if result.ndim:
+        result._mask = m
+    return result
+
+
+#................................................
+def rank(obj):
+    """Gets the rank of sequence a (the number of dimensions, not a matrix rank)
+The rank of a scalar is zero."""
+    return fromnumeric.rank(filled(obj))
+#
+def shape(obj):
+    """Returns the shape of `a` (as a function call which also works on nested sequences).
+    """
+    return fromnumeric.shape(filled(obj))
+#
+def size(obj, axis=None):
+    """Returns the number of elements in the array along the given axis,
+or in the sequence if `axis` is None.
+    """
+    return fromnumeric.size(filled(obj), axis)
+#................................................
+
+#####--------------------------------------------------------------------------
+#---- --- Extra functions ---
+#####--------------------------------------------------------------------------
+def where (condition, x, y):
+    """where(condition, x, y) is x where condition is nonzero, y otherwise.
+       condition must be convertible to an integer array.
+       Answer is always the shape of condition.
+       The type depends on x and y. It is integer if both x and y are
+       the value masked.
+    """
+    fc = filled(not_equal(condition, 0), 0)
+    xv = filled(x)
+    xm = getmask(x)
+    yv = filled(y)
+    ym = getmask(y)
+    d = numeric.choose(fc, (yv, xv))
+    md = numeric.choose(fc, (ym, xm))
+    m = getmask(condition)
+    m = make_mask(mask_or(m, md), copy=False, small_mask=True)
+    return masked_array(d, mask=m)
+
+def choose (indices, t, out=None, mode='raise'):
+    "Returns 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 = numeric.choose(c, a)
+    m = numeric.choose(c, masks)
+    m = make_mask(mask_or(m, getmask(indices)), copy=0, small_mask=1)
+    return masked_array(d, mask=m)
+
+def round_(a, decimals=0, out=None):
+    """Returns reference to result. Copies a and rounds to 'decimals' places.
+
+    Keyword arguments:
+        decimals -- number of decimals to round to (default 0). May be negative.
+        out -- existing array to use for output (default copy of a).
+
+    Return:
+        Reference to out, where None specifies a copy of the original array a.
+
+    Round to the specified number of decimals. 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."""
+    result = fromnumeric.round_(filled(a), decimals, out)
+    if isinstance(a,MaskedArray):
+        result = result.view(type(a))
+        result._mask = a._mask
+    else:
+        result = result.view(MaskedArray)
+    return result
+
+def arange(start, stop=None, step=1, dtype=None):
+    """Just like range() except it returns a array whose type can be specified
+    by the keyword argument dtype.
+    """
+    return array(numeric.arange(start, stop, step, dtype),mask=nomask)
+
+def inner(a, b):
+    """inner(a,b) returns the dot product of two arrays, which has
+    shape a.shape[:-1] + b.shape[:-1] with elements computed by summing the
+    product of the elements from the last dimensions of a and b.
+    Masked elements are replace by zeros.
+    """
+    fa = filled(a, 0)
+    fb = filled(b, 0)
+    if len(fa.shape) == 0:
+        fa.shape = (1,)
+    if len(fb.shape) == 0:
+        fb.shape = (1,)
+    return masked_array(numeric.inner(fa, fb))
+innerproduct = inner
+
+def outer(a, b):
+    """outer(a,b) = {a[i]*b[j]}, has shape (len(a),len(b))"""
+    fa = filled(a, 0).ravel()
+    fb = filled(b, 0).ravel()
+    d = numeric.outer(fa, fb)
+    ma = getmask(a)
+    mb = getmask(b)
+    if ma is nomask and mb is nomask:
+        return masked_array(d)
+    ma = getmaskarray(a)
+    mb = getmaskarray(b)
+    m = make_mask(1-numeric.outer(1-ma, 1-mb), copy=0)
+    return masked_array(d, mask=m)
+outerproduct = outer
+
+def allequal (a, b, fill_value=True):
+    """
+Returns `True` if all entries of  a and b are equal, using
+fill_value as a truth value where either or both are masked.
+    """
+    m = mask_or(getmask(a), getmask(b))
+    if m is nomask:
+        x = filled(a)
+        y = filled(b)
+        d = umath.equal(x, y)
+        return d.all()
+    elif fill_value:
+        x = filled(a)
+        y = filled(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):
+    """ Returns `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 = filled(a)
+    d2 = filled(b)
+    x = filled(array(d1, copy=0, mask=m), fill_value).astype(float)
+    y = filled(array(d2, copy=0, mask=m), 1).astype(float)
+    d = umath.less_equal(umath.absolute(x-y), atol + rtol * umath.absolute(y))
+    return fromnumeric.alltrue(fromnumeric.ravel(d))
+
+#..............................................................................
+def asarray(a, dtype=None):
+    """asarray(data, dtype) = array(data, dtype, copy=0)
+Returns `a` as an masked array.
+No copy is performed if `a` is already an array.
+Subclasses are converted to base class MaskedArray.
+    """
+    return masked_array(a, dtype=dtype, copy=False, keep_mask=True)
+
+def empty(new_shape, dtype=float):
+    """empty((d1,...,dn),dtype=float,order='C')
+Returns a new array of shape (d1,...,dn) and given type with all its
+entries uninitialized. This can be faster than zeros."""
+    return numeric.empty(new_shape, dtype).view(MaskedArray)
+
+def empty_like(a):
+    """empty_like(a)
+Returns an empty (uninitialized) array of the shape and typecode of a.
+Note that this does NOT initialize the returned array.
+If you require your array to be initialized, you should use zeros_like()."""
+    return numeric.empty_like(a).view(MaskedArray)
+
+def ones(new_shape, dtype=float):
+    """ones(shape, dtype=None)
+Returns an array of the given dimensions, initialized to all ones."""
+    return numeric.ones(new_shape, dtype).view(MaskedArray)
+
+def zeros(new_shape, dtype=float):
+    """zeros(new_shape, dtype=None)
+Returns an array of the given dimensions, initialized to all zeros."""
+    return numeric.zeros(new_shape, dtype).view(MaskedArray)
+
+#####--------------------------------------------------------------------------
+#---- --- Pickling ---
+#####--------------------------------------------------------------------------
+def dump(a,F):
+    """Pickles 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):
+    """Returns 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):
+    "Loads a pickle from the current string."""
+    return cPickle.loads(strg)
+
+
+################################################################################
+
+if __name__ == '__main__':
+    from testutils import assert_equal, assert_almost_equal
+    if 1:
+        narray = numpy.array
+        pi = numpy.pi
+        x = narray([1.,1.,1.,-2., pi/2., 4., 5., -10., 10., 1., 2., 3.])
+        y = narray([5.,0.,3., 2., -1.,  -4., 0., -10., 10., 1., 0., 3.])
+        a10 = 10.
+        m1 = [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]
+        m2 = [0, 0, 1, 0, 0, 1, 1, 0, 0, 0 ,0, 1]
+        xm = masked_array(x, mask=m1)
+        ym = masked_array(y, mask=m2)
+        
+    #
+    if 1:
+        n = [0,0,1,0,1]
+        m = make_mask(n)
+        m2 = make_mask(m)
+        assert(m is m2)
+        m3 = make_mask(m, copy=1)
+        assert(m is not m3)
+
+        x1 = numpy.arange(5)
+        y1 = array(x1, mask=m)
+        #assert( y1._data is x1)
+        assert_equal(y1._data.__array_interface__, x1.__array_interface__)
+        assert( allequal(x1,y1.raw_data()))
+        #assert( y1.mask is m)
+        assert_equal(y1._mask.__array_interface__, m.__array_interface__)
+
+        y1a = array(y1)
+        #assert( y1a.raw_data() is y1.raw_data())
+        assert( y1a._data.__array_interface__ == y1._data.__array_interface__)
+        assert( y1a.mask is y1.mask)
+
+        y2 = array(x1, mask=m)
+        #assert( y2.raw_data() is x1)
+        assert (y2._data.__array_interface__ == x1.__array_interface__)
+        #assert( y2.mask is m)
+        assert (y2._mask.__array_interface__ == m.__array_interface__)
+        assert( y2[2] is masked)
+        y2[2] = 9
+        assert( y2[2] is not masked)
+        #assert( y2.mask is not m)
+        assert (y2._mask.__array_interface__ != m.__array_interface__)
+        assert( allequal(y2.mask, 0))
\ No newline at end of file

Added: trunk/scipy/sandbox/maskedarray/bench.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/bench.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/bench.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -0,0 +1,199 @@
+#! python
+
+import timeit
+#import IPython.ipapi
+#ip = IPython.ipapi.get()
+#from IPython import ipmagic
+import numpy 
+import maskedarray
+from maskedarray import filled
+from maskedarray.testutils import assert_equal
+
+
+#####---------------------------------------------------------------------------
+#---- --- Global variables ---
+#####---------------------------------------------------------------------------
+
+# Small arrays ..................................
+xs = numpy.random.uniform(-1,1,6).reshape(2,3)
+ys = numpy.random.uniform(-1,1,6).reshape(2,3)
+zs = xs + 1j * ys
+m1 = [[True, False, False], [False, False, True]]
+m2 = [[True, False, True], [False, False, True]]
+nmxs = numpy.ma.array(xs, mask=m1)
+nmys = numpy.ma.array(ys, mask=m2)
+nmzs = numpy.ma.array(zs, mask=m1)
+mmxs = maskedarray.array(xs, mask=m1)
+mmys = maskedarray.array(ys, mask=m2)
+mmzs = maskedarray.array(zs, mask=m1)
+# Big arrays ....................................
+xl = numpy.random.uniform(-1,1,100*100).reshape(100,100)
+yl = numpy.random.uniform(-1,1,100*100).reshape(100,100)
+zl = xl + 1j * yl
+maskx = xl > 0.8
+masky = yl < -0.8
+nmxl = numpy.ma.array(xl, mask=maskx)
+nmyl = numpy.ma.array(yl, mask=masky)
+nmzl = numpy.ma.array(zl, mask=maskx)
+mmxl = maskedarray.array(xl, mask=maskx, shrink=True)
+mmyl = maskedarray.array(yl, mask=masky, shrink=True)
+mmzl = maskedarray.array(zl, mask=maskx, shrink=True)
+
+#####---------------------------------------------------------------------------
+#---- --- Functions ---
+#####---------------------------------------------------------------------------
+
+def timer(s, v='', nloop=500, nrep=3):
+    units = ["s", "ms", "\xb5s", "ns"]
+    scaling = [1, 1e3, 1e6, 1e9]
+    print "%s : %-50s : " % (v,s),
+    varnames = ["%ss,nm%ss,mm%ss,%sl,nm%sl,mm%sl" % tuple(x*6) for x in 'xyz']
+    setup = 'from __main__ import numpy, maskedarray, %s' % ','.join(varnames)
+    Timer = timeit.Timer(stmt=s, setup=setup)
+    best = min(Timer.repeat(nrep, nloop)) / nloop
+    if best > 0.0:
+        order = min(-int(numpy.floor(numpy.log10(best)) // 3), 3)
+    else:
+        order = 3
+    print "%d loops, best of %d: %.*g %s per loop" % (nloop, nrep,
+                                                      3,
+                                                      best * scaling[order],
+                                                      units[order])
+#    ip.magic('timeit -n%i %s' % (nloop,s))
+
+
+
+def compare_functions_1v(func, nloop=500, test=True,
+                       xs=xs, nmxs=nmxs, mmxs=mmxs,
+                       xl=xl, nmxl=nmxl, mmxl=mmxl):
+    funcname = func.__name__
+    print "-"*50
+    print "%s on small arrays" % funcname
+    if test:
+        assert_equal(filled(eval("numpy.ma.%s(nmxs)" % funcname),0),
+                     filled(eval("maskedarray.%s(mmxs)" % funcname),0))
+    for (module, data) in zip(("numpy", "numpy.ma","maskedarray","maskedarray._nfcore"),
+                              ("xs","nmxs","mmxs","mmxs")):
+        timer("%(module)s.%(funcname)s(%(data)s)" % locals())
+    #
+    print "%s on large arrays" % funcname
+    if test:
+        assert_equal(filled(eval("numpy.ma.%s(nmxl)" % funcname),0), 
+                     filled(eval("maskedarray.%s(mmxl)" % funcname),0))
+    for (module, data) in zip(("numpy", "numpy.ma","maskedarray","maskedarray._nfcore"),
+                              ("xl","nmxl","mmxl","mmxl")):
+        timer("%(module)s.%(funcname)s(%(data)s)" % locals())
+    return
+
+def compare_methods(methodname, args, vars='x', nloop=500, test=True,
+                    xs=xs, nmxs=nmxs, mmxs=mmxs,
+                    xl=xl, nmxl=nmxl, mmxl=mmxl):
+    print "-"*50
+    print "%s on small arrays" % methodname
+    if test:
+        assert_equal(filled(eval("nm%ss.%s(%s)" % (vars,methodname,args)),0), 
+                     filled(eval("mm%ss.%s(%s)" % (vars,methodname,args)),0))
+    for (data, ver) in zip(["nm%ss" % vars, "mm%ss" % vars], ('numpy.ma   ','maskedarray')):
+        timer("%(data)s.%(methodname)s(%(args)s)" % locals(), v=ver, nloop=nloop)
+    #
+    print "%s on large arrays" % methodname
+    if test:
+        assert_equal(filled(eval("nm%sl.%s(%s)" % (vars,methodname,args)),0), 
+                     filled(eval("mm%sl.%s(%s)" % (vars,methodname,args)),0))
+    for (data, ver) in zip(["nm%sl" % vars, "mm%sl" % vars], ('numpy.ma   ','maskedarray')):
+        timer("%(data)s.%(methodname)s(%(args)s)" % locals(), v=ver, nloop=nloop)
+    return
+
+def compare_functions_2v(func, nloop=500, test=True,
+                       xs=xs, nmxs=nmxs, mmxs=mmxs,
+                       ys=ys, nmys=nmys, mmys=mmys,
+                       xl=xl, nmxl=nmxl, mmxl=mmxl,
+                       yl=yl, nmyl=nmyl, mmyl=mmyl):
+    funcname = func.__name__
+    print "-"*50
+    print "%s on small arrays" % funcname
+    if test:
+        assert_equal(filled(eval("numpy.ma.%s(nmxs,nmys)" % funcname),0), 
+                     filled(eval("maskedarray.%s(mmxs,mmys)" % funcname),0))
+    for (module, data) in zip(("numpy", "numpy.ma","maskedarray","maskedarray._nfcore"),
+                              ("xs,ys","nmxs,nmys","mmxs,mmys","mmxs,mmys")):
+        timer("%(module)s.%(funcname)s(%(data)s)" % locals())
+    #
+    print "%s on large arrays" % funcname
+    if test:
+        assert_equal(filled(eval("numpy.ma.%s(nmxl, nmyl)" % funcname),0), 
+                     filled(eval("maskedarray.%s(mmxl, mmyl)" % funcname),0))
+    for (module, data) in zip(("numpy", "numpy.ma","maskedarray","maskedarray._nfcore"),
+                              ("xl,yl","nmxl,nmyl","mmxl,mmyl","mmxl,mmyl")):
+        timer("%(module)s.%(funcname)s(%(data)s)" % locals())
+    return
+
+
+###############################################################################
+
+
+################################################################################
+if __name__ == '__main__':
+#    # Small arrays ..................................
+#    xs = numpy.random.uniform(-1,1,6).reshape(2,3)
+#    ys = numpy.random.uniform(-1,1,6).reshape(2,3)
+#    zs = xs + 1j * ys
+#    m1 = [[True, False, False], [False, False, True]]
+#    m2 = [[True, False, True], [False, False, True]]
+#    nmxs = numpy.ma.array(xs, mask=m1)
+#    nmys = numpy.ma.array(ys, mask=m2)
+#    nmzs = numpy.ma.array(zs, mask=m1)
+#    mmxs = maskedarray.array(xs, mask=m1)
+#    mmys = maskedarray.array(ys, mask=m2)
+#    mmzs = maskedarray.array(zs, mask=m1)
+#    # Big arrays ....................................
+#    xl = numpy.random.uniform(-1,1,100*100).reshape(100,100)
+#    yl = numpy.random.uniform(-1,1,100*100).reshape(100,100)
+#    zl = xl + 1j * yl
+#    maskx = xl > 0.8
+#    masky = yl < -0.8
+#    nmxl = numpy.ma.array(xl, mask=maskx)
+#    nmyl = numpy.ma.array(yl, mask=masky)
+#    nmzl = numpy.ma.array(zl, mask=maskx)
+#    mmxl = maskedarray.array(xl, mask=maskx, shrink=True)
+#    mmyl = maskedarray.array(yl, mask=masky, shrink=True)
+#    mmzl = maskedarray.array(zl, mask=maskx, shrink=True)
+#
+    compare_functions_1v(numpy.sin)
+    compare_functions_1v(numpy.log)
+    compare_functions_1v(numpy.sqrt)
+    #....................................................................
+    compare_functions_2v(numpy.multiply)
+    compare_functions_2v(numpy.divide)
+    compare_functions_2v(numpy.power)
+    #....................................................................
+    compare_methods('ravel','', nloop=1000)
+    compare_methods('conjugate','','z', nloop=1000)
+    compare_methods('transpose','', nloop=1000)
+    compare_methods('compressed','', nloop=1000)
+    compare_methods('__getitem__','0', nloop=1000)
+    compare_methods('__getitem__','(0,0)', nloop=1000)
+    compare_methods('__getitem__','[0,-1]', nloop=1000)
+    compare_methods('__setitem__','0, 17', nloop=1000, test=False)
+    compare_methods('__setitem__','(0,0), 17', nloop=1000, test=False)
+    #....................................................................
+    print "-"*50
+    print "__setitem__ on small arrays"
+    timer('nmxs.__setitem__((-1,0),numpy.ma.masked)', 'numpy.ma   ',nloop=10000)
+    timer('mmxs.__setitem__((-1,0),maskedarray.masked)', 'maskedarray',nloop=10000)
+    print "-"*50
+    print "__setitem__ on large arrays"
+    timer('nmxl.__setitem__((-1,0),numpy.ma.masked)', 'numpy.ma   ',nloop=10000)
+    timer('mmxl.__setitem__((-1,0),maskedarray.masked)', 'maskedarray',nloop=10000)
+    #....................................................................
+    print "-"*50
+    print "where on small arrays"
+    assert_equal(eval("numpy.ma.where(nmxs>2,nmxs,nmys)"), 
+                 eval("maskedarray.where(mmxs>2, mmxs,mmys)"))
+    timer('numpy.ma.where(nmxs>2,nmxs,nmys)', 'numpy.ma   ',nloop=1000)
+    timer('maskedarray.where(mmxs>2, mmxs,mmys)', 'maskedarray',nloop=1000)
+    print "-"*50
+    print "where on large arrays"
+    timer('numpy.ma.where(nmxl>2,nmxl,nmyl)', 'numpy.ma   ',nloop=100)
+    timer('maskedarray.where(mmxl>2, mmxl,mmyl)', 'maskedarray',nloop=100)
+    
\ No newline at end of file

Modified: trunk/scipy/sandbox/maskedarray/core.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/core.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/core.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -27,13 +27,13 @@
                'amax', 'amin', 'anom', 'anomalies', 'any', 'arange',
                'arccos', 'arccosh', 'arcsin', 'arcsinh', 'arctan', 'arctan2',
                'arctanh', 'argmax', 'argmin', 'argsort', 'around',
-               'array', 'asarray',
+               'array', 'asarray','asanyarray',
            'bitwise_and', 'bitwise_or', 'bitwise_xor',
            'ceil', 'choose', 'compressed', 'concatenate', 'conjugate',
                'cos', 'cosh', 'count',
            'default_fill_value', 'diagonal', 'divide', 'dump', 'dumps',
            'empty', 'empty_like', 'equal', 'exp',
-           'fabs', 'fmod', 'filled', 'floor', 'floor_divide',
+           'fabs', 'fmod', 'filled', 'floor', 'floor_divide','fix_invalid',
            'getmask', 'getmaskarray', 'greater', 'greater_equal', 'hypot',
            'ids', 'inner', 'innerproduct',
                'isMA', 'isMaskedArray', 'is_mask', 'is_masked', 'isarray',
@@ -71,6 +71,7 @@
 import numpy.core.numerictypes as ntypes
 from numpy import bool_, dtype, typecodes, amax, amin, ndarray
 from numpy import expand_dims as n_expand_dims
+from numpy import array as narray
 import warnings
 
 
@@ -80,9 +81,8 @@
 divide_tolerance = 1.e-35
 numpy.seterr(all='ignore')
 
-# TODO: There's still a problem with N.add.reduce not working...
-# TODO: ...neither does N.add.accumulate
 
+
 #####--------------------------------------------------------------------------
 #---- --- Exceptions ---
 #####--------------------------------------------------------------------------
@@ -118,7 +118,6 @@
     max_filler.update([(numpy.float128,-numpy.inf)])
     min_filler.update([(numpy.float128, numpy.inf)])
 
-
 def default_fill_value(obj):
     "Calculates the default fill value for an object `obj`."
     if hasattr(obj,'dtype'):
@@ -198,26 +197,23 @@
         return t1
     return None
 
-#................................................
+#####--------------------------------------------------------------------------
 def filled(a, value = None):
-    """Returns `a` as an array with masked data replaced by `value`.
-If `value` is `None` or the special element `masked`, `get_fill_value(a)`
-is used instead.
+    """Returns 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 contiguous numeric array, `a` itself is returned.
-
-`filled(a)` can be used to be sure that the result is numeric when passing
-an object a to other software ignorant of MA, in particular to numpy itself.
+If a is already a ndarray, a itself is returned.
     """
     if hasattr(a, 'filled'):
         return a.filled(value)
     elif isinstance(a, ndarray): # and a.flags['CONTIGUOUS']:
         return a
     elif isinstance(a, dict):
-        return numeric.array(a, 'O')
+        return narray(a, 'O')
     else:
-        return numeric.array(a)
+        return narray(a)
 
+#####--------------------------------------------------------------------------
 def get_masked_subclass(*arrays):
     """Returns the youngest subclass of MaskedArray from a list of arrays,
  or MaskedArray. In case of siblings, the first takes over."""
@@ -237,21 +233,47 @@
                 rcls = cls
     return rcls
 
-def get_data(a, copy=False, subok=True):
-    """Return the ._data part of a (if any), or a as a ndarray."""
-    if hasattr(a,'_data'):
-        if copy:
-            if subok:
-                return a._data.copy()
-            return a._data.view(ndarray).copy()
-        elif subok:
-            return a._data
-        return a._data.view(ndarray)
-    return numpy.ndarray(a, copy=copy, subok=subok)
+#####--------------------------------------------------------------------------
+def get_data(a, subok=True):
+    """Returns the _data part of a (if any), or a as a ndarray.
+    
+:Parameters:
+    a : ndarray
+        A ndarray or a subclass of.
+    subok : boolean *[True]*
+        Whether to force a 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=False, fill_value=None):
+    """Returns (a copy of) a where invalid data (nan/inf) are masked and replaced
+    by fill_value.
+    If fill_value is None, a.fill_value is used instead.
+    
+:Parameters:
+    a : ndarray
+        A ndarray or a subclass of.
+    copy : boolean *[False]*
+        Whether to use a copy of a (True) or to fix a in place (False).
+    fill_value : var *[None]*
+        Value used for fixing invalid data. If None, use a default based on the
+        datatype.
+    """
+    a = masked_array(a, copy=copy, subok=True)
+    invalid = (numpy.isnan(a._data) | numpy.isinf(a._data))
+    a._mask |= invalid
+    if fill_value is None:
+        fill_value = a.fill_value
+    a._data[invalid] = fill_value
+    return a
+    
+    
 
-
-
 #####--------------------------------------------------------------------------
 #---- --- Ufuncs ---
 #####--------------------------------------------------------------------------
@@ -269,7 +291,7 @@
         self.b = b
 
     def __call__ (self, x):
-        "Execute the call behavior."
+        "Executes the call behavior."
         return umath.logical_or(umath.greater (x, self.b),
                                 umath.less(x, self.a))
 #............................
@@ -280,11 +302,11 @@
         "domain_tan(eps) = true where abs(cos(x)) < eps)"
         self.eps = eps
     def __call__ (self, x):
-        "Execute the call behavior."
+        "Executes the call behavior."
         return umath.less(umath.absolute(umath.cos(x)), self.eps)
 #............................
 class domain_safe_divide:
-    """defines a domain for safe division."""
+    """Defines a domain for safe division."""
     def __init__ (self, tolerance=divide_tolerance):
         self.tolerance = tolerance
     def __call__ (self, a, b):
@@ -297,7 +319,7 @@
         self.critical_value = critical_value
 
     def __call__ (self, x):
-        "Execute the call behavior."
+        "Executes the call behavior."
         return umath.less_equal(x, self.critical_value)
 #............................
 class domain_greater_equal:
@@ -307,17 +329,18 @@
         self.critical_value = critical_value
 
     def __call__ (self, x):
-        "Execute the call behavior."
+        "Executes the call behavior."
         return umath.less(x, self.critical_value)
+    
 #..............................................................................
 class masked_unary_operation:
     """Defines masked version of unary operations,
 where invalid values are pre-masked.
 
 :IVariables:
-    - `f` : function.
-    - `fill` : Default filling value *[0]*.
-    - `domain` : Default domain *[None]*.
+    f : function.
+    fill : Default filling value *[0]*.
+    domain : Default domain *[None]*.
     """
     def __init__ (self, mufunc, fill=0, domain=None):
         """ masked_unary_operation(aufunc, fill=0, domain=None)
@@ -334,16 +357,19 @@
         ufunc_fills[mufunc] = fill
     #
     def __call__ (self, a, *args, **kwargs):
-        "Execute the call behavior."
-# numeric tries to return scalars rather than arrays when given scalars.
+        "Executes the call behavior."
         m = getmask(a)
-        d1 = filled(a, self.fill)
+        d1 = get_data(a)
         if self.domain is not None:
-            m = mask_or(m, numeric.asarray(self.domain(d1)))
+            dm = narray(self.domain(d1), copy=False)
+            m = mask_or(m, narray(self.domain(d1)))
+            # The following two lines control the domain filling methods.
+            d1 = d1.copy()
+            numpy.putmask(d1, dm, self.fill)
         # Take care of the masked singletong first ...
-        if m.ndim == 0 and m:
+        if not m.ndim and m:
             return masked
-        # Get the result....
+        # Get the result..............................
         if isinstance(a, MaskedArray):
             result = self.f(d1, *args, **kwargs).view(type(a))
         else:
@@ -355,16 +381,17 @@
     #
     def __str__ (self):
         return "Masked version of %s. [Invalid values are masked]" % str(self.f)
+    
 #..............................................................................
 class masked_binary_operation:
     """Defines masked version of binary operations,
 where invalid values are pre-masked.
 
 :IVariables:
-    - `f` : function.
-    - `fillx` : Default filling value for first array*[0]*.
-    - `filly` : Default filling value for second array*[0]*.
-    - `domain` : Default domain *[None]*.
+    f : function.
+    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.
@@ -381,15 +408,14 @@
     def __call__ (self, a, b, *args, **kwargs):
         "Execute the call behavior."
         m = mask_or(getmask(a), getmask(b))
-        if (not m.ndim) and m:
+        (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
+        elif m:
             return masked
-        d1 = filled(a, self.fillx)
-        d2 = filled(b, self.filly)
-# CHECK : Do we really need to fill the arguments ? Pro'ly not        
-#        result = self.f(a, b, *args, **kwargs).view(get_masked_subclass(a,b))
-        result = self.f(d1, d2, *args, **kwargs).view(get_masked_subclass(a,b))
-        if result.ndim > 0:
-            result._mask = m
         return result
     #
     def reduce (self, target, axis=0, dtype=None):
@@ -409,9 +435,7 @@
             return self.f.reduce(t, axis).view(tclass)
         t = t.view(tclass)
         t._mask = m
-        # XXX: "or t.dtype" below is a workaround for what appears
-        # XXX: to be a bug in reduce.
-        tr = self.f.reduce(filled(t, self.filly), axis, dtype=dtype or t.dtype)
+        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:
@@ -434,7 +458,8 @@
         if (not m.ndim) and m:
             return masked
         rcls = get_masked_subclass(a,b)
-        d = self.f.outer(filled(a, self.fillx), filled(b, self.filly)).view(rcls)
+#        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
@@ -450,6 +475,7 @@
 
     def __str__ (self):
         return "Masked version of " + str(self.f)
+    
 #..............................................................................
 class domained_binary_operation:
     """Defines binary operations that have a domain, like divide.
@@ -458,10 +484,10 @@
 They have no reduce, outer or accumulate.
 
 :IVariables:
-    - `f` : function.
-    - `fillx` : Default filling value for first array*[0]*.
-    - `filly` : Default filling value for second array*[0]*.
-    - `domain` : Default domain *[None]*.
+    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.
@@ -480,13 +506,14 @@
         "Execute the call behavior."
         ma = getmask(a)
         mb = getmask(b)
-        d1 = filled(a, self.fillx)
-        d2 = filled(b, self.filly)
-        t = numeric.asarray(self.domain(d1, d2))
-
-        if fromnumeric.sometrue(t, None):
-            d2 = numeric.where(t, self.filly, d2)
+        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 two lines control the domain filling
+            d2 = d2.copy()
+            numpy.putmask(d2, t, self.filly)
         m = mask_or(ma, mb)
         if (not m.ndim) and m:
             return masked       
@@ -517,7 +544,7 @@
 ceil = masked_unary_operation(umath.ceil)
 around = masked_unary_operation(fromnumeric.round_)
 logical_not = masked_unary_operation(umath.logical_not)
-# Domained unary ufuncs
+# Domained unary ufuncs .......................................................
 sqrt = masked_unary_operation(umath.sqrt, 0.0, domain_greater_equal(0.0))
 log = masked_unary_operation(umath.log, 1.0, domain_greater(0.0))
 log10 = masked_unary_operation(umath.log10, 1.0, domain_greater(0.0))
@@ -529,7 +556,7 @@
 arccosh = masked_unary_operation(umath.arccosh, 1.0, domain_greater_equal(1.0))
 arctanh = masked_unary_operation(umath.arctanh, 0.0,
                                  domain_check_interval(-1.0+1e-15, 1.0-1e-15))
-# Binary ufuncs
+# Binary ufuncs ...............................................................
 add = masked_binary_operation(umath.add)
 subtract = masked_binary_operation(umath.subtract)
 multiply = masked_binary_operation(umath.multiply, 1, 1)
@@ -555,7 +582,7 @@
 bitwise_or = masked_binary_operation(umath.bitwise_or)
 bitwise_xor = masked_binary_operation(umath.bitwise_xor)
 hypot = masked_binary_operation(umath.hypot)
-# Domained binary ufuncs
+# Domained binary ufuncs ......................................................
 divide = domained_binary_operation(umath.divide, domain_safe_divide(), 0, 1)
 true_divide = domained_binary_operation(umath.true_divide,
                                         domain_safe_divide(), 0, 1)
@@ -569,27 +596,22 @@
 #####--------------------------------------------------------------------------
 #---- --- Mask creation functions ---
 #####--------------------------------------------------------------------------
-def getmask(a):
-    """Returns the mask of `a`, if any, or `nomask`.
-Returns `nomask` if `a` is not a masked array.
-To get an array for sure use getmaskarray."""
-    if hasattr(a, "_mask"):
-        return a._mask
-    else:
-        return nomask
+def get_mask(a):
+    """Returns 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):
-    """Returns the mask of `a`, if any.
-Otherwise, returns an array of `False`, with the same shape as `a`.
+    """Returns the mask of a, if any, or an array of the shape of a, full of False.
     """
     m = getmask(a)
     if m is nomask:
-        return make_mask_none(fromnumeric.shape(a))
-    else:
-        return m
+        m = make_mask_none(fromnumeric.shape(a))
+    return m
 
 def is_mask(m):
-    """Returns `True` if `m` is a legal mask.
+    """Returns True if m is a legal mask.
 Does not check contents, only type.
     """
     try:
@@ -597,78 +619,93 @@
     except AttributeError:
         return False
 #
-def make_mask(m, copy=False, small_mask=True, flag=None):
-    """make_mask(m, copy=0, small_mask=0)
-Returns `m` as a mask, creating a copy if necessary or requested.
-The function can accept any sequence of integers or `nomask`.
+def make_mask(m, copy=False, shrink=True, flag=None):
+    """make_mask(m, copy=0, shrink=0)
+Returns 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.
-If `small_mask=True`, returns `nomask` if `m` contains no true elements.
 
 :Parameters:
-    - `m` (ndarray) : Mask.
-    - `copy` (boolean, *[False]*) : Returns a copy of `m` if true.
-    - `small_mask` (boolean, *[False]*): Flattens mask to `nomask` if `m` is all false.
+    m : ndarray
+        Potential mask.
+    copy : boolean *[False]*
+        Whether to return a copy of m.
+    shrink : boolean *[True]*
+        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 'small_mask'!",
+        warnings.warn("The flag 'flag' is now called 'shrink'!",
                       DeprecationWarning)
-        small_mask = flag
+        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 = numeric.array(m, dtype=MaskType, copy=copy)
+                result = narray(m, dtype=MaskType, copy=copy)
             else:
                 result = m
         else:
-            result = numeric.array(m, dtype=MaskType)
+            result = narray(m, dtype=MaskType)
     else:
-        result = numeric.array(filled(m, True), dtype=MaskType)
+        result = narray(filled(m, True), dtype=MaskType)
     # Bas les masques !
-    if small_mask and not result.any():
+    if shrink and not result.any():
         return nomask
     else:
         return result
 
 def make_mask_none(s):
-    "Returns a mask of shape `s`, filled with `False`."
+    """Returns 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, small_mask=True):
-    """Returns the combination of two masks `m1` and `m2`.
-The masks are combined with the `logical_or` operator, treating `nomask` as false.
+def mask_or (m1, m2, copy=False, shrink=True):
+    """Returns 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:
-    - `m` (ndarray) : Mask.
-    - `copy` (boolean, *[False]*) : Returns a copy of `m` if true.
-    - `small_mask` (boolean, *[False]*): Flattens mask to `nomask` if `m` is all false.
+    m1 : ndarray
+        First mask.
+    m2 : ndarray
+        Second mask
+    copy : boolean *[False]*
+        Whether to return a copy.
+    shrink : boolean *[True]*
+        Whether to shrink m to nomask if all its values are False.
      """
     if m1 is nomask:
-        return make_mask(m2, copy=copy, small_mask=small_mask)
+        return make_mask(m2, copy=copy, shrink=shrink)
     if m2 is nomask:
-        return make_mask(m1, copy=copy, small_mask=small_mask)
+        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, small_mask=small_mask)
+    return make_mask(umath.logical_or(m1, m2), copy=copy, shrink=shrink)
 
 #####--------------------------------------------------------------------------
 #--- --- Masking functions ---
 #####--------------------------------------------------------------------------
 def masked_where(condition, a, copy=True):
-    """Returns `x` as an array masked where `condition` is true.
-Masked values of `x` or `condition` are kept.
+    """Returns a as an array masked where condition is true.
+Masked values of a or condition are kept.
 
 :Parameters:
-    - `condition` (ndarray) : Masking condition.
-    - `x` (ndarray) : Array to mask.
-    - `copy` (boolean, *[False]*) : Returns a copy of `m` if true.
+    condition : ndarray
+        Masking condition.
+    a : ndarray
+        Array to mask.
+    copy : boolean *[True]*
+        Whether to return a copy of a.
     """
     cond = filled(condition,1)
-    a = numeric.array(a, copy=copy, subok=True)
+    a = narray(a, copy=copy, subok=True)
     if hasattr(a, '_mask'):
         cond = mask_or(cond, a._mask)
         cls = type(a)
@@ -678,29 +715,29 @@
     result._mask = cond
     return result
 
-def masked_greater(x, value, copy=1):
-    "Shortcut to `masked_where`, with ``condition = (x > value)``."
+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=1):
-    "Shortcut to `masked_where`, with ``condition = (x >= value)``."
+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)``."
+    "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)``."
+    "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)``."
+    "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)``.
+    """Shortcut to masked_where, with condition = (x == value).
 For floating point, consider `masked_values(x, value)` instead.
     """
     return masked_where((x == value), x, copy=copy)
@@ -710,9 +747,9 @@
 #    return array(d, mask=m, 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.
+    """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.
     """
     if v2 < v1:
         (v1, v2) = (v2, v1)
@@ -721,9 +758,9 @@
     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.
+    """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.
     """
     if v2 < v1:
         (v1, v2) = (v2, v1)
@@ -733,46 +770,64 @@
 
 #
 def masked_object(x, value, copy=True):
-    """Masks 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.
+    """Masks 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.
 The mask is set to `nomask` if posible.
-
-:parameter copy (Boolean, *[True]*):  Returns a copy of `x` if true. """
+    """
     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, small_mask=True))
+    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):
-    """Masks the array `x` where the data are approximately equal to `value`
-(that is, ``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.
+    """Masks the array x where the data are approximately equal to value
+(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:
-    - `rtol` (Float, *[1e-5]*): Tolerance parameter.
-    - `atol` (Float, *[1e-8]*): Tolerance parameter.
-    - `copy` (boolean, *[False]*) : Returns a copy of `x` if True.
+    x : ndarray
+        Array to fill.
+    value : float
+        Masking value.
+    rtol : float *[1e-5]*
+        Tolerance parameter.
+    atol : float, *[1e-8]*
+        Tolerance parameter.
+    copy : boolean *[True]*
+        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))
-        try:
-            mask = x._mask
-        except AttributeError:
-            mask = nomask
+        mask = getattr(x, '_mask', nomask)
     else:
         condition = umath.equal(xnew, value)
         mask = nomask
-    mask = mask_or(mask, make_mask(condition, small_mask=True))
+    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):
+    """Masks the array for invalid values (nans or infs). 
+    Any preexisting mask is conserved."""
+    a = narray(a, copy=copy, subok=True)
+    condition = (numpy.isnan(a) | numpy.isinf(a))
+    if hasattr(a, '_mask'):
+        condition = mask_or(condition, a._mask)
+        cls = type(a)
+    else:
+        cls = MaskedArray
+    result = a.view(cls)
+    result._mask = cond
+    return result
+
+
 #####--------------------------------------------------------------------------
 #---- --- Printing options ---
 #####--------------------------------------------------------------------------
@@ -795,9 +850,9 @@
         "Is the use of the display value enabled?"
         return self._enabled
 
-    def enable(self, small_mask=1):
-        "Set the enabling small_mask to `small_mask`."
-        self._enabled = small_mask
+    def enable(self, shrink=1):
+        "Set the enabling shrink to `shrink`."
+        self._enabled = shrink
 
     def __str__ (self):
         return str(self._display)
@@ -810,90 +865,23 @@
 #####--------------------------------------------------------------------------
 #---- --- MaskedArray class ---
 #####--------------------------------------------------------------------------
-##def _getoptions(a_out, a_in):
-##    "Copies standards options of a_in to a_out."
-##    for att in [']
-#class _mathmethod(object):
-#    """Defines a wrapper for arithmetic methods.
-#Instead of directly calling a ufunc, the corresponding method of  the `array._data`
-#object is called instead.
-#    """
-#    def __init__ (self, methodname, fill_self=0, fill_other=0, domain=None):
-#        """
-#:Parameters:
-#    - `methodname` (String) : Method name.
-#    - `fill_self` (Float *[0]*) : Fill value for the instance.
-#    - `fill_other` (Float *[0]*) : Fill value for the target.
-#    - `domain` (Domain object *[None]*) : Domain of non-validity.
-#        """
-#        self.methodname = methodname
-#        self.fill_self = fill_self
-#        self.fill_other = fill_other
-#        self.domain = domain
-#        self.obj = None
-#        self.__doc__ = self.getdoc()
-#    #
-#    def getdoc(self):
-#        "Returns the doc of the function (from the doc of the method)."
-#        try:
-#            return getattr(MaskedArray, self.methodname).__doc__
-#        except:
-#            return getattr(ndarray, self.methodname).__doc__
-#    #
-#    def __get__(self, obj, objtype=None):
-#        self.obj = obj
-#        return self
-#    #
-#    def __call__ (self, other, *args):
-#        "Execute the call behavior."
-#        instance = self.obj
-#        m_self = instance._mask
-#        m_other = getmask(other)
-#        base = instance.filled(self.fill_self)
-#        target = filled(other, self.fill_other)
-#        if self.domain is not None:
-#            # We need to force the domain to a ndarray only.
-#            if self.fill_other > self.fill_self:
-#                domain = self.domain(base, target)
-#            else:
-#                domain = self.domain(target, base)
-#            if domain.any():
-#                #If `other` is a subclass of ndarray, `filled` must have the
-#                # same subclass, else we'll lose some info.
-#                #The easiest then is to fill `target` instead of creating
-#                # a pure ndarray.
-#                #Oh, and we better make a copy!
-#                if isinstance(other, ndarray):
-#                    # We don't want to modify other: let's copy target, then
-#                    target = target.copy()
-#                    target[fromnumeric.asarray(domain)] = self.fill_other
-#                else:
-#                    target = numeric.where(fromnumeric.asarray(domain),
-#                                           self.fill_other, target)
-#                m_other = mask_or(m_other, domain)
-#        m = mask_or(m_self, m_other)
-#        method = getattr(base, self.methodname)
-#        result = method(target, *args).view(type(instance))
-#        try:
-#            result._mask = m
-#        except AttributeError:
-#            if m:
-#                result = masked
-#        return result
+
 #...............................................................................
 class _arraymethod(object):
     """Defines 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`.
+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 calld on the initial mask.
-If `onmask` is False, the new mask is just a reference to the initial mask.
+If onmask is True, the new mask is the output of the method called on the initial mask.
+If onmask is False, the new mask is just a reference to the initial mask.
 
-:Parameters:
-    `funcname` : String
+:IVariables:
+    _name : String
         Name of the function to apply on data.
-    `onmask` : Boolean *[True]*
+    _onmask : Boolean *[True]*
         Whether the mask must be processed also (True) or left alone (False).
+    obj : Object
+        The object calling the arraymethod
     """
     def __init__(self, funcname, onmask=True):
         self._name = funcname
@@ -905,16 +893,8 @@
         "Returns the doc of the function (from the doc of the method)."
         methdoc = getattr(ndarray, self._name, None)
         methdoc = getattr(numpy, self._name, methdoc)
-#        methdoc = getattr(MaskedArray, self._name, methdoc)
         if methdoc is not None:
             return methdoc.__doc__
-#        try:
-#            return getattr(MaskedArray, self._name).__doc__
-#        except:
-#            try:
-#                return getattr(numpy, self._name).__doc__
-#            except:
-#                return getattr(ndarray, self._name).__doc
     #
     def __get__(self, obj, objtype=None):
         self.obj = obj
@@ -926,10 +906,11 @@
         mask = self.obj._mask
         cls = type(self.obj)
         result = getattr(data, methodname)(*args, **params).view(cls)
-        result._smallmask = self.obj._smallmask
+        result._update_from(self.obj)
+        #result._shrinkmask = self.obj._shrinkmask
         if result.ndim:
             if not self._onmask:
-                result._mask = mask
+                result.__setmask__(mask)
             elif mask is not nomask:
                 result.__setmask__(getattr(mask, methodname)(*args, **params))
         else:
@@ -969,54 +950,60 @@
 Masked values of True exclude the corresponding element from any computation.
 
 Construction:
-    x = array(data, dtype=None, copy=True, order=False,
-              mask = nomask, fill_value=None, small_mask=True)
+    x = MaskedArray(data, mask=nomask, dtype=None, copy=True, fill_value=None,
+              mask = nomask, fill_value=None, shrink=True)
 
-If copy=False, every effort is made not to copy the data:
-If `data` is a MaskedArray, and argument mask=nomask, then the candidate data
-is `data._data` and the mask used is `data._mask`.
-If `data` is a numeric array, it is used as the candidate raw data.
-If `dtype` is not None and is different from data.dtype.char then a data copy is required.
-Otherwise, the candidate is used.
-
-If a data copy is required, the raw (unmasked) data stored is the result of:
-numeric.array(data, dtype=dtype.char, copy=copy)
-
-If `mask` is `nomask` there are no masked values.
-Otherwise mask must be convertible to an array of booleans with the same shape as x.
-If `small_mask` is True, a mask consisting of zeros (False) only is compressed to `nomask`.
-Otherwise, the mask is not compressed.
-
-fill_value is used to fill in masked values when necessary, such as when
-printing and in method/function filled().
-The fill_value is not used for computation within this module.
+:Parameters:
+    data : var
+        Input data.
+    mask : sequence *[nomask]*
+        Mask. 
+        Must be convertible to an array of booleans with the same shape as data:
+        True indicates a masked (eg., invalid) data.
+    dtype : dtype *[None]*
+        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 : boolean *[False]*
+        Whether to copy the input data (True), or to use a reference instead.
+    fill_value : var *[None]*
+        Value used to fill in the masked values when necessary. If None, a default 
+        based on the datatype is used.
+    keep_mask : boolean *[True]*
+        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 : boolean *[False]*
+        Whether to use a hard mask or not. With a hard mask, masked values cannot
+        be unmasked.
+    subok : boolean *[True]*
+        Whether to return a subclass of MaskedArray (if possible) or a plain
+        MaskedArray.
     """
-    __array_priority__ = 10.1
+    
+    __array_priority__ = 15
     _defaultmask = nomask
     _defaulthardmask = False
     _baseclass =  numeric.ndarray
+    
     def __new__(cls, data=None, mask=nomask, dtype=None, copy=False, fill_value=None,
-                keep_mask=True, small_mask=True, hard_mask=False, flag=None,
+                keep_mask=True, hard_mask=False, flag=None,
                 subok=True, **options):
         """array(data, dtype=None, copy=True, mask=nomask, fill_value=None)
 
 If `data` is already a ndarray, its dtype becomes the default value of dtype.
         """
         if flag is not None:
-            warnings.warn("The flag 'flag' is now called 'small_mask'!",
+            warnings.warn("The flag 'flag' is now called 'shrink'!",
                           DeprecationWarning)
-            small_mask = flag
+            shrink = flag
         # Process data............
-        _data = numeric.array(data, dtype=dtype, copy=copy, subok=subok)
+        _data = narray(data, dtype=dtype, copy=copy, subok=True)
         _baseclass = getattr(data, '_baseclass', type(_data))
         _basedict = getattr(data, '_basedict', getattr(data, '__dict__', None))
-        if not isinstance(data, MaskedArray): 
+        if not isinstance(data, MaskedArray) or not subok: 
             _data = _data.view(cls)
-        elif not subok:
-            _data = data.view(cls)
         else:
             _data = _data.view(type(data))
-        # Backwards compat .......
+        # Backwards compatibility w/ numpy.core.ma .......
         if hasattr(data,'_mask') and not isinstance(data, ndarray):
             _data._mask = data._mask
             _sharedmask = True
@@ -1026,8 +1013,11 @@
                 _data._mask = nomask
             if copy:
                 _data._mask = _data._mask.copy()
+                _data._sharedmask = False
+            else:
+                _data._sharedmask = True
         else:
-            mask = numeric.array(mask, dtype=MaskType, copy=copy)
+            mask = narray(mask, dtype=MaskType, copy=copy)
             if mask.shape != _data.shape:
                 (nd, nm) = (_data.size, mask.size) 
                 if nm == 1:
@@ -1038,18 +1028,18 @@
                     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 = True
+                _data._sharedmask = not copy
             else:
-                # Make a copy of the mask to avoid propagation
-                _data._sharedmask = False
                 if not keep_mask:
                     _data._mask = mask
+                    _data._sharedmask = not copy
                 else:
-                    _data._mask = umath.logical_or(mask, _data._mask) 
+                    _data._mask = umath.logical_or(mask, _data._mask)
+                    _data._sharedmask = False                   
                     
-                    
         # Update fill_value.......
         if fill_value is None:
             _data._fill_value = getattr(data, '_fill_value', 
@@ -1058,66 +1048,80 @@
             _data._fill_value = fill_value
         # Process extra options ..
         _data._hardmask = hard_mask
-        _data._smallmask = small_mask
         _data._baseclass = _baseclass
         _data._basedict = _basedict
         return _data
+    #
+    def _update_from(self, obj):
+        """Copies some attributes of obj to self.
+        """
+        self._hardmask = getattr(obj, '_hardmask', self._defaulthardmask)
+        self._sharedmask = getattr(obj, '_sharedmask', False)
+        if obj is not None:
+            self._baseclass = getattr(obj, '_baseclass', type(obj))
+        else:
+            self._baseclass = ndarray
+        self._fill_value = getattr(obj, '_fill_value', None)
+        return
     #........................
     def __array_finalize__(self,obj):
         """Finalizes the masked array.
         """
-        # Finalize mask ...............
+        # Get main attributes .........
         self._mask = getattr(obj, '_mask', nomask)
-        if self._mask is not nomask:
-            self._mask.shape = self.shape
-        # Get the remaining options ...
-        self._hardmask = getattr(obj, '_hardmask', self._defaulthardmask)
-        self._smallmask = getattr(obj, '_smallmask', True)
-        self._sharedmask = True
-        self._baseclass = getattr(obj, '_baseclass', type(obj))
-        self._fill_value = getattr(obj, '_fill_value', None)
+        self._update_from(obj)
         # Update special attributes ...
         self._basedict = getattr(obj, '_basedict', getattr(obj, '__dict__', None))
         if self._basedict is not None:
             self.__dict__.update(self._basedict)
+        # 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.
         """
-        #TODO : Should we check for type result 
         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, [getmask(arg) for arg in args])
-            # Get domain mask
+            # 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
-            if not m.ndim and m:
-                if m:
-                    if result.shape == ():
-                        return masked
-                    result._mask = numeric.ones(result.shape, bool_)
+            # Make sure the mask has the proper size
+            if result.shape == () and m:
+                return masked
             else:
                 result._mask = m
+                result._sharedmask = False
         #....
-#        result._mask = m
-        result._fill_value = self._fill_value
-        result._hardmask = self._hardmask
-        result._smallmask = self._smallmask
-        result._baseclass = self._baseclass
         return result
     #.............................................
     def __getitem__(self, indx):
@@ -1128,20 +1132,26 @@
 #        if getmask(indx) is not nomask:
 #            msg = "Masked arrays must be filled before they can be used as indices!"
 #            raise IndexError, msg
-        # super() can't work here if the underlying data is a matrix...
-        dout = (self._data).__getitem__(indx)
+        dout = ndarray.__getitem__(self.view(ndarray), indx)
         m = self._mask
-        if hasattr(dout, 'shape') and len(dout.shape) > 0:
-            # Not a scalar: make sure that dout is a MA
+        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))
-            dout._smallmask = self._smallmask
-            dout._hardmask = self._hardmask
-            dout._fill_value = self._fill_value
+            # Inherit attributes from self
+            dout._update_from(self)
+            # Update the mask if needed
             if m is not nomask:
-                # use _set_mask to take care of the shape
-                dout.__setmask__(m[indx])
-        elif m is not nomask and m[indx]:
-            return masked
+                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):
@@ -1157,26 +1167,22 @@
         if value is masked:
             m = self._mask
             if m is nomask:
-                m = make_mask_none(self.shape)
-#            else:
-#                m = m.copy()
+                m = numpy.zeros(self.shape, dtype=MaskType)
             m[indx] = True
-            self.__setmask__(m)
+            self._mask = m
+            self._sharedmask = False
             return
         #....
-        dval = numeric.asarray(value).astype(self.dtype)
+        dval = getdata(value).astype(self.dtype)
         valmask = getmask(value)
         if self._mask is nomask:
             if valmask is not nomask:
-                self._mask = make_mask_none(self.shape)
+                self._mask = numpy.zeros(self.shape, dtype=MaskType)
                 self._mask[indx] = valmask
         elif not self._hardmask:
-            _mask = self._mask.copy()
-            if valmask is nomask:
-                _mask[indx] = False
-            else:
-                _mask[indx] = valmask
-            self._set_mask(_mask)
+            # Unshare the mask if necessary to avoid propagation
+            self.unshare_mask()
+            self._mask[indx] = valmask
         elif hasattr(indx, 'dtype') and (indx.dtype==bool_):
             indx = indx * umath.logical_not(self._mask)
         else:
@@ -1189,7 +1195,6 @@
             dval = dindx
             self._mask[indx] = mindx
         # Set data ..........
-        #dval = filled(value).astype(self.dtype)
         ndarray.__setitem__(self._data,indx,dval)
     #............................................
     def __getslice__(self, i, j):
@@ -1205,29 +1210,33 @@
         self.__setitem__(slice(i,j), value)
     #............................................
     def __setmask__(self, mask, copy=False):
-        newmask = make_mask(mask, copy=copy, small_mask=self._smallmask)
-#        self.unshare_mask()
+        """Sets the mask."""
+        if mask is not nomask:
+            mask = narray(mask, copy=copy, dtype=MaskType)
+#            if self._shrinkmask and not mask.any():
+#                mask = nomask
         if self._mask is nomask:
-            self._mask = newmask
+            self._mask = mask
         elif self._hardmask:
-            if newmask is not nomask:
-                self._mask.__ior__(newmask)
+            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 newmask is nomask:
+            if mask is nomask:
                 self._mask = nomask
             else:
-                self._mask.flat = newmask
+                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):
         """Returns the current mask."""
         return self._mask
-
+#        return self._mask.reshape(self.shape)
     mask = property(fget=_get_mask, fset=__setmask__, doc="Mask")
     #............................................
     def harden_mask(self):
@@ -1243,19 +1252,24 @@
         if self._sharedmask:
             self._mask = self._mask.copy()
             self._sharedmask = False
-        
     #............................................
     def _get_data(self):
         "Returns 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):
+        """Returns the `_data` part of the MaskedArray.
+You should really use `data` instead..."""
+        return self._data
     #............................................
     def _get_flat(self):
-        """Calculates the flat value.
-        """
+        "Returns a flat iterator."
         return flatiter(self)
     #
     def _set_flat (self, value):
+        "Sets a flattened version of self to value."
         "x.flat = value"
         y = self.ravel()
         y[:] = value
@@ -1269,8 +1283,8 @@
         return self._fill_value
 
     def set_fill_value(self, value=None):
-        """Sets the filling value to `value`.
-If None, uses the default, based on the data type."""
+        """Sets the filling value to value.
+If None, uses a default based on the data type."""
         if value is None:
             value = default_fill_value(self)
         self._fill_value = value
@@ -1279,11 +1293,18 @@
                           doc="Filling value")
 
     def filled(self, fill_value=None):
-        """Returns an array of the same class as `_data`,
- with masked values filled with `fill_value`.
-Subclassing is preserved.
-
-If `fill_value` is None, uses self.fill_value.
+        """Returns a copy of self._data, where masked values are filled with 
+    fill_value. If fill_value is None, self.fill_value is used instead.
+    Subclassing is preserved. 
+    Note : 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():
@@ -1298,9 +1319,8 @@
             result = self._data.copy()
             try:
                 numpy.putmask(result, m, fill_value)
-                #result[m] = fill_value
             except (TypeError, AttributeError):
-                fill_value = numeric.array(fill_value, dtype=object)
+                fill_value = narray(fill_value, dtype=object)
                 d = result.astype(object)
                 result = fromnumeric.choose(m, (d, fill_value))
             except IndexError:
@@ -1308,20 +1328,21 @@
                 if self._data.shape:
                     raise
                 elif m:
-                    result = numeric.array(fill_value, dtype=self.dtype)
+                    result = narray(fill_value, dtype=self.dtype)
                 else:
                     result = self._data
         return result
 
     def compressed(self):
-        "A 1-D array of all the non-masked data."
-        d = self.ravel()
-        if self._mask is nomask:
-            return d
-        elif not self._smallmask and not self._mask.any():
-            return d
-        else:
-            return d[numeric.logical_not(d._mask)]
+        """Returns a 1-D array of all the non-masked data."""
+        data = ndarray.ravel(self._data).view(type(self))
+        data._update_from(self)
+        if self._mask is not nomask:
+            data = data[numpy.logical_not(ndarray.ravel(self._mask))]
+#        if not self._shrinkmask:
+#            data._mask = numpy.zeros(data.shape, dtype=MaskType)
+        return data
+
     #............................................
     def __str__(self):
         """x.__str__() <==> str(x)
@@ -1384,9 +1405,34 @@
             'fill': str(self.fill_value),
             }
     #............................................
+    def __add__(self, other):
+        "Adds other to self, and returns a new masked array."
+        return add(self, other)
+    #
+    def __sub__(self, other):
+        "Subtracts other to self, and returns a new masked array."
+        return subtract(self, other)
+    #
+    def __mul__(self, other):
+        "Multiplies other to self, and returns a new masked array."
+        return multiply(self, other)
+    #
+    def __div__(self, other):
+        "Divides other to self, and returns a new masked array."
+        return divide(self, other)
+    #
+    def __truediv__(self, other):
+        "Divides other to self, and returns a new masked array."
+        return true_divide(self, other)
+    #
+    def __floordiv__(self, other):
+        "Divides other to self, and returns a new masked array."
+        return floor_divide(self, other)
+    
+    #............................................
     def __iadd__(self, other):
         "Adds other to self in place."
-        ndarray.__iadd__(self._data,other)
+        ndarray.__iadd__(self._data, getdata(other))
         m = getmask(other)
         if self._mask is nomask:
             self._mask = m
@@ -1396,7 +1442,7 @@
     #....
     def __isub__(self, other):
         "Subtracts other from self in place."
-        ndarray.__isub__(self._data,other)
+        ndarray.__isub__(self._data, getdata(other))
         m = getmask(other)
         if self._mask is nomask:
             self._mask = m
@@ -1406,7 +1452,7 @@
     #....
     def __imul__(self, other):
         "Multiplies self by other in place."
-        ndarray.__imul__(self._data,other)
+        ndarray.__imul__(self._data, getdata(other))
         m = getmask(other)
         if self._mask is nomask:
             self._mask = m
@@ -1416,10 +1462,15 @@
     #....
     def __idiv__(self, other):
         "Divides self by other in place."
-        dom_mask = domain_safe_divide().__call__(self, filled(other,1))
+        other_data = getdata(other)
+        dom_mask = domain_safe_divide().__call__(self._data, other_data)
         other_mask = getmask(other)
         new_mask = mask_or(other_mask, dom_mask)
-        ndarray.__idiv__(self._data, other)
+        # 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
     #............................................
@@ -1440,7 +1491,7 @@
     def count(self, axis=None):
         """Counts the non-masked elements of the array along a given axis,
 and returns a masked array where the mask is True where all data are masked.
-If `axis` is None, counts all the non-masked elements, and returns either a
+If axis is None, counts all the non-masked elements, and returns either a
 scalar or the masked singleton."""
         m = self._mask
         s = self.shape
@@ -1457,17 +1508,30 @@
                 t = list(s)
                 del t[axis]
                 return numeric.ones(t) * n
-        n1 = fromnumeric.size(m, axis)
+        n1 = numpy.size(m, axis)
         n2 = m.astype(int_).sum(axis)
         if axis is None:
             return (n1-n2)
         else:
             return masked_array(n1 - n2)
     #............................................
+    flatten = _arraymethod('flatten')
+#    ravel = _arraymethod('ravel')
+    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):
         """Reshapes the array to shape s.
-Returns a new masked array.
-If you want to modify the shape in place, please use `a.shape = s`"""
+    Returns a new masked array.
+    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:
@@ -1475,12 +1539,10 @@
             result._mask.shape = result.shape
         return result
     #
-    repeat = _arraymethod('repeat')
-    #
     def resize(self, newshape, refcheck=True, order=False):
         """Attempts to modify size and shape of self inplace.
-        The array must own its own memory and not be referenced by other arrays.
-        Returns None.
+    The array must own its own memory and not be referenced by other arrays.
+    Returns None.
         """
         try:
             self._data.resize(newshape, refcheck, order)
@@ -1492,22 +1554,19 @@
                              "Use the resize function.")
         return None
     #
-    flatten = _arraymethod('flatten')
-    #
     def put(self, indices, values, mode='raise'):
         """Sets storage-indexed locations to corresponding values.
-a.put(values, indices, mode) sets a.flat[n] = values[n] for each n in indices.
-`values` can be scalar or an array shorter than indices, and it will be repeated,
-if necessary.
-If `values` has some masked values, the initial mask is updated in consequence,
-else the corresponding values are unmasked.
+    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 = numeric.asarray(indices)
-            values = numeric.asanyarray(values)
+            indices = narray(indices, copy=False)
+            values = narray(values, copy=False, subok=True)
             values.resize(indices.shape)
             indices = indices[~mask]
             values = values[~mask]
@@ -1522,17 +1581,17 @@
                 m.put(indices, False, mode=mode)
             else:
                 m.put(indices, values._mask, mode=mode)
-            m = make_mask(m, copy=False, small_mask=True)
+            m = make_mask(m, copy=False, shrink=True)
         self._mask = m
     #............................................
     def ids (self):
-        """Return the address of the data and mask areas."""
+        """Returns the addresses of the data and mask areas."""
         return (self.ctypes.data, self._mask.ctypes.data)    
     #............................................
     def all(self, axis=None, out=None):
         """a.all(axis) returns True if all entries along the axis are True.
     Returns False otherwise. If axis is None, uses the flatten array.
-    Masked data are considered as True during computation.
+    Masked values are considered as True during computation.
     Outputs a masked array, where the mask is True if all data are masked along the axis.
     Note: the out argument is not really operational...
         """
@@ -1554,19 +1613,19 @@
         return d
     
     def nonzero(self):
-        """a.nonzero() returns a tuple of arrays
+        """a.nonzero() returns the indices of the elements of a that are not 
+    zero nor masked, as a tuple of arrays.
 
-    Returns a tuple of arrays, one for each dimension of a,
-    containing the indices of the non-zero elements in that
-    dimension.  The corresponding non-zero values can be obtained
-    with
+    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
         transpose(a.nonzero())
     instead. The result of this is always a 2d array, with a row for
     each non-zero element."""
-        return numeric.asarray(self.filled(0)).nonzero()
+        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)
@@ -1580,13 +1639,19 @@
             return result.astype(dtype)
         else:
             D = self.diagonal(offset=offset, axis1=axis1, axis2=axis2)
-            return D.astype(dtype).sum(axis=None)
+            return D.astype(dtype).filled(0).sum(axis=None)
     #............................................
     def sum(self, axis=None, dtype=None):
         """a.sum(axis=None, dtype=None)
-Sums the array `a` over the given axis `axis`.
-Masked values are set to 0.
-If `axis` is None, applies to a flattened version of the array.
+    Sums the array a over the given axis.
+    Masked elements are set to 0.
+    
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    dtype : dtype *[None]*
+        Datatype for the intermediary computation.
     """
         if self._mask is nomask:
             mask = nomask
@@ -1601,9 +1666,15 @@
 
     def cumsum(self, axis=None, dtype=None):
         """a.cumprod(axis=None, dtype=None)
-Returns the cumulative sum of the elements of array `a` along the given axis `axis`.
-Masked values are set to 0.
-If `axis` is None, applies to a flattened version of the array.
+    Returns the cumulative sum of the elements of a along the given axis.
+    Masked values are set to 0.
+    
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    dtype : dtype *[None]*
+        Datatype for the intermediary computation.
         """
         result = self.filled(0).cumsum(axis=axis, dtype=dtype).view(type(self))
         result.__setmask__(self.mask)
@@ -1611,9 +1682,15 @@
 
     def prod(self, axis=None, dtype=None):
         """a.prod(axis=None, dtype=None)
-Returns the product of the elements of array `a` along the given axis `axis`.
-Masked elements are set to 1.
-If `axis` is None, applies to a flattened version of the array.
+    Returns the product of the elements of a along the given axis.
+    Masked elements are set to 1.
+    
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    dtype : dtype *[None]*
+        Datatype for the intermediary computation.
         """
         if self._mask is nomask:
             mask = nomask
@@ -1629,9 +1706,15 @@
 
     def cumprod(self, axis=None, dtype=None):
         """a.cumprod(axis=None, dtype=None)
-Returns the cumulative product of ethe lements of array `a` along the given axis `axis`.
-Masked values are set to 1.
-If `axis` is None, applies to a flattened version of the array.
+    Returns the cumulative product of the elements of a along the given axis.
+    Masked values are set to 1.
+    
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    dtype : dtype *[None]*
+        Datatype for the intermediary computation.
         """
         result = self.filled(1).cumprod(axis=axis, dtype=dtype).view(type(self))
         result.__setmask__(self.mask)
@@ -1640,15 +1723,16 @@
     def mean(self, axis=None, dtype=None):
         """a.mean(axis=None, dtype=None)
 
-    Averages the array over the given axis.  If the axis is None,
-    averages over all dimensions of the array.  Equivalent to
+    Averages the array over the given axis.  Equivalent to
 
-      a.sum(axis, dtype) / size(a, axis).
-
-    The optional dtype argument is the data type for intermediate
-    calculations in the sum.
-
-    Returns a masked array, of the same class as a.
+      a.sum(axis, dtype) / size(a, axis).    
+    
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    dtype : dtype *[None]*
+        Datatype for the intermediary computation.
         """
         if self._mask is nomask:
             return super(MaskedArray, self).mean(axis=axis, dtype=dtype)
@@ -1660,6 +1744,13 @@
     def anom(self, axis=None, dtype=None):
         """a.anom(axis=None, dtype=None)
     Returns the anomalies, or deviation from the average.
+    
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    dtype : dtype *[None]*
+        Datatype for the intermediary computation.
             """
         m = self.mean(axis, dtype)
         if not axis:
@@ -1670,9 +1761,21 @@
     def var(self, axis=None, dtype=None):
         """a.var(axis=None, dtype=None)
 Returns 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((x - x.mean())**2).
+    
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    dtype : dtype *[None]*
+        Datatype for the intermediary computation.
+        
+
+Notes
+-----
+The value returned is a biased estimate of the true variance.
+For the more standard unbiased estimate, use varu.        
         """
         if self._mask is nomask:
             # TODO: Do we keep super, or var _data and take a view ?
@@ -1681,9 +1784,10 @@
             cnt = self.count(axis=axis)
             danom = self.anom(axis=axis, dtype=dtype)
             danom *= danom
-            dvar = numeric.array(danom.sum(axis) / cnt).view(type(self))
+            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):
@@ -1691,29 +1795,45 @@
 Returns 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((x - x.mean())**2)).
+deviations from the mean, i.e. std = sqrt(mean((x - x.mean())**2)).    
+
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    dtype : dtype *[None]*
+        Datatype for the intermediary computation.
+        
+
+Notes
+-----
+The value returned is a biased estimate of the true standard deviation.
+For the more standard unbiased estimate, use stdu.        
         """
         dvar = self.var(axis,dtype)
         if axis is not None or dvar is not masked:
             dvar = sqrt(dvar)
         return dvar
+        
     #............................................
     def argsort(self, axis=None, fill_value=None, kind='quicksort',
                 order=None):
-        """Returns an array of indices that sort 'a' along the specified axis.
-    Masked values are filled beforehand to `fill_value`.
-    If `fill_value` is None, uses the default for the data type.
+        """Returns a ndarray of indices that sort 'a' along the specified axis.
+    Masked values are filled beforehand to fill_value.
     Returns a numpy array.
 
-:Keywords:
-    `axis` : Integer *[None]*
-        Axis to be indirectly sorted (default -1)
-    `kind` : String *['quicksort']*
+:Parameters:
+    axis : integer *[None]*
+        Axis to be indirectly sorted.
+    fill_value : var *[None]*
+        Value used to fill in the masked values.
+        If None, use self.fill_value instead.
+    kind : String *['quicksort']*
         Sorting algorithm (default 'quicksort')
         Possible values: 'quicksort', 'mergesort', or 'heapsort'
 
-    Returns: array of indices that sort 'a' along the specified axis.
-
+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.
@@ -1741,17 +1861,16 @@
         return d.argsort(axis=axis, kind=kind, order=order)
     #........................
     def argmin(self, axis=None, fill_value=None):
-        """Returns a ndarray of indices for the minimum values of `a` along the
+        """Returns a 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`.
-    If `fill_value` is None, the default for the data type is used.
-    Returns a numpy array.
+    Masked values are treated as if they had the value fill_value.
 
-:Keywords:
-    `axis` : Integer *[None]*
-        Axis to be indirectly sorted (default -1)
-    `fill_value` : var *[None]*
-        Default filling value. If None, uses the minimum default for the data type.
+:Parameters:
+    axis : integer *[None]*
+        Axis to be indirectly sorted.
+    fill_value : var *[None]*
+        Value used to fill in the masked values.
+        If None, use the the output of minimum_fill_value().
         """
         if fill_value is None:
             fill_value = minimum_fill_value(self)
@@ -1762,14 +1881,13 @@
         """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`.
-    If `fill_value` is None, the maximum default for the data type is used.
-    Returns a numpy array.
 
-:Keywords:
-    `axis` : Integer *[None]*
-        Axis to be indirectly sorted (default -1)
-    `fill_value` : var *[None]*
-        Default filling value. If None, uses the data type default.
+:Parameters:
+    axis : integer *[None]*
+        Axis to be indirectly sorted.
+    fill_value : var *[None]*
+        Value used to fill in the masked values.
+        If None, use the the output of maximum_fill_value().
         """
         if fill_value is None:
             fill_value = maximum_fill_value(self._data)
@@ -1778,24 +1896,31 @@
 
     def sort(self, axis=-1, kind='quicksort', order=None, 
              endwith=True, fill_value=None):
-        """
-        Sort a along the given axis.
+        """Sort a along the given axis.
 
-    Keyword arguments:
+:Parameters:
+    axis : integer *[-1]*
+        Axis to be indirectly sorted.
+    kind : String *['quicksort']*
+        Sorting algorithm (default 'quicksort')
+        Possible values: 'quicksort', 'mergesort', or 'heapsort'.
+    order : var *[None]*
+        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 *[None]*
+        Value used to fill in the masked values.
+        If None, use the the output of minimum_fill_value().
+    endwith : boolean *[True]*
+        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).
 
-    axis  -- axis to be sorted (default -1)
-    kind  -- sorting algorithm (default 'quicksort')
-             Possible values: 'quicksort', 'mergesort', or 'heapsort'.
-    order -- 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.
-    endwith--Boolean flag indicating whether missing values (if any) should
-             be forced in the upper indices (at the end of the array) or
-             lower indices (at the beginning).
+:Returns:
+    When used as method, returns None.
+    When used as a function, returns an array.
 
-    Returns: None.
-
+Notes
+-----
     This method sorts 'a' in place along the given axis using the algorithm
     specified by the kind keyword.
 
@@ -1832,12 +1957,20 @@
             self.flat = tmp_data
             self._mask.flat = tmp_mask
         return
+            
     #............................................
     def min(self, axis=None, fill_value=None):
-        """Returns the minimum/a along the given axis.
-If `axis` is None, applies to the flattened array. Masked values are filled 
-with `fill_value` during processing. If `fill_value is None, it is set to the
-maximum_fill_value corresponding to the data type."""
+        """Returns the minimum of a along the given axis.
+    Masked values are filled with fill_value.
+        
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    fill_value : var *[None]*
+        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:
@@ -1860,9 +1993,16 @@
     #........................
     def max(self, axis=None, fill_value=None):
         """Returns the maximum/a along the given axis.
-If `axis` is None, applies to the flattened array. Masked values are filled 
-with `fill_value` during processing. If `fill_value is None, it is set to the
-maximum_fill_value corresponding to the data type."""
+    Masked values are filled with fill_value.
+        
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    fill_value : var *[None]*
+        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:
@@ -1885,17 +2025,25 @@
     #........................
     def ptp(self, axis=None, fill_value=None):
         """Returns the visible data range (max-min) along the given axis.
-If the axis is `None`, applies on a flattened array. Masked values are filled
-with `fill_value` for processing. If `fill_value` is None, the maximum is uses
-the maximum default, the minimum uses the minimum default."""
+        
+:Parameters:
+    axis : integer *[None]*
+        Axis along which to perform the operation.
+        If None, applies to a flattened version of the array.
+    fill_value : var *[None]*
+        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 ---------------------------------------
-    conj = conjugate = _arraymethod('conjugate')
+#    conj = conjugate = _arraymethod('conjugate')
     copy = _arraymethod('copy')
     diagonal = _arraymethod('diagonal')
     take = _arraymethod('take')
-    ravel = _arraymethod('ravel')
+#    ravel = _arraymethod('ravel')
     transpose = _arraymethod('transpose')
     T = property(fget=lambda self:self.transpose())
     swapaxes = _arraymethod('swapaxes')
@@ -1908,7 +2056,7 @@
         """Copies 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
+    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:
@@ -1931,37 +2079,25 @@
             
             
     #........................
-    def tostring(self, fill_value=None):
-        """a.tostring(order='C', fill_value=None) -> raw copy of array data as a Python string.
+    def tostring(self, fill_value=None, order='C'):
+        """a.tostring(order='C', fill_value=None)
+        
+    Returns a copy of array data as a Python string containing the 
+    raw bytes in the array. 
 
-    Keyword arguments:
-        order      : order of the data item in the copy {"C","F","A"} (default "C")
-        fill_value : value used in lieu of missing data 
-
-    Construct a Python string containing the raw bytes in the array. The order
-    of the data in arrays with ndim > 1 is specified by the 'order' keyword and
-    this keyword overrides the order of the array. The
-    choices are:
-
+:Parameters:
+    fill_value : var *[None]*
+        Value used to fill in the masked values.
+        If None, uses self.fill_value instead.
+    order : string *['C']*
+        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"
-    
-    Masked data are filled with fill_value. If fill_value is None, the data-type-
-    dependent default is used."""
-        return self.filled(fill_value).tostring()   
+    """
+        return self.filled(fill_value).tostring(order=order)   
     #--------------------------------------------
-    # Backwards Compatibility. Heck...
-    @property
-    def data(self):
-        """Returns the `_data` part of the MaskedArray."""
-        return self._data
-    def raw_data(self):
-        """Returns the `_data` part of the MaskedArray.
-You should really use `data` instead..."""
-        return self._data
-    #--------------------------------------------
     # Pickling
     def __getstate__(self):
         "Returns the internal state of the masked array, for pickling purposes."
@@ -2002,7 +2138,7 @@
 in a pickle."""
     _data = ndarray.__new__(baseclass, baseshape, basetype)
     _mask = ndarray.__new__(ndarray, baseshape, 'b1')
-    return subtype.__new__(subtype, _data, mask=_mask, dtype=basetype, small_mask=False)
+    return subtype.__new__(subtype, _data, mask=_mask, dtype=basetype, shrink=False)
 #MaskedArray.__dump__ = dump
 #MaskedArray.__dumps__ = dumps
     
@@ -2016,22 +2152,24 @@
     return isinstance(x, MaskedArray)
 isarray = isMaskedArray
 isMA = isMaskedArray  #backward compatibility
-#masked = MaskedArray(0, int, mask=1)
-masked_singleton = MaskedArray(0, dtype=int_, mask=True)
+# We define the masked singleton as a float for higher precedence...
+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, subok=True,
-          keep_mask=True, small_mask=True, hard_mask=None, fill_value=None):
+          keep_mask=True, hard_mask=False, fill_value=None):
     """array(data, dtype=None, copy=True, order=False, mask=nomask,
-             keep_mask=True, small_mask=True, fill_value=None)
+             keep_mask=True, shrink=True, fill_value=None)
 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, small_mask=small_mask,
-                       hard_mask=hard_mask, fill_value=fill_value)
+                       keep_mask=keep_mask, hard_mask=hard_mask, 
+                       fill_value=fill_value)
+array.__doc__ = masked_array.__doc__
 
 def is_masked(x):
     """Returns whether x has some masked values."""
@@ -2055,7 +2193,7 @@
         return where(self.compare(a, b), a, b)
     #.........
     def reduce(self, target, axis=None):
-        """Reduces target along the given axis."""
+        "Reduces target along the given axis."
         m = getmask(target)
         if axis is not None:
             kargs = { 'axis' : axis }
@@ -2089,6 +2227,7 @@
         result = self.ufunc.outer(filled(a), filled(b))
         result._mask = m
         return result
+        
 #............................
 class _minimum_operation(_extrema_operation):
     "Object to calculate minima"
@@ -2100,6 +2239,7 @@
         self.afunc = amin
         self.compare = less
         self.fill_value_func = minimum_fill_value
+        
 #............................
 class _maximum_operation(_extrema_operation):
     "Object to calculate maxima"
@@ -2111,6 +2251,7 @@
         self.afunc = amax
         self.compare = greater
         self.fill_value_func = maximum_fill_value
+        
 #..........................................................
 def min(array, axis=None, out=None):
     """Returns the minima along the given axis.
@@ -2121,16 +2262,16 @@
         return minimum(array)
     else:
         return minimum.reduce(array, axis)
+min.__doc__ = MaskedArray.min.__doc__
 #............................
 def max(obj, axis=None, out=None):
-    """Returns the maxima 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 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)"""
@@ -2138,6 +2279,7 @@
         return obj.max(axis)-obj.min(axis)
     except AttributeError:
         return max(obj, axis=axis) - min(obj, axis=axis)
+ptp.__doc__ = MaskedArray.ptp.__doc__
 
 
 #####---------------------------------------------------------------------------
@@ -2165,7 +2307,7 @@
         #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(fromnumeric.asarray(a), self._methodname)
+        method = getattr(narray(a, copy=False), self._methodname)
         try:
             return method(*args, **params)
         except SystemError:
@@ -2200,93 +2342,40 @@
     ma = getmask(a)
     mb = getmask(b)
     m = mask_or(ma, mb)
-    fa = filled(a, 1)
-    fb = filled(b, 1)
+    fa = getdata(a)
+    fb = getdata(b)
     if fb.dtype.char in typecodes["Integer"]:
         return masked_array(umath.power(fa, fb), m)
-    md = make_mask((fa < 0), small_mask=1)
+    md = make_mask((fa < 0), shrink=True)
     m = mask_or(m, md)
     if m is nomask:
         return masked_array(umath.power(fa, fb))
     else:
-        fa[m] = 1
+        fa = fa.copy()
+        fa[(fa < 0)] = 1
         return masked_array(umath.power(fa, fb), m)
 
 #..............................................................................
 def argsort(a, axis=None, kind='quicksort', order=None, fill_value=None):
-    """Returns an array of indices that sort 'a' along the specified axis.
-    Masked values are filled beforehand to `fill_value`.
-    If `fill_value` is None, uses the default for the data type.
-    Returns a numpy array.
-
-:Keywords:
-    `axis` : Integer *[None]*
-        Axis to be indirectly sorted (default -1)
-    `kind` : String *['quicksort']*
-        Sorting algorithm (default 'quicksort')
-        Possible values: 'quicksort', 'mergesort', or 'heapsort'
-
-    Returns: array of indices that sort 'a' along the specified axis.
-
-    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.
-    """
+    "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):
-    """Returns the array of indices for the minimum values of `a` along the
-    specified axis.
-    Masked values are treated as if they had the value `fill_value`.
-    If `fill_value` is None, the default for the data type is used.
-    Returns a numpy array.
-
-:Keywords:
-    `axis` : Integer *[None]*
-        Axis to be indirectly sorted (default -1)
-    `fill_value` : var *[None]*
-        Default filling value. If None, uses the data type default.
-    """
+    "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):
-    """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`.
-    If `fill_value` is None, the default for the data type is used.
-    Returns a numpy array.
-
-:Keywords:
-    `axis` : Integer *[None]*
-        Axis to be indirectly sorted (default -1)
-    `fill_value` : var *[None]*
-        Default filling value. If None, uses the data type default.
-    """
+    "Function version of the eponymous method."
     if fill_value is None:
         fill_value = default_fill_value(a)
         try:
@@ -2295,49 +2384,11 @@
             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):
-    """
-    Sort a along the given axis.
-
-Keyword arguments:
-
-axis  -- axis to be sorted (default -1)
-kind  -- sorting algorithm (default 'quicksort')
-         Possible values: 'quicksort', 'mergesort', or 'heapsort'.
-order -- 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.
-endwith--Boolean flag indicating whether missing values (if any) should
-         be forced in the upper indices (at the end of the array) or
-         lower indices (at the beginning).
-
-Returns: None.
-
-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 with 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  |
-|------------------------------------------------------|
-
-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.
-
-"""
-    a = numeric.asanyarray(a)
+    "Function version of the eponymous method."
+    a = narray(a, copy=False, subok=True)
     if fill_value is None:
         if endwith:
             filler = minimum_fill_value(a)
@@ -2349,37 +2400,42 @@
     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):
-    """Returns a compressed version of a masked array (or just the array if it
-    wasn't masked first)."""
-    if getmask(x) is None:
+    """Returns a 1-D array of all the non-masked data."""
+    if getmask(x) is nomask:
         return x
     else:
         return x.compressed()
 
-def count(a, axis = None):
-    "Count of the non-masked elements in a, or along a certain axis."
-    a = masked_array(a)
-    return a.count(axis)
-
 def concatenate(arrays, axis=0):
     "Concatenates the arrays along the given axis"
-    d = numeric.concatenate([filled(a) for a in arrays], 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
+        return data
+    # OK, so we have to concatenate the masks
+    dm = numpy.concatenate([getmaskarray(a) for a in arrays], axis)
+    shrink = numpy.logical_or.reduce([getattr(a,'_shrinkmask',True) for a in arrays])
+    if shrink and not dm.any():
+        data._mask = nomask
     else:
-        return data
-    dm = numeric.concatenate([getmaskarray(a) for a in arrays], axis)
-    dm = make_mask(dm, copy=False, small_mask=True)
-    data._mask = dm
+        data._mask = dm.reshape(d.shape)
     return data
 
+def count(a, axis = None):
+    "Count of the non-masked elements in a, or along a certain axis."
+    return masked_array(a, copy=False).count(axis)
+
+
 def expand_dims(x,axis):
-    """Expand the shape of a by including newaxis before given axis."""
+    """Expands the shape of a by including newaxis before given axis."""
     result = n_expand_dims(x,axis)
     if isinstance(x, MaskedArray):
         new_shape = result.shape
@@ -2409,6 +2465,7 @@
     else:
         d = umath.right_shift(filled(a, 0), n)
         return masked_array(d, mask=m)
+    
 #......................................
 def put(a, indices, values, mode='raise'):
     """Sets storage-indexed locations to corresponding values.
@@ -2417,7 +2474,7 @@
     try:
         return a.put(indices, values, mode=mode)
     except AttributeError:
-        return fromnumeric.asarray(a).put(indices, values, mode=mode)
+        return narray(a, copy=False).put(indices, values, mode=mode)
 
 def putmask(a, mask, values): #, mode='raise'):
     """`putmask(a, mask, v)` results in `a = v` for all places where `mask` is true.
@@ -2427,7 +2484,7 @@
     try:
         return a.putmask(values, mask)
     except AttributeError:
-        return fromnumeric.asarray(a).putmask(values, mask)
+        return narray(a, copy=False).putmask(values, mask)
 
 def transpose(a,axes=None):
     """Returns a view of the array with dimensions permuted according to axes.
@@ -2437,15 +2494,15 @@
     try:
         return a.transpose(axes)
     except AttributeError:
-        return fromnumeric.asarray(a).transpose(axes)
+        return narray(a, copy=False).transpose(axes)
 
 def reshape(a, new_shape):
-    """Changes the shape of the array `a` to `new_shape`."""
+    """Changes 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 fromnumeric.asarray(a).reshape(new_shape)
+        return narray(a, copy=False).reshape(new_shape)
 
 def resize(x, new_shape):
     """resize(a,new_shape) returns a new array with the specified shape.
@@ -2456,8 +2513,8 @@
     # We can't use _frommethods here, as N.resize is notoriously whiny.
     m = getmask(x)
     if m is not nomask:
-        m = fromnumeric.resize(m, new_shape)
-    result = fromnumeric.resize(x, new_shape).view(get_masked_subclass(x))
+        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
@@ -2467,46 +2524,82 @@
 def rank(obj):
     """Gets the rank of sequence a (the number of dimensions, not a matrix rank)
 The rank of a scalar is zero."""
-    return fromnumeric.rank(filled(obj))
+    return fromnumeric.rank(getdata(obj))
 #
 def shape(obj):
     """Returns the shape of `a` (as a function call which also works on nested sequences).
     """
-    return fromnumeric.shape(filled(obj))
+    return fromnumeric.shape(getdata(obj))
 #
 def size(obj, axis=None):
     """Returns the number of elements in the array along the given axis,
 or in the sequence if `axis` is None.
     """
-    return fromnumeric.size(filled(obj), axis)
+    return fromnumeric.size(getdata(obj), axis)
 #................................................
 
 #####--------------------------------------------------------------------------
 #---- --- Extra functions ---
 #####--------------------------------------------------------------------------
-def where (condition, x, y):
-    """where(condition, x, y) is x where condition is nonzero, y otherwise.
-       condition must be convertible to an integer array.
-       Answer is always the shape of condition.
-       The type depends on x and y. It is integer if both x and y are
-       the value masked.
+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.
+    condition must be convertible to an integer array.
+    If neither x nor y are given, returns a tuple of indices where condition is
+    True (a la condition.nonzero()).  
     """
-    fc = filled(not_equal(condition, 0), 0)
-    xv = filled(x)
-    xm = getmask(x)
-    yv = filled(y)
-    ym = getmask(y)
-    d = numeric.choose(fc, (yv, xv))
-    md = numeric.choose(fc, (ym, xm))
-    m = getmask(condition)
-    m = make_mask(mask_or(m, md), copy=False, small_mask=True)
-    return masked_array(d, mask=m)
+    if x is None and y is None:
+        return filled(condition, 0).nonzero()
+    elif x is None or y is None:
+        raise ValueError, "Either bioth 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
+#    # Get the data as a (subclass of) MaskedArray
+#    xv = getdata(x)
+#    yv = getdata(y)
+#    d = numpy.choose(fc, (yv, xv)).view(MaskedArray)
+#    # Get the mask ....................
+#    xm = getmask(x)
+#    ym = getmask(y)
+#    d.mask = numpy.choose(fc, (ym, xm)) | getmask(condition)
+#    # Fix the dtype if one of the values was masked, to prevent an upload to float
+#    if y is masked:
+#        ndtype = xv.dtype
+#    elif x is masked:
+#        ndtype = yv.dtype
+#    else:
+#        ndtype = d.dtype
+#    return d.astype(ndtype)
 
 def choose (indices, t, out=None, mode='raise'):
     "Returns 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``."
+        "Returns the filled array, or True if masked."
         if x is masked:
             return 1
         return filled(x)
@@ -2521,27 +2614,28 @@
     c = filled(indices, 0)
     masks = [nmask(x) for x in t]
     a = [fmask(x) for x in t]
-    d = numeric.choose(c, a)
-    m = numeric.choose(c, masks)
-    m = make_mask(mask_or(m, getmask(indices)), copy=0, small_mask=1)
+    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):
-    """Returns reference to result. Copies a and rounds to 'decimals' places.
+    """Returns a copy of a rounded to 'decimals' places.
 
-    Keyword arguments:
-        decimals -- number of decimals to round to (default 0). May be negative.
-        out -- existing array to use for output (default copy of a).
+:Parameters:
+    decimals : integer *[0]*
+        Number of decimals to round to. May be negative. 
+    out : ndarray
+        Existing array to use for output (default copy of a).
 
-    Return:
-        Reference to out, where None specifies a copy of the original array a.
-
-    Round to the specified number of decimals. 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.
+Notes
+-----
+    Rounds to the specified number of decimals. 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."""
-    result = fromnumeric.round_(filled(a), decimals, out)
+    result = fromnumeric.round_(getdata(a), decimals, out)
     if isinstance(a,MaskedArray):
         result = result.view(type(a))
         result._mask = a._mask
@@ -2549,29 +2643,25 @@
         result = result.view(MaskedArray)
     return result
 
-def arange(start, stop=None, step=1, dtype=None):
-    """Just like range() except it returns a array whose type can be specified
-    by the keyword argument dtype.
-    """
-    return array(numeric.arange(start, stop, step, dtype),mask=nomask)
+def arange(stop, start=None, step=1, dtype=None):
+    "maskedarray version of the numpy function."
+    return numpy.arange(start, stop, step, dtype).view(MaskedArray)
+arange.__doc__ = numpy.arange.__doc__
 
 def inner(a, b):
-    """inner(a,b) returns the dot product of two arrays, which has
-    shape a.shape[:-1] + b.shape[:-1] with elements computed by summing the
-    product of the elements from the last dimensions of a and b.
-    Masked elements are replace by zeros.
-    """
+    "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 masked_array(numeric.inner(fa, fb))
+    return numpy.inner(fa, fb).view(MaskedArray)
+inner.__doc__ = numpy.inner.__doc__ + "\nMasked values are replaced by 0."
 innerproduct = inner
 
 def outer(a, b):
-    """outer(a,b) = {a[i]*b[j]}, has shape (len(a),len(b))"""
+    "maskedarray version of the numpy function."
     fa = filled(a, 0).ravel()
     fb = filled(b, 0).ravel()
     d = numeric.outer(fa, fb)
@@ -2583,22 +2673,22 @@
     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__ + "\nMasked values are replaced by 0."
 outerproduct = outer
 
 def allequal (a, b, fill_value=True):
+    """Returns True if all entries of  a and b are equal, using fill_value 
+    as a truth value where either or both are masked.
     """
-Returns `True` if all entries of  a and b are equal, using
-fill_value as a truth value where either or both are masked.
-    """
     m = mask_or(getmask(a), getmask(b))
     if m is nomask:
-        x = filled(a)
-        y = filled(b)
+        x = getdata(a)
+        y = getdata(b)
         d = umath.equal(x, y)
         return d.all()
     elif fill_value:
-        x = filled(a)
-        y = filled(b)
+        x = getdata(a)
+        y = getdata(b)
         d = umath.equal(x, y)
         dm = array(d, mask=m, copy=False)
         return dm.filled(True).all(None)
@@ -2606,16 +2696,16 @@
         return False
 
 def allclose (a, b, fill_value=True, rtol=1.e-5, atol=1.e-8):
-    """ Returns `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.
+    """ Returns 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.
+The absolute error atol comes into play for those elements of b that are very small 
+or zero; it says how small `a` must be also.
     """
     m = mask_or(getmask(a), getmask(b))
-    d1 = filled(a)
-    d2 = filled(b)
+    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))
@@ -2623,35 +2713,41 @@
 
 #..............................................................................
 def asarray(a, dtype=None):
-    """asarray(data, dtype) = array(data, dtype, copy=0)
-Returns `a` as an masked array.
-No copy is performed if `a` is already an array.
-Subclasses are converted to base class MaskedArray.
+    """asarray(data, dtype) = array(data, dtype, copy=0, subok=0)
+Returns a as a MaskedArray object.
+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)
+    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)
+Returns a as an masked array.
+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):
-    """empty((d1,...,dn),dtype=float,order='C')
-Returns a new array of shape (d1,...,dn) and given type with all its
-entries uninitialized. This can be faster than zeros."""
-    return numeric.empty(new_shape, dtype).view(MaskedArray)
+    "maskedarray version of the numpy function."
+    return numpy.empty(new_shape, dtype).view(MaskedArray)
+empty.__doc__ = numpy.empty.__doc__
 
 def empty_like(a):
-    """empty_like(a)
-Returns an empty (uninitialized) array of the shape and typecode of a.
-Note that this does NOT initialize the returned array.
-If you require your array to be initialized, you should use zeros_like()."""
-    return numeric.empty_like(a).view(MaskedArray)
+    "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):
-    """ones(shape, dtype=None)
-Returns an array of the given dimensions, initialized to all ones."""
-    return numeric.ones(new_shape, dtype).view(MaskedArray)
+    "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):
-    """zeros(new_shape, dtype=None)
-Returns an array of the given dimensions, initialized to all zeros."""
-    return numeric.zeros(new_shape, dtype).view(MaskedArray)
+    "maskedarray version of the numpy function."
+    return numpy.zeros(new_shape, dtype).view(MaskedArray)
+zeros.__doc__ = numpy.zeros.__doc__
 
 #####--------------------------------------------------------------------------
 #---- --- Pickling ---
@@ -2680,82 +2776,16 @@
     return cPickle.loads(strg)
 
 
-################################################################################
+###############################################################################
 
-if __name__ == '__main__':
-    from testutils import assert_equal, assert_almost_equal
-    if 0:
-        x = arange(10)
-        assert(x.ctypes.data == x.filled().ctypes.data)
-    if 0:
-        a = array([1,2,3,4],mask=[0,0,0,0],small_mask=True)
-        a[1] = masked
-        a[1] = 1
-        assert(a.ravel()._mask, [0,0,0,0])
-        assert(a.compressed(), a)
-        a[0] = masked
-        assert(a.compressed()._mask, [0,0,0])
-    if 0:
-        x = array(0, mask=0)
-        I = x.ctypes.data
-        J = x.filled().ctypes.data
-        print (I,J)
-        x = array([0,0], mask=0)
-        (I,J) = (x.ctypes.data, x.filled().ctypes.data)
-        print (I,J)
-    if 0:
-        x = array(numpy.arange(12))
-        x[[1,-2]] = masked
-        xlist = x.tolist()
-        assert(xlist[1] is None)
-        assert(xlist[-2] is None)
-        #
-        x.shape = (3,4) 
-        xlist = x.tolist()
-        #
-        assert_equal(xlist[0],[0,None,2,3])
-        assert_equal(xlist[1],[4,5,6,7])
-        assert_equal(xlist[2],[8,9,None,11])
-        
-    if 0:
-        xl = numpy.random.rand(100,100)
-        yl = numpy.random.rand(100,100)
-        maskx = xl > 0.8
-        masky = yl < 0.2
-        mxl = array(xl, mask=maskx)
-        myl = array(yl, mask=masky)
-        
-        zz = mxl + myl
+#if __name__ == '__main__':
+    #from maskedarray.testutils import assert_equal, assert_almost_equal
     
-    if 0:
-        print "x is ndarray"
-        x = array(numpy.random.rand(50,50))
-        print "set x._mask"
-        x[x > 0.8] = masked
-        print "set y"
-        y = array(numpy.random.rand(50,50))
-        print "set y._mask"
-        ymask = y._data < 0.2
-        print "set y._mask"
-        y.__setmask__(ymask)
-        print "add x + y"
-        z = x + y
-        
-        r.__setmask__()
-        
-    if 1:
-        "Check that we don't lose the fill_value"
-        data = masked_array([1,2,3],fill_value=-999)
-        series = data[[0,2,1]]
-        assert_equal(series._fill_value, data._fill_value)
-        
-    if 1:
-        "Check squeeze"
-        data = masked_array([[1,2,3]])
-        assert_equal(data.squeeze(), [1,2,3])
-        data = masked_array([[1,2,3]], mask=[[1,1,1]])
-        assert_equal(data.squeeze(), [1,2,3])
-        assert_equal(data.squeeze()._mask, [1,1,1])
-        data = masked_array([[1]], mask=True)
-        assert(data.squeeze() is masked)
-        
+    #xm = array(numpy.random.uniform(-1,1,25))
+    #xm[xm>0.5] = masked
+    #xm.fill_value = -999
+    ##
+    #z = 3//where(xm.mask,0,xm)
+    #assert_equal(z._mask, numpy.logical_or(xm==0,xm._mask))
+    #assert_equal(z._data[xm._mask], 1)
+

Modified: trunk/scipy/sandbox/maskedarray/extras.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/extras.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/extras.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -672,12 +672,7 @@
     import numpy as N
     from maskedarray.testutils import assert_equal
     if 1:
-        n = N.arange(1,7)
-        #
-        m = [1,0,0,0,0,0]
-        a = masked_array(n, mask=m).reshape(2,3)
-        b = masked_array(n, mask=m).reshape(3,2)
-        c = dot(a,b, True)
-        assert_equal(c.mask, [[1,1],[1,0]])
-        c = dot(a,b,False)
-        assert_equal(c, N.dot(a.filled(0), b.filled(0)))
\ No newline at end of file
+        b = ones(5)
+        m = [1,0,0,0,0]
+        d = masked_array(b,mask=m)
+        c = mr_[d,0,0,d]
\ No newline at end of file

Modified: trunk/scipy/sandbox/maskedarray/morestats.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/morestats.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/morestats.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -365,6 +365,37 @@
 
 ###############################################################################
 if __name__ == '__main__':
-    data = numpy.arange(100).reshape(4,25)
-#    tmp = hdquantiles(data, prob=[0.25,0.75,0.5], axis=1, var=False)
     
+    if 0:
+        from maskedarray.testutils import assert_almost_equal
+        data = [0.706560797,0.727229578,0.990399276,0.927065621,0.158953014,
+            0.887764025,0.239407086,0.349638551,0.972791145,0.149789972,
+            0.936947700,0.132359948,0.046041972,0.641675031,0.945530547,
+            0.224218684,0.771450991,0.820257774,0.336458052,0.589113496,
+            0.509736129,0.696838829,0.491323573,0.622767425,0.775189248,
+            0.641461450,0.118455200,0.773029450,0.319280007,0.752229111,
+            0.047841438,0.466295911,0.583850781,0.840581845,0.550086491,
+            0.466470062,0.504765074,0.226855960,0.362641207,0.891620942,
+            0.127898691,0.490094097,0.044882048,0.041441695,0.317976349,
+            0.504135618,0.567353033,0.434617473,0.636243375,0.231803616,
+            0.230154113,0.160011327,0.819464108,0.854706985,0.438809221,
+            0.487427267,0.786907310,0.408367937,0.405534192,0.250444460,
+            0.995309248,0.144389588,0.739947527,0.953543606,0.680051621,
+            0.388382017,0.863530727,0.006514031,0.118007779,0.924024803,
+            0.384236354,0.893687694,0.626534881,0.473051932,0.750134705,
+            0.241843555,0.432947602,0.689538104,0.136934797,0.150206859,
+            0.474335206,0.907775349,0.525869295,0.189184225,0.854284286,
+            0.831089744,0.251637345,0.587038213,0.254475554,0.237781276,
+            0.827928620,0.480283781,0.594514455,0.213641488,0.024194386,
+            0.536668589,0.699497811,0.892804071,0.093835427,0.731107772]
+        #
+        assert_almost_equal(hdquantiles(data,[0., 1.]),
+                            [0.006514031, 0.995309248])
+        hdq = hdquantiles(data,[0.25, 0.5, 0.75])
+        assert_almost_equal(hdq, [0.253210762, 0.512847491, 0.762232442,])
+        hdq = hdquantiles_sd(data,[0.25, 0.5, 0.75])
+        assert_almost_equal(hdq, [0.03786954, 0.03805389, 0.03800152,], 4)
+        #
+        data = numpy.array(data).reshape(10,10)
+        hdq = hdquantiles(data,[0.25,0.5,0.75],axis=0)
+    

Modified: trunk/scipy/sandbox/maskedarray/mrecords.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/mrecords.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/mrecords.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -15,6 +15,7 @@
 
 import numpy
 from numpy import bool_, complex_, float_, int_, str_, object_
+from numpy import array as narray
 import numpy.core.numeric as numeric
 import numpy.core.numerictypes as ntypes
 from numpy.core.defchararray import chararray
@@ -27,10 +28,9 @@
 _byteorderconv = numpy.core.records._byteorderconv
 _typestr = ntypes._typestr
 
-import maskedarray as MA
-from maskedarray import masked, nomask, mask_or, filled, getmask, getmaskarray, \
-    masked_array, make_mask
-from maskedarray import MaskedArray
+import maskedarray
+from maskedarray import MaskedArray, masked, nomask, masked_array,\
+    make_mask, mask_or, getmask, getmaskarray, filled     
 from maskedarray.core import default_fill_value, masked_print_option
 
 import warnings
@@ -194,9 +194,12 @@
         if attr in _names:
             _data = self._data
             _mask = self._fieldmask
-            obj = numeric.asarray(_data.__getattribute__(attr)).view(MaskedArray)
-            obj.__setmask__(_mask.__getattribute__(attr))
-            if (obj.ndim == 0) and obj._mask:
+#            obj = masked_array(_data.__getattribute__(attr), copy=False,
+#                               mask=_mask.__getattribute__(attr))
+            # Use a view in order to avoid the copy of the mask in MaskedArray.__new__
+            obj = narray(_data.__getattribute__(attr), copy=False).view(MaskedArray)
+            obj._mask = _mask.__getattribute__(attr)
+            if not obj.ndim and obj._mask:
                 return masked
             return obj
         raise AttributeError,"No attribute '%s' !" % attr
@@ -213,7 +216,7 @@
                 exctype, value = sys.exc_info()[:2]
                 raise exctype, value
         else:
-            if attr not in list(self.dtype.names) + ['_mask']:
+            if attr not in list(self.dtype.names) + ['_mask','mask']:
                 return ret
             if newattr:         # We just added this one
                 try:            #  or this setattr worked on an internal
@@ -245,6 +248,9 @@
         if isinstance(indx, str):           
             obj = _data[indx].view(MaskedArray)
             obj._set_mask(_localdict['_fieldmask'][indx])
+            # Force to nomask if the mask is empty
+            if not obj._mask.any():
+                obj._mask = nomask
             return obj
         # We want some elements ..
         # First, the data ........
@@ -309,6 +315,7 @@
         else:
             for n in names:
                 fmask[n].flat = newmask
+        return
                 
     def _getmask(self):
         """Returns the mask of the mrecord: a record is masked when all the fields
@@ -435,7 +442,7 @@
        
 
     """
-    arraylist = [MA.asarray(x) for x in arraylist]
+    arraylist = [masked_array(x) for x in arraylist]
     # Define/check the shape.....................
     if shape is None or shape == 0:
         shape = arraylist[0].shape
@@ -611,7 +618,7 @@
     if varnames is None:
         varnames = _varnames
     # Get the data ..............................
-    _variables = MA.asarray([line.strip().split(delimitor) for line in f
+    _variables = masked_array([line.strip().split(delimitor) for line in f
                                   if line[0] != commentchar and len(line) > 1])
     (_, nfields) = _variables.shape
     # Try to guess the dtype ....................
@@ -643,7 +650,7 @@
     _mask = mrecord._fieldmask
     if newfieldname is None or newfieldname in reserved_fields:
         newfieldname = 'f%i' % len(_data.dtype)
-    newfield = MA.asarray(newfield)
+    newfield = masked_array(newfield)
     # Get the new data ............
     # Create a new empty recarray
     newdtype = numeric.dtype(_data.dtype.descr + \
@@ -674,58 +681,23 @@
     from maskedarray.testutils import assert_equal
     if 1:
         d = N.arange(5)
-        m = MA.make_mask([1,0,0,1,1])
+        m = maskedarray.make_mask([1,0,0,1,1])
         base_d = N.r_[d,d[::-1]].reshape(2,-1).T
         base_m = N.r_[[m, m[::-1]]].T
-        base = MA.array(base_d, mask=base_m).T
+        base = masked_array(base_d, mask=base_m).T
         mrecord = fromarrays(base,dtype=[('a',N.float_),('b',N.float_)])
         mrec = MaskedRecords(mrecord.copy())
         #
-        mrec.a[3:] = 5
-        assert_equal(mrec.a, [0,1,2,5,5])
-        assert_equal(mrec.a._mask, [1,0,0,0,0])
-        #
-        mrec.b[3:] = masked
-        assert_equal(mrec.b, [4,3,2,1,0])
-        assert_equal(mrec.b._mask, [1,1,0,1,1])
-        #
-        mrec[:2] = masked
-        assert_equal(mrec._mask, [1,1,0,0,0])
-        mrec[-1] = masked
-        assert_equal(mrec._mask, [1,1,0,0,1])
-    if 1:
-        nrec = N.core.records.fromarrays(N.r_[[d,d[::-1]]],
-                                         dtype=[('a',N.float_),('b',N.float_)])
-        mrec = mrecord
-        #....................
-        mrecfr = fromrecords(nrec)
-        assert_equal(mrecfr.a, mrec.a)
-        assert_equal(mrecfr.dtype, mrec.dtype)
-        #....................
-        tmp = mrec[::-1] #.tolist()
-        mrecfr = fromrecords(tmp)
-        assert_equal(mrecfr.a, mrec.a[::-1])
-        #....................        
-        mrecfr = fromrecords(nrec.tolist(), names=nrec.dtype.names)
-        assert_equal(mrecfr.a, mrec.a)
-        assert_equal(mrecfr.dtype, mrec.dtype)
-    if 0:
-        assert_equal(mrec.a, MA.array(d,mask=m))
-        assert_equal(mrec.b, MA.array(d[::-1],mask=m[::-1]))
-        assert((mrec._fieldmask == N.core.records.fromarrays([m, m[::-1]])).all())
+    if 1:   
+        mrec = mrec.copy()
+        mrec.harden_mask()
+        assert(mrec._hardmask)
+        mrec._mask = nomask
         assert_equal(mrec._mask, N.r_[[m,m[::-1]]].all(0))
-        assert_equal(mrec.a[1], mrec[1].a)
-
-    if 0:
-        x = [(1.,10.,'a'),(2.,20,'b'),(3.14,30,'c'),(5.55,40,'d')]
-        desc = [('ffloat', N.float_), ('fint', N.int_), ('fstr', 'S10')] 
-        mr = MaskedRecords(x,dtype=desc)
-        mr[0] = masked
-        mr.ffloat[-1] = masked
-        #
-        mrlast = mr[-1]
-        assert(isinstance(mrlast,MaskedRecords))
-        assert(hasattr(mrlast,'ffloat'))
-        assert_equal(mrlast.ffloat, masked)
-        
+        mrec.soften_mask()
+        assert(not mrec._hardmask)
+        mrec.mask = nomask
+        tmp = mrec['b']._mask
+        assert(mrec['b']._mask is nomask)
+        assert_equal(mrec['a']._mask,mrec['b']._mask)   
         
\ No newline at end of file

Modified: trunk/scipy/sandbox/maskedarray/mstats.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/mstats.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/mstats.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -19,7 +19,7 @@
 import numpy.core.numeric as numeric
 from numpy.core.numeric import concatenate
 
-import maskedarray as MA
+import maskedarray
 from maskedarray.core import masked, nomask, MaskedArray, masked_array
 from maskedarray.extras import apply_along_axis, dot
 
@@ -391,4 +391,14 @@
     nlo = (data[:,None] < points[None,:] - h).sum(0)
     return (nhi-nlo) / (2.*n*h)
 
-    
\ No newline at end of file
+################################################################################
+if __name__ == '__main__':
+    from maskedarray.testutils import assert_almost_equal
+    if 1:
+        a = maskedarray.arange(1,101)
+        a[1::2] = masked
+        b = maskedarray.resize(a, (100,100))
+        assert_almost_equal(mquantiles(b), [25., 50., 75.])
+        assert_almost_equal(mquantiles(b, axis=0), maskedarray.resize(a,(3,100)))
+        assert_almost_equal(mquantiles(b, axis=1), 
+                            maskedarray.resize([24.9, 50., 75.1], (100,3)))        
\ No newline at end of file

Modified: trunk/scipy/sandbox/maskedarray/tests/test_core.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/tests/test_core.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/tests/test_core.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -229,6 +229,17 @@
         assert_equal(x,y/a)
         assert_equal(xm,y/a)
         assert_equal(xm.mask, mask_or(mask_or(m,a.mask), (a==0))) 
+        #
+        (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
+        z = xm/ym
+        assert_equal(z._mask, [1,1,1,0,0,1,1,0,0,0,1,1])
+        assert_equal(z._data, [0.2,1.,1./3.,-1.,-pi/2.,-1.,5.,1.,1.,1.,2.,1.])
+        xm = xm.copy()
+        xm /= ym
+        assert_equal(xm._mask, [1,1,1,0,0,1,1,0,0,0,1,1])
+        assert_equal(xm._data, [1/5.,1.,1./3.,-1.,-pi/2.,-1.,5.,1.,1.,1.,2.,1.])
+        
+        
     #..........................
     def check_scalararithmetic(self):
         "Tests some scalar arithmetics on MaskedArrays."
@@ -479,6 +490,46 @@
 #        assert_not_equal(id(y._mask), id(x._mask))
         assert_not_equal(y._data.ctypes.data, x._data.ctypes.data)
         assert_not_equal(y._mask.ctypes.data, x._mask.ctypes.data)
+    #........................
+    def check_where(self):
+        "Test the where function"
+        (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
+        d = where(xm>2,xm,-9)
+        assert_equal(d, [-9.,-9.,-9.,-9., -9., 4., -9., -9., 10., -9., -9., 3.])
+        assert_equal(d._mask, xm._mask)
+        d = where(xm>2,-9,ym)
+        assert_equal(d, [5.,0.,3., 2., -1.,-9.,-9., -10., -9., 1., 0., -9.])
+        assert_equal(d._mask, [1,0,1,0,0,0,1,0,0,0,0,0])
+        d = where(xm>2, xm, masked)
+        assert_equal(d, [-9.,-9.,-9.,-9., -9., 4., -9., -9., 10., -9., -9., 3.])
+        tmp = xm._mask.copy()
+        tmp[(xm<=2).filled(True)] = True
+        assert_equal(d._mask, tmp)
+        #
+        ixm = xm.astype(int_)
+        d = where(ixm>2, ixm, masked)
+        assert_equal(d, [-9,-9,-9,-9, -9, 4, -9, -9, 10, -9, -9, 3])
+        assert_equal(d.dtype, ixm.dtype)
+        #
+        x = arange(10)
+        x[3] = masked
+        c = x >= 8
+        z = where(c , x, masked)
+        assert z.dtype is x.dtype
+        assert z[3] is masked
+        assert z[4] is masked
+        assert z[7] is masked
+        assert z[8] is not masked
+        assert z[9] is not masked
+        assert_equal(x,z)
+        #
+        z = where(c , masked, x)
+        assert z.dtype is x.dtype
+        assert z[3] is masked
+        assert z[4] is not masked
+        assert z[7] is not masked
+        assert z[8] is masked
+        assert z[9] is masked
         
     #........................
     def check_oddfeatures_1(self):
@@ -500,23 +551,6 @@
         assert count(where(c,masked,masked)) == 0
         assert shape(where(c,masked,masked)) == c.shape
         #
-        z = where(c , x, masked)
-        assert z.dtype is x.dtype
-        assert z[3] is masked
-        assert z[4] is masked
-        assert z[7] is masked
-        assert z[8] is not masked
-        assert z[9] is not masked
-        assert_equal(x,z)
-        #
-        z = where(c , masked, x)
-        assert z.dtype is x.dtype
-        assert z[3] is masked
-        assert z[4] is not masked
-        assert z[7] is not masked
-        assert z[8] is masked
-        assert z[9] is masked
-        #
         z = masked_where(c, x)
         assert z.dtype is x.dtype
         assert z[3] is masked
@@ -631,6 +665,8 @@
         y = x * masked
         assert_equal(y.shape, x.shape)
         assert_equal(y._mask, [True, True])
+        y = x[0] * masked
+        assert y is masked
         y = x + masked
         assert_equal(y.shape, x.shape)
         assert_equal(y._mask, [True, True])
@@ -734,6 +770,19 @@
         data = masked_array([1,2,3],fill_value=-999)
         series = data[[0,2,1]]
         assert_equal(series._fill_value, data._fill_value)
+    #
+    def check_asarray(self):
+        (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
+        xmm = asarray(xm)
+        assert_equal(xmm._data, xm._data)
+        assert_equal(xmm._mask, xm._mask)
+    #
+    def check_fix_invalid(self):
+        "Checks fix_invalid."
+        data = masked_array(numpy.sqrt([-1., 0., 1.]), mask=[0,0,1])
+        data_fixed = fix_invalid(data)
+        assert_equal(data_fixed._data, [data.fill_value, 0., 1.])
+        assert_equal(data_fixed._mask, [1., 0., 1.])    
         
 #...............................................................................
         
@@ -784,6 +833,7 @@
         assert(sometrue(a,axis=0))
         assert_equal(sum(a[:3],axis=0), 0)
         assert_equal(product(a,axis=0), 0)
+        assert_equal(add.reduce(a), pi)
     #........................
     def test_minmax(self):
         "Tests extrema on MaskedArrays."
@@ -929,7 +979,7 @@
         assert( x[3] is masked)
         assert( x[4] is masked)
         x[[1,4]] = [10,40]
-        assert( x.mask is not m)
+#        assert( x.mask is not m)
         assert( x[3] is masked)
         assert( x[4] is not masked)
         assert_equal(x, [0,10,2,-1,40])        
@@ -1202,7 +1252,7 @@
         assert_equal(a.shape,(1,5))
         assert_equal(a._mask.shape, a.shape)
         # Checs that small_mask is preserved
-        a = array([1,2,3,4],mask=[0,0,0,0],small_mask=False)
+        a = array([1,2,3,4],mask=[0,0,0,0],shrink=False)
         assert_equal(a.ravel()._mask, [0,0,0,0])
         
     def check_reshape(self):
@@ -1217,17 +1267,13 @@
         
     def check_compressed(self):
         "Tests compressed"
-        a = array([1,2,3,4],mask=[0,0,0,0],small_mask=False)
+        a = array([1,2,3,4],mask=[0,0,0,0])
         b = a.compressed()
         assert_equal(b, a)
-        assert_equal(b._mask, a._mask)
+        assert_equal(b._mask, nomask)
         a[0] = masked
         b = a.compressed()
         assert_equal(b._data, [2,3,4])
-        assert_equal(b._mask, [0,0,0])
-        a._smallmask = True
-        b = a.compressed()
-        assert_equal(b._data, [2,3,4])
         assert_equal(b._mask, nomask)
         
     def check_tolist(self):

Modified: trunk/scipy/sandbox/maskedarray/tests/test_mrecords.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/tests/test_mrecords.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/tests/test_mrecords.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -20,7 +20,9 @@
 import maskedarray.testutils
 from maskedarray.testutils import *
 
-import maskedarray.core as MA
+import maskedarray
+from maskedarray import masked_array, masked, nomask
+
 #import maskedarray.mrecords
 #from maskedarray.mrecords import mrecarray, fromarrays, fromtextfile, fromrecords
 import maskedarray.mrecords
@@ -37,10 +39,10 @@
     def setup(self):       
         "Generic setup" 
         d = N.arange(5)
-        m = MA.make_mask([1,0,0,1,1])
+        m = maskedarray.make_mask([1,0,0,1,1])
         base_d = N.r_[d,d[::-1]].reshape(2,-1).T
         base_m = N.r_[[m, m[::-1]]].T
-        base = MA.array(base_d, mask=base_m)    
+        base = masked_array(base_d, mask=base_m)    
         mrecord = fromarrays(base.T, dtype=[('a',N.float_),('b',N.float_)])
         self.data = [d, m, mrecord]
         
@@ -48,8 +50,8 @@
         "Tests fields retrieval"
         [d, m, mrec] = self.data
         mrec = mrec.copy()
-        assert_equal(mrec.a, MA.array(d,mask=m))
-        assert_equal(mrec.b, MA.array(d[::-1],mask=m[::-1]))
+        assert_equal(mrec.a, masked_array(d,mask=m))
+        assert_equal(mrec.b, masked_array(d[::-1],mask=m[::-1]))
         assert((mrec._fieldmask == N.core.records.fromarrays([m, m[::-1]], dtype=mrec._fieldmask.dtype)).all())
         assert_equal(mrec._mask, N.r_[[m,m[::-1]]].all(0))
         assert_equal(mrec.a[1], mrec[1].a)
@@ -65,13 +67,13 @@
         mrecord.a = 1
         assert_equal(mrecord['a']._data, [1]*5)
         assert_equal(getmaskarray(mrecord['a']), [0]*5)
-        mrecord.b = MA.masked
+        mrecord.b = masked
         assert_equal(mrecord.b.mask, [1]*5)
         assert_equal(getmaskarray(mrecord['b']), [1]*5)
-        mrecord._mask = MA.masked
+        mrecord._mask = masked
         assert_equal(getmaskarray(mrecord['b']), [1]*5)
         assert_equal(mrecord['a']._mask, mrecord['b']._mask)
-        mrecord._mask = MA.nomask
+        mrecord._mask = nomask
         assert_equal(getmaskarray(mrecord['b']), [0]*5)
         assert_equal(mrecord['a']._mask, mrecord['b']._mask)   
         #

Modified: trunk/scipy/sandbox/maskedarray/tests/test_subclassing.py
===================================================================
--- trunk/scipy/sandbox/maskedarray/tests/test_subclassing.py	2007-09-26 19:48:46 UTC (rev 3370)
+++ trunk/scipy/sandbox/maskedarray/tests/test_subclassing.py	2007-09-27 03:34:44 UTC (rev 3371)
@@ -45,8 +45,8 @@
         _data.info = subarr.info
         return _data
     def __array_finalize__(self,obj):
-        SubArray.__array_finalize__(self, obj)
-        MaskedArray.__array_finalize__(self,obj)    
+        MaskedArray.__array_finalize__(self,obj)   
+        SubArray.__array_finalize__(self, obj) 
         return
     def _get_series(self):
         return self.view(MaskedArray)
@@ -136,15 +136,38 @@
         assert isinstance(mxsub, MaskedArray)
         assert_equal(mxsub._mask, m)
         #
+        mxsub = asarray(xsub)
+        assert not isinstance(mxsub, MSubArray)
+        assert isinstance(mxsub, MaskedArray)
+        assert_equal(mxsub._mask, m)
+        #
         mxsub = masked_array(xsub, subok=True)
         assert isinstance(mxsub, MSubArray)
         assert_equal(mxsub.info, xsub.info)
         assert_equal(mxsub._mask, xsub._mask)
+        #
+        mxsub = asanyarray(xsub)
+        assert isinstance(mxsub, MSubArray)
+        assert_equal(mxsub.info, xsub.info)
+        assert_equal(mxsub._mask, m)
         
         
 ################################################################################
 if __name__ == '__main__':
     NumpyTest().run()
+    #
+    if 0:
+        x = array(arange(5), mask=[0]+[1]*4)
+        my = masked_array(subarray(x))
+        ym = msubarray(x)
+        #
+        z = (my+1)
+        assert isinstance(z,MaskedArray)
+        assert not isinstance(z, MSubArray)
+        assert isinstance(z._data, SubArray)
+        assert_equal(z._data.info, {})
+        #
+        z = ym+1
 
                      
 



More information about the Scipy-svn mailing list