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

scipy-svn at scipy.org scipy-svn at scipy.org
Fri Dec 22 19:00:49 EST 2006


Author: pierregm
Date: 2006-12-22 18:00:43 -0600 (Fri, 22 Dec 2006)
New Revision: 2464

Added:
   trunk/Lib/sandbox/maskedarray/timer_comparison.py
Modified:
   trunk/Lib/sandbox/maskedarray/CHANGELOG
   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
Log:
cf changelog

Modified: trunk/Lib/sandbox/maskedarray/CHANGELOG
===================================================================
--- trunk/Lib/sandbox/maskedarray/CHANGELOG	2006-12-22 14:39:32 UTC (rev 2463)
+++ trunk/Lib/sandbox/maskedarray/CHANGELOG	2006-12-23 00:00:43 UTC (rev 2464)
@@ -1,3 +1,13 @@
+#2006-12-22 : Core
+#           : - Optimized(?) default_/maximum_/minimum_fill_value
+#           : - Force __new__ to not return a MaskedArray, in order to ...
+#           :   ... optimize __array_finalize__
+#           : - Add the hard_mask flag to __new__ (*[False]*)
+#2006-12-19 : Core
+#			: - Fixed a problem on _set_mask which prevented to set a mask to nomask
+#           : Extras
+#           : - Renamed compress2d to compress_rowcols
+#           : - Added dot
 #2006-12-18 : Extras
 #           : - Added compress2d and mask_rowcols
 #2006-12-13 : - moved 'average' to 'extras'

Modified: trunk/Lib/sandbox/maskedarray/core.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/core.py	2006-12-22 14:39:32 UTC (rev 2463)
+++ trunk/Lib/sandbox/maskedarray/core.py	2006-12-23 00:00:43 UTC (rev 2464)
@@ -55,13 +55,12 @@
            'var', 'where',
            'zeros']
 
-
 import sys
 import types
 import cPickle
 #
 import numpy
-from numpy import bool_, complex_, float_, int_, object_
+from numpy import bool_, complex_, float_, int_, object_, str_
 
 import numpy.core.umath as umath 
 import numpy.core.fromnumeric  as fromnumeric
@@ -73,11 +72,16 @@
 from numpy.lib.shape_base import expand_dims as n_expand_dims
 import warnings
 
+import logging
+logging.basicConfig(level=logging.WARNING,
+                    format='%(levelname)s %(message)s',)
 
+
 MaskType = bool_
 nomask = MaskType(0)
 
 divide_tolerance = 1.e-35
+numpy.seterr(all='ignore')
 
 #####--------------------------------------------------------------------------
 #---- --- Helper functions ---
@@ -128,70 +132,67 @@
 #####--------------------------------------------------------------------------
 #---- --- Filling options ---
 #####--------------------------------------------------------------------------
-# Use single element arrays or scalars.
-default_real_fill_value = 1.e20
-default_complex_fill_value = 1.e20 + 0.0j
-default_character_fill_value = '-'
-default_integer_fill_value = 999999
-default_object_fill_value = '?'
+# 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',
+                  }
+max_filler = {'b': False,
+              'f' : -numeric.inf,
+              'i' : -sys.maxint,
+              }
+min_filler = {'b' : True,
+              'f' : numeric.inf,
+              'i' : sys.maxint,
+              }
 
+
 def default_fill_value (obj):
     "Calculates the default fill value for an object `obj`."
-    if isinstance(obj, types.FloatType):
-        return default_real_fill_value
-    elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType):
-        return default_integer_fill_value
-    elif isinstance(obj, types.StringType):
-        return default_character_fill_value
-    elif isinstance(obj, types.ComplexType):
-        return default_complex_fill_value
-    elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray):
-        x = obj.dtype.char
-        if x in typecodes['Float']:
-            return default_real_fill_value
-        if x in typecodes['Integer']:
-            return default_integer_fill_value
-        if x in typecodes['Complex']:
-            return default_complex_fill_value
-        if x in typecodes['Character']:
-            return default_character_fill_value
-        if x in typecodes['UnsignedInteger']:
-            return umath.absolute(default_integer_fill_value)
-        return default_object_fill_value
+    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']
     else:
-        return default_object_fill_value
+        return default_filler['O']
 
 def minimum_fill_value (obj):
     "Calculates the default fill value suitable for taking the minimum of `obj`."
-    if isinstance(obj, types.FloatType):
-        return numeric.inf
-    elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType):
-        return sys.maxint
-    elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray):
-        x = obj.dtype.char
-        if x in typecodes['Float']:
-            return numeric.inf
-        if x in typecodes['Integer']:
-            return sys.maxint
-        if x in typecodes['UnsignedInteger']:
-            return sys.maxint
+    if hasattr(obj, 'dtype'):
+        objtype = obj.dtype.kind
+        try:
+            return min_filler[objtype]
+        except KeyError:
+            raise TypeError, 'Unsuitable type for calculating minimum.'
+    elif isinstance(obj, float):
+        return min_filler['f']
+    elif isinstance(obj, int) or isinstance(obj, long):
+        return min_filler['i']
     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 isinstance(obj, types.FloatType):
-        return -numeric.inf
-    elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType):
-        return -sys.maxint
-    elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray):
-        x = obj.dtype.char
-        if x in typecodes['Float']:
-            return -numeric.inf
-        if x in typecodes['Integer']:
-            return -sys.maxint
-        if x in typecodes['UnsignedInteger']:
-            return 0
+    if hasattr(obj, 'dtype'):
+        objtype = obj.dtype.kind
+        try:
+            return max_filler[objtype]
+        except KeyError:
+            raise TypeError, 'Unsuitable type for calculating maximum.'
+    elif isinstance(obj, float):
+        return max_filler['f']
+    elif isinstance(obj, int) or isinstance(obj, long):
+        #TODO: Check what happens to 'UnisgnedInteger'!
+        return max_filler['i']
     else:
         raise TypeError, 'Unsuitable type for calculating maximum.'
 
@@ -567,17 +568,17 @@
     except AttributeError:
         return False
 #
-def make_mask(m, copy=False, flag=False):
-    """make_mask(m, copy=0, flag=0)
+def make_mask(m, copy=False, small_mask=False):
+    """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 `flag=True`, returns `nomask` if `m` contains no true elements.
+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.
-    - `flag` (boolean, *[False]*): Flattens mask to `nomask` if `m` is all false.
+    - `small_mask` (boolean, *[False]*): Flattens mask to `nomask` if `m` is all false.
     """
     if m is nomask:
         return nomask
@@ -593,7 +594,7 @@
     else:
         result = numeric.array(filled(m, True), dtype=MaskType)
     # Bas les masques !
-    if flag and not result.any():
+    if small_mask and not result.any():
         return nomask
     else:
         return result
@@ -603,7 +604,7 @@
     result = numeric.zeros(s, dtype=MaskType)
     return result
 
-def mask_or (m1, m2, copy=False, flag=True):
+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.
@@ -611,15 +612,15 @@
 :Parameters:
     - `m` (ndarray) : Mask.
     - `copy` (boolean, *[False]*) : Returns a copy of `m` if true.
-    - `flag` (boolean, *[False]*): Flattens mask to `nomask` if `m` is all false.
+    - `small_mask` (boolean, *[False]*): Flattens mask to `nomask` if `m` is all false.
      """
     if m1 is nomask: 
-        return make_mask(m2, copy=copy, flag=flag)
+        return make_mask(m2, copy=copy, small_mask=small_mask)
     if m2 is nomask: 
-        return make_mask(m1, copy=copy, flag=flag)
+        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, flag=flag)
+    return make_mask(umath.logical_or(m1, m2), copy=copy, small_mask=small_mask)
 
 #####--------------------------------------------------------------------------
 #--- --- Masking functions ---
@@ -707,7 +708,7 @@
     else:
         condition = umath.equal(fromnumeric.asarray(x), value)
         mask = nomask
-    mask = mask_or(mask, make_mask(condition, flag=True))
+    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):
@@ -732,7 +733,7 @@
     else:
         condition = umath.equal(xnew, value)
         mask = nomask
-    mask = mask_or(mask, make_mask(condition, flag=True))
+    mask = mask_or(mask, make_mask(condition, small_mask=True))
     return masked_array(xnew, mask=mask, copy=copy, fill_value=value)
 
 #####--------------------------------------------------------------------------
@@ -757,9 +758,9 @@
         "Is the use of the display value enabled?"
         return self._enabled
 
-    def enable(self, flag=1):
-        "Set the enabling flag to `flag`."
-        self._enabled = flag
+    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)
@@ -778,7 +779,7 @@
 
 Construction:
     x = array(data, dtype=None, copy=True, order=False,
-              mask = nomask, fill_value=None, flag=True)
+              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 
@@ -792,7 +793,7 @@
 
 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 `flag` is True, a mask consisting of zeros (False) only is compressed to `nomask`.
+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 
@@ -800,108 +801,95 @@
 The fill_value is not used for computation within this module.
     """
     __array_priority__ = 10.1
-    
+    #TODO: There some reorganization to do round here
     def __new__(cls, data, mask=nomask, dtype=None, copy=False, fill_value=None,
-                flag=True, keep_mask=True):
+                keep_mask=True, small_mask=True, hard_mask=False):
         """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 dtype is not None:
-            dtype = numeric.dtype(dtype)
+#        logging.debug("__new__ received %s" % type(data))
         # 1. Argument is MA ...........
         if isinstance(data, MaskedArray) or\
            (hasattr(data,"_mask") and hasattr(data,"_data")) :
             if keep_mask:
                 if mask is nomask:
-                    cls._basemask = data._mask
+                    cls.__defaultmask = data._mask
                 else:
-                    cls._basemask = mask_or(data._mask, mask)
+                    cls.__defaultmask = mask_or(data._mask, mask, 
+                                                copy=copy, small_mask=small_mask)
             else:
-                # Force copy of mask if it changes
-                cls._basemask = make_mask(mask, copy=copy, flag=flag)
+                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
-            return numeric.array(data._data, dtype=dtype, copy=copy).view(cls)
+            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 != dtype:
+            if dtype is not None and data.dtype != numeric.dtype(dtype):
                 _data = data.astype(dtype)
             elif copy:
                 _data = data.copy()
             else:
                 _data = data
         else:
-            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._basemask = getmask(_data)
-                    return _data.view(cls)
+            _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=False, flag=flag)
+        if not is_mask(mask):
+            mask = make_mask(mask, small_mask=small_mask)
         #....Check shapes compatibility
-        if _mask is not nomask:
-            (nd, nm) = (_data.size, _mask.size)
+        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)
+                    mask = fromnumeric.resize(mask, _data.shape)
                 elif nd == 1:
-                    _data = fromnumeric.resize(_data, _mask.shape)
+                    _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).copy()
+            elif (mask.shape != _data.shape):
+                mask = mask.reshape(_data.shape)
+#                mask = _data.shape
         #....
         cls._fill_value = fill_value
-        cls._basemask = _mask
+        cls.__defaulthardmask = hard_mask
+        cls.__defaultmask = mask
+#        logging.debug("__new__ returned %s as %s" % (type(_data), cls))
         return numeric.asanyarray(_data).view(cls)
     #..................................
-    def __array__ (self, t=None, context=None):
-        "Special hook for numeric. Converts to numeric if possible."
-        # Er... Do we really need __array__ ?
-        if self._mask is not nomask:
-            if fromnumeric.ravel(self._mask).any():
-                if context is None:
-                    # Hardliner stand: raise an exception
-                    # We may wanna use warnings.warn instead
-                    raise MAError,\
-                          "Cannot automatically convert masked array to "\
-                          "numeric because data\n    is masked in one or "\
-                          "more locations."
-                    #return self._data
-                else:
-                    func, args, i = context
-                    fills = ufunc_fills.get(func)
-                    if fills is None:
-                        raise MAError, "%s not known to ma" % func
-                    return self.filled(fills[i])
-            else:  # Mask is all false
-                   # Optimize to avoid future invocations of this section.
-                self._mask = nomask
-                self._shared_mask = 0
-        if t:
-            return self._data.astype(t)
-        else:
-            return self._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__
+#        mclass = self.__class__
         #..........
+#        logging.debug("__wrap__ received %s" % type(obj))
         if context is None:
-            print "DEBUG _wrap_: no context"
-            return mclass(obj, mask=self._mask, copy=False)
+#            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])
@@ -918,83 +906,122 @@
             else:
                 if m.shape != dshape:
                     m = reduce(mask_or, [getmaskarray(arg) for arg in args])
-        return mclass(obj, copy=False, mask=m)
+#        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 not hasattr(self, "_data"):
-            try:
-                self._data = obj._data
-            except AttributeError:
-                self._data = obj
-        # 
-        self.fill_value = self._fill_value
-        #
-        if not hasattr(self, '_mask'):
-            self._mask = self._basemask
-        #
+#        logging.debug("__finalize__ received %s" % type(obj))
+        if isMaskedArray(obj):
+            self._data = obj._data
+            self._mask = obj._mask
+            self._hardmask = obj._hardmask
+            self._fill_value = obj._fill_value
+        else:
+            self._data = obj
+            self._mask = self.__defaultmask
+            self._hardmask = self.__defaulthardmask
+            self.fill_value = self._fill_value
+#        #
+#        logging.debug("__finalize__ returned %s" % type(self))       
         return    
     #............................................
-    def __getitem__(self, i):
+    def __getitem__(self, index):
         """x.__getitem__(y) <==> x[y]
 Returns the item described by i. Not a copy as in previous versions.
         """
-        dout = self._data[i]
-        if self._mask is nomask:
-            if numeric.size(dout)==1:
+        if getmask(index) is not nomask:
+            msg = "Masked arrays must be filled before they can be used as indices!"
+            raise IndexError, msg
+        dout = self._data[index]
+        m = self._mask
+        scalardout = (len(numeric.shape(dout))==0)
+        # 
+        if m is nomask:
+            if scalardout:
                 return dout
             else:
-                return self.__class__(dout, mask=nomask, 
-                                      fill_value=self._fill_value,
-                                      dtype = self.dtype,)
+                return self.__class__(dout, mask=nomask, keep_mask=True,
+                                      fill_value=self._fill_value)
         #....
-#        m = self._mask.copy()   
-        m = self._mask
-        mi = m[i]
+        mi = m[index]
         if mi.size == 1:
             if mi:
                 return masked
-            else:
-                return dout
+            return dout
         else:
-            return self.__class__(dout, mask=mi, dtype = self.dtype,
-                                  fill_value=self._fill_value)
+            return self.__class__(dout, mask=mi, fill_value=self._fill_value)
     #........................
     def __setitem__(self, index, value):
         """x.__setitem__(i, y) <==> x[i]=y
 Sets item described by index. If value is masked, masks those locations.
         """
-        d = self._data
         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 self._mask is nomask:
-                _mask = make_mask_none(d.shape)
+            if m is nomask:
+                m = make_mask_none(d.shape)
             else:
-                _mask = self._mask.copy()
-            _mask[index] = True
-            self._mask = _mask
+                m = m.copy()
+            m[index] = True
+            self._mask = m
             return
-        #....
-        m = getmask(value)
-        value = filled(value).astype(d.dtype)
-        d[index] = value
+        #....   
         if m is nomask:
-            if self._mask is not nomask:
-                _mask = self._mask.copy()
-                _mask[index] = False
+            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:
-                _mask = nomask
+                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
+            else:
+                dindx = d[index]
+                numeric.putmask(dindx, ~valmask, valdata)
+                d[index] = dindx
+                numeric.putmask(mindx, valmask, True)
+                m[index] = mindx
+#        else:
+#            mindx = m[index]
+#            value = masked_array(value, mask=mindx, keep_mask=True)
+#            valmask = getmask(value)
+#            if valmask is nomask:
+#                d[vindx] = filled(value)
+#            else:
+#                d[vindx[~valmask]] = filled(value[~valmask])
+#                m[vindx[valmask]] = True
+        #.....    
+        if not m.any():
+            self._mask = nomask
         else:
-            if self._mask is nomask:
-                _mask = make_mask_none(d.shape)
-            else:
-                _mask = self._mask.copy()
-            _mask[index] = m
-        self._mask = _mask
+            self._mask = m
     #............................................
     def __getslice__(self, i, j):
         """x.__getslice__(i, j) <==> x[i:j]
@@ -1003,46 +1030,81 @@
         m = self._mask
         dout = self._data[i:j]
         if m is nomask:
-            return self.__class__(dout, dtype = self.dtype,
-                                  fill_value=self._fill_value)
+            return self.__class__(dout, fill_value=self._fill_value)
         else:
-            return self.__class__(dout, mask=m[i:j], dtype = self.dtype, 
-                                  fill_value=self._fill_value)
+            return self.__class__(dout, mask=m[i:j], fill_value=self._fill_value)
     #........................
     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."""
-        d = self._data
         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 self._mask is nomask:
-                _mask = make_mask_none(d.shape)
-            else:
-                _mask = self._mask.copy()
-            _mask[i:j] = True
-            self._mask = _mask
+            if m is nomask:
+                m = make_mask_none(d.shape)
+            m[i:j] = True
             return
         #....
-        m = getmask(value)
-        value = filled(value).astype(d.dtype)
-        d[i:j] = value
+#        valmask = getmask(value)
+#        valdata = filled(value).astype(d.dtype)
+#        if valmask is nomask:
+#            if m is nomask:
+#                d[i:j] = valdata
+#            elif self._hardmask:
+#                d[(~m)[i:j]] = valdata
+#            elif not self._hardmask:
+#                d[i:j] = valdata
+#                m[i:j] = False
+#                if not m.any():
+#                    self._mask = nomask
+#            return
+#        else:
+#            if m is nomask:
+#                m = make_mask_none(d.shape)
+#                m[i:j] = valmask
+#            elif self._hardmask:
+#                m[i:j] = mask_or(m[i:j], valmask)
+#            else:
+#                m[i:j] = valmask
+#            d[(~m)[i:j]] = valdata
+#            if not m.any():
+#                self._mask = nomask
+#            return
+        #....   
         if m is nomask:
-            if self._mask is not nomask:
-                _mask = self._mask.copy()
-                _mask[i:j] = False
+            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:
-                _mask = nomask
+                m[i:j] = valmask
         else:
-            if self._mask is nomask:
-                _mask = make_mask_none(d.shape)
+            mindx = m[i:j]
+            value = masked_array(value, mask=mindx, keep_mask=True)
+            valmask = getmask(value)
+            if valmask is None:
+                d[i:j][~mindx] = filled(value)
             else:
-                _mask = self._mask.copy()
-            _mask[i:j] = m
-        self._mask = make_mask(_mask, flag=True)
+                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):
@@ -1124,7 +1186,7 @@
 the absolute of the inital `_data`.
         """
         return self.__class__(self._data.__abs__(), mask=self._mask,
-                             dtype = self.dtype,)
+                              fill_value = self._fill_value,)
     #
     def __neg__(self):
         """x.__abs__() <==> neg(x)
@@ -1132,53 +1194,41 @@
 the negative of the inital `_data`."""
         try:
             return self.__class__(self._data.__neg__(), mask=self._mask,
-                                  dtype = self.dtype,)
+                                  fill_value = self._fill_value,)
         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._data += f
-            m = getmask(other)
             self._mask = m
-            ###self._shared_mask = m is not nomask
-        else: 
-            tmp = masked_array(f, mask=getmask(other))
-            self._data += tmp._data
-            self._mask = mask_or(self._mask, tmp._mask)
-            ###self._shared_mask = 1
+        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._data -= f
-            m = getmask(other)
             self._mask = m
-            ###self._shared_mask = m is not nomask
-        else:
-            tmp = masked_array(f, mask=getmask(other))
-            self._data -= tmp._data
-            self._mask = mask_or(self._mask, tmp._mask)
-            ###self._shared_mask = 1
+        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, 0), self._data.dtype.char) 
+        f = convert_typecode(filled(other, 1), self._data.dtype.char) 
+        m = getmask(other)
+        self._data *= f
         if self._mask is nomask:
-            self._data *= f
-            m = getmask(other)
             self._mask = m
-            ####self._shared_mask = m is not nomask
-        else:
-            tmp = masked_array(f, mask=getmask(other))
-            self._data *= tmp._data
-            self._mask = mask_or(self._mask, tmp._mask)
-            ###self._shared_mask = 1
+        elif m is not nomask:
+            self._mask += m
         return self
     #
     def __idiv__(self, other):
@@ -1215,7 +1265,7 @@
     def __float__(self):
         "Converts self to float."
         if self._mask is not nomask:
-            print "Warning: converting a masked element to nan."
+            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())
@@ -1236,10 +1286,15 @@
 Subclassing is preserved."""
         if tc == self._data.dtype:
             return self
-        d = self._data.astype(tc)
-#        print "DEBUG: _astype: d", d
-#        print "DEBUG: _astype: m", self._mask        
+        d = self._data.astype(tc) 
         return self.__class__(d, mask=self._mask, dtype=tc)
+    #............................................
+    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):
@@ -1247,10 +1302,10 @@
         """
         if self._mask is nomask:
             return masked_array(self._data.ravel(), mask=nomask, copy=False,
-                                fill_value = self.fill_value)
+                                fill_value = self._fill_value)
         else:
             return masked_array(self._data.ravel(), mask=self._mask.ravel(),
-                                copy=False, fill_value = self.fill_value)
+                                copy=False, fill_value = self._fill_value)
     #
     def _set_flat (self, value):
         "x.flat = value"
@@ -1263,7 +1318,7 @@
     def _get_real(self):
         "Returns the real part of a complex array."
         return masked_array(self._data.real, mask=self.mask,
-                            fill_value = self.fill_value)
+                            fill_value = self._fill_value)
 #        if self.mask is nomask:
 #            return masked_array(self._data.real, mask=nomask,
 #                            fill_value = self.fill_value)
@@ -1281,7 +1336,7 @@
     def _get_imaginary(self):
         "Returns the imaginary part of a complex array."
         return masked_array(self._data.imag, mask=nomask,
-                            fill_value = self.fill_value)
+                            fill_value = self._fill_value)
 
     def _set_imaginary (self, value):
         "Sets the imaginary part of a complex array to `value`."
@@ -1295,17 +1350,17 @@
     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, flag=True)
+        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
-    
+            self._mask = mask  
+        else:
+            self._mask = nomask
     mask = property(fget=_get_mask, fset=_set_mask, doc="Mask")
     #............................................
     def get_fill_value(self):
@@ -1332,11 +1387,10 @@
         d = self._data
         m = self._mask
         if m is nomask:
-#            return fromnumeric.asarray(d)
             return d
         #
         if fill_value is None:
-            value = self._fill_value
+            value = self.fill_value
         else:
             value = fill_value
         #
@@ -1367,11 +1421,8 @@
         d = self._data.ravel()
         if self._mask is nomask:
             return d
-#            return numeric.asarray(d)
         else:
-#            m = 1 - self._mask.ravel()
-#            return numeric.asarray(d.compress(m))
-            return d.compress(-self._mask.ravel())
+            return d[~self._mask.ravel()]
     #............................................
     def count(self, axis=None):
         """Counts the non-masked elements of the array along a given axis,
@@ -1433,11 +1484,9 @@
         if self._mask is not nomask:
             return self.__class__(self._data.reshape(*s), 
                                   mask=self._mask.reshape(*s),
-                                  dtype = self.dtype,
                                   fill_value=self.fill_value)
         else:
             return self.__class__(self._data.reshape(*s),
-                                  dtype = self.dtype,
                                   fill_value=self.fill_value)
     #
     def repeat(self, repeats, axis=None):
@@ -1458,25 +1507,18 @@
         if m is not nomask:
             m = fromnumeric.repeat(m, repeats, axis)
         d = fromnumeric.repeat(f, repeats, axis)
-        return self.__class__(d, mask=m, dtype = self.dtype,
-                              fill_value=self.fill_value)
+        return self.__class__(d, mask=m, fill_value=self.fill_value)
     #
     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.
         """
-        raiseit = False
         try:
             self._data.resize(newshape,)
+            if self.mask is not nomask:
+                self._mask.resize(newshape,)
         except ValueError:
-            raiseit = True
-        if self.mask is not nomask:
-            try:
-                self._mask.resize(newshape,)
-            except ValueError:
-                raiseit = True
-        if raiseit:
             msg = "Cannot resize an array that has been referenced or "+\
                   "is referencing another array in this way.\n"+\
                   "Use the resize function."
@@ -1511,7 +1553,7 @@
                 m.put(ind, values._mask, mode=mode)
             else:
                 m.put(ind, False, mode=mode)
-            self._mask = make_mask(m, copy=False, flag=True)
+            self._mask = make_mask(m, copy=False, small_mask=True)
     #............................................
     def ids (self):
         """Return the ids of the data and mask areas."""
@@ -1576,14 +1618,12 @@
 #            if axis is None:
 #                return self._data.sum(None, dtype=dtype)
             return self.__class__(self._data.sum(axis, dtype=dtype),
-                                  mask=nomask, dtype = self.dtype,
-                                  fill_value=self.fill_value)
+                                  mask=nomask, fill_value=self.fill_value)
         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),
-                                  dtype = self.dtype,
                                   fill_value=self.fill_value)
             
     def cumsum(self, axis=None, dtype=None):
@@ -1596,15 +1636,12 @@
 #            if axis is None:
 #                return self._data.cumsum(None, dtype=dtype)
             return self.__class__(self._data.cumsum(axis=axis, dtype=dtype),
-                                  dtype = self.dtype,
                                   fill_value=self.fill_value)
         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,
-                                  dtype = self.dtype,
-                                  fill_value=self.fill_value)
+                                  mask=self._mask, fill_value=self.fill_value)
         
     def prod(self, axis=None, dtype=None):
         """a.prod(axis=None, dtype=None)
@@ -1616,16 +1653,13 @@
 #            if axis is None:
 #                return self._data.prod(None, dtype=dtype)
             return self.__class__(self._data.prod(axis, dtype=dtype),
-                                  mask=nomask,
-                                  dtype = self.dtype,
-                                  fill_value=self.fill_value)
+                                  mask=nomask, fill_value=self.fill_value)
 #            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),
-                                  dtype = self.dtype,
                                   fill_value=self.fill_value)
     product = prod
             
@@ -1639,15 +1673,12 @@
 #            if axis is None:
 #                return self._data.cumprod(None, dtype=dtype)
             return self.__class__(self._data.cumprod(axis=axis, dtype=dtype),
-                                  mask=nomask,
-                                  dtype = self.dtype,
-                                  fill_value=self.fill_value,)
+                                  mask=nomask, fill_value=self.fill_value,)
         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,
-                                  dtype = self.dtype,
                                   fill_value=self.fill_value)        
             
     def mean(self, axis=None, dtype=None):
@@ -1667,8 +1698,7 @@
 #            if axis is None:
 #                return self._data.mean(axis=None, dtype=dtype)
             return self.__class__(self._data.mean(axis=axis, dtype=dtype),
-                                  mask=nomask, dtype = self.dtype,
-                                  fill_value=self.fill_value)
+                                  mask=nomask, fill_value=self.fill_value)
         else:
             dsum = fromnumeric.sum(self.filled(0), axis=axis, dtype=dtype)
             cnt = self.count(axis=axis)
@@ -1676,7 +1706,6 @@
             if axis is None and mask:
                 return masked
             return self.__class__(dsum*1./cnt, mask=mask,
-                                  dtype = self.dtype,
                                   fill_value=self.fill_value,)
     
     def anom(self, axis=None, dtype=None):
@@ -1701,7 +1730,6 @@
 #                return self._data.var(axis=None, dtype=dtype)
             return self.__class__(self._data.var(axis=axis, dtype=dtype),
                                   mask=nomask,
-                                  dtype = self.dtype,
                                   fill_value=self.fill_value)
         else:
             cnt = self.count(axis=axis)
@@ -1713,7 +1741,6 @@
                 return dvar
             return self.__class__(dvar, 
                                   mask=mask_or(self._mask.all(axis), (cnt==1)),
-                                  dtype = self.dtype,
                                   fill_value=self.fill_value)
             
     def std(self, axis=None, dtype=None):
@@ -1970,13 +1997,13 @@
     def __call__ (self, other, *args):
         "Execute the call behavior."
         instance = self.obj
-        m = mask_or(instance._mask, getmask(other), flag=False)
+        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,
-                                  dtype = instance.dtype,
-                                  fill_value=instance.fill_value)  
+#                                  fill_value=instance.fill_value
+                                  )  
 #..........................................................
 MaskedArray.__add__ = _arithmethods('__add__')
 MaskedArray.__radd__ = _arithmethods('__add__')
@@ -2016,15 +2043,16 @@
 masked = masked_singleton
 
 masked_array = MaskedArray
-def array(data, dtype=None, copy=False, order=False,
-          mask=nomask, keep_mask=True, flag=True, fill_value=None):
+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, flag=True, fill_value=None)
+             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, flag=flag, fill_value=fill_value)
+                       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."""
@@ -2344,7 +2372,7 @@
     fb = filled(b, 1)
     if fb.dtype.char in typecodes["Integer"]:
         return masked_array(umath.power(fa, fb), m)
-    md = make_mask((fa < 0), flag=1)
+    md = make_mask((fa < 0), small_mask=1)
     m = mask_or(m, md)
     if m is nomask:
         return masked_array(umath.power(fa, fb))
@@ -2469,7 +2497,7 @@
     dm = []
     for x in arrays:
         dm.append(getmaskarray(x))
-    dm = make_mask(numeric.concatenate(dm, axis), copy=False, flag=True)
+    dm = make_mask(numeric.concatenate(dm, axis), copy=False, small_mask=True)
     return masked_array(d, mask=dm)
 
 def expand_dims(x,axis):
@@ -2598,7 +2626,7 @@
     d = numeric.choose(fc, (yv, xv))
     md = numeric.choose(fc, (ym, xm))
     m = getmask(condition)
-    m = make_mask(mask_or(m, md), copy=False, flag=True)
+    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'):
@@ -2622,7 +2650,7 @@
     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, flag=1)
+    m = make_mask(mask_or(m, getmask(indices)), copy=0, small_mask=1)
     return masked_array(d, mask=m)
 
 def sort (x, axis=-1, fill_value=None, kind='quicksort'):
@@ -2780,7 +2808,7 @@
 in a pickle."""
     _data = ndarray.__new__(ndarray, baseshape, basetype)
     _mask = ndarray.__new__(ndarray, baseshape, basetype)
-    return MaskedArray.__new__(subtype, _data, mask=_mask, dtype=basetype, flag=False)
+    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."

Modified: trunk/Lib/sandbox/maskedarray/extras.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/extras.py	2006-12-22 14:39:32 UTC (rev 2463)
+++ trunk/Lib/sandbox/maskedarray/extras.py	2006-12-23 00:00:43 UTC (rev 2464)
@@ -14,8 +14,10 @@
 __all__ = ['apply_along_axis', 'atleast_1d', 'atleast_2d', 'atleast_3d',
                'average',
            'vstack', 'hstack', 'dstack', 'row_stack', 'column_stack',
-           'compress2d', 'count_masked', 
-           'mask_rowcols','masked_all', 'masked_all_like', 'mr_',
+           'compress_rowcols', 'compress_rows', 'compress_cols', 'count_masked', 
+           'dot',
+           'mask_rowcols', 'mask_rows', 'mask_cols', 'masked_all', 
+               'masked_all_like', 'mr_',
            'notmasked_edges', 'notmasked_contiguous',
            'stdu', 'varu',
            ]
@@ -229,7 +231,6 @@
             outarr[tuple(i.tolist())] = res
             dtypes.append(asarray(res).dtype)
             k += 1
-    print dtypes
     if not hasattr(arr, '_mask'):
         return numeric.asarray(outarr, dtype=max(dtypes))
     else:
@@ -349,13 +350,12 @@
             if isinstance(d, ndarray) and (not d.shape == result.shape):
                 d = ones(result.shape, dtype=float_) * d
     if returned:
-        print type(result)
         return result, d
     else:
         return result
     
 #..............................................................................
-def compress2d(x, axis=None):
+def compress_rowcols(x, axis=None):
     """Suppresses the rows and/or columns of a 2D array that contains masked values.
     The suppression behavior is selected with the `axis`parameter.
         - If axis is None, rows and columns are suppressed. 
@@ -369,7 +369,7 @@
     m = getmask(x)
     # Nothing is masked: return x
     if m is nomask or not m.any():
-        return nxasarray(x)
+        return x._data
     # All is masked: return empty
     if m.all():
         return nxarray([])
@@ -382,11 +382,19 @@
     if axis in [None, 1, -1]:
         for j in function_base.unique(masked[1]):
             idxc.remove(j)
-    return nxasarray(x[idxr][:,idxc])    
+    return x._data[idxr][:,idxc]
 
+def compress_rows(a):
+    """Suppresses whole rows of a 2D array that contain masked values."""
+    return compress_rowcols(a,0)
+
+def compress_cols(a):
+    """Suppresses whole columnss of a 2D array that contain masked values."""
+    return compress_rowcols(a,1)
+
 def mask_rowcols(a, axis=None):
-    """Suppresses the rows and/or columns of a 2D array that contains masked values.
-    The suppression behavior is selected with the `axis`parameter.
+    """Masks whole rows and/or columns of a 2D array that contain masked values.
+    The masking behavior is selected with the `axis`parameter.
         - If axis is None, rows and columns are suppressed. 
         - If axis is 0, only rows are suppressed. 
         - If axis is 1 or -1, only columns are suppressed.
@@ -406,7 +414,37 @@
         a[:,function_base.unique(maskedval[1])] = masked
     return a
 
+def mask_rows(a, axis=None):
+    """Masks whole rows of a 2D array that contain masked values."""
+    return mask_rowcols(a, 0)
 
+def mask_cols(a, axis=None):
+    """Masks whole columns of a 2D array that contain masked values."""
+    return mask_rowcols(a, 1)
+
+        
+def dot(a,b):
+    """Returns the dot product of two 2D masked arrays a and b.
+    Like the generic numpy equivalent the product sum is over
+    the last dimension of a and the second-to-last dimension of b.
+    
+    Masked values are propagated: if a masked value appears in a row or column,
+    the whole row or column is considered masked.
+    
+    NB: The first argument is not conjugated.
+    """
+    #TODO: Works only with 2D arrays. There should be a way to get it to run with higher dimension
+    a = mask_rows(a)
+    b = mask_cols(b)
+    #
+    d = numpy.dot(a.filled(0), b.filled(0))
+    #
+    am = (~getmaskarray(a))
+    bm = (~getmaskarray(b))
+    m = ~numpy.dot(am,bm)
+    return masked_array(d, mask=m)
+
+
 #####--------------------------------------------------------------------------
 #---- --- Concatenation helpers ---
 #####--------------------------------------------------------------------------

Modified: trunk/Lib/sandbox/maskedarray/tests/test_core.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/tests/test_core.py	2006-12-22 14:39:32 UTC (rev 2463)
+++ trunk/Lib/sandbox/maskedarray/tests/test_core.py	2006-12-23 00:00:43 UTC (rev 2464)
@@ -950,6 +950,63 @@
         #We default to true
         mx = masked_array(x, mask=[0,1,0])
         assert_equal(mx.mask, [1,1,0])
+    
+    def check_hardmask(self):        
+        "Test hard_mask"
+        d = arange(5)
+        n = [0,0,0,1,1]
+        m = make_mask(n)
+        xh = array(d, mask = m, hard_mask=True)
+        # We need to copy, to avoid updating d in xh!
+        xs = array(d, mask = m, hard_mask=False, copy=True)
+        xh[[1,4]] = [10,40]
+        xs[[1,4]] = [10,40]
+        assert_equal(xh._data, [0,10,2,3,4])
+        assert_equal(xs._data, [0,10,2,3,40])
+        assert(xh.mask is m)
+        assert_equal(xs.mask, [0,0,0,1,0])
+        assert(xh._hardmask)
+        assert(not xs._hardmask)
+        xh[1:4] = [10,20,30]
+        xs[1:4] = [10,20,30]
+        assert_equal(xh._data, [0,10,20,3,4])
+        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
+        assert_equal(xh.mask, [1,0,0,1,1])
+        assert_equal(xs.mask, [1,0,0,0,0])
+        xh[:] = 1
+        xs[:] = 1
+        assert_equal(xh._data, [0,1,1,3,4])
+        assert_equal(xs._data, [1,1,1,1,1])
+        assert_equal(xh.mask, [1,0,0,1,1])
+        assert_equal(xs.mask, nomask)
+        # Switch to soft mask
+        xh.soften_mask()
+        xh[:] = arange(5)
+        assert_equal(xh._data, [0,1,2,3,4])
+        assert_equal(xh.mask, nomask)
+        # Switch back to hard mask
+        xh.harden_mask()
+        xh[xh<3] = maskedarray.core.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
+        assert_equal(xh._data, [0,1,2,5,5])
+        assert_equal(xh._mask, [1,1,1,0,0])
+        #
+        xh = array([[1,2],[3,4]], mask = [[1,0],[0,0]], hard_mask=True)
+        xh[0] = 0
+        assert_equal(xh._data, [[1,0],[3,4]])
+        assert_equal(xh._mask, [[1,0],[0,0]])
+        xh[-1,-1] = 5
+        assert_equal(xh._data, [[1,0],[3,5]])
+        assert_equal(xh._mask, [[1,0],[0,0]])
+        xh[filled(xh<5,False)] = 2
+        assert_equal(xh._data, [[1,2],[2,5]])
+        assert_equal(xh._mask, [[1,0],[0,0]])
 #..............................................................................
 
 #..............................................................................

Modified: trunk/Lib/sandbox/maskedarray/tests/test_extras.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/tests/test_extras.py	2006-12-22 14:39:32 UTC (rev 2463)
+++ trunk/Lib/sandbox/maskedarray/tests/test_extras.py	2006-12-23 00:00:43 UTC (rev 2464)
@@ -173,26 +173,26 @@
         assert(tmp[1] is None)
         assert_equal(tmp[2][-1], (6, (0,5)))
         
-class test_compress2d(NumpyTestCase):
-    "Tests compress2d and mask_row_columns."
+class test_2dfunctions(NumpyTestCase):
+    "Tests 2D functions"
     def check_compress2d(self):
         "Tests compress2d"
         x = array(N.arange(9).reshape(3,3), mask=[[1,0,0],[0,0,0],[0,0,0]])
-        assert_equal(compress2d(x), [[4,5],[7,8]] )
-        assert_equal(compress2d(x,0), [[3,4,5],[6,7,8]] )
-        assert_equal(compress2d(x,1), [[1,2],[4,5],[7,8]] )
+        assert_equal(compress_rowcols(x), [[4,5],[7,8]] )
+        assert_equal(compress_rowcols(x,0), [[3,4,5],[6,7,8]] )
+        assert_equal(compress_rowcols(x,1), [[1,2],[4,5],[7,8]] )
         x = array(x._data, mask=[[0,0,0],[0,1,0],[0,0,0]])
-        assert_equal(compress2d(x), [[0,2],[6,8]] )
-        assert_equal(compress2d(x,0), [[0,1,2],[6,7,8]] )
-        assert_equal(compress2d(x,1), [[0,2],[3,5],[6,8]] )
+        assert_equal(compress_rowcols(x), [[0,2],[6,8]] )
+        assert_equal(compress_rowcols(x,0), [[0,1,2],[6,7,8]] )
+        assert_equal(compress_rowcols(x,1), [[0,2],[3,5],[6,8]] )
         x = array(x._data, mask=[[1,0,0],[0,1,0],[0,0,0]])
-        assert_equal(compress2d(x), [[8]] )
-        assert_equal(compress2d(x,0), [[6,7,8]] )
-        assert_equal(compress2d(x,1,), [[2],[5],[8]] )
+        assert_equal(compress_rowcols(x), [[8]] )
+        assert_equal(compress_rowcols(x,0), [[6,7,8]] )
+        assert_equal(compress_rowcols(x,1,), [[2],[5],[8]] )
         x = array(x._data, mask=[[1,0,0],[0,1,0],[0,0,1]])
-        assert_equal(compress2d(x).size, 0 )
-        assert_equal(compress2d(x,0).size, 0 )
-        assert_equal(compress2d(x,1).size, 0 )
+        assert_equal(compress_rowcols(x).size, 0 )
+        assert_equal(compress_rowcols(x,0).size, 0 )
+        assert_equal(compress_rowcols(x,1).size, 0 )
     #
     def check_mask_rowcols(self):
         "Tests mask_rowcols."
@@ -212,6 +212,55 @@
         assert(mask_rowcols(x).all())
         assert(mask_rowcols(x,0).all())
         assert(mask_rowcols(x,1).all())
+    #
+    def test_dot(self):        
+        "Tests dot product"
+        n = N.arange(1,7)
+        #
+        m = [1,0,0,0,0,0]
+        a = masked_array(n, mask=m).reshape(2,3)
+        b = masked_array(n, mask=m).reshape(3,2)
+        c = dot(a,b)
+        assert_equal(c.mask, [[1,1],[1,0]])
+        c = dot(b,a)
+        assert_equal(c.mask, [[1,1,1],[1,0,0],[1,0,0]])
+        #        
+        m = [0,0,0,0,0,1]
+        a = masked_array(n, mask=m).reshape(2,3)
+        b = masked_array(n, mask=m).reshape(3,2)
+        c = dot(a,b)
+        assert_equal(c.mask,[[0,1],[1,1]])        
+        c = dot(b,a)
+        assert_equal(c.mask, [[0,0,1],[0,0,1],[1,1,1]])
+        #        
+        m = [0,0,0,0,0,0]
+        a = masked_array(n, mask=m).reshape(2,3)
+        b = masked_array(n, mask=m).reshape(3,2)
+        c = dot(a,b)
+        assert_equal(c.mask,nomask)
+        c = dot(b,a)
+        assert_equal(c.mask,nomask)
+        #        
+        a = masked_array(n, mask=[1,0,0,0,0,0]).reshape(2,3)
+        b = masked_array(n, mask=[0,0,0,0,0,0]).reshape(3,2)
+        c = dot(a,b)
+        assert_equal(c.mask,[[1,1],[0,0]])
+        c = dot(b,a)
+        assert_equal(c.mask,[[1,0,0],[1,0,0],[1,0,0]])
+        #        
+        a = masked_array(n, mask=[0,0,0,0,0,1]).reshape(2,3)
+        b = masked_array(n, mask=[0,0,0,0,0,0]).reshape(3,2)
+        c = dot(a,b)
+        assert_equal(c.mask,[[0,0],[1,1]])
+        c = dot(b,a)
+        assert_equal(c.mask,[[0,0,1],[0,0,1],[0,0,1]])
+        #        
+        a = masked_array(n, mask=[0,0,0,0,0,1]).reshape(2,3)
+        b = masked_array(n, mask=[0,0,1,0,0,0]).reshape(3,2)
+        c = dot(a,b)
+        assert_equal(c.mask,[[1,0],[1,1]])
+        c = dot(b,a)
+        assert_equal(c.mask,[[0,0,1],[1,1,1],[0,0,1]])
 
 ###############################################################################
 #------------------------------------------------------------------------------

Modified: trunk/Lib/sandbox/maskedarray/testutils.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/testutils.py	2006-12-22 14:39:32 UTC (rev 2463)
+++ trunk/Lib/sandbox/maskedarray/testutils.py	2006-12-23 00:00:43 UTC (rev 2464)
@@ -17,7 +17,7 @@
 from numpy.testing.utils import build_err_msg, rand
 
 import core
-reload(core)
+#reload(core)
 from core import mask_or, getmask, getmaskarray, masked_array, nomask
 from core import filled, equal, less
 

Added: trunk/Lib/sandbox/maskedarray/timer_comparison.py
===================================================================
--- trunk/Lib/sandbox/maskedarray/timer_comparison.py	2006-12-22 14:39:32 UTC (rev 2463)
+++ trunk/Lib/sandbox/maskedarray/timer_comparison.py	2006-12-23 00:00:43 UTC (rev 2464)
@@ -0,0 +1,364 @@
+
+import timeit
+
+import numpy 
+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
+
+
+numpy.seterr(all='ignore')
+
+pi = numpy.pi
+
+class moduletester:
+    #-----------------------------------
+    def __init__(self, module):
+        self.module = module
+        self.allequal = module.allequal
+        self.arange = module.arange
+        self.array = module.array
+        self.average =  module.average
+        self.concatenate = module.concatenate
+        self.count = module.count
+        self.filled = module.filled
+        self.getmask = module.getmask 
+        self.id = id
+        self.inner = module.inner
+        self.make_mask = module.make_mask
+        self.masked = module.masked
+        self.masked_array = module.masked_array
+        self.masked_values = module.masked_values
+        self.mask_or = module.mask_or
+        self.ones = module.ones
+        self.outer = module.outer
+        self.repeat = module.repeat
+        self.resize = module.resize
+        self.sort = module.sort
+        self.take = module.take
+        self.transpose = module.transpose
+        self.zeros = module.zeros
+        self.MaskType = module.MaskType
+        try:
+            self.umath = module.umath
+        except AttributeError:
+            self.umath = module.core.umath
+    #----------------------------------  
+    def test_1(self):
+        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.
+        m1 = [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]
+        m2 = [0, 0, 1, 0, 0, 1, 1, 0, 0, 0 ,0, 1]
+        xm = self.masked_array(x, mask=m1)
+        ym = self.masked_array(y, mask=m2)
+        z = numpy.array([-.5, 0., .5, .8])
+        zm = self.masked_array(z, mask=[0,1,0,0])
+        xf = numpy.where(m1, 1.e+20, x)
+        xm.set_fill_value(1.e+20)          
+        #.....
+        assert((xm-ym).filled(0).any())
+        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))
+        #.....
+        for s in [(4,3), (6,2)]:
+            x.shape = s
+            y.shape = s
+            xm.shape = s
+            ym.shape = s
+            xf.shape = s
+    
+            assert_equal(self.count(xm) , len(m1) - reduce(lambda x,y:x+y, m1))
+    #----------------------------------        
+    def test_2(self):        
+        "Tests conversions and indexing"
+        x1 = numpy.array([1,2,4,3])
+        x2 = self.array(x1, mask=[1,0,0,0])
+        x3 = self.array(x1, mask=[0,1,0,1])
+        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))
+    # 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:])
+        x1[2] = 9
+        x2[2] = 9
+        assert_equal(x1,x2)
+        x1[1:3] = 99
+        x2[1:3] = 99
+        assert_equal(x1,x2)
+        x2[1] = self.masked
+        assert_equal(x1,x2)
+        x2[1:3] = self.masked
+        assert_equal(x1,x2)
+        x2[:] = x1
+        x2[1] = self.masked
+        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]))
+        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]))
+        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)
+        x1 = self.array([1,'hello',2,3],object)
+        x2 = numpy.array([1,'hello',2,3],object)
+        s1 = x1[1]
+        s2 = x2[1]
+        assert x1[1:1].shape == (0,)
+        # Tests copy-size
+        n = [0,0,1,0,0]
+        m = self.make_mask(n)
+        m2 = self.make_mask(m)
+        assert(m is m2)
+        m3 = self.make_mask(m, copy=1)
+        assert(m is not m3)
+    
+        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])
+        y5 = self.repeat(x4, (2,2,2,2), axis=0)
+        assert_equal(y5, [0,0,1,1,2,2,3,3])
+        y6 = self.repeat(x4, 2, axis=0)
+        assert_equal(y5, y6)
+        y7 = x4.repeat((2,2,2,2), axis=0)
+        assert_equal(y5,y7)
+        y8 = x4.repeat(2,0)
+        assert_equal(y5,y8)
+        
+        #"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)),
+                            self.inner(x, y))
+        assert_equal(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
+        t = self.take(y,[0,3,4])
+        assert t[0] == 'abc'
+        assert t[1] == 2
+        assert t[2] == 3
+        # Tests in place
+        y = self.arange(10)
+    
+        x = self.arange(10)
+        xm = self.arange(10)
+        xm[2] = self.masked
+        x += 1
+        assert_equal(x, y+1)
+        xm += 1
+        assert_equal(xm, y+1)
+    
+        x = self.arange(10)
+        xm = self.arange(10)
+        xm[2] = self.masked
+        x -= 1
+        assert_equal(x, y-1)
+        xm -= 1
+        assert_equal(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)
+        xm *= 2.0
+        assert_equal(xm, y*2)
+    
+        x = self.arange(10)*2
+        xm = self.arange(10)*2
+        xm[2] = self.masked
+        x /= 2
+        assert_equal(x, y)
+        xm /= 2
+        assert_equal(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)
+        xm /= self.arange(10)
+        assert_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.)
+        
+        x = self.arange(10, dtype=float_)
+        xm = self.arange(10, dtype=float_)
+        xm[2] = self.masked
+        m = xm.mask
+        a = self.arange(10, dtype=float_)
+        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))
+        
+        x = self.arange(10, dtype=float_)
+        xm = self.arange(10, dtype=float_)
+        xm[2] = self.masked
+        m = xm.mask
+        a = self.arange(10, dtype=float_)
+        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))        
+        
+        x = self.arange(10, dtype=float_)
+        xm = self.arange(10, dtype=float_)
+        xm[2] = self.masked
+        m = xm.mask
+        a = self.arange(10, dtype=float_)
+        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))        
+                
+        x = self.arange(10, dtype=float_)
+        xm = self.arange(10, dtype=float_)
+        xm[2] = self.masked
+        m = xm.mask
+        a = self.arange(10, dtype=float_)
+        a[-1] = self.masked
+        x /= a
+        xm /= a
+    #----------------------------------    
+    def test_3(self):
+        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',
+                  ]:
+            #print f
+            try:
+                uf = getattr(self.umath, f)
+            except AttributeError:
+                uf = getattr(fromnumeric, f)
+            mf = getattr(self.module, f)
+            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)
+        # 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.]))
+        result, wts = self.average(ott, weights=[1.,1.,2.,1.], returned=1)
+        assert_equal(2.0, result)
+        assert(wts == 4.0)
+        ott[:] = self.masked
+    #    assert(average(ott,axis=0) is 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))
+        result, wts = self.average(ott, axis=0, returned=1)
+        assert_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)
+        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])
+        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)
+    #    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)
+        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])
+
+
+################################################################################    
+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
+    
+    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




More information about the Scipy-svn mailing list