[Scipy-svn] r2700 - in trunk/Lib/sandbox/maskedarray: . tests

scipy-svn at scipy.org scipy-svn at scipy.org
Sun Feb 11 18:24:28 EST 2007


Author: pierregm
Date: 2007-02-11 17:24:17 -0600 (Sun, 11 Feb 2007)
New Revision: 2700

Added:
   trunk/Lib/sandbox/maskedarray/core_ini.py
   trunk/Lib/sandbox/maskedarray/tests/test_subclassing.py
Modified:
   trunk/Lib/sandbox/maskedarray/core.py
   trunk/Lib/sandbox/maskedarray/extras.py
   trunk/Lib/sandbox/maskedarray/tests/test_core.py
   trunk/Lib/sandbox/maskedarray/tests/test_extras.py
   trunk/Lib/sandbox/maskedarray/testutils.py
   trunk/Lib/sandbox/maskedarray/timer_comparison.py
Log:
major update !
core     : MaskedArray now really behave like ndarrays.
         : got rid of (most) class default.
         : now using views instead of calls to self.__class__.
         : added a _smallmask attribute to prevent the automatic compression of masks.
         : added a _baseclass attribute, to record the class type of the input data.
         :  _data is now simply a view of self as _baseclass.
extras   : corrected a pb w/ apply_along_axis.
core_ini : previous version of core

Modified: trunk/Lib/sandbox/maskedarray/core.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/core.py	2007-02-11 00:33:13 UTC (rev 2699)
+++ trunk/Lib/sandbox/maskedarray/core.py	2007-02-11 23:24:17 UTC (rev 2700)
@@ -1,3 +1,4 @@
+# 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.
@@ -22,36 +23,36 @@
 __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',  
+               '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', 
+           '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', 
+           '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', 
+           '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', 
+           '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_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', 
+               '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', 
+           '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', 
+           'take', 'tan', 'tanh', 'transpose', 'true_divide',
            'var', 'where',
            'zeros']
 
@@ -62,7 +63,7 @@
 import numpy
 from numpy import bool_, complex_, float_, int_, object_, str_
 
-import numpy.core.umath as umath 
+import numpy.core.umath as umath
 import numpy.core.fromnumeric  as fromnumeric
 from numpy.core.numeric import ndarray
 from numpy.core.fromnumeric import amax, amin
@@ -81,38 +82,38 @@
 divide_tolerance = 1.e-35
 numpy.seterr(all='ignore')
 
-#####--------------------------------------------------------------------------
-#---- --- Helper functions ---
-#####--------------------------------------------------------------------------
-def convert_typecode(f,dtchar):
-    """Converts the type of `f` to a type compatible with `dtchar`, for inline operations."""
-    ftype = f.dtype.char
-    if dtchar == ftype:
-        return f
-    elif dtchar in typecodes['Integer']:
-        if ftype in typecodes['Integer']:
-            f = f.astype(dtchar)
-        else:
-            raise TypeError, 'Incorrect type for in-place operation.'
-    elif dtchar in typecodes['Float']:
-        if ftype in typecodes['Integer']:
-            f = f.astype(dtchar)
-        elif ftype in typecodes['Float']:
-            f = f.astype(dtchar)
-        else:
-            raise TypeError, 'Incorrect type for in-place operation.'
-    elif dtchar in typecodes['Complex']:
-        if ftype in typecodes['Integer']:
-            f = f.astype(dtchar)
-        elif ftype in typecodes['Float']:
-            f = f.astype(dtchar)
-        elif ftype in typecodes['Complex']:
-            f = f.astype(dtchar)
-        else:
-            raise TypeError, 'Incorrect type for in-place operation.'
-    else:
-        raise TypeError, 'Incorrect type for in-place operation.'
-    return f
+######--------------------------------------------------------------------------
+##---- --- Helper functions ---
+######--------------------------------------------------------------------------
+#def convert_typecode(f,dtchar):
+#    """Converts the type of `f` to a type compatible with `dtchar`, for inline operations."""
+#    ftype = f.dtype.char
+#    if dtchar == ftype:
+#        return f
+#    elif dtchar in typecodes['Integer']:
+#        if ftype in typecodes['Integer']:
+#            f = f.astype(dtchar)
+#        else:
+#            raise TypeError, 'Incorrect type for in-place operation.'
+#    elif dtchar in typecodes['Float']:
+#        if ftype in typecodes['Integer']:
+#            f = f.astype(dtchar)
+#        elif ftype in typecodes['Float']:
+#            f = f.astype(dtchar)
+#        else:
+#            raise TypeError, 'Incorrect type for in-place operation.'
+#    elif dtchar in typecodes['Complex']:
+#        if ftype in typecodes['Integer']:
+#            f = f.astype(dtchar)
+#        elif ftype in typecodes['Float']:
+#            f = f.astype(dtchar)
+#        elif ftype in typecodes['Complex']:
+#            f = f.astype(dtchar)
+#        else:
+#            raise TypeError, 'Incorrect type for in-place operation.'
+#    else:
+#        raise TypeError, 'Incorrect type for in-place operation.'
+#    return f
 
 #####--------------------------------------------------------------------------
 #---- --- Exceptions ---
@@ -139,56 +140,36 @@
                   'O' : '?',
                   'S' : 'N/A',
                   'u' : 999999,
-                  'V' : '???',        
+                  'V' : '???',
                   }
-#{0: <type 'numpy.bool_'>,
-# 1: <type 'numpy.int8'>,
-# 2: <type 'numpy.uint8'>,
-# 3: <type 'numpy.int16'>,
-# 4: <type 'numpy.uint16'>,
-# 5: <type 'numpy.int32'>,
-# 6: <type 'numpy.uint32'>,
-# 7: <type 'numpy.int64'>,
-# 8: <type 'numpy.uint64'>,
-# 9: <type 'numpy.int64'>,
-# 10: <type 'numpy.uint64'>,
-# 11: <type 'numpy.float32'>,
-# 12: <type 'numpy.float64'>,
-# 13: <type 'numpy.float128'>,
-# 14: <type 'numpy.complex64'>,
-# 15: <type 'numpy.complex128'>,
-# 16: <type 'numpy.complex256'>,
-# 17: <type 'numpy.object_'>,
-# 18: <type 'numpy.string_'>,
-# 19: <type 'numpy.unicode_'>,
-# 20: <type 'numpy.void'>,
-max_filler = ntypes._minvals   
-max_filler.update([(k,-numeric.inf) for k in [numpy.float32, numpy.float64]])         
+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]])    
+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)])       
+    min_filler.update([(numpy.float128, numeric.inf)])
 
 
-def default_fill_value (obj):
+def default_fill_value(obj):
     "Calculates the default fill value for an object `obj`."
     if hasattr(obj,'dtype'):
-        return default_filler[obj.dtype.kind]
+        defval = default_filler[obj.dtype.kind]
+    elif isinstance(obj, numeric.dtype):
+        defval = default_filler[obj.kind]
     elif isinstance(obj, float):
-        return default_filler['f']
+        defval = default_filler['f']
     elif isinstance(obj, int) or isinstance(obj, long):
-        return default_filler['i']
+        defval = default_filler['i']
     elif isinstance(obj, str):
-        return default_filler['S']
+        defval = default_filler['S']
     elif isinstance(obj, complex):
-        return default_filler['c']
-    elif isinstance(obj, numeric.dtype):
-        return default_filler[obj.kind]
+        defval = default_filler['c']
     else:
-        return default_filler['O']
+        defval = default_filler['O']
+    return defval
 
-def minimum_fill_value (obj):
+def minimum_fill_value(obj):
     "Calculates the default fill value suitable for taking the minimum of `obj`."
     if hasattr(obj, 'dtype'):
         objtype = obj.dtype
@@ -207,7 +188,7 @@
     else:
         raise TypeError, 'Unsuitable type for calculating minimum.'
 
-def maximum_fill_value (obj):
+def maximum_fill_value(obj):
     "Calculates the default fill value suitable for taking the maximum of `obj`."
     if hasattr(obj, 'dtype'):
         objtype = obj.dtype
@@ -226,12 +207,12 @@
     else:
         raise TypeError, 'Unsuitable type for calculating minimum.'
 
-def set_fill_value (a, fill_value):
+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):
+def get_fill_value(a):
     """Returns the fill value of `a`, if any.
     Otherwise, returns the default fill value for that type.
     """
@@ -241,11 +222,11 @@
         result = default_fill_value(a)
     return result
 
-def common_fill_value (a, b):
+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: 
+    if t1 == t2:
         return t1
     return None
 
@@ -257,18 +238,37 @@
 
 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 
+`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, types.DictType):
+    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 ---
 #####--------------------------------------------------------------------------
@@ -276,7 +276,7 @@
 ufunc_fills = {}
 
 class domain_check_interval:
-    """Defines a valid 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"
@@ -312,7 +312,7 @@
     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)
@@ -329,12 +329,12 @@
 #..............................................................................
 class masked_unary_operation:
     """Defines masked version of unary operations,
-where invalid values are pre-masked. 
+where invalid values are pre-masked.
 
 :IVariables:
     - `f` : function.
     - `fill` : Default filling value *[0]*.
-    - `domain` : Default domain *[None]*. 
+    - `domain` : Default domain *[None]*.
     """
     def __init__ (self, mufunc, fill=0, domain=None):
         """ masked_unary_operation(aufunc, fill=0, domain=None)
@@ -357,24 +357,31 @@
         d1 = filled(a, self.fill)
         if self.domain is not None:
             m = mask_or(m, numeric.asarray(self.domain(d1)))
-        result = self.f(d1, *args, **kwargs)
-        #
-        if isinstance(result, MaskedArray):
-            return result.__class__(result, mask=m)
-        return masked_array(result, mask=m)
+        # 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. 
+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]*. 
+    - `domain` : Default domain *[None]*.
     """
     def __init__ (self, mbfunc, fillx=0, filly=0):
         """abfunc(fillx, filly) must be defined.
@@ -391,21 +398,19 @@
     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)
-        result = self.f(d1, d2, *args, **kwargs)
-#        if isinstance(result, ndarray) \
-#               and m.ndim != 0 \
-#               and m.shape != result.shape:
-#            m = mask_or(getmaskarray(a), getmaskarray(b))
-        if isinstance(result, MaskedArray):
-            return result.__class__(result, mask=m)
-        return masked_array(result, mask=m)
+        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 = target.__class__
+            tclass = type(target)
         else:
             tclass = MaskedArray
         m = getmask(target)
@@ -416,19 +421,20 @@
                 m = make_mask(m, copy=1)
                 m.shape = (1,)
         if m is nomask:
-            return tclass(self.f.reduce (t, axis))
-        else:
-            t = tclass(t, mask=m)
-            # XXX: "or t.dtype" below is a workaround for what appears
-            # XXX: to be a bug in reduce.
-            t = self.f.reduce(filled(t, self.filly), axis, dtype=dtype or t.dtype)
-            m = umath.logical_and.reduce(m, axis)
-            if isinstance(t, ndarray):
-                return tclass(t, mask=m, fill_value=get_fill_value(target))
-            elif m:
-                return masked
-            else:
-                return t
+            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."
@@ -440,34 +446,37 @@
             ma = getmaskarray(a)
             mb = getmaskarray(b)
             m = umath.logical_or.outer(ma, mb)
-        d = self.f.outer(filled(a, self.fillx), filled(b, self.filly))
-        if isinstance(d, MaskedArray):
-            return d.__class__(d, mask=m)
-        return masked_array(d, mask=m)
+        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 = target.__class__
+            tclass = type(target)
         else:
             tclass = masked_array
         t = filled(target, self.filly)
-        return tclass(self.f.accumulate(t, axis))
-    
+        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. 
+    """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]*. 
+    - `domain` : Default domain *[None]*.
     """
     def __init__ (self, dbfunc, domain, fillx=0, filly=0):
         """abfunc(fillx, filly) must be defined.
@@ -494,10 +503,12 @@
             d2 = numeric.where(t, self.filly, d2)
             mb = mask_or(mb, t)
         m = mask_or(ma, mb)
-        result =  self.f(d1, d2)
-        if isinstance(result, MaskedArray):
-            return result.__class__(result, mask=m) 
-        return masked_array(result, mask=m)
+        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)
@@ -526,12 +537,12 @@
 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, 
+arcsin = masked_unary_operation(umath.arcsin, 0.0,
                                 domain_check_interval(-1.0, 1.0))
-arccos = masked_unary_operation(umath.arccos, 0.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, 
+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)
@@ -561,11 +572,11 @@
 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, 
+true_divide = domained_binary_operation(umath.true_divide,
                                         domain_safe_divide(), 0, 1)
-floor_divide = domained_binary_operation(umath.floor_divide, 
+floor_divide = domained_binary_operation(umath.floor_divide,
                                          domain_safe_divide(), 0, 1)
-remainder = domained_binary_operation(umath.remainder, 
+remainder = domained_binary_operation(umath.remainder,
                                       domain_safe_divide(), 0, 1)
 fmod = domained_binary_operation(umath.fmod, domain_safe_divide(), 0, 1)
 
@@ -604,7 +615,7 @@
 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`. 
+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.
 
@@ -651,11 +662,11 @@
     - `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: 
+    if m1 is nomask:
         return make_mask(m2, copy=copy, small_mask=small_mask)
-    if m2 is nomask: 
+    if m2 is nomask:
         return make_mask(m1, copy=copy, small_mask=small_mask)
-    if m1 is m2 and is_mask(m1): 
+    if m1 is m2 and is_mask(m1):
         return m1
     return make_mask(umath.logical_or(m1, m2), copy=copy, small_mask=small_mask)
 
@@ -710,7 +721,7 @@
 #    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 
+    """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.
     """
@@ -721,7 +732,7 @@
     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 
+    """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.
     """
@@ -734,7 +745,7 @@
 #
 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, 
+This function is suitable only for `object` arrays: for floating point,
 please use `masked_values` instead.
 The mask is set to `nomask` if posible.
 
@@ -810,6 +821,156 @@
 #####--------------------------------------------------------------------------
 #---- --- 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)
+        if result.ndim:
+            if not self._onmask:
+                result._mask = mask
+            elif mask is not nomask:
+                result._set_mask(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.
@@ -819,8 +980,8 @@
               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 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.
@@ -828,81 +989,56 @@
 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. 
+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 
+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
-    #TODO: There some reorganization to do round here
+    _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,
                 **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
-        # 1. Argument is MA ...........
-        if isinstance(data, MaskedArray) or\
-           (hasattr(data,"_mask") and hasattr(data,"_data")) :
+        # Process data...........................
+        _data = numeric.array(data, dtype=dtype, copy=copy, subok=True)
+        _baseclass = getattr(_data, '_baseclass', type(_data))
+#        _data = numeric.ndarray.__new__(cls, shape=_data.shape,
+#                                        dtype=_data.dtype, buffer=_data)
+        _data = _data.view(cls)
+        # Process mask ..........................
+        if hasattr(data,"_mask"):
             if data is masked:
                 return masked.view(cls)
+            # Keep the original mask if needed ..
             if keep_mask:
                 if mask is nomask:
                     if copy:
-                        cls._defaultmask = data._mask.copy()
+                        mask = data._mask.copy()
                     else:
-                        cls._defaultmask = data._mask
+                        mask = data._mask
                 else:
-                    cls._defaultmask = mask_or(data._mask, mask, 
-                                                copy=copy, small_mask=small_mask)
+                    mask = mask_or(data._mask, mask,
+                                   copy=copy, small_mask=small_mask)
+            # Create a new mask from scratch ....
             else:
-                cls._defaultmask = make_mask(mask, copy=copy, small_mask=small_mask)
-            # Update fille_value
-            if fill_value is None:
-                cls._fill_value = data._fill_value
-            else:
-                cls._fill_value = fill_value
-            cls._defaulthardmask = hard_mask
-            _data = data._data
-            if dtype is not None and _data.dtype != numeric.dtype(dtype):
-                return _data.astype(dtype).view(cls)
-            elif copy:
-                return _data.copy().view(cls)
-            else:
-                return _data.view(cls)
-        # 2. Argument is not MA .......
-        if isinstance(data, ndarray):
-            if dtype is not None and data.dtype != numeric.dtype(dtype):
-                _data = data.astype(dtype)
-            elif copy:
-                _data = data.copy()
-            else:
-                _data = data
+                mask = make_mask(mask, copy=copy, small_mask=small_mask)
         else:
-            _data = numeric.array(data, dtype=dtype, copy=copy)
-#            try:
-#                _data = numeric.array(data, dtype=dtype, copy=copy)
-#            except TypeError:
-#                _data = empty(len(data), dtype=dtype)
-#                for (k,v) in enumerate(data):
-#                    _data[k] = v
-#                if mask is nomask:
-#                    cls.__defaultmask = getmask(_data)
-#                    return _data.view(cls)
-        # Define mask .................
-        mask = make_mask(mask, copy=copy, small_mask=small_mask)
-        #....Check shapes compatibility
+            mask = make_mask(mask, copy=copy, small_mask=small_mask)
+        # Check shape compatibility between mask and data
         if mask is not nomask:
             (nd, nm) = (_data.size, mask.size)
             if (nm != nd):
@@ -917,160 +1053,146 @@
                     raise MAError, msg % (nd, nm)
             elif (mask.shape != _data.shape):
                 mask = mask.reshape(_data.shape)
-#                mask = _data.shape
-        #....
-        cls._fill_value = fill_value
-        cls._defaulthardmask = hard_mask
-        cls._defaultmask = mask
-        cls._defaultoptions = options
-        return numeric.asanyarray(_data).view(cls)
+
+        # Update fille_value
+        _data._fill_value = getattr(data, '_fill_value', fill_value)
+        if _data._fill_value is None:
+            _data._fill_value = default_fill_value(_data)
+        _data._hardmask = hard_mask
+        _data._smallmask = small_mask
+        _data._baseclass = _baseclass
+        _data._mask = mask
+        return _data
     #..................................
     def __array_wrap__(self, obj, context=None):
         """Special hook for ufuncs.
 Wraps the numpy array and sets the mask according to context.
         """
-#        mclass = self.__class__
+        result = obj.view(type(self))
         #..........
         if context is None:
-#            return mclass(obj, mask=self._mask, copy=False)
-            return MaskedArray(obj, mask=self._mask, copy=False,
-                               dtype=obj.dtype,
-                               fill_value=self.fill_value, )
+            m = self._mask
         #..........
-        (func, args) = context[:2]
-        m = reduce(mask_or, [getmask(arg) for arg in args])
-        # Get domain mask
-        domain = ufunc_domain.get(func, None)
-        if domain is not None:
-            m = mask_or(m, domain(*[getattr(arg, '_data', arg) for arg in args]))
-        # Update mask
-        if m is not nomask:
-            try:
-                dshape = obj.shape
-            except AttributeError:
-                pass
-            else:
-                if m.shape != dshape:
+        else:
+            (func, args, out_index) = 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)
+                m = mask_or(m, d)
+            # Update mask
+            if m is not nomask and hasattr(obj,'shape'):
+                if m.shape != obj.shape:
                     m = reduce(mask_or, [getmaskarray(arg) for arg in args])
-#        return mclass(obj, copy=False, mask=m)
-        return MaskedArray(obj, copy=False, mask=m,)
-#                           dtype=obj.dtype, fill_value=self._fill_value)
+        #....
+#        result = obj.view(type(self))
+        result._mask = m
+        result._fill_value = self._fill_value
+        result._hardmask = self._hardmask
+        result._smallmask = self._smallmask
+        result._baseclass = self._baseclass
+        return result
     #........................
-    #TODO: there should be some reorganization to do round here.
     def __array_finalize__(self,obj):
         """Finalizes the masked array.
         """
-        #
-        if isinstance(obj, MaskedArray):
-            # We came here from a MaskedArray
-            self._data = obj._data
-            self._mask = obj._mask
-            self._hardmask = obj._hardmask
-            self._fill_value = obj._fill_value
-            self.options = obj.options
-        else:
-            # We came here from a .view()
-            if hasattr(obj,'_data') and hasattr(obj, '_mask'):
-                # obj is an old masked array or a smart record
-                self._data = obj._data
-                self._mask = obj._mask
-            else:
-                # obj is anything but...
-                self._data = obj
-                self._mask = self._defaultmask
-            # Set the instance default
-            self._hardmask = self._defaulthardmask
-            self.fill_value = self._fill_value
-            self.options = self._defaultoptions
-            # Reset the class default
-            MaskedArray._defaultmask = nomask
-            MaskedArray._defaulthardmask = False
-            MaskedArray._fill_value = None
-#        #
-        return    
+        # Finalize mask ...............
+        self._mask = getattr(obj, '_mask', self._defaultmask)
+        if self._mask is not nomask:
+            self._mask.shape = self.shape
+#        if mask is not nomask and mask.shape != self.shape:
+#            # try and set the shape if we're coming from a assignment to shape
+#            # __array_interface__['data'] is REALLY faster than ctypes.data
+#            if self.__array_interface__['data'] == obj.__array_interface__['data']:
+#                try:
+#                    mask.shape = self.shape
+#                except:
+#                    mask = nomask
+#            else:
+#                mask = nomask
+#        self._mask = mask
+        # Get the remaining options ...
+        self._hardmask = getattr(obj, '_hardmask', self._defaulthardmask)
+        self._smallmask = getattr(obj, '_smallmask', True)
+        self._baseclass = getattr(obj, '_baseclass', type(obj))
+        self._fill_value = getattr(obj, '_fill_value', None)
+        return
     #............................................
     def __getitem__(self, indx):
         """x.__getitem__(y) <==> x[y]
 Returns the item described by i. Not a copy as in previous versions.
         """
-        if getmask(indx) is not nomask:
-            msg = "Masked arrays must be filled before they can be used as indices!"
-            raise IndexError, msg
-        dout = self._data[indx]
+        # 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
-        scalardout = (len(numeric.shape(dout))==0)
-        # 
-        if m is nomask:
-            if scalardout:
-                return dout
-            else:
-                return self.__class__(dout, mask=nomask, keep_mask=True,
-                                      fill_value=self._fill_value,
-                                      **self.options)
-        #....
-        mi = m[indx]
-        if mi.size == 1:
-            if mi:
-                return masked
-            return dout
-        else:
-            return self.__class__(dout, mask=mi, fill_value=self._fill_value,
-                                  **self.options)
+        if hasattr(dout, 'shape') and len(dout.shape) > 0:
+            # Not a scalar: make sure that dout is a MA
+            dout = dout.view(type(self))
+            if m is not nomask:
+                # use _set_mask to take care of the shape
+                dout._set_mask(m[indx])
+        elif m is not nomask and m[indx]:
+            dout = masked
+        return dout
     #........................
-    def __setitem__(self, index, value):
+    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(index) is not nomask:
-            msg = "Masked arrays must be filled before they can be used as indices!"
-            raise IndexError, msg
+#        if getmask(indx) is not nomask:
+#            msg = "Masked arrays must be filled before they can be used as indices!"
+#            raise IndexError, msg
         #....
-        (d, m) = (self._data, self._mask)
+        m = self._mask
         #....
         if value is masked:
             if m is nomask:
-                m = make_mask_none(d.shape)
+                m = make_mask_none(self.shape)
             else:
                 m = m.copy()
-            m[index] = True
+            m[indx] = True
             self._mask = m
             return
-        #....   
+        #....
+        dval = numeric.asarray(value).astype(self.dtype)
+        valmask = getmask(value)
         if m is nomask:
-            d[index] = filled(value)
-            valmask = getmask(value)
             if valmask is not nomask:
-                m = make_mask_none(d.shape)
-                m[index] = valmask
+                m = make_mask_none(self.shape)
+                m[indx] = valmask
         elif not self._hardmask:
-            d[index] = filled(value)
-            valmask = getmask(value)
             m = m.copy()
             if valmask is nomask:
-                m[index] = False
+                m[indx] = False
             else:
-                m[index] = valmask
-        elif hasattr(index, 'dtype') and (index.dtype==bool_):
-            index *= ~m
-            d[index] = filled(value)
+                m[indx] = valmask
+        elif hasattr(indx, 'dtype') and (indx.dtype==bool_):
+            indx = indx * ~m
 #        elif isinstance(index, int):
-        else:            
-            mindx = m[index]
-            value = masked_array(value, mask=mindx, keep_mask=True)
-            valdata = filled(value)
-            valmask = getmask(value)
-            if valmask is nomask:
-                d[index] = valdata
-            elif valmask.size > 1:
-                dindx = d[index]
-                numeric.putmask(dindx, ~valmask, valdata)
-                d[index] = dindx
-                numeric.putmask(mindx, valmask, True)
-                m[index] = mindx
-        #.....    
-        if not m.any():
+        else:
+            mindx = mask_or(m[indx], valmask, copy=True)
+            dindx = self._data[indx]
+            if dindx.size > 1:
+                dindx[~mindx] = dval
+            elif mindx is nomask:
+                dindx = dval
+            dval = dindx
+            m[indx] = mindx
+        # Set data ..........
+        #dval = filled(value).astype(self.dtype)
+        ndarray.__setitem__(self._data,indx,dval)
+        #.....
+        if m is nomask or not m.any():
             self._mask = nomask
         else:
             self._mask = m
@@ -1079,66 +1201,18 @@
         """x.__getslice__(i, j) <==> x[i:j]
 Returns the slice described by i, j.
 The use of negative indices is not supported."""
-        m = self._mask
-        dout = self._data[i:j]
-        if m is nomask:
-            return self.__class__(dout, fill_value=self._fill_value,
-                                  **self.options)
-        else:
-            return self.__class__(dout, mask=m[i:j], fill_value=self._fill_value,
-                                  **self.options)
+        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`. 
+Sets a slice i:j to `value`.
 If `value` is masked, masks those locations."""
-        if self is masked:
-            #TODO: Well, maybe we could/should
-            raise MAError, "Cannot alter the 'masked' object."
-        #....
-        (d, m) = (self._data, self._mask)
-        #....
-        if value is masked:
-            if m is nomask:
-                m = make_mask_none(d.shape)
-            m[i:j] = True
-            self._mask = m
-            return
-        #....   
-        if m is nomask:
-            valmask = getmask(value)
-            valdata = filled(value)
-            d[i:j] = valdata
-            if valmask is not nomask:
-                m = make_mask_none(d.shape)
-                m[i:j] = valmask
-        elif not self._hardmask:
-            valmask = getmask(value)
-            valdata = filled(value)
-            d[i:j] = valdata
-            if valmask is nomask:
-                m[i:j] = False
-            else:
-                m[i:j] = valmask
-        else:
-            mindx = m[i:j]
-            value = masked_array(value, mask=mindx, keep_mask=True)
-            valmask = value._mask
-            if valmask is nomask:
-                d[i:j][~mindx] = filled(value)
-            elif valmask.size > 1:
-                d[i:j][~mindx] = value[~valmask]
-                m[i:j][valmask] = True
-        #.....    
-        if not m.any():
-            self._mask = nomask
-        else:
-            self._mask = m
+        self.__setitem__(slice(i,j), value)
     #............................................
     # If we don't want to crash the performance, we better leave __getattribute__ alone...
 #    def __getattribute__(self, name):
 #        """x.__getattribute__('name') = x.name
-#Returns the chosen attribute. 
+#Returns the chosen attribute.
 #If the attribute cannot be directly accessed, checks the _data section.
 #        """
 #        try:
@@ -1152,7 +1226,7 @@
     #............................................
     def __str__(self):
         """x.__str__() <==> str(x)
-Calculates the string representation, using masked for fill if it is enabled. 
+Calculates the string representation, using masked for fill if it is enabled.
 Otherwise, fills with fill value.
         """
         if masked_print_option.enabled():
@@ -1170,14 +1244,14 @@
 #                res = numeric.empty(self._data.shape, object_)
 #                numeric.putmask(res,~m,self._data)
                 res = self._data.astype("|O8")
-                res[self._mask] = f
+                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. 
+Calculates the repr representation, using masked for fill if it is enabled.
 Otherwise fill with fill value.
         """
         with_mask = """\
@@ -1208,87 +1282,65 @@
             'fill': str(self.fill_value),
             }
     #............................................
-    def __abs__(self):
-        """x.__abs__() <==> abs(x)
-Returns a masked array of the current subclass, with the new `_data` 
-the absolute of the inital `_data`.
-        """
-        return self.__class__(self._data.__abs__(), mask=self._mask,
-                              fill_value = self._fill_value, **self.options)
-    #
-    def __neg__(self):
-        """x.__abs__() <==> neg(x)
-Returns a masked array of the current subclass, with the new `_data` 
-the negative of the inital `_data`."""
-        try:
-            return self.__class__(self._data.__neg__(), mask=self._mask,
-                                  fill_value = self._fill_value, **self.options)
-        except MAError:
-            return negative(self)
-    #
     def __iadd__(self, other):
         "Adds other to self in place."
-        f = convert_typecode(filled(other, 0), self._data.dtype.char) 
+        ndarray.__iadd__(self._data,other)
         m = getmask(other)
-        self._data += f
         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."
-        f = convert_typecode(filled(other, 0), self._data.dtype.char) 
+        ndarray.__isub__(self._data,other)
         m = getmask(other)
-        self._data -= f
         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."
-        f = convert_typecode(filled(other, 1), self._data.dtype.char) 
+        ndarray.__imul__(self._data,other)
         m = getmask(other)
-        self._data *= f
         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."
-        f = convert_typecode(filled(other, 0), self._data.dtype.char) 
-        mo = getmask(other)
-        result = divide(self, masked_array(f, mask=mo))
-        self._data = result._data
-        dm = result._mask
-        if dm is not self._mask:
-            self._mask = dm
-        return self    
+        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
+    #....
+    __add__ = _mathmethod('__add__')
+    __radd__ = _mathmethod('__add__')
+    __sub__ = _mathmethod('__sub__')
+    __rsub__ = _mathmethod('__rsub__')
+    __pow__ = _mathmethod('__pow__')
+    __mul__ = _mathmethod('__mul__', 1, 1)
+    __rmul__ = _mathmethod('__mul__', 1, 1)
+    __div__ = _mathmethod('__div__', 0, 1, domain_safe_divide())
+    __rdiv__ = _mathmethod('__rdiv__', 1, 0, domain_safe_divide())
+    __truediv__ = _mathmethod('__truediv__', 0, 1, domain_safe_divide())
+    __rtruediv__ = _mathmethod('__rtruediv__', 1, 0, domain_safe_divide())
+    __floordiv__ = _mathmethod('__floordiv__', 0, 1, domain_safe_divide())
+    __rfloordiv__ = _mathmethod('__rfloordiv__', 1, 0, domain_safe_divide())
+    __eq__ = _mathmethod('__eq__')
+    __ne__ = _mathmethod('__ne__')
+    __le__ = _mathmethod('__le__')
+    __lt__ = _mathmethod('__lt__')
+    __ge__ = _mathmethod('__ge__')
+    __gt__ = _mathmethod('__gt__')
 
-#    #
-#    def __eq__(self, other):
-#        return equal(self,other)
-#
-#    def __ne__(self, other):
-#        return not_equal(self,other)
-#
-#    def __lt__(self, other):
-#        return less(self,other)
-#
-#    def __le__(self, other):
-#        return less_equal(self,other)
-#
-#    def __gt__(self, other):
-#        return greater(self,other)
-#
-#    def __ge__(self, other):
-#        return greater_equal(self,other)    
-
     #............................................
     def __float__(self):
         "Converts self to float."
@@ -1296,33 +1348,13 @@
             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._data.item())
+        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._data.item())    
-    
-    @property
-    def dtype(self):
-        """returns the data type of `_data`."""
-        return self._data.dtype
-    
-    def astype (self, tc):
-        """Returns self as an array of given type.
-Subclassing is preserved."""
-        if tc == self._data.dtype:
-            return self
-        try:
-            return self.__class__(self, mask=self._mask, dtype=tc, copy=True,
-                                  **self.options)
-        except:
-#            d = self._data.astype(tc) 
-            return self.__class__(self._data.astype(tc), mask=self._mask, 
-                                  dtype=tc, **self.options)
-#        
-#        
+        return int(self.item())
     #............................................
     def harden_mask(self):
         "Forces the mask to hard"
@@ -1330,19 +1362,11 @@
     def soften_mask(self):
         "Forces the mask to soft"
         self._hardmask = False
-    #............................................        
-    #TODO: FIX THAT: THAT"S NOT A REAL FLATITER
+    #............................................
     def _get_flat(self):
         """Calculates the flat value.
         """
-        if self._mask is nomask:
-            return masked_array(self._data.ravel(), mask=nomask, copy=False,
-                                fill_value = self._fill_value,
-                                **self.options)
-        else:
-            return masked_array(self._data.ravel(), mask=self._mask.ravel(),
-                                copy=False, fill_value = self._fill_value,
-                                **self.options)
+        return flatiter(self)
     #
     def _set_flat (self, value):
         "x.flat = value"
@@ -1350,107 +1374,85 @@
         y[:] = value
     #
     flat = property(fget=_get_flat, fset=_set_flat, doc="Flat version")
-    #   
     #............................................
-    def _get_real(self):
-        "Returns the real part of a complex array."
-        return self.__class__(self._data.real, mask=self.mask,
-                              fill_value = self._fill_value, **self.options)
-    def _set_real (self, value):
-        "Sets the real part of a complex array to `value`."
-        y = self.real
-        y[...] = value
-        
-    real = property(fget=_get_real, fset=_set_real, doc="Get it real!")
-
-    def _get_imaginary(self):
-        "Returns the imaginary part of a complex array."
-        return self.__class__(self._data.imag, mask=self.mask,
-                              fill_value = self._fill_value, **self.options)
-
-    def _set_imaginary (self, value):
-        "Sets the imaginary part of a complex array to `value`."
-        y = self.imaginary
-        y[...] = value
-        
-    imag = property(fget=_get_imaginary, fset=_set_imaginary, 
-                    doc="Imaginary part.")
-    imaginary = imag
-    #............................................
     def _get_mask(self):
         """Returns the current mask."""
         return self._mask
     def _set_mask(self, mask):
         """Sets the mask to `mask`."""
-        mask = make_mask(mask, copy=False, small_mask=True)
+        mask = make_mask(mask, copy=False, small_mask=self._smallmask)
         if mask is not nomask:
-            if mask.size != self._data.size:
+            if mask.size != self.size:
                 raise ValueError, "Inconsistent shape between data and mask!"
-            if mask.shape != self._data.shape:
-                mask.shape = self._data.shape
-            self._mask = mask  
+            if mask.shape != self.shape:
+                mask.shape = self.shape
+            self._mask = mask
         else:
             self._mask = nomask
     mask = property(fget=_get_mask, fset=_set_mask, doc="Mask")
     #............................................
+    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_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`. 
+        """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._data)
+            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.
         """
-        d = self._data
         m = self._mask
         if m is nomask:
-            return d
+            return self._data
         #
         if fill_value is None:
-            value = self.fill_value
-        else:
-            value = fill_value
+            fill_value = self.fill_value
         #
         if self is masked_singleton:
-            result = numeric.asanyarray(value)
+            result = numeric.asanyarray(fill_value)
         else:
-            result = d.copy()
+            result = self._data.copy()
             try:
-                result[m] = value
+                result[m] = fill_value
             except (TypeError, AttributeError):
-                value = numeric.array(value, dtype=object)
-                d = d.astype(object)
-                result = fromnumeric.choose(m, (d, value))
+                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 d.shape:
+                if self._data.shape:
                     raise
                 elif m:
-                    result = numeric.array(value, dtype=d.dtype)
+                    result = numeric.array(fill_value, dtype=self.dtype)
                 else:
-                    result = d
+                    result = self._data
         return result
-    
+
     def compressed(self):
         "A 1-D array of all the non-masked data."
-        d = self._data.ravel()
+        d = self.ravel()
         if self._mask is nomask:
             return d
         else:
-            return d[~self._mask.ravel()]
+            return d[numeric.logical_not(d._mask)]
     #............................................
     def count(self, axis=None):
         """Counts the non-masked elements of the array along a given axis,
@@ -1458,7 +1460,7 @@
 If `axis` is None, counts all the non-masked elements, and returns either a
 scalar or the masked singleton."""
         m = self._mask
-        s = self._data.shape
+        s = self.shape
         ls = len(s)
         if m is nomask:
             if ls == 0:
@@ -1466,12 +1468,12 @@
             if ls == 1:
                 return s[0]
             if axis is None:
-                return self._data.size
+                return self.size
             else:
                 n = s[axis]
                 t = list(s)
                 del t[axis]
-                return numeric.ones(t) * n      
+                return numeric.ones(t) * n
         n1 = fromnumeric.size(m, axis)
         n2 = m.astype(int_).sum(axis)
         if axis is None:
@@ -1479,135 +1481,108 @@
         else:
             return masked_array(n1 - n2)
     #............................................
-    def _get_shape(self):
-        "Returns the current shape."
-        return self._data.shape
+#    def _get_shape(self):
+#        "Returns the current shape."
+#        return self._data.shape
+#    #
+#    def _set_shape (self, newshape):
+#        "Sets the array's shape."
+#        self._data.shape = newshape
+#        if self._mask is not nomask:
+#            #self._mask = self._mask.copy()
+#            self._mask.shape = newshape
+#    #
+#    shape = property(fget=_get_shape, fset=_set_shape,
+#                     doc="Shape of the array, as a tuple.")
     #
-    def _set_shape (self, newshape):
-        "Sets the array's shape."
-        self._data.shape = newshape
-        if self._mask is not nomask:
-            #self._mask = self._mask.copy()
-            self._mask.shape = newshape
-    #
-    shape = property(fget=_get_shape, fset=_set_shape, 
-                     doc="Shape of the array, as a tuple.")
-    #
-    def _get_size(self):
-        "Returns the current size."
-        return self._data.size
-    size = property(fget=_get_size,
-                    doc="Size (number of elements) of the array.")
-    #
-    def _get_ndim(self):
-        "Returns the number of dimensions."
-        return self._data.ndim
-    ndim = property(fget=_get_ndim,
-                    doc="Number of dimensions of the array.")
-    #
     def reshape (self, *s):
         """Reshapes the array to shape s.
-Returns a new masked array. 
+Returns a new masked array.
 If you want to modify the shape in place, please use `a.shape = s`"""
+        # TODO: Do we keep super, or reshape _data and take a view ?
+        result = super(MaskedArray, self).reshape(*s)
         if self._mask is not nomask:
-            return self.__class__(self._data.reshape(*s), 
-                                  mask=self._mask.reshape(*s),
-                                  fill_value=self.fill_value, **self.options)
-        else:
-            return self.__class__(self._data.reshape(*s),
-                                  fill_value=self.fill_value, **self.options)
+            result._mask = self._mask.reshape(*s)
+        return result
     #
-    def repeat(self, repeats, axis=None):
-        """Repeat elements of `a` `repeats` times along `axis`.
-`repeats` is a sequence of length `a.shape[axis]` telling how many times 
-each element should be repeated.
-The mask is repeated accordingly.
-        """
-        f = self.filled()
-        if isinstance(repeats, types.IntType):
-            if axis is None:
-                num = f.size
-            else:
-                num = f.shape[axis]
-            repeats = tuple([repeats]*num)
-    
-        m = self._mask
-        if m is not nomask:
-            m = fromnumeric.repeat(m, repeats, axis)
-        d = fromnumeric.repeat(f, repeats, axis)
-        return self.__class__(d, mask=m, fill_value=self.fill_value, 
-                              **self.options)
+    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.    
+        """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,)
+            self._data.resize(newshape, refcheck, order)
             if self.mask is not nomask:
-                self._mask.resize(newshape,)
+                self._mask.resize(newshape, refcheck, order)
         except ValueError:
-            msg = "Cannot resize an array that has been referenced or "+\
-                  "is referencing another array in this way.\n"+\
-                  "Use the resize function."
-            raise ValueError, msg
+            raise ValueError("Cannot resize an array that has been referenced "
+                             "or is referencing another array in this way.\n"
+                             "Use the resize function.")
         return None
     #
-    def flatten(self):
-        """Flattens the array in place.
-        """
-        flatsize = self.size
-        self._data.resize((flatsize,))
-        if self.mask is not nomask:
-            self._mask.resize((flatsize,))
-        return self
-
+    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 repeat,
+`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.
         """
-        #TODO: Check that
-        (d, m) = (self._data, self._mask)
-        ind = filled(indices)
-        v = filled(values)
-        d.put(ind, v, mode=mode)
-        if m is not nomask:
-            if getmask(values) is not nomask:
-                m.put(ind, values._mask, mode=mode)
+        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(ind, False, mode=mode)
-            self._mask = make_mask(m, copy=False, small_mask=True)
+                m.put(indices, values._mask, mode=mode)
+            m = make_mask(m, copy=False, small_mask=True)
+        self._mask = m
     #............................................
     def ids (self):
-        """Return the ids of the data and mask areas."""
-        return (id(self._data), id(self._mask))
+        """Return the address of the data and mask areas."""
+        return (self.ctypes.data, self._mask.ctypes.data)    
     #............................................
-    def all(self, axis=None):
+    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 = filled(self, True).all(axis)
-        m = self._mask.all(axis)
-        return self.__class__(d, mask=m, dtype=bool_,
-                              fill_value=self._fill_value, **self.options)
-    def any(self, axis=None):
+        d = self.filled(True).all(axis=axis, out=out).view(type(self))
+        if d.ndim > 0:
+            d._set_mask(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 = filled(self, False).any(axis)
-        m = self._mask.all(axis)
-        return self.__class__(d, mask=m, dtype=bool_,
-                              fill_value=self._fill_value, **self.options) 
+        d = self.filled(False).any(axis=axis, out=out).view(type(self))
+        if d.ndim > 0:
+            d._set_mask(self._mask.all(axis))
+        return d
+    
     def nonzero(self):
         """a.nonzero() returns a tuple of arrays
 
@@ -1627,97 +1602,74 @@
         """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`?
-        (d,m) = (self._data, self._mask)
+        # TODO: What are we doing with `out`?
+        m = self._mask
         if m is nomask:
-            return d.trace(offset=offset, axis1=axis1, axis2=axis2, 
-                           out=out).astype(dtype)
+            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, 
-                              ).astype(dtype)
-            return D.sum(axis=None)
+            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) 
+        """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:
-#            if axis is None:
-#                return self._data.sum(None, dtype=dtype)
-            return self.__class__(self._data.sum(axis, dtype=dtype),
-                                  mask=nomask, fill_value=self.fill_value,
-                                  **self.options)
+            mask = nomask
         else:
-#            if axis is None:
-#                return self.filled(0).sum(None, dtype=dtype)
-            return self.__class__(self.filled(0).sum(axis, dtype=dtype),
-                                  mask=self._mask.all(axis),
-                                  fill_value=self.fill_value, **self.options)
-            
+            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._set_mask(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`. 
+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.
         """
-        if self._mask is nomask:
-#            if axis is None:
-#                return self._data.cumsum(None, dtype=dtype)
-            return self.__class__(self._data.cumsum(axis=axis, dtype=dtype),
-                                  fill_value=self.fill_value, **self.options)
-        else:
-#            if axis is None:
-#                return self.filled(0).cumsum(None, dtype=dtype)
-            return self.__class__(self.filled(0).cumsum(axis=axis, dtype=dtype),
-                                  mask=self._mask, fill_value=self.fill_value,
-                                  **self.options)
-        
+        result = self.filled(0).cumsum(axis=axis, dtype=dtype).view(type(self))
+        result._set_mask(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`. 
+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:
-#            if axis is None:
-#                return self._data.prod(None, dtype=dtype)
-            return self.__class__(self._data.prod(axis, dtype=dtype),
-                                  mask=nomask, fill_value=self.fill_value,
-                                  **self.options)
-#            return self.__class__(self._data.prod(axis=axis, dtype=dtype))
+            mask = nomask
         else:
-#            if axis is None:
-#                return self.filled(1).prod(None, dtype=dtype)
-            return self.__class__(self.filled(1).prod(axis=axis, dtype=dtype),
-                                  mask=self._mask.all(axis),
-                                  fill_value=self.fill_value,
-                                  **self.options)
+            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._set_mask(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`. 
+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.
         """
-        if self._mask is nomask:
-#            if axis is None:
-#                return self._data.cumprod(None, dtype=dtype)
-            return self.__class__(self._data.cumprod(axis=axis, dtype=dtype),
-                                  mask=nomask, fill_value=self.fill_value,
-                                  **self.options)
-        else:
-#            if axis is None:
-#                return self.filled(1).cumprod(None, dtype=dtype)
-            return self.__class__(self.filled(1).cumprod(axis=axis, dtype=dtype),
-                                  mask=self._mask,
-                                  fill_value=self.fill_value, **self.options)        
-            
+        result = self.filled(1).cumprod(axis=axis, dtype=dtype).view(type(self))
+        result._set_mask(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
 
@@ -1729,30 +1681,22 @@
     Returns a masked array, of the same class as a.
         """
         if self._mask is nomask:
-#            if axis is None:
-#                return self._data.mean(axis=None, dtype=dtype)
-            return self.__class__(self._data.mean(axis=axis, dtype=dtype),
-                                  mask=nomask, fill_value=self.fill_value,
-                                  **self.options)
+            return super(MaskedArray, self).mean(axis=axis, dtype=dtype)
         else:
-            dsum = fromnumeric.sum(self.filled(0), axis=axis, dtype=dtype)
+            dsum = self.sum(axis=axis, dtype=dtype)
             cnt = self.count(axis=axis)
-            mask = self._mask.all(axis)
-            if axis is None and mask:
-                return masked
-            return self.__class__(dsum*1./cnt, mask=mask,
-                                  fill_value=self.fill_value, **self.options)
-    
+            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.
@@ -1761,23 +1705,18 @@
 i.e. var = mean((x - x.mean())**2).
         """
         if self._mask is nomask:
-#            if axis is None:
-#                return self._data.var(axis=None, dtype=dtype)
-            return self.__class__(self._data.var(axis=axis, dtype=dtype),
-                                  mask=nomask,
-                                  fill_value=self.fill_value, **self.options)
+            # 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 = danom.sum(axis) / cnt
 #            dvar /= cnt
-            if axis is None:
-                return dvar
-            return self.__class__(dvar, 
-                                  mask=mask_or(self._mask.all(axis), (cnt==1)),
-                                  fill_value=self.fill_value, **self.options)
-            
+            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.
@@ -1786,20 +1725,14 @@
 deviations from the mean, i.e. std = sqrt(mean((x - x.mean())**2)).
         """
         dvar = self.var(axis,dtype)
-        if axis is None:
-            if dvar is masked:
-                return masked
-            else:
-                # Should we use umath.sqrt instead ?
-                return sqrt(dvar)
-        return self.__class__(sqrt(dvar._data), mask=dvar._mask,
-                              dtype = self.dtype,
-                              fill_value=self.fill_value, **self.options)
+        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`.        
+    Masked values are filled beforehand to `fill_value`.
     If `fill_value` is None, uses the default for the data type.
     Returns a numpy array.
 
@@ -1834,16 +1767,14 @@
     less space than sorts along other axis.
         """
         if fill_value is None:
-            fill_value = default_fill_value(self._data)
+            fill_value = default_fill_value(self)
         d = self.filled(fill_value)
-        if axis is None:
-            return d.argsort(kind=kind, order=order)
-        return d.argsort(axis, kind=kind, order=order)
-
+        return d.argsort(axis=axis, kind=kind, order=order)
+    #........................
     def argmin(self, axis=None, fill_value=None):
-        """Returns the array 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`.        
+    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.
 
@@ -1854,16 +1785,14 @@
         Default filling value. If None, uses the data type default.
         """
         if fill_value is None:
-            fill_value = default_fill_value(self._data)
+            fill_value = minimum_fill_value(self)
         d = self.filled(fill_value)
-        if axis is None:
-            return d.argmin()
         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 
+        """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`.        
+    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.
 
@@ -1874,17 +1803,12 @@
         Default filling value. If None, uses the data type default.
         """
         if fill_value is None:
-            fill_value = default_fill_value(self._data)
-            try:
-                fill_value = - fill_value
-            except:
-                pass
+            fill_value = maximum_fill_value(self._data)
         d = self.filled(fill_value)
-        if axis is None:
-            return d.argmax()
         return d.argmax(axis)
-    
-    def sort(self, axis=-1, kind='quicksort', order=None, endwith=True):
+
+    def sort(self, axis=-1, kind='quicksort', order=None, 
+             endwith=True, fill_value=None):
         """
         Sort a along the given axis.
 
@@ -1925,26 +1849,123 @@
     and use less space than sorts along other axis.
 
     """
-        
-        if endwith:
-            filler = minimum_fill_value(self.dtype)
+        if fill_value is None:
+            if endwith:
+                filler = minimum_fill_value(self)
+            else:
+                filler = maximum_fill_value(self)
         else:
-            filler = maximum_fill_value(self.dtype)
+            filler = fill_value
         indx = self.filled(filler).argsort(axis=axis,kind=kind,order=order)
-        self._data = self._data[indx]
-        m = self._mask
-        if m is not nomask:
-            self._mask = m[indx]
-        return 
+        self[:] = self[indx]
+        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 teh 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 filled with `fill_value`"""
+        return self.filled(fill_value).tolist()
+    #........................
+    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. 
+        """Returns the `_data` part of the MaskedArray.
 You should really use `_data` instead..."""
         return self._data
     def raw_data(self):
-        """Returns the `_data` part of the MaskedArray. 
+        """Returns the `_data` part of the MaskedArray.
 You should really use `_data` instead..."""
         return self._data
 
@@ -1954,7 +1975,7 @@
 
 #class _arithmethods:
 #    """Defines a wrapper for arithmetic methods.
-#Instead of directly calling a ufunc, the corresponding method of  the `array._data` 
+#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):
@@ -1963,7 +1984,7 @@
 #    - `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. 
+#    - `domain` (Domain object *[None]*) : Domain of non-validity.
 #        """
 #        self.methodname = methodname
 #        self.fill_self = fill_self
@@ -1983,7 +2004,7 @@
 #            else:
 #                domain = self.domain(target, base)
 #            if domain.any():
-#                #If `other` is a subclass of ndarray, `filled` must have the 
+#                #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.
@@ -1992,155 +2013,25 @@
 #                    if target is other:
 #                        # We don't want to modify other: let's copy target, then
 #                        target = target.copy()
-#                    target[:] = numeric.where(fromnumeric.asarray(domain), 
+#                    target[:] = numeric.where(fromnumeric.asarray(domain),
 #                                              self.fill_other, target)
 #                else:
-#                    target = numeric.where(fromnumeric.asarray(domain), 
+#                    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)       
+#        m = mask_or(m_self, m_other)
+#        method = getattr(base, self.methodname)
 #        return instance.__class__(method(target, *args), mask=m)
 #    #
 #    def patch(self):
 #        """Applies the method `func` from class `method` to MaskedArray"""
 #        return types.MethodType(self,None,MaskedArray)
 #..............................................................................
-class _arithmethods(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(numpy, 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 = filled(instance,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):
-                    if target is other or target is base:
-                        # We don't want to modify other: let's copy target, then
-                        # Same if target us base, instead...
-                        target = target.copy()
-                    target[:] = numeric.where(fromnumeric.asarray(domain), 
-                                              self.fill_other, target)
-                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)       
-        return instance.__class__(method(target, *args), mask=m,
-                                  fill_value=instance.fill_value,
-                                  **instance.options)
-#......................................
-class _compamethods(object):
-    """Defines comparison methods (eq, ge, gt...).
-Instead of calling a ufunc, the method of the masked object is called.
-    """
-    def __init__ (self, methodname, fill_self=0, fill_other=0):
-        """
-: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.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(numpy, 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 = mask_or(instance._mask, getmask(other), small_mask=False)
-        base = instance.filled(self.fill_self)
-        target = filled(other, self.fill_other)
-        method = getattr(base, self.methodname)      
-        return instance.__class__(method(target, *args), mask=m,
-                                  **instance.options)  
-#..........................................................
-MaskedArray.__add__ = _arithmethods('__add__')
-MaskedArray.__radd__ = _arithmethods('__add__')
-MaskedArray.__sub__ = _arithmethods('__sub__')
-MaskedArray.__rsub__ = _arithmethods('__rsub__')
-MaskedArray.__pow__ = _arithmethods('__pow__')
-MaskedArray.__mul__ = _arithmethods('__mul__', 1, 1)
-MaskedArray.__rmul__ = _arithmethods('__mul__', 1, 1)
-MaskedArray.__div__ = _arithmethods('__div__', 0, 1, 
-                                    domain_safe_divide())
-MaskedArray.__rdiv__ = _arithmethods('__rdiv__', 1, 0, 
-                                     domain_safe_divide())
-MaskedArray.__truediv__ = _arithmethods('__truediv__', 0, 1, 
-                                        domain_safe_divide())
-MaskedArray.__rtruediv__ = _arithmethods('__rtruediv__', 1, 0, 
-                                         domain_safe_divide())
-MaskedArray.__floordiv__ = _arithmethods('__floordiv__', 0, 1, 
-                                         domain_safe_divide())
-MaskedArray.__rfloordiv__ = _arithmethods('__rfloordiv__', 1, 0, 
-                                          domain_safe_divide())
-MaskedArray.__eq__ = _compamethods('__eq__')
-MaskedArray.__ne__ = _compamethods('__ne__')
-MaskedArray.__le__ = _compamethods('__le__')
-MaskedArray.__lt__ = _compamethods('__lt__')
-MaskedArray.__ge__ = _compamethods('__ge__')
-MaskedArray.__gt__ = _compamethods('__gt__')
+
 #####--------------------------------------------------------------------------
 #---- --- Shortcuts ---
 #####---------------------------------------------------------------------------
-def isMaskedArray (x):
+def isMaskedArray(x):
     "Is x a masked array, that is, an instance of MaskedArray?"
     return isinstance(x, MaskedArray)
 isarray = isMaskedArray
@@ -2150,15 +2041,16 @@
 masked = masked_singleton
 
 masked_array = MaskedArray
-def array(data, dtype=None, copy=False, order=False, mask=nomask, 
+def array(data, dtype=None, copy=False, order=False, mask=nomask,
           keep_mask=True, small_mask=True, hard_mask=None, fill_value=None):
-    """array(data, dtype=None, copy=True, order=False, mask=nomask, 
+    """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...
     """
-    return MaskedArray(data, mask=mask, dtype=dtype, copy=copy, 
-                       keep_mask=keep_mask, small_mask=small_mask, 
+    #TODO: we should try to put 'order' somwehere
+    return MaskedArray(data, mask=mask, dtype=dtype, copy=copy,
+                       keep_mask=keep_mask, small_mask=small_mask,
                        hard_mask=hard_mask, fill_value=fill_value)
 
 def is_masked(x):
@@ -2181,10 +2073,10 @@
 #
 #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. 
+#        Name of the function to apply on data.
 #    `onmask` : Boolean *[True]*
 #        Whether the mask must be processed also (True) or left alone (False).
 #    """
@@ -2202,13 +2094,13 @@
 #            return C(getattr(d,methodname).__call__(*args, **params),
 #                     mask=getattr(m,methodname)(*args, **params) )
 #        else:
-#            return C(getattr(d,methodname).__call__(*args, **params), mask=m) 
-#        
+#            return C(getattr(d,methodname).__call__(*args, **params), mask=m)
+#
 #    def patch(self):
 #        "Adds the new method to MaskedArray."
 #        return types.MethodType(self, None, MaskedArray)
 ##......................................
-#MaskedArray.conj = MaskedArray.conjugate = _arraymethod('conjugate').patch()  
+#MaskedArray.conj = MaskedArray.conjugate = _arraymethod('conjugate').patch()
 #MaskedArray.diagonal = _arraymethod('diagonal').patch()
 #MaskedArray.take = _arraymethod('take').patch()
 #MaskedArray.ravel = _arraymethod('ravel').patch()
@@ -2220,112 +2112,39 @@
 #MaskedArray.resize = _arraymethod('resize').patch()
 #MaskedArray.copy = _arraymethod('copy').patch()
 
-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)."
-        try:
-            return getattr(MaskedArray, self._name).__doc__
-        except:
-            return getattr(numpy, self._name).__doc__
-    #
-    def __get__(self, obj, objtype=None):
-        self.obj = obj
-        return self
-    #
-    def __call__(self, *args, **params):
-        methodname = self._name
-        obj = self.obj
-        (d, m) = (obj._data, obj._mask)
-        (t, f) = (obj.dtype, obj._fill_value)
-        C = self.obj.__class__
-        if m is nomask:
-            return C(getattr(d,methodname).__call__(*args, **params),
-                     dtype=t, fill_value=f)
-        elif self._onmask:
-            return C(getattr(d,methodname).__call__(*args, **params),
-                     mask=getattr(m,methodname)(*args, **params),
-                     dtype=t, fill_value=f, **obj.options)
-        else:
-            return C(getattr(d,methodname).__call__(*args, **params), mask=m,
-                     dtype=t, fill_value=f, **obj.options) 
-#......................................
-MaskedArray.conj = MaskedArray.conjugate = _arraymethod('conjugate') 
-MaskedArray.copy = _arraymethod('copy') 
-MaskedArray.diagonal = _arraymethod('diagonal')
-MaskedArray.take = _arraymethod('take')
-MaskedArray.ravel = _arraymethod('ravel')
-MaskedArray.transpose = _arraymethod('transpose')
-MaskedArray.T = property(fget=lambda self:self.transpose())
-MaskedArray.swapaxes = _arraymethod('swapaxes')
-MaskedArray.clip = _arraymethod('clip', onmask=False)
-MaskedArray.compress = _arraymethod('compress')
-MaskedArray.copy = _arraymethod('copy')
-MaskedArray.squeeze = _arraymethod('squeeze')
 
-#####--------------------------------------------------------------------------
+#####---------------------------------------------------------------------------
 #---- --- Extrema functions ---
-#####--------------------------------------------------------------------------
-class _minimum_operation:
-    "Object to calculate minima"
-    def __init__ (self):
-        """minimum(a, b) or minimum(a)
-In one argument case, returns the scalar minimum.
-        """
-        pass
-    #.........
-    def __call__ (self, a, b=None):
-        "Execute the call behavior."
+#####---------------------------------------------------------------------------
+class _extrema_operation(object):
+    def __call__(self, a, b=None):
+        "Executes the call behavior."
         if b is None:
-            m = getmask(a)
-            if m is nomask:
-                d = amin(filled(a).ravel())
-                return d
-            ac = a.compressed()
-            if len(ac) == 0:
-                return masked
-            else:
-                return amin(ac)
-        else:
-            return where(less(a, b), a, b)
+            return self.reduce(a)
+        return where(self.compare(a, b), a, b)
     #.........
-    def reduce(self, target, axis=0):
-        """Reduces `target` along the given `axis`."""
+    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 = {}
+
         if m is nomask:
-            t = filled(target)
-            return masked_array (umath.minimum.reduce (t, axis))
+            t = self.ufunc.reduce(target, **kargs)
         else:
-            t = umath.minimum.reduce(filled(target, minimum_fill_value(target)), 
-                                     axis)
-            m = umath.logical_and.reduce(m, axis)
-#            return masked_array(t, mask=m, fill_value=get_fill_value(target))
-            try:
-                return target.__class__(t, mask=m, dtype=t.dtype,
-                                        fill_value=get_fill_value(target))
-            except AttributeError:
-                return masked_array(t, mask=m, dtype=t.dtype,
-                                    fill_value=get_fill_value(target))
+            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):
+    def outer (self, a, b):
         "Returns the function applied to the outer product of a and b."
         ma = getmask(a)
         mb = getmask(b)
@@ -2335,9 +2154,32 @@
             ma = getmaskarray(a)
             mb = getmaskarray(b)
             m = logical_or.outer(ma, mb)
-        d = umath.minimum.outer(filled(a), filled(b))
-        return masked_array(d, mask=m)
-
+        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."""
@@ -2347,60 +2189,7 @@
         return minimum(array)
     else:
         return minimum.reduce(array, axis)
-#................................................
-class _maximum_operation:
-    "Object to calculate maxima"
-    def __init__ (self):
-        """maximum(a, b) or maximum(a)
-           In one argument case returns the scalar maximum.
-        """
-        pass
-    #.........
-    def __call__ (self, a, b=None):
-        "Executes the call behavior."
-        if b is None:
-            m = getmask(a)
-            if m is nomask:
-                d = amax(filled(a).ravel())
-                return d
-            ac = a.compressed()
-            if len(ac) == 0:
-                return masked
-            else:
-                return amax(ac)
-        else:
-            return where(greater(a, b), a, b)
-    #.........
-    def reduce (self, target, axis=0):
-        """Reduces target along the given axis."""
-        m = getmask(target)
-        if m is nomask:
-            t = filled(target)
-            return masked_array(umath.maximum.reduce (t, axis))
-        else:
-            t = umath.maximum.reduce(filled(target, maximum_fill_value(target)), 
-                                     axis)
-            m = umath.logical_and.reduce(m, axis)
-            try:
-                return target.__class__(t, mask=m, dtype=t.dtype,
-                                        fill_value=get_fill_value(target))
-            except AttributeError:
-                return masked_array(t, mask=m, dtype=t.dtype,
-                                    fill_value=get_fill_value(target))
-    #.........
-    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)
-        d = umath.maximum.outer(filled(a), filled(b))
-        return masked_array(d, mask=m)
-
+#............................
 def max(obj, axis=None, out=None):
     """Returns the maxima along the given axis.
 If `axis` is None, applies to the flattened array."""
@@ -2410,18 +2199,15 @@
         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)
-#................................................
-MaskedArray.min = min
-MaskedArray.max = max
-MaskedArray.ptp = ptp
 
+
 #####---------------------------------------------------------------------------
 #---- --- Definition of functions from the corresponding methods ---
 #####---------------------------------------------------------------------------
@@ -2438,19 +2224,20 @@
             return getattr(MaskedArray, self._methodname).__doc__
         except:
             return getattr(numpy, self._methodname).__doc__
-    def __call__(self, x, *args, **params):
-        if isinstance(x, MaskedArray):
-            return getattr(x, self._methodname).__call__(*args, **params)
-        #FIXME: As x is not a MaskedArray, we transform it to a ndarray with asarray
-        #FIXME: ... and call the corresponding method.
-        #FIXME: Except that sometimes it doesn't work (try reshape([1,2,3,4],(2,2)))
-        #FIXME: we end up with a "SystemError: NULL result without error in PyObject_Call"
-        #FIXME: A dirty trick is then to call the initial numpy function...
-        method = getattr(fromnumeric.asarray(x), self._methodname)
+    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__(x, *args, **params)
+            return getattr(numpy,self._methodname).__call__(a, *args, **params)
 
 all = _frommethod('all')
 anomalies = anom = _frommethod('anom')
@@ -2494,9 +2281,9 @@
         return masked_array(umath.power(fa, fb), m)
 
 #..............................................................................
-def argsort(a, axis=None, kind='quicksort', fill_value=None):
+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`.        
+    Masked values are filled beforehand to `fill_value`.
     If `fill_value` is None, uses the default for the data type.
     Returns a numpy array.
 
@@ -2534,13 +2321,13 @@
         fill_value = default_fill_value(a)
     d = filled(a, fill_value)
     if axis is None:
-        return d.argsort(kind=kind)
-    return d.argsort(axis, kind)
+        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 
+    """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`.        
+    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.
 
@@ -2553,14 +2340,12 @@
     if fill_value is None:
         fill_value = default_fill_value(a)
     d = filled(a, fill_value)
-    if axis is None:
-        return d.argmin(axis=None)
     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 
+    """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`.        
+    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.
 
@@ -2577,11 +2362,9 @@
         except:
             pass
     d = filled(a, fill_value)
-    if axis is None:
-        return d.argmax(axis=None)
     return d.argmax(axis=axis)
 
-def sort(a, axis=-1, kind='quicksort', order=None, endwith=True):
+def sort(a, axis=-1, kind='quicksort', order=None, endwith=True, fill_value=None):
     """
     Sort a along the given axis.
 
@@ -2623,14 +2406,17 @@
 
 """
     a = numeric.asanyarray(a)
-    if endwith:
-        filler = minimum_fill_value(a)
+    if fill_value is None:
+        if endwith:
+            filler = minimum_fill_value(a)
+        else:
+            filler = maximum_fill_value(a)
     else:
-        filler = maximum_fill_value(a)
+        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]   
+    return a[indx]
 
 def compressed(x):
     """Returns a compressed version of a masked array (or just the array if it
@@ -2647,36 +2433,29 @@
 
 def concatenate(arrays, axis=0):
     "Concatenates the arrays along the given axis"
-    #TODO: We lose the subclass, here! We should keep track of the classes...
-    #TODO: ...and find the max ? the lowest according to MRO?
-    d = []
+    d = numeric.concatenate([filled(a) for a in arrays], axis)
+    rcls = get_masked_subclass(*arrays)
+    data = d.view(rcls)
     for x in arrays:
-        d.append(filled(x))
-    d = numeric.concatenate(d, axis)
-    for x in arrays:
-        if getmask(x) is not nomask: 
+        if getmask(x) is not nomask:
             break
     else:
-        return masked_array(d)
-    dm = []
-    for x in arrays:
-        dm.append(getmaskarray(x))
-    dm = make_mask(numeric.concatenate(dm, axis), copy=False, small_mask=True)
-    return masked_array(d, mask=dm)
+        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):
-        (d,m) = (x._data, x._mask)
-        if m is nomask:
-            return masked_array(n_expand_dims(d,axis),                                
-                                dtype=d.dtype, fill_value=x._fill_value)
-        else:
-            return masked_array(n_expand_dims(d,axis), 
-                                mask=n_expand_dims(m,axis),
-                                dtype=d.dtype, fill_value=x._fill_value)
-    else:
-        return n_expand_dims(x,axis)
+        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):
@@ -2719,7 +2498,7 @@
         return fromnumeric.asarray(a).putmask(values, mask)
 
 def transpose(a,axes=None):
-    """Returns a view of the array with dimensions permuted according to axes.  
+    """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
@@ -2727,7 +2506,7 @@
         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.
@@ -2738,7 +2517,7 @@
 
 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 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.
     """
@@ -2746,11 +2525,9 @@
     m = getmask(x)
     if m is not nomask:
         m = fromnumeric.resize(m, new_shape)
-    if isinstance(x, MaskedArray):
-        result = x.__class__(fromnumeric.resize(filled(x), new_shape), mask=m)
-    else:
-        result = masked_array(fromnumeric.resize(filled(x), new_shape), mask=m)
-    result.set_fill_value(get_fill_value(x))
+    result = fromnumeric.resize(x, new_shape).view(get_masked_subclass(x))
+    if result.ndim:
+        result._mask = m
     return result
 
 
@@ -2766,10 +2543,10 @@
     return fromnumeric.shape(filled(obj))
 #
 def size(obj, axis=None):
-    """Returns the number of elements in the array along the given axis, 
+    """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(filled(obj), axis)
 #................................................
 
 #####--------------------------------------------------------------------------
@@ -2798,15 +2575,15 @@
     #TODO: implement options `out` and `mode`, if possible.
     def fmask (x):
         "Returns the filled array, or True if ``masked``."
-        if x is 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: 
+        if x is masked:
             return 1
         m = getmask(x)
-        if m is nomask: 
+        if m is nomask:
             return 0
         return m
     c = filled(indices, 0)
@@ -2832,22 +2609,20 @@
     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."""
-    if not hasattr(a, "_mask"):
-        mask = nomask
+    result = fromnumeric.round_(filled(a), decimals, out)
+    if isinstance(a,MaskedArray):
+        result = result.view(type(a))
+        result._mask = a._mask
     else:
-        mask = a._mask
-    if out is None:
-        return a.__class__(fromnumeric.round_(a, decimals, None), mask=mask)
-    else:
-        out = a.__class__(fromnumeric.round_(a, decimals, out), mask=mask)
-        return out
+        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)
-    
+    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
@@ -2856,9 +2631,9 @@
     """
     fa = filled(a, 0)
     fb = filled(b, 0)
-    if len(fa.shape) == 0: 
+    if len(fa.shape) == 0:
         fa.shape = (1,)
-    if len(fb.shape) == 0: 
+    if len(fb.shape) == 0:
         fb.shape = (1,)
     return masked_array(numeric.inner(fa, fb))
 innerproduct = inner
@@ -2876,7 +2651,7 @@
     mb = getmaskarray(b)
     m = make_mask(1-numeric.outer(1-ma, 1-mb), copy=0)
     return masked_array(d, mask=m)
-outerproduct = outer    
+outerproduct = outer
 
 def allequal (a, b, fill_value=True):
     """
@@ -2897,13 +2672,13 @@
         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` 
+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))
@@ -2927,24 +2702,24 @@
     """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 masked_array(numeric.empty(new_shape, dtype), mask=nomask)
+    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.  
+Note that this does NOT initialize the returned array.
 If you require your array to be initialized, you should use zeros_like()."""
-    return masked_array(numeric.empty_like(a), mask=nomask)
+    return numeric.empty_like(a).view(MaskedArray)
 
 def ones(new_shape, dtype=float):
-    """ones(shape, dtype=None) 
+    """ones(shape, dtype=None)
 Returns an array of the given dimensions, initialized to all ones."""
-    return masked_array(numeric.ones(new_shape, dtype), mask=nomask)
+    return numeric.ones(new_shape, dtype).view(MaskedArray)
 
 def zeros(new_shape, dtype=float):
-    """zeros(new_shape, dtype=None) 
+    """zeros(new_shape, dtype=None)
 Returns an array of the given dimensions, initialized to all zeros."""
-    return masked_array(numeric.zeros(new_shape, dtype), mask=nomask)
+    return numeric.zeros(new_shape, dtype).view(MaskedArray)
 
 #####--------------------------------------------------------------------------
 #---- --- Pickling ---
@@ -2960,13 +2735,13 @@
 def _getstate(a):
     "Returns the internal state of the masked array, for pickling purposes."
     state = (1,
-             a.shape, 
+             a.shape,
              a.dtype,
              a.flags.fnc,
-             (a._data).__reduce__()[-1][-1],
-             getmaskarray(a).__reduce__()[-1][-1])
+             a.tostring(),
+             getmaskarray(a).tostring())
     return state
-    
+
 def _setstate(a, 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:
@@ -2978,15 +2753,15 @@
     - a binary string for the mask.
         """
     (ver, shp, typ, isf, raw, msk) = state
-    (a._data).__setstate__((shp, typ, isf, raw))
+    super(MaskedArray, a).__setstate__((shp, typ, isf, raw))
     (a._mask).__setstate__((shp, dtype('|b1'), isf, msk))
-        
+
 def _reduce(a):
     """Returns a 3-tuple for pickling a MaskedArray."""
     return (_mareconstruct,
             (a.__class__, (0,), 'b', ),
             a.__getstate__())
-    
+
 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.
@@ -2994,15 +2769,15 @@
     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)  
+    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'): 
+    if not hasattr(F, 'readline'):
         F = open(F,'r')
     return cPickle.load(F)
 
@@ -3019,14 +2794,6 @@
 ################################################################################
 
 if __name__ == '__main__':
-    if 1:
-        a = masked_array([0,0])
-        b = a/a
-        assert (b._mask == [1,1]).all()
-        assert (a._data == [0,0]).all()
-    if 1:
-        a = arange(4)
-        a[1:-1] = masked
-        b = a[:-5]
-    if 1:
-        assert(masked_array(masked) is masked)
+    import numpy as N
+    from maskedarray.testutils import assert_equal, assert_array_equal
+    pi = N.pi

Added: trunk/Lib/sandbox/maskedarray/core_ini.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/core_ini.py	2007-02-11 00:33:13 UTC (rev 2699)
+++ trunk/Lib/sandbox/maskedarray/core_ini.py	2007-02-11 23:24:17 UTC (rev 2700)
@@ -0,0 +1,3033 @@
+"""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
+
+:author: Pierre Gerard-Marchant
+:contact: pierregm_at_uga_dot_edu
+:version: $Id$
+"""
+__author__ = "Pierre GF Gerard-Marchant ($Author$)"
+__version__ = '1.0'
+__revision__ = "$Revision$"
+__date__     = '$Date$'
+
+__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 numpy
+from numpy import bool_, complex_, float_, int_, object_, str_
+
+import numpy.core.umath as umath 
+import numpy.core.fromnumeric  as fromnumeric
+from numpy.core.numeric import ndarray
+from numpy.core.fromnumeric import amax, amin
+import numpy.core.numerictypes as ntypes
+from numpy.core.numerictypes import bool_, typecodes
+from numpy.core.multiarray import dtype
+import numpy.core.numeric as numeric
+from numpy.lib.shape_base import expand_dims as n_expand_dims
+import warnings
+
+
+
+MaskType = bool_
+nomask = MaskType(0)
+
+divide_tolerance = 1.e-35
+numpy.seterr(all='ignore')
+
+#####--------------------------------------------------------------------------
+#---- --- Helper functions ---
+#####--------------------------------------------------------------------------
+def convert_typecode(f,dtchar):
+    """Converts the type of `f` to a type compatible with `dtchar`, for inline operations."""
+    ftype = f.dtype.char
+    if dtchar == ftype:
+        return f
+    elif dtchar in typecodes['Integer']:
+        if ftype in typecodes['Integer']:
+            f = f.astype(dtchar)
+        else:
+            raise TypeError, 'Incorrect type for in-place operation.'
+    elif dtchar in typecodes['Float']:
+        if ftype in typecodes['Integer']:
+            f = f.astype(dtchar)
+        elif ftype in typecodes['Float']:
+            f = f.astype(dtchar)
+        else:
+            raise TypeError, 'Incorrect type for in-place operation.'
+    elif dtchar in typecodes['Complex']:
+        if ftype in typecodes['Integer']:
+            f = f.astype(dtchar)
+        elif ftype in typecodes['Float']:
+            f = f.astype(dtchar)
+        elif ftype in typecodes['Complex']:
+            f = f.astype(dtchar)
+        else:
+            raise TypeError, 'Incorrect type for in-place operation.'
+    else:
+        raise TypeError, 'Incorrect type for in-place operation.'
+    return f
+
+#####--------------------------------------------------------------------------
+#---- --- 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' : '???',        
+                  }
+#{0: <type 'numpy.bool_'>,
+# 1: <type 'numpy.int8'>,
+# 2: <type 'numpy.uint8'>,
+# 3: <type 'numpy.int16'>,
+# 4: <type 'numpy.uint16'>,
+# 5: <type 'numpy.int32'>,
+# 6: <type 'numpy.uint32'>,
+# 7: <type 'numpy.int64'>,
+# 8: <type 'numpy.uint64'>,
+# 9: <type 'numpy.int64'>,
+# 10: <type 'numpy.uint64'>,
+# 11: <type 'numpy.float32'>,
+# 12: <type 'numpy.float64'>,
+# 13: <type 'numpy.float128'>,
+# 14: <type 'numpy.complex64'>,
+# 15: <type 'numpy.complex128'>,
+# 16: <type 'numpy.complex256'>,
+# 17: <type 'numpy.object_'>,
+# 18: <type 'numpy.string_'>,
+# 19: <type 'numpy.unicode_'>,
+# 20: <type 'numpy.void'>,
+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'):
+        return default_filler[obj.dtype.kind]
+    elif isinstance(obj, float):
+        return default_filler['f']
+    elif isinstance(obj, int) or isinstance(obj, long):
+        return default_filler['i']
+    elif isinstance(obj, str):
+        return default_filler['S']
+    elif isinstance(obj, complex):
+        return default_filler['c']
+    elif isinstance(obj, numeric.dtype):
+        return default_filler[obj.kind]
+    else:
+        return default_filler['O']
+
+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, types.DictType):
+        return numeric.array(a, 'O')
+    else:
+        return numeric.array(a)
+        
+#####--------------------------------------------------------------------------
+#---- --- 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)))
+        result = self.f(d1, *args, **kwargs)
+        #
+        if isinstance(result, MaskedArray):
+            return result.__class__(result, mask=m)
+        return masked_array(result, mask=m)
+    #
+    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))
+        d1 = filled(a, self.fillx)
+        d2 = filled(b, self.filly)
+        result = self.f(d1, d2, *args, **kwargs)
+#        if isinstance(result, ndarray) \
+#               and m.ndim != 0 \
+#               and m.shape != result.shape:
+#            m = mask_or(getmaskarray(a), getmaskarray(b))
+        if isinstance(result, MaskedArray):
+            return result.__class__(result, mask=m)
+        return masked_array(result, mask=m)
+    #
+    def reduce (self, target, axis=0, dtype=None):
+        """Reduces `target` along the given `axis`."""
+        if isinstance(target, MaskedArray):
+            tclass = target.__class__
+        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 tclass(self.f.reduce (t, axis))
+        else:
+            t = tclass(t, mask=m)
+            # XXX: "or t.dtype" below is a workaround for what appears
+            # XXX: to be a bug in reduce.
+            t = self.f.reduce(filled(t, self.filly), axis, dtype=dtype or t.dtype)
+            m = umath.logical_and.reduce(m, axis)
+            if isinstance(t, ndarray):
+                return tclass(t, mask=m, fill_value=get_fill_value(target))
+            elif m:
+                return masked
+            else:
+                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 = umath.logical_or.outer(ma, mb)
+        d = self.f.outer(filled(a, self.fillx), filled(b, self.filly))
+        if isinstance(d, MaskedArray):
+            return d.__class__(d, mask=m)
+        return masked_array(d, mask=m)
+
+    def accumulate (self, target, axis=0):
+        """Accumulates `target` along `axis` after filling with y fill value."""
+        if isinstance(target, MaskedArray):
+            tclass = target.__class__
+        else:
+            tclass = masked_array
+        t = filled(target, self.filly)
+        return tclass(self.f.accumulate(t, axis))
+    
+    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)
+        result =  self.f(d1, d2)
+        if isinstance(result, MaskedArray):
+            return result.__class__(result, mask=m) 
+        return masked_array(result, mask=m)
+
+    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, x, 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.
+    """
+    cm = filled(condition,1)
+    if isinstance(x,MaskedArray):
+        m = mask_or(x._mask, cm)
+        return x.__class__(x._data, mask=m, copy=copy)
+    else:
+        return MaskedArray(fromnumeric.asarray(x), copy=copy, mask=cm)
+
+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 ---
+#####--------------------------------------------------------------------------
+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
+    #TODO: There some reorganization to do round here
+    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,
+                **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
+        # 1. Argument is MA ...........
+        if isinstance(data, MaskedArray) or\
+           (hasattr(data,"_mask") and hasattr(data,"_data")) :
+            if data is masked:
+                return masked.view(cls)
+            if keep_mask:
+                if mask is nomask:
+                    if copy:
+                        cls._defaultmask = data._mask.copy()
+                    else:
+                        cls._defaultmask = data._mask
+                else:
+                    cls._defaultmask = mask_or(data._mask, mask, 
+                                                copy=copy, small_mask=small_mask)
+            else:
+                cls._defaultmask = make_mask(mask, copy=copy, small_mask=small_mask)
+            # Update fille_value
+            if fill_value is None:
+                cls._fill_value = data._fill_value
+            else:
+                cls._fill_value = fill_value
+            cls._defaulthardmask = hard_mask
+            _data = data._data
+            if dtype is not None and _data.dtype != numeric.dtype(dtype):
+                return _data.astype(dtype).view(cls)
+            elif copy:
+                return _data.copy().view(cls)
+            else:
+                return _data.view(cls)
+        # 2. Argument is not MA .......
+        if isinstance(data, ndarray):
+            if dtype is not None and data.dtype != numeric.dtype(dtype):
+                _data = data.astype(dtype)
+            elif copy:
+                _data = data.copy()
+            else:
+                _data = data
+        else:
+            _data = numeric.array(data, dtype=dtype, copy=copy)
+#            try:
+#                _data = numeric.array(data, dtype=dtype, copy=copy)
+#            except TypeError:
+#                _data = empty(len(data), dtype=dtype)
+#                for (k,v) in enumerate(data):
+#                    _data[k] = v
+#                if mask is nomask:
+#                    cls.__defaultmask = getmask(_data)
+#                    return _data.view(cls)
+        # Define mask .................
+        mask = make_mask(mask, copy=copy, small_mask=small_mask)
+        #....Check shapes compatibility
+        if mask is not nomask:
+            (nd, nm) = (_data.size, mask.size)
+            if (nm != nd):
+                # We need to resize w/ a function, in case _data is only a reference
+                if nm == 1:
+                    mask = fromnumeric.resize(mask, _data.shape)
+                elif nd == 1:
+                    _data = fromnumeric.resize(_data, mask.shape)
+                else:
+                    msg = "Mask and data not compatible: data size is %i, "+\
+                          "mask size is %i."
+                    raise MAError, msg % (nd, nm)
+            elif (mask.shape != _data.shape):
+                mask = mask.reshape(_data.shape)
+#                mask = _data.shape
+        #....
+        cls._fill_value = fill_value
+        cls._defaulthardmask = hard_mask
+        cls._defaultmask = mask
+        cls._defaultoptions = options
+        return numeric.asanyarray(_data).view(cls)
+    #..................................
+    def __array_wrap__(self, obj, context=None):
+        """Special hook for ufuncs.
+Wraps the numpy array and sets the mask according to context.
+        """
+#        mclass = self.__class__
+        #..........
+        if context is None:
+#            return mclass(obj, mask=self._mask, copy=False)
+            return MaskedArray(obj, mask=self._mask, copy=False,
+                               dtype=obj.dtype,
+                               fill_value=self.fill_value, )
+        #..........
+        (func, args) = context[:2]
+        m = reduce(mask_or, [getmask(arg) for arg in args])
+        # Get domain mask
+        domain = ufunc_domain.get(func, None)
+        if domain is not None:
+            m = mask_or(m, domain(*[getattr(arg, '_data', arg) for arg in args]))
+        # Update mask
+        if m is not nomask:
+            try:
+                dshape = obj.shape
+            except AttributeError:
+                pass
+            else:
+                if m.shape != dshape:
+                    m = reduce(mask_or, [getmaskarray(arg) for arg in args])
+#        return mclass(obj, copy=False, mask=m)
+        return MaskedArray(obj, copy=False, mask=m,)
+#                           dtype=obj.dtype, fill_value=self._fill_value)
+    #........................
+    #TODO: there should be some reorganization to do round here.
+    def __array_finalize__(self,obj):
+        """Finalizes the masked array.
+        """
+        #
+        if isinstance(obj, MaskedArray):
+            # We came here from a MaskedArray
+            self._data = obj._data
+            self._mask = obj._mask
+            self._hardmask = obj._hardmask
+            self._fill_value = obj._fill_value
+            self.options = obj.options
+        else:
+            # We came here from a .view()
+            if hasattr(obj,'_data') and hasattr(obj, '_mask'):
+                # obj is an old masked array or a smart record
+                self._data = obj._data
+                self._mask = obj._mask
+            else:
+                # obj is anything but...
+                self._data = obj
+                self._mask = self._defaultmask
+            # Set the instance default
+            self._hardmask = self._defaulthardmask
+            self.fill_value = self._fill_value
+            self.options = self._defaultoptions
+            # Reset the class default
+            MaskedArray._defaultmask = nomask
+            MaskedArray._defaulthardmask = False
+            MaskedArray._fill_value = None
+#        #
+        return    
+    #............................................
+    def __getitem__(self, indx):
+        """x.__getitem__(y) <==> x[y]
+Returns the item described by i. Not a copy as in previous versions.
+        """
+        if getmask(indx) is not nomask:
+            msg = "Masked arrays must be filled before they can be used as indices!"
+            raise IndexError, msg
+        dout = self._data[indx]
+        m = self._mask
+        scalardout = (len(numeric.shape(dout))==0)
+        # 
+        if m is nomask:
+            if scalardout:
+                return dout
+            else:
+                return self.__class__(dout, mask=nomask, keep_mask=True,
+                                      fill_value=self._fill_value,
+                                      **self.options)
+        #....
+        mi = m[indx]
+        if mi.size == 1:
+            if mi:
+                return masked
+            return dout
+        else:
+            return self.__class__(dout, mask=mi, fill_value=self._fill_value,
+                                  **self.options)
+    #........................
+    def __setitem__(self, index, 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(index) is not nomask:
+            msg = "Masked arrays must be filled before they can be used as indices!"
+            raise IndexError, msg
+        #....
+        (d, m) = (self._data, self._mask)
+        #....
+        if value is masked:
+            if m is nomask:
+                m = make_mask_none(d.shape)
+            else:
+                m = m.copy()
+            m[index] = True
+            self._mask = m
+            return
+        #....   
+        if m is nomask:
+            d[index] = filled(value)
+            valmask = getmask(value)
+            if valmask is not nomask:
+                m = make_mask_none(d.shape)
+                m[index] = valmask
+        elif not self._hardmask:
+            d[index] = filled(value)
+            valmask = getmask(value)
+            m = m.copy()
+            if valmask is nomask:
+                m[index] = False
+            else:
+                m[index] = valmask
+        elif hasattr(index, 'dtype') and (index.dtype==bool_):
+            index *= ~m
+            d[index] = filled(value)
+#        elif isinstance(index, int):
+        else:            
+            mindx = m[index]
+            value = masked_array(value, mask=mindx, keep_mask=True)
+            valdata = filled(value)
+            valmask = getmask(value)
+            if valmask is nomask:
+                d[index] = valdata
+            elif valmask.size > 1:
+                dindx = d[index]
+                numeric.putmask(dindx, ~valmask, valdata)
+                d[index] = dindx
+                numeric.putmask(mindx, valmask, True)
+                m[index] = mindx
+        #.....    
+        if not m.any():
+            self._mask = nomask
+        else:
+            self._mask = m
+    #............................................
+    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."""
+        m = self._mask
+        dout = self._data[i:j]
+        if m is nomask:
+            return self.__class__(dout, fill_value=self._fill_value,
+                                  **self.options)
+        else:
+            return self.__class__(dout, mask=m[i:j], fill_value=self._fill_value,
+                                  **self.options)
+    #........................
+    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."""
+        if self is masked:
+            #TODO: Well, maybe we could/should
+            raise MAError, "Cannot alter the 'masked' object."
+        #....
+        (d, m) = (self._data, self._mask)
+        #....
+        if value is masked:
+            if m is nomask:
+                m = make_mask_none(d.shape)
+            m[i:j] = True
+            self._mask = m
+            return
+        #....   
+        if m is nomask:
+            valmask = getmask(value)
+            valdata = filled(value)
+            d[i:j] = valdata
+            if valmask is not nomask:
+                m = make_mask_none(d.shape)
+                m[i:j] = valmask
+        elif not self._hardmask:
+            valmask = getmask(value)
+            valdata = filled(value)
+            d[i:j] = valdata
+            if valmask is nomask:
+                m[i:j] = False
+            else:
+                m[i:j] = valmask
+        else:
+            mindx = m[i:j]
+            value = masked_array(value, mask=mindx, keep_mask=True)
+            valmask = value._mask
+            if valmask is nomask:
+                d[i:j][~mindx] = filled(value)
+            elif valmask.size > 1:
+                d[i:j][~mindx] = value[~valmask]
+                m[i:j][valmask] = True
+        #.....    
+        if not m.any():
+            self._mask = nomask
+        else:
+            self._mask = m
+    #............................................
+    # If we don't want to crash the performance, we better leave __getattribute__ alone...
+#    def __getattribute__(self, name):
+#        """x.__getattribute__('name') = x.name
+#Returns the chosen attribute. 
+#If the attribute cannot be directly accessed, checks the _data section.
+#        """
+#        try:
+#            return ndarray.__getattribute__(self, name)
+#        except AttributeError:
+#            pass
+#        try:
+#            return self._data.__getattribute__(name)
+#        except AttributeError:
+#            raise AttributeError
+    #............................................
+    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 == () and m:
+                    return str(f)
+                # 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[self._mask] = 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 __abs__(self):
+        """x.__abs__() <==> abs(x)
+Returns a masked array of the current subclass, with the new `_data` 
+the absolute of the inital `_data`.
+        """
+        return self.__class__(self._data.__abs__(), mask=self._mask,
+                              fill_value = self._fill_value, **self.options)
+    #
+    def __neg__(self):
+        """x.__abs__() <==> neg(x)
+Returns a masked array of the current subclass, with the new `_data` 
+the negative of the inital `_data`."""
+        try:
+            return self.__class__(self._data.__neg__(), mask=self._mask,
+                                  fill_value = self._fill_value, **self.options)
+        except MAError:
+            return negative(self)
+    #
+    def __iadd__(self, other):
+        "Adds other to self in place."
+        f = convert_typecode(filled(other, 0), self._data.dtype.char) 
+        m = getmask(other)
+        self._data += f
+        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."
+        f = convert_typecode(filled(other, 0), self._data.dtype.char) 
+        m = getmask(other)
+        self._data -= f
+        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."
+        f = convert_typecode(filled(other, 1), self._data.dtype.char) 
+        m = getmask(other)
+        self._data *= f
+        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."
+        f = convert_typecode(filled(other, 0), self._data.dtype.char) 
+        mo = getmask(other)
+        result = divide(self, masked_array(f, mask=mo))
+        self._data = result._data
+        dm = result._mask
+        if dm is not self._mask:
+            self._mask = dm
+        return self    
+
+#    #
+#    def __eq__(self, other):
+#        return equal(self,other)
+#
+#    def __ne__(self, other):
+#        return not_equal(self,other)
+#
+#    def __lt__(self, other):
+#        return less(self,other)
+#
+#    def __le__(self, other):
+#        return less_equal(self,other)
+#
+#    def __gt__(self, other):
+#        return greater(self,other)
+#
+#    def __ge__(self, other):
+#        return greater_equal(self,other)    
+
+    #............................................
+    def __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._data.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._data.item())    
+    
+    @property
+    def dtype(self):
+        """returns the data type of `_data`."""
+        return self._data.dtype
+    
+    def astype (self, tc):
+        """Returns self as an array of given type.
+Subclassing is preserved."""
+        if tc == self._data.dtype:
+            return self
+        try:
+            return self.__class__(self, mask=self._mask, dtype=tc, copy=True,
+                                  **self.options)
+        except:
+#            d = self._data.astype(tc) 
+            return self.__class__(self._data.astype(tc), mask=self._mask, 
+                                  dtype=tc, **self.options)
+#        
+#        
+    #............................................
+    def harden_mask(self):
+        "Forces the mask to hard"
+        self._hardmask = True
+    def soften_mask(self):
+        "Forces the mask to soft"
+        self._hardmask = False
+    #............................................        
+    #TODO: FIX THAT: THAT"S NOT A REAL FLATITER
+    def _get_flat(self):
+        """Calculates the flat value.
+        """
+        if self._mask is nomask:
+            return masked_array(self._data.ravel(), mask=nomask, copy=False,
+                                fill_value = self._fill_value,
+                                **self.options)
+        else:
+            return masked_array(self._data.ravel(), mask=self._mask.ravel(),
+                                copy=False, fill_value = self._fill_value,
+                                **self.options)
+    #
+    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_real(self):
+        "Returns the real part of a complex array."
+        return self.__class__(self._data.real, mask=self.mask,
+                              fill_value = self._fill_value, **self.options)
+    def _set_real (self, value):
+        "Sets the real part of a complex array to `value`."
+        y = self.real
+        y[...] = value
+        
+    real = property(fget=_get_real, fset=_set_real, doc="Get it real!")
+
+    def _get_imaginary(self):
+        "Returns the imaginary part of a complex array."
+        return self.__class__(self._data.imag, mask=self.mask,
+                              fill_value = self._fill_value, **self.options)
+
+    def _set_imaginary (self, value):
+        "Sets the imaginary part of a complex array to `value`."
+        y = self.imaginary
+        y[...] = value
+        
+    imag = property(fget=_get_imaginary, fset=_set_imaginary, 
+                    doc="Imaginary part.")
+    imaginary = imag
+    #............................................
+    def _get_mask(self):
+        """Returns the current mask."""
+        return self._mask
+    def _set_mask(self, mask):
+        """Sets the mask to `mask`."""
+        mask = make_mask(mask, copy=False, small_mask=True)
+        if mask is not nomask:
+            if mask.size != self._data.size:
+                raise ValueError, "Inconsistent shape between data and mask!"
+            if mask.shape != self._data.shape:
+                mask.shape = self._data.shape
+            self._mask = mask  
+        else:
+            self._mask = nomask
+    mask = property(fget=_get_mask, fset=_set_mask, doc="Mask")
+    #............................................
+    def get_fill_value(self):
+        "Returns the filling value."
+        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._data)
+        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.
+        """
+        d = self._data
+        m = self._mask
+        if m is nomask:
+            return d
+        #
+        if fill_value is None:
+            value = self.fill_value
+        else:
+            value = fill_value
+        #
+        if self is masked_singleton:
+            result = numeric.asanyarray(value)
+        else:
+            result = d.copy()
+            try:
+                result[m] = value
+            except (TypeError, AttributeError):
+                value = numeric.array(value, dtype=object)
+                d = d.astype(object)
+                result = fromnumeric.choose(m, (d, value))
+            except IndexError:
+                #ok, if scalar
+                if d.shape:
+                    raise
+                elif m:
+                    result = numeric.array(value, dtype=d.dtype)
+                else:
+                    result = d
+        return result
+    
+    def compressed(self):
+        "A 1-D array of all the non-masked data."
+        d = self._data.ravel()
+        if self._mask is nomask:
+            return d
+        else:
+            return d[~self._mask.ravel()]
+    #............................................
+    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._data.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._data.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 _get_shape(self):
+        "Returns the current shape."
+        return self._data.shape
+    #
+    def _set_shape (self, newshape):
+        "Sets the array's shape."
+        self._data.shape = newshape
+        if self._mask is not nomask:
+            #self._mask = self._mask.copy()
+            self._mask.shape = newshape
+    #
+    shape = property(fget=_get_shape, fset=_set_shape, 
+                     doc="Shape of the array, as a tuple.")
+    #
+    def _get_size(self):
+        "Returns the current size."
+        return self._data.size
+    size = property(fget=_get_size,
+                    doc="Size (number of elements) of the array.")
+    #
+    def _get_ndim(self):
+        "Returns the number of dimensions."
+        return self._data.ndim
+    ndim = property(fget=_get_ndim,
+                    doc="Number of dimensions of the array.")
+    #
+    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`"""
+        if self._mask is not nomask:
+            return self.__class__(self._data.reshape(*s), 
+                                  mask=self._mask.reshape(*s),
+                                  fill_value=self.fill_value, **self.options)
+        else:
+            return self.__class__(self._data.reshape(*s),
+                                  fill_value=self.fill_value, **self.options)
+    #
+    def repeat(self, repeats, axis=None):
+        """Repeat elements of `a` `repeats` times along `axis`.
+`repeats` is a sequence of length `a.shape[axis]` telling how many times 
+each element should be repeated.
+The mask is repeated accordingly.
+        """
+        f = self.filled()
+        if isinstance(repeats, types.IntType):
+            if axis is None:
+                num = f.size
+            else:
+                num = f.shape[axis]
+            repeats = tuple([repeats]*num)
+    
+        m = self._mask
+        if m is not nomask:
+            m = fromnumeric.repeat(m, repeats, axis)
+        d = fromnumeric.repeat(f, repeats, axis)
+        return self.__class__(d, mask=m, fill_value=self.fill_value, 
+                              **self.options)
+    #
+    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,)
+            if self.mask is not nomask:
+                self._mask.resize(newshape,)
+        except ValueError:
+            msg = "Cannot resize an array that has been referenced or "+\
+                  "is referencing another array in this way.\n"+\
+                  "Use the resize function."
+            raise ValueError, msg
+        return None
+    #
+    def flatten(self):
+        """Flattens the array in place.
+        """
+        flatsize = self.size
+        self._data.resize((flatsize,))
+        if self.mask is not nomask:
+            self._mask.resize((flatsize,))
+        return self
+
+    #
+    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 repeat,
+if necessary.
+If `values` has some masked values, the initial mask is updated in consequence,
+else the corresponding values are unmasked.
+        """
+        #TODO: Check that
+        (d, m) = (self._data, self._mask)
+        ind = filled(indices)
+        v = filled(values)
+        d.put(ind, v, mode=mode)
+        if m is not nomask:
+            if getmask(values) is not nomask:
+                m.put(ind, values._mask, mode=mode)
+            else:
+                m.put(ind, False, mode=mode)
+            self._mask = make_mask(m, copy=False, small_mask=True)
+    #............................................
+    def ids (self):
+        """Return the ids of the data and mask areas."""
+        return (id(self._data), id(self._mask))
+    #............................................
+    def all(self, axis=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.
+        """
+        d = filled(self, True).all(axis)
+        m = self._mask.all(axis)
+        return self.__class__(d, mask=m, dtype=bool_,
+                              fill_value=self._fill_value, **self.options)
+    def any(self, axis=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.
+        """
+        d = filled(self, False).any(axis)
+        m = self._mask.all(axis)
+        return self.__class__(d, mask=m, dtype=bool_,
+                              fill_value=self._fill_value, **self.options) 
+    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`?
+        (d,m) = (self._data, self._mask)
+        if m is nomask:
+            return d.trace(offset=offset, axis1=axis1, axis2=axis2, 
+                           out=out).astype(dtype)
+        else:
+            D = self.diagonal(offset=offset, axis1=axis1, axis2=axis2, 
+                              ).astype(dtype)
+            return D.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:
+#            if axis is None:
+#                return self._data.sum(None, dtype=dtype)
+            return self.__class__(self._data.sum(axis, dtype=dtype),
+                                  mask=nomask, fill_value=self.fill_value,
+                                  **self.options)
+        else:
+#            if axis is None:
+#                return self.filled(0).sum(None, dtype=dtype)
+            return self.__class__(self.filled(0).sum(axis, dtype=dtype),
+                                  mask=self._mask.all(axis),
+                                  fill_value=self.fill_value, **self.options)
+            
+    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.
+        """
+        if self._mask is nomask:
+#            if axis is None:
+#                return self._data.cumsum(None, dtype=dtype)
+            return self.__class__(self._data.cumsum(axis=axis, dtype=dtype),
+                                  fill_value=self.fill_value, **self.options)
+        else:
+#            if axis is None:
+#                return self.filled(0).cumsum(None, dtype=dtype)
+            return self.__class__(self.filled(0).cumsum(axis=axis, dtype=dtype),
+                                  mask=self._mask, fill_value=self.fill_value,
+                                  **self.options)
+        
+    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:
+#            if axis is None:
+#                return self._data.prod(None, dtype=dtype)
+            return self.__class__(self._data.prod(axis, dtype=dtype),
+                                  mask=nomask, fill_value=self.fill_value,
+                                  **self.options)
+#            return self.__class__(self._data.prod(axis=axis, dtype=dtype))
+        else:
+#            if axis is None:
+#                return self.filled(1).prod(None, dtype=dtype)
+            return self.__class__(self.filled(1).prod(axis=axis, dtype=dtype),
+                                  mask=self._mask.all(axis),
+                                  fill_value=self.fill_value,
+                                  **self.options)
+    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.
+        """
+        if self._mask is nomask:
+#            if axis is None:
+#                return self._data.cumprod(None, dtype=dtype)
+            return self.__class__(self._data.cumprod(axis=axis, dtype=dtype),
+                                  mask=nomask, fill_value=self.fill_value,
+                                  **self.options)
+        else:
+#            if axis is None:
+#                return self.filled(1).cumprod(None, dtype=dtype)
+            return self.__class__(self.filled(1).cumprod(axis=axis, dtype=dtype),
+                                  mask=self._mask,
+                                  fill_value=self.fill_value, **self.options)        
+            
+    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:
+#            if axis is None:
+#                return self._data.mean(axis=None, dtype=dtype)
+            return self.__class__(self._data.mean(axis=axis, dtype=dtype),
+                                  mask=nomask, fill_value=self.fill_value,
+                                  **self.options)
+        else:
+            dsum = fromnumeric.sum(self.filled(0), axis=axis, dtype=dtype)
+            cnt = self.count(axis=axis)
+            mask = self._mask.all(axis)
+            if axis is None and mask:
+                return masked
+            return self.__class__(dsum*1./cnt, mask=mask,
+                                  fill_value=self.fill_value, **self.options)
+    
+    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:
+#            if axis is None:
+#                return self._data.var(axis=None, dtype=dtype)
+            return self.__class__(self._data.var(axis=axis, dtype=dtype),
+                                  mask=nomask,
+                                  fill_value=self.fill_value, **self.options)
+        else:
+            cnt = self.count(axis=axis)
+            danom = self.anom(axis=axis, dtype=dtype)
+            danom *= danom
+            dvar = danom.sum(axis) / cnt
+#            dvar /= cnt
+            if axis is None:
+                return dvar
+            return self.__class__(dvar, 
+                                  mask=mask_or(self._mask.all(axis), (cnt==1)),
+                                  fill_value=self.fill_value, **self.options)
+            
+    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 None:
+            if dvar is masked:
+                return masked
+            else:
+                # Should we use umath.sqrt instead ?
+                return sqrt(dvar)
+        return self.__class__(sqrt(dvar._data), mask=dvar._mask,
+                              dtype = self.dtype,
+                              fill_value=self.fill_value, **self.options)
+    #............................................
+    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._data)
+        d = self.filled(fill_value)
+        if axis is None:
+            return d.argsort(kind=kind, order=order)
+        return d.argsort(axis, kind=kind, order=order)
+
+    def argmin(self, 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(self._data)
+        d = self.filled(fill_value)
+        if axis is None:
+            return d.argmin()
+        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 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(self._data)
+            try:
+                fill_value = - fill_value
+            except:
+                pass
+        d = self.filled(fill_value)
+        if axis is None:
+            return d.argmax()
+        return d.argmax(axis)
+    
+    def sort(self, axis=-1, kind='quicksort', order=None, endwith=True):
+        """
+        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.
+
+    """
+        
+        if endwith:
+            filler = minimum_fill_value(self.dtype)
+        else:
+            filler = maximum_fill_value(self.dtype)
+        indx = self.filled(filler).argsort(axis=axis,kind=kind,order=order)
+        self._data = self._data[indx]
+        m = self._mask
+        if m is not nomask:
+            self._mask = m[indx]
+        return 
+    #............................................
+    # Backwards Compatibility. Heck...
+    @property
+    def data(self):
+        """Returns the `_data` part of the MaskedArray. 
+You should really use `_data` instead..."""
+        return self._data
+    def raw_data(self):
+        """Returns the `_data` part of the MaskedArray. 
+You should really use `_data` instead..."""
+        return self._data
+
+##..............................................................................
+
+
+
+#class _arithmethods:
+#    """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
+#    #
+#    def __call__ (self, instance, other, *args):
+#        "Execute the call behavior."
+#        m_self = instance._mask
+#        m_other = getmask(other)
+#        base = filled(instance,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):
+#                    if target is other:
+#                        # We don't want to modify other: let's copy target, then
+#                        target = target.copy()
+#                    target[:] = numeric.where(fromnumeric.asarray(domain), 
+#                                              self.fill_other, target)
+#                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)       
+#        return instance.__class__(method(target, *args), mask=m)
+#    #
+#    def patch(self):
+#        """Applies the method `func` from class `method` to MaskedArray"""
+#        return types.MethodType(self,None,MaskedArray)
+#..............................................................................
+class _arithmethods(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(numpy, 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 = filled(instance,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):
+                    if target is other or target is base:
+                        # We don't want to modify other: let's copy target, then
+                        # Same if target us base, instead...
+                        target = target.copy()
+                    target[:] = numeric.where(fromnumeric.asarray(domain), 
+                                              self.fill_other, target)
+                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)       
+        return instance.__class__(method(target, *args), mask=m,
+                                  fill_value=instance.fill_value,
+                                  **instance.options)
+#......................................
+class _compamethods(object):
+    """Defines comparison methods (eq, ge, gt...).
+Instead of calling a ufunc, the method of the masked object is called.
+    """
+    def __init__ (self, methodname, fill_self=0, fill_other=0):
+        """
+: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.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(numpy, 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 = mask_or(instance._mask, getmask(other), small_mask=False)
+        base = instance.filled(self.fill_self)
+        target = filled(other, self.fill_other)
+        method = getattr(base, self.methodname)      
+        return instance.__class__(method(target, *args), mask=m,
+                                  **instance.options)  
+#..........................................................
+MaskedArray.__add__ = _arithmethods('__add__')
+MaskedArray.__radd__ = _arithmethods('__add__')
+MaskedArray.__sub__ = _arithmethods('__sub__')
+MaskedArray.__rsub__ = _arithmethods('__rsub__')
+MaskedArray.__pow__ = _arithmethods('__pow__')
+MaskedArray.__mul__ = _arithmethods('__mul__', 1, 1)
+MaskedArray.__rmul__ = _arithmethods('__mul__', 1, 1)
+MaskedArray.__div__ = _arithmethods('__div__', 0, 1, 
+                                    domain_safe_divide())
+MaskedArray.__rdiv__ = _arithmethods('__rdiv__', 1, 0, 
+                                     domain_safe_divide())
+MaskedArray.__truediv__ = _arithmethods('__truediv__', 0, 1, 
+                                        domain_safe_divide())
+MaskedArray.__rtruediv__ = _arithmethods('__rtruediv__', 1, 0, 
+                                         domain_safe_divide())
+MaskedArray.__floordiv__ = _arithmethods('__floordiv__', 0, 1, 
+                                         domain_safe_divide())
+MaskedArray.__rfloordiv__ = _arithmethods('__rfloordiv__', 1, 0, 
+                                          domain_safe_divide())
+MaskedArray.__eq__ = _compamethods('__eq__')
+MaskedArray.__ne__ = _compamethods('__ne__')
+MaskedArray.__le__ = _compamethods('__le__')
+MaskedArray.__lt__ = _compamethods('__lt__')
+MaskedArray.__ge__ = _compamethods('__ge__')
+MaskedArray.__gt__ = _compamethods('__gt__')
+#####--------------------------------------------------------------------------
+#---- --- 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, 
+          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...
+    """
+    return MaskedArray(data, mask=mask, dtype=dtype, copy=copy, 
+                       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
+
+
+#####--------------------------------------------------------------------------
+#---- --- Patch methods ---
+#####--------------------------------------------------------------------------
+#class _arraymethod:
+#    """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.__doc__ = getattr(ndarray, self._name).__doc__
+#    def __call__(self, instance, *args, **params):
+#        methodname = self._name
+#        (d,m) = (instance._data, instance._mask)
+#        C = instance.__class__
+#        if m is nomask:
+#            return C(getattr(d,methodname).__call__(*args, **params))
+#        elif self._onmask:
+#            return C(getattr(d,methodname).__call__(*args, **params),
+#                     mask=getattr(m,methodname)(*args, **params) )
+#        else:
+#            return C(getattr(d,methodname).__call__(*args, **params), mask=m) 
+#        
+#    def patch(self):
+#        "Adds the new method to MaskedArray."
+#        return types.MethodType(self, None, MaskedArray)
+##......................................
+#MaskedArray.conj = MaskedArray.conjugate = _arraymethod('conjugate').patch()  
+#MaskedArray.diagonal = _arraymethod('diagonal').patch()
+#MaskedArray.take = _arraymethod('take').patch()
+#MaskedArray.ravel = _arraymethod('ravel').patch()
+#MaskedArray.transpose = _arraymethod('transpose').patch()
+#MaskedArray.T = _arraymethod('transpose').patch()
+#MaskedArray.swapaxes = _arraymethod('swapaxes').patch()
+#MaskedArray.clip = _arraymethod('clip', onmask=False).patch()
+#MaskedArray.compress = _arraymethod('compress').patch()
+#MaskedArray.resize = _arraymethod('resize').patch()
+#MaskedArray.copy = _arraymethod('copy').patch()
+
+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)."
+        try:
+            return getattr(MaskedArray, self._name).__doc__
+        except:
+            return getattr(numpy, self._name).__doc__
+    #
+    def __get__(self, obj, objtype=None):
+        self.obj = obj
+        return self
+    #
+    def __call__(self, *args, **params):
+        methodname = self._name
+        obj = self.obj
+        (d, m) = (obj._data, obj._mask)
+        (t, f) = (obj.dtype, obj._fill_value)
+        C = self.obj.__class__
+        if m is nomask:
+            return C(getattr(d,methodname).__call__(*args, **params),
+                     dtype=t, fill_value=f)
+        elif self._onmask:
+            return C(getattr(d,methodname).__call__(*args, **params),
+                     mask=getattr(m,methodname)(*args, **params),
+                     dtype=t, fill_value=f, **obj.options)
+        else:
+            return C(getattr(d,methodname).__call__(*args, **params), mask=m,
+                     dtype=t, fill_value=f, **obj.options) 
+#......................................
+MaskedArray.conj = MaskedArray.conjugate = _arraymethod('conjugate') 
+MaskedArray.copy = _arraymethod('copy') 
+MaskedArray.diagonal = _arraymethod('diagonal')
+MaskedArray.take = _arraymethod('take')
+MaskedArray.ravel = _arraymethod('ravel')
+MaskedArray.transpose = _arraymethod('transpose')
+MaskedArray.T = property(fget=lambda self:self.transpose())
+MaskedArray.swapaxes = _arraymethod('swapaxes')
+MaskedArray.clip = _arraymethod('clip', onmask=False)
+MaskedArray.compress = _arraymethod('compress')
+MaskedArray.copy = _arraymethod('copy')
+MaskedArray.squeeze = _arraymethod('squeeze')
+
+#####--------------------------------------------------------------------------
+#---- --- Extrema functions ---
+#####--------------------------------------------------------------------------
+class _minimum_operation:
+    "Object to calculate minima"
+    def __init__ (self):
+        """minimum(a, b) or minimum(a)
+In one argument case, returns the scalar minimum.
+        """
+        pass
+    #.........
+    def __call__ (self, a, b=None):
+        "Execute the call behavior."
+        if b is None:
+            m = getmask(a)
+            if m is nomask:
+                d = amin(filled(a).ravel())
+                return d
+            ac = a.compressed()
+            if len(ac) == 0:
+                return masked
+            else:
+                return amin(ac)
+        else:
+            return where(less(a, b), a, b)
+    #.........
+    def reduce(self, target, axis=0):
+        """Reduces `target` along the given `axis`."""
+        m = getmask(target)
+        if m is nomask:
+            t = filled(target)
+            return masked_array (umath.minimum.reduce (t, axis))
+        else:
+            t = umath.minimum.reduce(filled(target, minimum_fill_value(target)), 
+                                     axis)
+            m = umath.logical_and.reduce(m, axis)
+#            return masked_array(t, mask=m, fill_value=get_fill_value(target))
+            try:
+                return target.__class__(t, mask=m, dtype=t.dtype,
+                                        fill_value=get_fill_value(target))
+            except AttributeError:
+                return masked_array(t, mask=m, dtype=t.dtype,
+                                    fill_value=get_fill_value(target))
+    #.........
+    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)
+        d = umath.minimum.outer(filled(a), filled(b))
+        return masked_array(d, mask=m)
+
+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)
+#................................................
+class _maximum_operation:
+    "Object to calculate maxima"
+    def __init__ (self):
+        """maximum(a, b) or maximum(a)
+           In one argument case returns the scalar maximum.
+        """
+        pass
+    #.........
+    def __call__ (self, a, b=None):
+        "Executes the call behavior."
+        if b is None:
+            m = getmask(a)
+            if m is nomask:
+                d = amax(filled(a).ravel())
+                return d
+            ac = a.compressed()
+            if len(ac) == 0:
+                return masked
+            else:
+                return amax(ac)
+        else:
+            return where(greater(a, b), a, b)
+    #.........
+    def reduce (self, target, axis=0):
+        """Reduces target along the given axis."""
+        m = getmask(target)
+        if m is nomask:
+            t = filled(target)
+            return masked_array(umath.maximum.reduce (t, axis))
+        else:
+            t = umath.maximum.reduce(filled(target, maximum_fill_value(target)), 
+                                     axis)
+            m = umath.logical_and.reduce(m, axis)
+            try:
+                return target.__class__(t, mask=m, dtype=t.dtype,
+                                        fill_value=get_fill_value(target))
+            except AttributeError:
+                return masked_array(t, mask=m, dtype=t.dtype,
+                                    fill_value=get_fill_value(target))
+    #.........
+    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)
+        d = umath.maximum.outer(filled(a), filled(b))
+        return masked_array(d, mask=m)
+
+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)
+#................................................
+MaskedArray.min = min
+MaskedArray.max = max
+MaskedArray.ptp = ptp
+
+#####---------------------------------------------------------------------------
+#---- --- 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, x, *args, **params):
+        if isinstance(x, MaskedArray):
+            return getattr(x, self._methodname).__call__(*args, **params)
+        #FIXME: As x is not a MaskedArray, we transform it to a ndarray with asarray
+        #FIXME: ... and call the corresponding method.
+        #FIXME: Except that sometimes it doesn't work (try reshape([1,2,3,4],(2,2)))
+        #FIXME: we end up with a "SystemError: NULL result without error in PyObject_Call"
+        #FIXME: A dirty trick is then to call the initial numpy function...
+        method = getattr(fromnumeric.asarray(x), self._methodname)
+        try:
+            return method(*args, **params)
+        except SystemError:
+            return getattr(numpy,self._methodname).__call__(x, *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', 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)
+    return d.argsort(axis, kind)
+
+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)
+    if axis is None:
+        return d.argmin(axis=None)
+    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)
+    if axis is None:
+        return d.argmax(axis=None)
+    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"
+    #TODO: We lose the subclass, here! We should keep track of the classes...
+    #TODO: ...and find the max ? the lowest according to MRO?
+    d = []
+    for x in arrays:
+        d.append(filled(x))
+    d = numeric.concatenate(d, axis)
+    for x in arrays:
+        if getmask(x) is not nomask: 
+            break
+    else:
+        return masked_array(d)
+    dm = []
+    for x in arrays:
+        dm.append(getmaskarray(x))
+    dm = make_mask(numeric.concatenate(dm, axis), copy=False, small_mask=True)
+    return masked_array(d, mask=dm)
+
+def expand_dims(x,axis):
+    """Expand the shape of a by including newaxis before given axis."""
+    if isinstance(x, MaskedArray):
+        (d,m) = (x._data, x._mask)
+        if m is nomask:
+            return masked_array(n_expand_dims(d,axis),                                
+                                dtype=d.dtype, fill_value=x._fill_value)
+        else:
+            return masked_array(n_expand_dims(d,axis), 
+                                mask=n_expand_dims(m,axis),
+                                dtype=d.dtype, fill_value=x._fill_value)
+    else:
+        return n_expand_dims(x,axis)
+
+#......................................
+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)
+    if isinstance(x, MaskedArray):
+        result = x.__class__(fromnumeric.resize(filled(x), new_shape), mask=m)
+    else:
+        result = masked_array(fromnumeric.resize(filled(x), new_shape), mask=m)
+    result.set_fill_value(get_fill_value(x))
+    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."""
+    if not hasattr(a, "_mask"):
+        mask = nomask
+    else:
+        mask = a._mask
+    if out is None:
+        return a.__class__(fromnumeric.round_(a, decimals, None), mask=mask)
+    else:
+        out = a.__class__(fromnumeric.round_(a, decimals, out), mask=mask)
+        return out
+
+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 masked_array(numeric.empty(new_shape, dtype), mask=nomask)
+
+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 masked_array(numeric.empty_like(a), mask=nomask)
+
+def ones(new_shape, dtype=float):
+    """ones(shape, dtype=None) 
+Returns an array of the given dimensions, initialized to all ones."""
+    return masked_array(numeric.ones(new_shape, dtype), mask=nomask)
+
+def zeros(new_shape, dtype=float):
+    """zeros(new_shape, dtype=None) 
+Returns an array of the given dimensions, initialized to all zeros."""
+    return masked_array(numeric.zeros(new_shape, dtype), mask=nomask)
+
+#####--------------------------------------------------------------------------
+#---- --- Pickling ---
+#####--------------------------------------------------------------------------
+#FIXME: We're kinda stuck with forcing the mask to have the same shape as the data
+def _mareconstruct(subtype, baseshape, basetype,):
+    """Internal function that builds a new MaskedArray from the information stored
+in a pickle."""
+    _data = ndarray.__new__(ndarray, baseshape, basetype)
+    _mask = ndarray.__new__(ndarray, baseshape, basetype)
+    return MaskedArray.__new__(subtype, _data, mask=_mask, dtype=basetype, small_mask=False)
+
+def _getstate(a):
+    "Returns the internal state of the masked array, for pickling purposes."
+    state = (1,
+             a.shape, 
+             a.dtype,
+             a.flags.fnc,
+             (a._data).__reduce__()[-1][-1],
+             getmaskarray(a).__reduce__()[-1][-1])
+    return state
+    
+def _setstate(a, 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) = state
+    (a._data).__setstate__((shp, typ, isf, raw))
+    (a._mask).__setstate__((shp, dtype('|b1'), isf, msk))
+        
+def _reduce(a):
+    """Returns a 3-tuple for pickling a MaskedArray."""
+    return (_mareconstruct,
+            (a.__class__, (0,), 'b', ),
+            a.__getstate__())
+    
+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)
+
+MaskedArray.__getstate__ = _getstate
+MaskedArray.__setstate__ = _setstate
+MaskedArray.__reduce__ = _reduce
+MaskedArray.__dump__ = dump
+MaskedArray.__dumps__ = dumps
+
+################################################################################
+
+if __name__ == '__main__':
+    import numpy as N
+    if 1:
+        x = N.array([[ 0.13,  0.26,  0.90],
+                     [ 0.28,  0.33,  0.63],
+                     [ 0.31,  0.87,  0.70]])
+        m = N.array([[ True, False, False],
+                     [False, False, False],
+                     [True,  True, False]], dtype=N.bool_)
+        X = N.asmatrix(x)
+        mX = masked_array(X, mask=m)
\ No newline at end of file


Property changes on: trunk/Lib/sandbox/maskedarray/core_ini.py
___________________________________________________________________
Name: svn:keywords
   + Date 
Author 
Revision
Id

Modified: trunk/Lib/sandbox/maskedarray/extras.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/extras.py	2007-02-11 00:33:13 UTC (rev 2699)
+++ trunk/Lib/sandbox/maskedarray/extras.py	2007-02-11 23:24:17 UTC (rev 2700)
@@ -193,7 +193,7 @@
     if asscalar:
         dtypes.append(numeric.asarray(res).dtype)
         outarr = zeros(outshape, object_)
-        outarr[ind] = res
+        outarr[tuple(ind)] = res
         Ntot = numeric.product(outshape)
         k = 1
         while k < Ntot:
@@ -206,7 +206,7 @@
                 n -= 1
             i.put(indlist,ind)
             res = func1d(arr[tuple(i.tolist())],*args)
-            outarr[ind] = res
+            outarr[tuple(ind)] = res
             dtypes.append(asarray(res).dtype)
             k += 1
     else:

Modified: trunk/Lib/sandbox/maskedarray/tests/test_core.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/tests/test_core.py	2007-02-11 00:33:13 UTC (rev 2699)
+++ trunk/Lib/sandbox/maskedarray/tests/test_core.py	2007-02-11 23:24:17 UTC (rev 2700)
@@ -1,6 +1,5 @@
 # pylint: disable-msg=W0611, W0612, W0511,R0201
-"""Tests suite for MaskedArray.
-Adapted from the original test_ma by Pierre Gerard-Marchant
+"""Tests suite for MaskedArray & subclassing.
 
 :author: Pierre Gerard-Marchant
 :contact: pierregm_at_uga_dot_edu
@@ -22,8 +21,8 @@
 reload(maskedarray.testutils)
 from maskedarray.testutils import *
 
-import maskedarray.core
-reload(maskedarray.core)
+import maskedarray.core as coremodule
+reload(coremodule)
 from maskedarray.core import *
 
 pi = N.pi
@@ -50,7 +49,7 @@
         xm.set_fill_value(1.e+20)  
         self.d = (x, y, a10, m1, m2, xm, ym, z, zm, xf)
     #........................
-    def check_testBasic1d(self):
+    def check_basic1d(self):
         "Test of basic array creation and properties in 1 dimension."
         (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
         assert(not isMaskedArray(x))
@@ -68,7 +67,7 @@
         assert_array_equal(filled(xm, 1.e20), xf)
         assert_array_equal(x, xm)   
     #........................
-    def check_testBasic2d(self):
+    def check_basic2d(self):
         "Test of basic array creation and properties in 2 dimensions."
         (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
         for s in [(4,3), (6,2)]:
@@ -88,7 +87,7 @@
             assert_equal(filled(xm, 1.e20), xf)
             assert_equal(x, xm)
     #........................
-    def check_testArithmetic (self):
+    def check_basic_arithmetic (self):
         "Test of basic arithmetic."
         (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
         a2d = array([[1,2],[0,4]])
@@ -123,14 +122,131 @@
             assert_equal(N.multiply(x,y), multiply(xm, ym))
             assert_equal(N.divide(x,y), divide(xm, ym))
     #........................
-    def check_testMixedArithmetic(self):
+    def check_mixed_arithmetic(self):
         "Tests mixed arithmetics."
         na = N.array([1])
         ma = array([1])
         self.failUnless(isinstance(na + ma, MaskedArray))
-        self.failUnless(isinstance(ma + na, MaskedArray))    
+        self.failUnless(isinstance(ma + na, MaskedArray))   
+    #........................
+    def check_inplace_arithmetic(self):
+        """Test of inplace operations and rich comparisons"""
+        # addition
+        x = arange(10)
+        y = arange(10)
+        xm = arange(10)
+        xm[2] = masked
+        x += 1
+        assert_equal(x, y+1)
+        xm += 1
+        assert_equal(xm, y+1)
+        # subtraction
+        x = arange(10)
+        xm = arange(10)
+        xm[2] = masked
+        x -= 1
+        assert_equal(x, y-1)
+        xm -= 1
+        assert_equal(xm, y-1)
+        # multiplication
+        x = arange(10)*1.0
+        xm = arange(10)*1.0
+        xm[2] = masked
+        x *= 2.0
+        assert_equal(x, y*2)
+        xm *= 2.0
+        assert_equal(xm, y*2)
+        # division
+        x = arange(10)*2
+        xm = arange(10)*2
+        xm[2] = masked
+        x /= 2
+        assert_equal(x, y)
+        xm /= 2
+        assert_equal(xm, y)
+        # division, pt 2
+        x = arange(10)*1.0
+        xm = arange(10)*1.0
+        xm[2] = masked
+        x /= 2.0
+        assert_equal(x, y/2.0)
+        xm /= arange(10)
+        assert_equal(xm, ones((10,)))
+
+        x = arange(10).astype(float_)
+        xm = arange(10)
+        xm[2] = masked
+#        id1 = id(x.raw_data())
+        id1 = x.raw_data().ctypes.data
+        x += 1.
+#        assert id1 == id(x.raw_data())
+        assert (id1 == x.raw_data().ctypes.data)
+        assert_equal(x, y+1.)
+        # addition w/ array
+        x = arange(10, dtype=float_)
+        xm = arange(10, dtype=float_)
+        xm[2] = masked
+        m = xm.mask
+        a = arange(10, dtype=float_)
+        a[-1] = masked
+        x += a
+        xm += a
+        assert_equal(x,y+a)
+        assert_equal(xm,y+a)
+        assert_equal(xm.mask, mask_or(m,a.mask))
+        # subtraction w/ array
+        x = arange(10, dtype=float_)
+        xm = arange(10, dtype=float_)
+        xm[2] = masked
+        m = xm.mask
+        a = arange(10, dtype=float_)
+        a[-1] = masked
+        x -= a
+        xm -= a
+        assert_equal(x,y-a)
+        assert_equal(xm,y-a)
+        assert_equal(xm.mask, mask_or(m,a.mask))        
+        # multiplication w/ array
+        x = arange(10, dtype=float_)
+        xm = arange(10, dtype=float_)
+        xm[2] = masked
+        m = xm.mask
+        a = arange(10, dtype=float_)
+        a[-1] = masked
+        x *= a
+        xm *= a
+        assert_equal(x,y*a)
+        assert_equal(xm,y*a)
+        assert_equal(xm.mask, mask_or(m,a.mask))        
+        # division w/ array      
+        x = arange(10, dtype=float_)
+        xm = arange(10, dtype=float_)
+        xm[2] = masked
+        m = xm.mask
+        a = arange(10, dtype=float_)
+        a[-1] = masked
+        x /= a
+        xm /= a
+        assert_equal(x,y/a)
+        assert_equal(xm,y/a)
+        assert_equal(xm.mask, mask_or(mask_or(m,a.mask), (a==0))) 
+    #..........................
+    def check_scalararithmetic(self):
+        "Tests some scalar arithmetics on MaskedArrays."
+        xm = array(0, mask=1)
+        assert((1/array(0)).mask)
+        assert((1 + xm).mask)
+        assert((-xm).mask)
+        assert((-xm).mask)
+        assert(maximum(xm, xm).mask)
+        assert(minimum(xm, xm).mask)
+        assert(xm.filled().dtype is xm.data.dtype)
+        x = array(0, mask=0)
+#        assert(x.filled() is x.data)
+        assert_equal(x.filled().ctypes.data, x.ctypes.data)
+        assert_equal(str(xm), str(masked_print_option))
     #.........................
-    def check_testUfuncs1 (self):
+    def check_basic_ufuncs (self):
         "Test various functions such as sin, cos."
         (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
         assert_equal(N.cos(x), cos(xm))
@@ -155,15 +271,11 @@
         assert_equal(N.less_equal(x,y), less_equal(xm, ym))
         assert_equal(N.greater_equal(x,y), greater_equal(xm, ym))
         assert_equal(N.conjugate(x), conjugate(xm))
-        assert_equal(N.concatenate((x,y)), concatenate((xm,ym)))
-        assert_equal(N.concatenate((x,y)), concatenate((x,y)))
-        assert_equal(N.concatenate((x,y)), concatenate((xm,y)))
-        assert_equal(N.concatenate((x,y,x)), concatenate((x,ym,x)))
     #........................
-    def check_xtestCount (self):
+    def check_count_func (self):
         "Tests count"
         ott = array([0.,1.,2.,3.], mask=[1,0,0,0])
-        assert( isinstance(count(ott), types.IntType))
+        assert( isinstance(count(ott), int))
         assert_equal(3, count(ott))
         assert_equal(1, count(1))
         assert_equal(0, array(1,mask=[1]))
@@ -174,7 +286,7 @@
         assert getmask(count(ott,0)) is nomask
         assert_equal([1,2],count(ott,0))
     #........................
-    def check_testMinMax (self):
+    def check_minmax_func (self):
         "Tests minimum and maximum."
         (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
         xr = N.ravel(x) #max doesn't work if shaped
@@ -192,8 +304,25 @@
         assert_equal(maximum(x,y), where(greater(x,y), x, y))
         assert minimum(x) == 0
         assert maximum(x) == 4
+        
+    def check_minmax_methods(self):        
+        "Additional tests on max/min"
+        (_, _, _, _, _, xm, _, _, _, _) = self.d
+        xm.shape = (xm.size,)
+        assert_equal(xm.max(), 10)
+        assert(xm[0].max() is masked)
+        assert(xm[0].max(0) is masked)
+        assert(xm[0].max(-1) is masked)
+        assert_equal(xm.min(), -10.)
+        assert(xm[0].min() is masked)
+        assert(xm[0].min(0) is masked)
+        assert(xm[0].min(-1) is masked)
+        assert_equal(xm.ptp(), 20.)
+        assert(xm[0].ptp() is masked)
+        assert(xm[0].ptp(0) is masked)
+        assert(xm[0].ptp(-1) is masked)
     #........................
-    def check_testAddSumProd (self):
+    def check_addsumprod (self):
         "Tests add, sum, product."
         (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
         assert_equal(N.add.reduce(x), add.reduce(x))
@@ -217,6 +346,12 @@
     def check_concat(self):
         "Tests concatenations."
         (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
+        # basic concatenation
+        assert_equal(N.concatenate((x,y)), concatenate((xm,ym)))
+        assert_equal(N.concatenate((x,y)), concatenate((x,y)))
+        assert_equal(N.concatenate((x,y)), concatenate((xm,y)))
+        assert_equal(N.concatenate((x,y,x)), concatenate((x,ym,x)))
+        # Concatenation along an axis
         s = (3,4)
         x.shape = y.shape = xm.shape = ym.shape = s
         assert_equal(xm.mask, N.reshape(m1, s))
@@ -225,7 +360,7 @@
         assert_equal(N.concatenate((x,y),1), xmym)
         assert_equal(N.concatenate((xm.mask,ym.mask),1), xmym._mask)
     #........................
-    def check_testCI(self):
+    def check_indexing(self):
         "Tests conversions and indexing"
         x1 = N.array([1,2,4,3])
         x2 = array(x1, mask=[1,0,0,0])
@@ -275,7 +410,7 @@
         assert_equal(s1, s2)
         assert x1[1:1].shape == (0,)
     #........................
-    def check_testCopySize(self):
+    def check_copy(self):
         "Tests of some subtle points of copying and sizing."
         n = [0,0,1,0,0]
         m = make_mask(n)
@@ -286,21 +421,27 @@
 
         x1 = N.arange(5)
         y1 = array(x1, mask=m)
-        assert( y1._data is x1)
+        #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( 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.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.mask is 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 is not m)
+        assert (y2._mask.__array_interface__ != m.__array_interface__)
         assert( allequal(y2.mask, 0))
 
         y3 = array(x1*1.0, mask=m)
@@ -327,14 +468,18 @@
         x = masked_array([1,2,3], mask=[0,1,0])
         # Copy is False by default
         y = masked_array(x)
-        assert_equal(id(y._data), id(x._data))
-        assert_equal(id(y._mask), id(x._mask))
+#        assert_equal(id(y._data), id(x._data))
+#        assert_equal(id(y._mask), id(x._mask))
+        assert_equal(y._data.ctypes.data, x._data.ctypes.data)
+        assert_equal(y._mask.ctypes.data, x._mask.ctypes.data)
         y = masked_array(x, copy=True)
-        assert_not_equal(id(y._data), id(x._data))
-        assert_not_equal(id(y._mask), id(x._mask))
+#        assert_not_equal(id(y._data), id(x._data))
+#        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_testOddFeatures_1(self):
+    def check_oddfeatures_1(self):
         "Test of other odd features"
         x = arange(20)
         x = x.reshape(4,5)
@@ -380,7 +525,7 @@
         assert_equal(x,z)
         #
     #........................
-    def check_testOddFeatures_2(self):
+    def check_oddfeatures_2(self):
         "Tests some more features."
         x = array([1.,2.,3.,4.,5.])
         c = array([1,1,1,0,0])
@@ -410,7 +555,7 @@
         z = where(c, 1, masked)
         assert_equal(z, [99, 1, 1, 99, 99, 99])
     #........................
-    def check_testOddFeatures_3(self):
+    def check_oddfeatures_3(self):
         """Tests some generic features."""
         atest = ones((10,10,10), dtype=float_)
         btest = zeros(atest.shape, MaskType)
@@ -436,7 +581,7 @@
         assert_equal(masked_not_equal(array([2,2,1,2,1], mask=[1,0,0,0,0]), 2).mask, [1,0,1,0,1])
         assert_equal(masked_where([1,1,0,0,0], [1,2,3,4,5]), [99,99,3,4,5])
     #........................
-    def check_testTakeTransposeInnerOuter(self):
+    def check_TakeTransposeInnerOuter(self):
         "Test of take, transpose, inner, outer products"
         x = arange(24)
         y = N.arange(24)
@@ -456,107 +601,7 @@
         assert t[1] == 2
         assert t[2] == 3
     #........................
-    def check_testInplace(self):
-        """Test of inplace operations and rich comparisons"""
-        y = arange(10)
-
-        x = arange(10)
-        xm = arange(10)
-        xm[2] = masked
-        x += 1
-        assert_equal(x, y+1)
-        xm += 1
-        assert_equal(xm, y+1)
-
-        x = arange(10)
-        xm = arange(10)
-        xm[2] = masked
-        x -= 1
-        assert_equal(x, y-1)
-        xm -= 1
-        assert_equal(xm, y-1)
-
-        x = arange(10)*1.0
-        xm = arange(10)*1.0
-        xm[2] = masked
-        x *= 2.0
-        assert_equal(x, y*2)
-        xm *= 2.0
-        assert_equal(xm, y*2)
-
-        x = arange(10)*2
-        xm = arange(10)*2
-        xm[2] = masked
-        x /= 2
-        assert_equal(x, y)
-        xm /= 2
-        assert_equal(xm, y)
-
-        x = arange(10)*1.0
-        xm = arange(10)*1.0
-        xm[2] = masked
-        x /= 2.0
-        assert_equal(x, y/2.0)
-        xm /= arange(10)
-        assert_equal(xm, ones((10,)))
-
-        x = arange(10).astype(float_)
-        xm = arange(10)
-        xm[2] = masked
-        id1 = id(x.raw_data())
-        x += 1.
-        assert id1 == id(x.raw_data())
-        assert_equal(x, y+1.)
-        
-        x = arange(10, dtype=float_)
-        xm = arange(10, dtype=float_)
-        xm[2] = masked
-        m = xm.mask
-        a = arange(10, dtype=float_)
-        a[-1] = masked
-        x += a
-        xm += a
-        assert_equal(x,y+a)
-        assert_equal(xm,y+a)
-        assert_equal(xm.mask, mask_or(m,a.mask))
-        
-        x = arange(10, dtype=float_)
-        xm = arange(10, dtype=float_)
-        xm[2] = masked
-        m = xm.mask
-        a = arange(10, dtype=float_)
-        a[-1] = masked
-        x -= a
-        xm -= a
-        assert_equal(x,y-a)
-        assert_equal(xm,y-a)
-        assert_equal(xm.mask, mask_or(m,a.mask))        
-        
-        x = arange(10, dtype=float_)
-        xm = arange(10, dtype=float_)
-        xm[2] = masked
-        m = xm.mask
-        a = arange(10, dtype=float_)
-        a[-1] = masked
-        x *= a
-        xm *= a
-        assert_equal(x,y*a)
-        assert_equal(xm,y*a)
-        assert_equal(xm.mask, mask_or(m,a.mask))        
-                
-        x = arange(10, dtype=float_)
-        xm = arange(10, dtype=float_)
-        xm[2] = masked
-        m = xm.mask
-        a = arange(10, dtype=float_)
-        a[-1] = masked
-        x /= a
-        xm /= a
-        assert_equal(x,y/a)
-        assert_equal(xm,y/a)
-        assert_equal(xm.mask, mask_or(mask_or(m,a.mask), (a==0)))
-    #........................
-    def check_testPickle(self):
+    def check_pickling(self):
         "Test of pickling"
         import pickle
         x = arange(12)
@@ -566,7 +611,7 @@
         y = pickle.loads(s)
         assert_equal(x,y)
     #.......................
-    def check_testMasked(self):
+    def check_maskedelement(self):
         "Test of masked element"
         x = arange(6)
         x[1] = masked
@@ -579,7 +624,7 @@
         #self.failUnlessRaises(Exception, lambda x,y: x+y, masked, xx)
         #self.failUnlessRaises(Exception, lambda x,y: x+y, xx, masked)
     #........................
-    def check_testToPython(self):
+    def check_topython(self):
         "Tests some communication issues with Python."
         assert_equal(1, int(array(1)))
         assert_equal(1.0, float(array(1)))
@@ -592,22 +637,8 @@
 #TODO:        self.failUnless(bool(array([0,0],mask=[0,1])))
 #TODO:        self.failIf(bool(array([0,0])))
 #TODO:        self.failIf(bool(array([0,0],mask=[0,0])))
-    #..........................
-    def check_testScalarArithmetic(self):
-        "Tests some scalar arithmetics on MaskedArrays."
-        xm = array(0, mask=1)
-        assert((1/array(0)).mask)
-        assert((1 + xm).mask)
-        assert((-xm).mask)
-        assert((-xm).mask)
-        assert(maximum(xm, xm).mask)
-        assert(minimum(xm, xm).mask)
-        assert(xm.filled().dtype is xm.data.dtype)
-        x = array(0, mask=0)
-        assert(x.filled() is x.data)
-        assert_equal(str(xm), str(masked_print_option))
     #........................
-    def check_testArrayMethods(self):
+    def check_arraymethods(self):
         "Tests some MaskedArray methods."
         a = array([1,3,2])
         b = array([1,3,2], mask=[1,0,1])
@@ -626,7 +657,7 @@
         assert_equal(a.take([1,2]), a.data.take([1,2]))
         assert_equal(m.transpose(), m.data.transpose())
     #........................
-    def check_testArrayAttributes(self):
+    def check_basicattributes(self):
         "Tests some basic array attributes."
         a = array([1,3,2])
         b = array([1,3,2], mask=[1,0,1])
@@ -637,7 +668,7 @@
         assert_equal(a.shape, (3,))
         assert_equal(b.shape, (3,))
     #........................
-    def check_testSingleElementSubscript(self):
+    def check_single_element_subscript(self):
         "Tests single element subscripts of Maskedarrays."
         a = array([1,3,2])
         b = array([1,3,2], mask=[1,0,1])
@@ -669,7 +700,7 @@
         assert_equal(X._mask, x.mask)
         assert_equal(getmask(x), [0,0,1,0,0])
         
-#..............................................................................
+#...............................................................................
         
 class test_ufuncs(NumpyTestCase):
     "Test class for the application of ufuncs on MaskedArrays."
@@ -704,7 +735,7 @@
                 uf = getattr(umath, f)
             except AttributeError:
                 uf = getattr(fromnumeric, f)
-            mf = getattr(maskedarray, f)
+            mf = getattr(coremodule, f)
             args = self.d[:uf.nin]
             ur = uf(*args)
             mr = mf(*args)
@@ -730,7 +761,7 @@
         assert(amask.max(1)[0].mask)
         assert(amask.min(1)[0].mask)
 
-#..............................................................................
+#...............................................................................
 
 class test_array_methods(NumpyTestCase):
     "Test class for miscellaneous MaskedArrays methods."
@@ -767,7 +798,7 @@
         self.d =  (x,X,XX,m,mx,mX,mXX,m2x,m2X,m2XX)
 
     #------------------------------------------------------
-    def test_trace(self):
+    def check_trace(self):
         "Tests trace on MaskedArrays."
         (x,X,XX,m,mx,mX,mXX,m2x,m2X,m2XX) = self.d
         mXdiag = mX.diagonal()
@@ -775,7 +806,7 @@
         assert_almost_equal(mX.trace(),
                             X.trace() - sum(mXdiag.mask*X.diagonal(),axis=0))
 
-    def test_clip(self):
+    def check_clip(self):
         "Tests clip on MaskedArrays."
         (x,X,XX,m,mx,mX,mXX,m2x,m2X,m2XX) = self.d
         clipped = mx.clip(2,8)
@@ -783,7 +814,7 @@
         assert_equal(clipped.data,x.clip(2,8))
         assert_equal(clipped.data,mx.data.clip(2,8))
 
-    def test_ptp(self):
+    def check_ptp(self):
         "Tests ptp on MaskedArrays."
         (x,X,XX,m,mx,mX,mXX,m2x,m2X,m2XX) = self.d
         (n,m) = X.shape
@@ -797,7 +828,7 @@
         assert_equal(mX.ptp(0),cols)
         assert_equal(mX.ptp(1),rows)
 
-    def test_swapaxes(self):
+    def check_swapaxes(self):
         "Tests swapaxes on MaskedArrays."
         (x,X,XX,m,mx,mX,mXX,m2x,m2X,m2XX) = self.d
         mXswapped = mX.swapaxes(0,1)
@@ -805,7 +836,7 @@
         mXXswapped = mXX.swapaxes(0,2)
         assert_equal(mXXswapped.shape,(2,2,3,3))
 
-    def test_cumsumprod(self):
+    def check_cumsumprod(self):
         "Tests cumsum & cumprod on MaskedArrays."
         (x,X,XX,m,mx,mX,mXX,m2x,m2X,m2XX) = self.d
         mXcp = mX.cumsum(0)
@@ -818,7 +849,7 @@
         mXcp = mX.cumprod(1)
         assert_equal(mXcp.data,mX.filled(1).cumprod(1))
                                                       
-    def test_varstd(self):
+    def check_varstd(self):
         "Tests var & std on MaskedArrays."
         (x,X,XX,m,mx,mX,mXX,m2x,m2X,m2XX) = self.d
         assert_almost_equal(mX.var(axis=None),mX.compressed().var())
@@ -831,7 +862,7 @@
             assert_almost_equal(mXvar0[k],mX[:,k].compressed().var())
             assert_almost_equal(N.sqrt(mXvar0[k]), mX[:,k].compressed().std())
     
-    def test_argmin(self):
+    def check_argmin(self):
         "Tests argmin & argmax on MaskedArrays."
         (x,X,XX,m,mx,mX,mXX,m2x,m2X,m2XX) = self.d
         #
@@ -884,6 +915,15 @@
         put(x, i, masked_array([0,2,4,6],[1,0,1,0]))
         assert_array_equal(x, [0,1,2,3,4,5,6,7,8,9,])
         assert_equal(x.mask, [1,0,0,0,1,1,0,0,0,0])  
+        
+    def check_put_hardmask(self):
+        "Tests put on hardmask"        
+        d = arange(5)
+        n = [0,0,0,1,1]
+        m = make_mask(n)
+        xh = array(d+1, mask = m, hard_mask=True, copy=True)
+        xh.put([4,2,0,1,3],[1,2,3,4,5])
+        assert_equal(xh._data, [3,4,2,4,5]) 
     
     def check_take(self):
         "Tests take"
@@ -982,8 +1022,8 @@
         assert_equal(xs._data, [0,10,20,30,40])
         assert(xh.mask is m)
         assert_equal(xs.mask, nomask)
-        xh[0] = maskedarray.core.masked
-        xs[0] = maskedarray.core.masked
+        xh[0] = masked
+        xs[0] = masked
         assert_equal(xh.mask, [1,0,0,1,1])
         assert_equal(xs.mask, [1,0,0,0,0])
         xh[:] = 1
@@ -999,7 +1039,7 @@
         assert_equal(xh.mask, nomask)
         # Switch back to hard mask
         xh.harden_mask()
-        xh[xh<3] = maskedarray.core.masked
+        xh[xh<3] = masked
         assert_equal(xh._data, [0,1,2,3,4])
         assert_equal(xh._mask, [1,1,1,0,0])
         xh[filled(xh>1,False)] = 5
@@ -1060,35 +1100,21 @@
         sortedx = sort(x, endwith=False)
         assert_equal(sortedx._data, [1,2,-2,-1,0])
         assert_equal(sortedx._mask, [1,1,0,0,0])
+    
+    def check_ravel(self):
+        a = array([[1,2,3,4,5]], mask=[[0,1,0,0,0]])
+        aravel = a.ravel()
+        assert_equal(a._mask.shape, a.shape)
+        a = array([0,0], mask=[1,1])
+        aravel = a.ravel()
+        assert_equal(a._mask.shape, a.shape)
+        a = array(N.matrix([1,2,3,4,5]), mask=[[0,1,0,0,0]])
+        aravel = a.ravel()
+        assert_equal(a.shape,(1,5))
+        assert_equal(a._mask.shape, a.shape)
+        
 #..............................................................................
 
-#..............................................................................
-
-class test_subclassing(NumpyTestCase):
-    """Test suite for masked subclasses of ndarray."""
-    
-    class SubArray(N.ndarray):
-        """Defines a generic N.ndarray subclass, that stores some metadata
-        in the  dictionary `info`."""
-        def __new__(cls,arr,info={}):
-            x = N.array(arr).view(cls)
-            x.info = info
-            return x
-        def __array_finalize__(self, obj):
-            if hasattr(obj,'info'):
-                self.info = obj.info
-            return
-            
-    def check_subclassing(self):
-        "Tests whether the subclass is kept."
-        x = N.arange(5)
-        m = [0,0,1,0,0]
-        xsub = test_subclassing.SubArray(x)
-        xmsub = masked_array(xsub, mask=m)
-        assert isinstance(xmsub, MaskedArray)
-        assert_equal(xmsub._data, xsub)
-        assert isinstance(xmsub._data, test_subclassing.SubArray)
-
 ###############################################################################
 #------------------------------------------------------------------------------
 if __name__ == "__main__":

Modified: trunk/Lib/sandbox/maskedarray/tests/test_extras.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/tests/test_extras.py	2007-02-11 00:33:13 UTC (rev 2699)
+++ trunk/Lib/sandbox/maskedarray/tests/test_extras.py	2007-02-11 23:24:17 UTC (rev 2700)
@@ -296,6 +296,15 @@
         assert_equal(dx._data, N.r_[0,difx_d,0])
         assert_equal(dx._mask, N.r_[1,0,0,0,0,1])
 
+class test_apply_along_axis(NumpyTestCase):
+    "Tests 2D functions"
+    def check_3d(self):
+        a = arange(12.).reshape(2,2,3)
+        def myfunc(b):
+            return b[1]
+        xa = apply_along_axis(myfunc,2,a)
+        assert_equal(xa,[[1,4],[7,10]])
+
 ###############################################################################
 #------------------------------------------------------------------------------
 if __name__ == "__main__":

Added: trunk/Lib/sandbox/maskedarray/tests/test_subclassing.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/tests/test_subclassing.py	2007-02-11 00:33:13 UTC (rev 2699)
+++ trunk/Lib/sandbox/maskedarray/tests/test_subclassing.py	2007-02-11 23:24:17 UTC (rev 2700)
@@ -0,0 +1,115 @@
+import numpy as N
+import numpy.core.numeric as numeric
+
+from numpy.testing import NumpyTest, NumpyTestCase
+
+import maskedarray.testutils
+reload(maskedarray.testutils)
+from maskedarray.testutils import *
+
+import maskedarray.core as coremodule
+reload(coremodule)
+from maskedarray.core import *
+
+
+
+class SubArray(N.ndarray):
+    """Defines a generic N.ndarray subclass, that stores some metadata
+    in the  dictionary `info`."""
+    def __new__(cls,arr,info={}):
+        x = N.asanyarray(arr).view(cls)
+        x.info = info
+        return x
+    def __array_finalize__(self, obj):
+        self.info = getattr(obj,'info',{})
+        return
+    def __add__(self, other):
+        result = N.ndarray.__add__(self, other)
+        result.info.update({'added':result.info.pop('added',0)+1})
+        return result
+subarray = SubArray
+
+class MSubArray(SubArray,MaskedArray):
+    def __new__(cls, data, info=None, mask=nomask):
+        subarr = SubArray(data, info)
+        _data = MaskedArray.__new__(cls, data=subarr, mask=mask)
+        _data.info = subarr.info
+        return _data
+    def __array_finalize__(self,obj):
+        SubArray.__array_finalize__(self, obj)
+        MaskedArray.__array_finalize__(self,obj)
+        return
+msubarray = MSubArray
+
+class MMatrix(MaskedArray, N.matrix,):
+    def __new__(cls, data, mask=nomask):
+        mat = N.matrix(data)
+        _data = MaskedArray.__new__(cls, data=mat, mask=mask)
+        return _data
+    def __array_finalize__(self,obj):
+        N.matrix.__array_finalize__(self, obj)
+        MaskedArray.__array_finalize__(self,obj)
+        return     
+mmatrix = MMatrix 
+        
+
+
+class test_subclassing(NumpyTestCase):
+    """Test suite for masked subclasses of ndarray."""
+                
+    def check_data_subclassing(self):
+        "Tests whether the subclass is kept."
+        x = N.arange(5)
+        m = [0,0,1,0,0]
+        xsub = SubArray(x)
+        xmsub = masked_array(xsub, mask=m)
+        assert isinstance(xmsub, MaskedArray)
+        assert_equal(xmsub._data, xsub)
+        assert isinstance(xmsub._data, SubArray)
+        
+    def check_maskedarray_subclassing(self):
+        "Tests subclassing MaskedArray"
+        x = N.arange(5)
+        mx = mmatrix(x,mask=[0,1,0,0,0])
+        assert isinstance(mx._data, N.matrix)
+        "Tests masked_unary_operation"
+        assert isinstance(add(mx,mx), mmatrix)
+        assert isinstance(add(mx,x), mmatrix)
+        assert_equal(add(mx,x), mx+x)
+        assert isinstance(add(mx,mx)._data, N.matrix)
+        assert isinstance(add.outer(mx,mx), mmatrix)
+        "Tests masked_binary_operation"
+        assert isinstance(hypot(mx,mx), mmatrix)
+        assert isinstance(hypot(mx,x), mmatrix)
+
+    def check_attributepropagation(self):
+        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)
+        assert isinstance(z, MaskedArray)
+        assert isinstance(z, MSubArray)
+        assert isinstance(z._data, SubArray)
+        assert z._data.info['added'] > 0
+        
+################################################################################
+if __name__ == '__main__':
+    NumpyTest().run()
+    if 1:
+        x = N.arange(5)
+        mx = mmatrix(x,mask=[0,1,0,0,0])
+        
+        c = mx.ravel()
+        
+        a = add(mx,x)
+        b = mx+x
+        assert_equal(a,b)
+
+


Property changes on: trunk/Lib/sandbox/maskedarray/tests/test_subclassing.py
___________________________________________________________________
Name: svn:keywords
   + Date 
Author 
Revision
Id

Modified: trunk/Lib/sandbox/maskedarray/testutils.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/testutils.py	2007-02-11 00:33:13 UTC (rev 2699)
+++ trunk/Lib/sandbox/maskedarray/testutils.py	2007-02-11 23:24:17 UTC (rev 2700)
@@ -11,6 +11,7 @@
 
 
 import numpy as N
+from numpy.core import ndarray
 from numpy.core.numerictypes import float_
 import numpy.core.umath as umath
 from numpy.testing import NumpyTest, NumpyTestCase
@@ -54,7 +55,6 @@
         for k in range(len(desired)):
             assert_equal(actual[k], desired[k], 'item=%r\n%s' % (k,err_msg))
         return
-    from numpy.core import ndarray
     if isinstance(actual, ndarray) or isinstance(desired, ndarray):
         return assert_array_equal(actual, desired, err_msg)
     msg = build_err_msg([actual, desired], err_msg,)
@@ -98,8 +98,8 @@
     yf = filled(y)
     m = mask_or(getmask(x), getmask(y))
     
-    x = filled(masked_array(xf, copy=False, mask=m), fill_value)
-    y = filled(masked_array(yf, copy=False, mask=m), fill_value)
+    x = masked_array(xf, copy=False, mask=m).filled(fill_value)
+    y = masked_array(yf, copy=False, mask=m).filled(fill_value)
     if (x.dtype.char != "O"):
         x = x.astype(float_)
         if isinstance(x, N.ndarray) and x.size > 1:

Modified: trunk/Lib/sandbox/maskedarray/timer_comparison.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/timer_comparison.py	2007-02-11 00:33:13 UTC (rev 2699)
+++ trunk/Lib/sandbox/maskedarray/timer_comparison.py	2007-02-11 23:24:17 UTC (rev 2700)
@@ -5,8 +5,7 @@
 from numpy import int_, float_, bool_
 import numpy.core.fromnumeric as fromnumeric
 
-from maskedarray.testutils import assert_equal,assert_array_equal, \
-    fail_if_equal, assert_mask_equal
+from numpy.testing.utils import build_err_msg, rand
 
 
 numpy.seterr(all='ignore')
@@ -20,9 +19,10 @@
         self.allequal = module.allequal
         self.arange = module.arange
         self.array = module.array
-        self.average =  module.average
+#        self.average =  module.average
         self.concatenate = module.concatenate
         self.count = module.count
+        self.equal = module.equal
         self.filled = module.filled
         self.getmask = module.getmask 
         self.id = id
@@ -32,6 +32,7 @@
         self.masked_array = module.masked_array
         self.masked_values = module.masked_values
         self.mask_or = module.mask_or
+        self.nomask = module.nomask
         self.ones = module.ones
         self.outer = module.outer
         self.repeat = module.repeat
@@ -45,8 +46,76 @@
             self.umath = module.umath
         except AttributeError:
             self.umath = module.core.umath
+        self.testnames = []
+    #........................
+    def assert_array_compare(self, comparison, x, y, err_msg='', header='', 
+                         fill_value=True):
+        """Asserts that a comparison relation between two masked arrays is satisfied
+        elementwise."""
+        xf = self.filled(x)
+        yf = self.filled(y)
+        m = self.mask_or(self.getmask(x), self.getmask(y))
+        
+        x = self.filled(self.masked_array(xf, mask=m), fill_value)
+        y = self.filled(self.masked_array(yf, mask=m), fill_value)
+        if (x.dtype.char != "O"):
+            x = x.astype(float_)
+            if isinstance(x, numpy.ndarray) and x.size > 1:
+                x[numpy.isnan(x)] = 0
+            elif numpy.isnan(x):
+                x = 0
+        if (y.dtype.char != "O"):
+            y = y.astype(float_)
+            if isinstance(y, numpy.ndarray) and y.size > 1:
+                y[numpy.isnan(y)] = 0
+            elif numpy.isnan(y):
+                y = 0
+        try:
+            cond = (x.shape==() or y.shape==()) or x.shape == y.shape
+            if not cond:
+                msg = build_err_msg([x, y],
+                                    err_msg
+                                    + '\n(shapes %s, %s mismatch)' % (x.shape,
+                                                                      y.shape),
+                                    header=header,
+                                    names=('x', 'y'))
+                assert cond, msg
+            val = comparison(x,y)
+            if m is not self.nomask and fill_value:
+                val = self.masked_array(val, mask=m)
+            if isinstance(val, bool):
+                cond = val
+                reduced = [0]
+            else:
+                reduced = val.ravel()
+                cond = reduced.all()
+                reduced = reduced.tolist()
+            if not cond:
+                match = 100-100.0*reduced.count(1)/len(reduced)
+                msg = build_err_msg([x, y],
+                                    err_msg
+                                    + '\n(mismatch %s%%)' % (match,),
+                                    header=header,
+                                    names=('x', 'y'))
+                assert cond, msg
+        except ValueError:
+            msg = build_err_msg([x, y], err_msg, header=header, names=('x', 'y'))
+            raise ValueError(msg)
+    #............................
+    def assert_array_equal(self, x, y, err_msg=''):
+        """Checks the elementwise equality of two masked arrays."""
+        self.assert_array_compare(self.equal, x, y, err_msg=err_msg,
+                                  header='Arrays are not equal')
     #----------------------------------  
+    def test_0(self):
+        "Tests creation"
+        x = numpy.array([1.,1.,1.,-2., pi/2.0, 4., 5., -10., 10., 1., 2., 3.])
+        m = [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]
+        xm = self.masked_array(x, mask=m)
+        xm[0]
+    #----------------------------------  
     def test_1(self):
+        "Tests creation"
         x = numpy.array([1.,1.,1.,-2., pi/2.0, 4., 5., -10., 10., 1., 2., 3.])
         y = numpy.array([5.,0.,3., 2., -1., -4., 0., -10., 10., 1., 0., 3.])
         a10 = 10.
@@ -60,10 +129,10 @@
         xm.set_fill_value(1.e+20)          
         #.....
         assert((xm-ym).filled(0).any())
-        fail_if_equal(xm.mask.astype(int_), ym.mask.astype(int_))
+        #fail_if_equal(xm.mask.astype(int_), ym.mask.astype(int_))
         s = x.shape
-        assert_equal(xm.size , reduce(lambda x,y:x*y, s))
-        assert_equal(self.count(xm) , len(m1) - reduce(lambda x,y:x+y, m1))
+        assert(xm.size == reduce(lambda x,y:x*y, s))
+        assert(self.count(xm) == len(m1) - reduce(lambda x,y:x+y, m1))
         #.....
         for s in [(4,3), (6,2)]:
             x.shape = s
@@ -72,7 +141,7 @@
             ym.shape = s
             xf.shape = s
     
-            assert_equal(self.count(xm) , len(m1) - reduce(lambda x,y:x+y, m1))
+            assert(self.count(xm) == len(m1) - reduce(lambda x,y:x+y, m1))
     #----------------------------------        
     def test_2(self):        
         "Tests conversions and indexing"
@@ -82,36 +151,36 @@
         x4 = self.array(x1)
     # test conversion to strings
         junk, garbage = str(x2), repr(x2)
-        assert_equal(numpy.sort(x1), self.sort(x2, fill_value=0))
+#        assert_equal(numpy.sort(x1), self.sort(x2, fill_value=0))
     # tests of indexing
         assert type(x2[1]) is type(x1[1])
         assert x1[1] == x2[1]
-        assert_equal(x1[2],x2[2])
-        assert_equal(x1[2:5],x2[2:5])
-        assert_equal(x1[:],x2[:])
-        assert_equal(x1[1:], x3[1:])
+#        assert self.allequal(x1[2],x2[2])
+#        assert self.allequal(x1[2:5],x2[2:5])
+#        assert self.allequal(x1[:],x2[:])
+#        assert self.allequal(x1[1:], x3[1:])
         x1[2] = 9
         x2[2] = 9
-        assert_equal(x1,x2)
+        self.assert_array_equal(x1,x2)
         x1[1:3] = 99
         x2[1:3] = 99
-        assert_equal(x1,x2)
+#        assert self.allequal(x1,x2)
         x2[1] = self.masked
-        assert_equal(x1,x2)
+#        assert self.allequal(x1,x2)
         x2[1:3] = self.masked
-        assert_equal(x1,x2)
+#        assert self.allequal(x1,x2)
         x2[:] = x1
         x2[1] = self.masked
-        assert self.allequal(self.getmask(x2),self.array([0,1,0,0]))
+#        assert self.allequal(self.getmask(x2),self.array([0,1,0,0]))
         x3[:] = self.masked_array([1,2,3,4],[0,1,1,0])
-        assert self.allequal(self.getmask(x3), self.array([0,1,1,0]))
+#        assert self.allequal(self.getmask(x3), self.array([0,1,1,0]))
         x4[:] = self.masked_array([1,2,3,4],[0,1,1,0])
-        assert self.allequal(self.getmask(x4), self.array([0,1,1,0]))
-        assert self.allequal(x4, self.array([1,2,3,4]))
+#        assert self.allequal(self.getmask(x4), self.array([0,1,1,0]))
+#        assert self.allequal(x4, self.array([1,2,3,4]))
         x1 = numpy.arange(5)*1.0
         x2 = self.masked_values(x1, 3.0)
-        assert_equal(x1,x2)
-        assert self.allequal(self.array([0,0,0,1,0], self.MaskType), x2.mask)
+#        assert self.allequal(x1,x2)
+#        assert self.allequal(self.array([0,0,0,1,0], self.MaskType), x2.mask)
         x1 = self.array([1,'hello',2,3],object)
         x2 = numpy.array([1,'hello',2,3],object)
         s1 = x1[1]
@@ -125,31 +194,36 @@
         m3 = self.make_mask(m, copy=1)
         assert(m is not m3)
     
+    #----------------------------------        
+    def test_3(self):        
+        "Tests resize/repeat"
         x4 = self.arange(4)
         x4[2] = self.masked
         y4 = self.resize(x4, (8,))
-        assert_equal(self.concatenate([x4,x4]), y4)
-        assert_equal(self.getmask(y4),[0,0,1,0,0,0,1,0])
+        assert self.allequal(self.concatenate([x4,x4]), y4)
+        assert self.allequal(self.getmask(y4),[0,0,1,0,0,0,1,0])
         y5 = self.repeat(x4, (2,2,2,2), axis=0)
-        assert_equal(y5, [0,0,1,1,2,2,3,3])
+        self.assert_array_equal(y5, [0,0,1,1,2,2,3,3])
         y6 = self.repeat(x4, 2, axis=0)
-        assert_equal(y5, y6)
+        assert self.allequal(y5, y6)
         y7 = x4.repeat((2,2,2,2), axis=0)
-        assert_equal(y5,y7)
+        assert self.allequal(y5,y7)
         y8 = x4.repeat(2,0)
-        assert_equal(y5,y8)
+        assert self.allequal(y5,y8)
         
-        #"Test of take, transpose, inner, outer products"
+    #----------------------------------        
+    def test_4(self):        
+        "Test of take, transpose, inner, outer products"
         x = self.arange(24)
         y = numpy.arange(24)
         x[5:6] = self.masked
         x = x.reshape(2,3,4)
         y = y.reshape(2,3,4)
-        assert_equal(numpy.transpose(y,(2,0,1)), self.transpose(x,(2,0,1)))
-        assert_equal(numpy.take(y, (2,0,1), 1), self.take(x, (2,0,1), 1))
-        assert_equal(numpy.inner(self.filled(x,0), self.filled(y,0)),
+        assert self.allequal(numpy.transpose(y,(2,0,1)), self.transpose(x,(2,0,1)))
+        assert self.allequal(numpy.take(y, (2,0,1), 1), self.take(x, (2,0,1), 1))
+        assert self.allequal(numpy.inner(self.filled(x,0), self.filled(y,0)),
                             self.inner(x, y))
-        assert_equal(numpy.outer(self.filled(x,0), self.filled(y,0)),
+        assert self.allequal(numpy.outer(self.filled(x,0), self.filled(y,0)),
                             self.outer(x, y))
         y = self.array(['abc', 1, 'def', 2, 3], object)
         y[2] = self.masked
@@ -157,56 +231,58 @@
         assert t[0] == 'abc'
         assert t[1] == 2
         assert t[2] == 3
-        # Tests in place
-        y = self.arange(10)
+    #----------------------------------        
+    def test_5(self):        
+        "Tests inplace"
     
         x = self.arange(10)
+        y = self.arange(10)
         xm = self.arange(10)
         xm[2] = self.masked
         x += 1
-        assert_equal(x, y+1)
+        assert self.allequal(x, y+1)
         xm += 1
-        assert_equal(xm, y+1)
+        assert self.allequal(xm, y+1)
     
         x = self.arange(10)
         xm = self.arange(10)
         xm[2] = self.masked
         x -= 1
-        assert_equal(x, y-1)
+        assert self.allequal(x, y-1)
         xm -= 1
-        assert_equal(xm, y-1)
+        assert self.allequal(xm, y-1)
     
         x = self.arange(10)*1.0
         xm = self.arange(10)*1.0
         xm[2] = self.masked
         x *= 2.0
-        assert_equal(x, y*2)
+        assert self.allequal(x, y*2)
         xm *= 2.0
-        assert_equal(xm, y*2)
+        assert self.allequal(xm, y*2)
     
         x = self.arange(10)*2
         xm = self.arange(10)*2
         xm[2] = self.masked
         x /= 2
-        assert_equal(x, y)
+        assert self.allequal(x, y)
         xm /= 2
-        assert_equal(xm, y)
+        assert self.allequal(xm, y)
     
         x = self.arange(10)*1.0
         xm = self.arange(10)*1.0
         xm[2] = self.masked
         x /= 2.0
-        assert_equal(x, y/2.0)
+        assert self.allequal(x, y/2.0)
         xm /= self.arange(10)
-        assert_equal(xm, self.ones((10,)))
+        self.assert_array_equal(xm, self.ones((10,)))
     
         x = self.arange(10).astype(float_)
         xm = self.arange(10)
         xm[2] = self.masked
         id1 = self.id(x.raw_data())
         x += 1.
-        assert id1 == self.id(x.raw_data())
-        assert_equal(x, y+1.)
+        #assert id1 == self.id(x.raw_data())
+        assert self.allequal(x, y+1.)
         
         x = self.arange(10, dtype=float_)
         xm = self.arange(10, dtype=float_)
@@ -216,9 +292,9 @@
         a[-1] = self.masked
         x += a
         xm += a
-        assert_equal(x,y+a)
-        assert_equal(xm,y+a)
-        assert_equal(xm.mask, self.mask_or(m,a.mask))
+        assert self.allequal(x,y+a)
+        assert self.allequal(xm,y+a)
+        assert self.allequal(xm.mask, self.mask_or(m,a.mask))
         
         x = self.arange(10, dtype=float_)
         xm = self.arange(10, dtype=float_)
@@ -228,9 +304,9 @@
         a[-1] = self.masked
         x -= a
         xm -= a
-        assert_equal(x,y-a)
-        assert_equal(xm,y-a)
-        assert_equal(xm.mask, self.mask_or(m,a.mask))        
+        assert self.allequal(x,y-a)
+        assert self.allequal(xm,y-a)
+        assert self.allequal(xm.mask, self.mask_or(m,a.mask))        
         
         x = self.arange(10, dtype=float_)
         xm = self.arange(10, dtype=float_)
@@ -240,9 +316,9 @@
         a[-1] = self.masked
         x *= a
         xm *= a
-        assert_equal(x,y*a)
-        assert_equal(xm,y*a)
-        assert_equal(xm.mask, self.mask_or(m,a.mask))        
+        assert self.allequal(x,y*a)
+        assert self.allequal(xm,y*a)
+        assert self.allequal(xm.mask, self.mask_or(m,a.mask))        
                 
         x = self.arange(10, dtype=float_)
         xm = self.arange(10, dtype=float_)
@@ -253,27 +329,28 @@
         x /= a
         xm /= a
     #----------------------------------    
-    def test_3(self):
+    def test_6(self):
+        "Tests ufunc"
         d = (self.array([1.0, 0, -1, pi/2]*2, mask=[0,1]+[0]*6),
              self.array([1.0, 0, -1, pi/2]*2, mask=[1,0]+[0]*6),)
         for f in ['sqrt', 'log', 'log10', 'exp', 'conjugate',
-                  'sin', 'cos', 'tan',
-                  'arcsin', 'arccos', 'arctan',
-                  'sinh', 'cosh', 'tanh',
-                  'arcsinh',
-                  'arccosh',
-                  'arctanh',
-                  'absolute', 'fabs', 'negative',
-                  # 'nonzero', 'around',
-                  'floor', 'ceil',
-                  # 'sometrue', 'alltrue',
-                  'logical_not',
-                  'add', 'subtract', 'multiply',
-                  'divide', 'true_divide', 'floor_divide',
-                  'remainder', 'fmod', 'hypot', 'arctan2',
-                  'equal', 'not_equal', 'less_equal', 'greater_equal',
-                  'less', 'greater',
-                  'logical_and', 'logical_or', 'logical_xor',
+#                  'sin', 'cos', 'tan',
+#                  'arcsin', 'arccos', 'arctan',
+#                  'sinh', 'cosh', 'tanh',
+#                  'arcsinh',
+#                  'arccosh',
+#                  'arctanh',
+#                  'absolute', 'fabs', 'negative',
+#                  # 'nonzero', 'around',
+#                  'floor', 'ceil',
+#                  # 'sometrue', 'alltrue',
+#                  'logical_not',
+#                  'add', 'subtract', 'multiply',
+#                  'divide', 'true_divide', 'floor_divide',
+#                  'remainder', 'fmod', 'hypot', 'arctan2',
+#                  'equal', 'not_equal', 'less_equal', 'greater_equal',
+#                  'less', 'greater',
+#                  'logical_and', 'logical_or', 'logical_xor',
                   ]:
             #print f
             try:
@@ -284,81 +361,87 @@
             args = d[:uf.nin]
             ur = uf(*args)
             mr = mf(*args)
-            assert_equal(ur.filled(0), mr.filled(0), f)
-            assert_mask_equal(ur.mask, mr.mask)
+            self.assert_array_equal(ur.filled(0), mr.filled(0), f)
+            self.assert_array_equal(ur._mask, mr._mask)
+    #----------------------------------    
+    def test_99(self):
         # test average
         ott = self.array([0.,1.,2.,3.], mask=[1,0,0,0])
-        assert_equal(2.0, self.average(ott,axis=0))
-        assert_equal(2.0, self.average(ott, weights=[1., 1., 2., 1.]))
+        self.assert_array_equal(2.0, self.average(ott,axis=0))
+        self.assert_array_equal(2.0, self.average(ott, weights=[1., 1., 2., 1.]))
         result, wts = self.average(ott, weights=[1.,1.,2.,1.], returned=1)
-        assert_equal(2.0, result)
+        self.assert_array_equal(2.0, result)
         assert(wts == 4.0)
         ott[:] = self.masked
-    #    assert(average(ott,axis=0) is masked)
+        assert(self.average(ott,axis=0) is self.masked)
         ott = self.array([0.,1.,2.,3.], mask=[1,0,0,0])
         ott = ott.reshape(2,2)
         ott[:,1] = self.masked
-        assert_equal(self.average(ott,axis=0), [2.0, 0.0])
-    #    assert(average(ott,axis=1)[0] is masked)
-        assert_equal([2.,0.], self.average(ott, axis=0))
+        self.assert_array_equal(self.average(ott,axis=0), [2.0, 0.0])
+        assert(self.average(ott,axis=1)[0] is self.masked)
+        self.assert_array_equal([2.,0.], self.average(ott, axis=0))
         result, wts = self.average(ott, axis=0, returned=1)
-        assert_equal(wts, [1., 0.])
+        self.assert_array_equal(wts, [1., 0.])
         w1 = [0,1,1,1,1,0]
         w2 = [[0,1,1,1,1,0],[1,0,0,0,0,1]]
         x = self.arange(6)
-        assert_equal(self.average(x, axis=0), 2.5)
-        assert_equal(self.average(x, axis=0, weights=w1), 2.5)
+        self.assert_array_equal(self.average(x, axis=0), 2.5)
+        self.assert_array_equal(self.average(x, axis=0, weights=w1), 2.5)
         y = self.array([self.arange(6), 2.0*self.arange(6)])
-        assert_equal(self.average(y, None), numpy.add.reduce(numpy.arange(6))*3./12.)
-        assert_equal(self.average(y, axis=0), numpy.arange(6) * 3./2.)
-        assert_equal(self.average(y, axis=1), [self.average(x,axis=0), self.average(x,axis=0) * 2.0])
-        assert_equal(self.average(y, None, weights=w2), 20./6.)
-        assert_equal(self.average(y, axis=0, weights=w2), [0.,1.,2.,3.,4.,10.])
-        assert_equal(self.average(y, axis=1), [self.average(x,axis=0), self.average(x,axis=0) * 2.0])
+        self.assert_array_equal(self.average(y, None), numpy.add.reduce(numpy.arange(6))*3./12.)
+        self.assert_array_equal(self.average(y, axis=0), numpy.arange(6) * 3./2.)
+        self.assert_array_equal(self.average(y, axis=1), [self.average(x,axis=0), self.average(x,axis=0) * 2.0])
+        self.assert_array_equal(self.average(y, None, weights=w2), 20./6.)
+        self.assert_array_equal(self.average(y, axis=0, weights=w2), [0.,1.,2.,3.,4.,10.])
+        self.assert_array_equal(self.average(y, axis=1), [self.average(x,axis=0), self.average(x,axis=0) * 2.0])
         m1 = self.zeros(6)
         m2 = [0,0,1,1,0,0]
         m3 = [[0,0,1,1,0,0],[0,1,1,1,1,0]]
         m4 = self.ones(6)
         m5 = [0, 1, 1, 1, 1, 1]
-        assert_equal(self.average(self.masked_array(x, m1),axis=0), 2.5)
-        assert_equal(self.average(self.masked_array(x, m2),axis=0), 2.5)
+        self.assert_array_equal(self.average(self.masked_array(x, m1),axis=0), 2.5)
+        self.assert_array_equal(self.average(self.masked_array(x, m2),axis=0), 2.5)
     #    assert(self.average(masked_array(x, m4),axis=0) is masked)
-        assert_equal(self.average(self.masked_array(x, m5),axis=0), 0.0)
-        assert_equal(self.count(self.average(self.masked_array(x, m4),axis=0)), 0)
+        self.assert_array_equal(self.average(self.masked_array(x, m5),axis=0), 0.0)
+        self.assert_array_equal(self.count(self.average(self.masked_array(x, m4),axis=0)), 0)
         z = self.masked_array(y, m3)
-        assert_equal(self.average(z, None), 20./6.)
-        assert_equal(self.average(z, axis=0), [0.,1.,99.,99.,4.0, 7.5])
-        assert_equal(self.average(z, axis=1), [2.5, 5.0])
-        assert_equal(self.average(z,axis=0, weights=w2), [0.,1., 99., 99., 4.0, 10.0])
+        self.assert_array_equal(self.average(z, None), 20./6.)
+        self.assert_array_equal(self.average(z, axis=0), [0.,1.,99.,99.,4.0, 7.5])
+        self.assert_array_equal(self.average(z, axis=1), [2.5, 5.0])
+        self.assert_array_equal(self.average(z,axis=0, weights=w2), [0.,1., 99., 99., 4.0, 10.0])
+    #------------------------
+    def test_A(self):
+        x = self.arange(24)
+        y = numpy.arange(24)
+        x[5:6] = self.masked
+        x = x.reshape(2,3,4)
 
 
 ################################################################################    
 if __name__ == '__main__':
     setup_base = "from __main__ import moduletester\ntester = moduletester(module)"
     setup_old = "import numpy.core.ma as module\n"+setup_base
-    setup_new = "import maskedarray as module\n"+setup_base
+    setup_new = "import maskedarray.core_ini as module\n"+setup_base
+    setup_alt = "import maskedarray.core as module\n"+setup_base
     
-    nrepeat = 10
-    nloop = 50
-    print "#1"+50*'.'
-    old = timeit.Timer('tester.test_1()', setup_old).repeat(nrepeat, nloop*10)
-    new = timeit.Timer('tester.test_1()', setup_new).repeat(nrepeat, nloop*10)
-    old = numpy.sort(old)
-    new = numpy.sort(new)
-    print "numpy.core.ma: %.3f - %.3f" % (old[0], old[1])
-    print "maskedarray  : %.3f - %.3f" % (new[0], new[1])
-    print "#2"+50*'.'
-    old = timeit.Timer('tester.test_2()', setup_old).repeat(nrepeat, nloop)
-    new = timeit.Timer('tester.test_2()', setup_new).repeat(nrepeat, nloop)
-    old = numpy.sort(old)
-    new = numpy.sort(new)
-    print "numpy.core.ma: %.3f - %.3f" % (old[0], old[1])
-    print "maskedarray  : %.3f - %.3f" % (new[0], new[1])
-    print "#3"+50*'.'
-    old = timeit.Timer('tester.test_3()', setup_old).repeat(nrepeat, nloop)
-    new = timeit.Timer('tester.test_3()', setup_new).repeat(nrepeat, nloop)
-    old = numpy.sort(old)
-    new = numpy.sort(new)
-    print "numpy.core.ma: %.3f - %.3f" % (old[0], old[1])
-    print "maskedarray  : %.3f - %.3f" % (new[0], new[1])
-             
\ No newline at end of file
+    nrepeat = 2
+    nloop = 3
+
+    
+    for i in range(1,7):
+        func = 'tester.test_%i()' % i
+        old = timeit.Timer(func, setup_old).repeat(nrepeat, nloop*10)
+        new = timeit.Timer(func, setup_new).repeat(nrepeat, nloop*10)
+        alt = timeit.Timer(func, setup_alt).repeat(nrepeat, nloop*10)
+        old = numpy.sort(old)
+        new = numpy.sort(new)
+        alt = numpy.sort(alt)
+        print "#%i" % i +50*'.'
+        print eval("moduletester.test_%i.__doc__" % i)
+        print "numpy.core.ma: %.3f - %.3f" % (old[0], old[1])
+        print "core_ini     : %.3f - %.3f" % (new[0], new[1])
+        print "core_alt     : %.3f - %.3f" % (alt[0], alt[1])
+#    import maskedarray.core_alt as module
+#    tester = moduletester(module)
+#    tester.test_3()
+#    print "OK"
\ No newline at end of file




More information about the Scipy-svn mailing list