Changes between Version 1 and Version 2 of MaskedArrayAlternative

Timestamp:

08/25/07 19:00:52 (6 years ago)

Author:

pierregm

Comment:

--

MaskedArrayAlternative

@@ -3 +3 @@
 '''Note: the new implementation of MaskedArray is now available in the scipy sandbox. '''
 
-=== History ===
+== History ==
 
 As a regular user of MaskedArray, I (Pierre G.F. Gerard-Marchant) became increasingly frustrated with the subclassing of masked arrays (even if I can only blame my inexperience). I needed to develop a class of arrays that could store some additional information along with numerical values, while keeping the possibility for missing data (picture storing a series of dates along with measurements, what would later become the {{{TimeSeries}}} package). 
@@ -14 +14 @@
 Note that if the subclass has some special methods and attributes, they are not propagated to the masked version: this would require a modification of the {{{__getattribute__}}} method (first trying {{{ndarray.__getattribute__}}}, then trying {{{self._data.__getattribute__}}} if an exception is raised in the first place), which really slows things down. 
 
-=== Main differences ===
+== Main differences ==
   * The {{{_data}}} part of the masked array can be any subclass of ndarray (but not recarray, cf below).
   * {{{fill_value}}} is now a property, not a function.
@@ -27 +27 @@
   * {{{cumprod}}} works as if the {{{_data}}} array was filled with 1. The mask is preserved, but not updated.
 
-=== New features ===
+== New features ==
+  This list is non-exhaustive...
   * the {{{mr_}}} function mimics {{{r_}}} for masked arrays.
   * the {{{anom}}} method returns the anomalies (deviations from the average)
   * the {{{stdu}}} and {{{varu}}} return unbiased estimates of the standard deviation and variance, respectively.
 
-=== Using the new package with numpy.core.ma ===
+== Using the new package with numpy.core.ma ==
   I tried to make sure that the new package can understand old masked arrays. Unfortunately, there's no upward compatibility.
 For example:
@@ -64 +65 @@
 }}}
   
-=== Using maskedarray with matplotlib ===
+== Using maskedarray with matplotlib ==
 By default matplotlib still uses numpy.ma, but there is an rcParams setting that you can use to select maskedarray instead.  In the matplotlibrc file you will find:
 
@@ -79 +80 @@
 }}}
 
-=== Masked records ===
+== Masked records ==
   Like {{{numpy.core.ma}}}, the {{{ndarray}}}-based implementation of {{{MaskedArray}}} is limited when working with records: you can mask any record of the array, but not a field in a record. If you need this feature, you may want to give the {{{mrecords}}} package a try (available in the {{{maskedarray}}} directory in the scipy sandbox). This module defines a new class, {{{MaskedRecord}}}. An instance of this class accepts a {{{recarray}}} as data, and uses two masks: the {{{fieldmask}}} has as many entries as records in the array, each entry with the same fields as a record, but of boolean types: they indicate whether the field is masked or not; a record entry is flagged as masked in the {{{mask}}} array if all the fields are masked. A few examples in the file should give you an idea of what can be done. Note that {{{mrecords}}} is still experimental...
 
+== Optimizing maskedarray ==
 
-=== Thanks ===
+=== Should masked arrays be filled before processing or not ? ===
+  In the current implementation, most operations on masked arrays involve the following steps:
+* the input arrays are filled
+* the operation is performed on the filled arrays
+* the mask is set for the results, from the combination of the input masks and the mask corresponding to the domain of the operation.
+
+For example, consider the division of two masked arrays:
+{{{
+#!python
+import numpy
+import maskedarray as ma
+x = ma.array([1,2,3,4],mask=[1,0,0,0], dtype=numpy.float_)
+y = ma.array([-1,0,1,2], mask=[0,0,0,1], dtype=numpy.float_)
+}}}
+
+The division of x by y is then computed as
+{{{
+#!python
+d1 = x.filled(0) # d1 = array([0., 2., 3., 4.])
+d2 = y.filled(1) # array([-1.,  0.,  1.,  1.])
+m = ma.mask_or(ma.getmask(x), ma.getmask(y)) # m = array([True,False,False,True])
+dm = ma.divide.domain(d1,d2) # array([False,  True, False, False])
+result = (d1/d2).view(MaskedArray) # masked_array([-0. inf, 3., 4.])
+result._mask = logical_or(m, dm)
+}}}
+
+Note that a division by zero takes place. To avoid it, we can consider to fill the input arrays, taking the domain mask into account, so that:
+{{{
+#!python
+d1 = x._data.copy() # d1 = array([1., 2., 3., 4.])
+d2 = y._data.copy() # array([-1.,  0.,  1.,  2.])
+dm = ma.divide.domain(d1,d2) # array([False,  True, False, False])
+numpy.putmask(d2, dm, 1) # d2 = array([-1.,  1.,  1.,  2.])
+m = ma.mask_or(ma.getmask(x), ma.getmask(y)) # m = array([True,False,False,True])
+result = (d1/d2).view(MaskedArray) # masked_array([-1. 0., 3., 2.])
+result._mask = logical_or(m, dm)
+}}}
+Note that the {{{.copy()}}} is required to avoid updating the inputs with {{{putmask}}}.
+In the previous version, the {{{.filled}}} methods involved a {{{.copy()}}}.
+
+A third possibility consists in avoid filling the arrays:
+{{{
+#!python
+d1 = x._data # d1 = array([1., 2., 3., 4.])
+d2 = y._data # array([-1.,  0.,  1.,  2.])
+dm = ma.divide.domain(d1,d2) # array([False,  True, False, False])
+m = ma.mask_or(ma.getmask(x), ma.getmask(y)) # m = array([True,False,False,True])
+result = (d1/d2).view(MaskedArray) # masked_array([-1. inf, 3., 2.])
+result._mask = logical_or(m, dm)
+}}}
+Note that here again the division by zero takes place.
+
+A quick benchmark gives the following results:
+  * {{{numpy.ma.divide}}}  : 2.84 ms per loop
+  * classical division     : 2.99 ms per loop
+  * division w/ prefilling : 2.20 ms per loop
+  * division w/o filling   : 1.54 ms per loop
+
+So, is it worth filling the arrays beforehand ? Yes, if we are interested in avoiding floating-point exceptions that may fill the result with infs and nans. No, if we are only interested into speed...
+
+
+
+
+== Thanks ==
   I'd like to thank Paul Dubois, Travis Oliphant and Sasha for the original masked array package: without you, I would never have started that (it might be argued that I shouldn't have anyway, but that's another story...).
   I also wish to extend this thanks to Reggie Dugard and Eric Firing for their suggestions and numerous improvements.