"""

numpy.ma : a package to handle missing or invalid values.

This package was initially written for numarray by Paul F. Dubois

at Lawrence Livermore National Laboratory.

In 2006, the package was completely rewritten by Pierre Gerard-Marchant

(University of Georgia) to make the MaskedArray class a subclass of ndarray,

and to improve support of structured arrays.

Copyright 1999, 2000, 2001 Regents of the University of California.

Released for unlimited redistribution.

* Adapted for numpy_core 2005 by Travis Oliphant and (mainly) Paul Dubois.

* Subclassing of the base `ndarray` 2006 by Pierre Gerard-Marchant

(pgmdevlist_AT_gmail_DOT_com)

* Improvements suggested by Reggie Dugard (reggie_AT_merfinllc_DOT_com)

.. moduleauthor:: Pierre Gerard-Marchant

"""

# pylint: disable-msg=E1002

from __future__ import division, absolute_import, print_function

import sys

import warnings

from functools import reduce

import numpy as np

import numpy.core.umath as umath

import numpy.core.numerictypes as ntypes

from numpy import ndarray, amax, amin, iscomplexobj, bool_

from numpy import array as narray

from numpy.lib.function_base import angle

from numpy.compat import getargspec, formatargspec, long, basestring

from numpy import expand_dims as n_expand_dims

if sys.version_info[0] >= 3:

import pickle

else:

import cPickle as pickle

__author__ = "Pierre GF Gerard-Marchant"

__docformat__ = "restructuredtext en"

__all__ = ['MAError', 'MaskError', 'MaskType', 'MaskedArray',

'bool_',

'abs', 'absolute', 'add', 'all', 'allclose', 'allequal', 'alltrue',

'amax', 'amin', 'angle', 'anom', 'anomalies', 'any', 'arange',

'arccos', 'arccosh', 'arcsin', 'arcsinh', 'arctan', 'arctan2',

'arctanh', 'argmax', 'argmin', 'argsort', 'around',

'array', 'asarray', 'asanyarray',

'bitwise_and', 'bitwise_or', 'bitwise_xor',

'ceil', 'choose', 'clip', 'common_fill_value', 'compress',

'compressed', 'concatenate', 'conjugate', 'copy', 'cos', 'cosh',

'count', 'cumprod', 'cumsum',

'default_fill_value', 'diag', 'diagonal', 'diff', 'divide', 'dump',

'dumps',

'empty', 'empty_like', 'equal', 'exp', 'expand_dims',

'fabs', 'flatten_mask', 'fmod', 'filled', 'floor', 'floor_divide',

'fix_invalid', 'flatten_structured_array', 'frombuffer', 'fromflex',

'fromfunction',

'getdata', 'getmask', 'getmaskarray', 'greater', 'greater_equal',

'harden_mask', 'hypot',

'identity', 'ids', 'indices', 'inner', 'innerproduct',

'isMA', 'isMaskedArray', 'is_mask', 'is_masked', 'isarray',

'left_shift', 'less', 'less_equal', 'load', 'loads', 'log', 'log2',

'log10', 'logical_and', 'logical_not', 'logical_or', 'logical_xor',

'make_mask', 'make_mask_descr', 'make_mask_none', 'mask_or',

'masked', 'masked_array', 'masked_equal', 'masked_greater',

'masked_greater_equal', 'masked_inside', 'masked_invalid',

'masked_less', 'masked_less_equal', 'masked_not_equal',

'masked_object', 'masked_outside', 'masked_print_option',

'masked_singleton', 'masked_values', 'masked_where', 'max', 'maximum',

'maximum_fill_value', 'mean', 'min', 'minimum', 'minimum_fill_value',

'mod', 'multiply', 'mvoid',

'negative', 'nomask', 'nonzero', 'not_equal',

'ones', 'outer', 'outerproduct',

'power', 'prod', 'product', 'ptp', 'put', 'putmask',

'rank', 'ravel', 'remainder', 'repeat', 'reshape', 'resize',

'right_shift', 'round_', 'round',

'set_fill_value', 'shape', 'sin', 'sinh', 'size', 'sometrue',

'sort', 'soften_mask', 'sqrt', 'squeeze', 'std', 'subtract', 'sum',

'swapaxes',

'take', 'tan', 'tanh', 'trace', 'transpose', 'true_divide',

'var', 'where',

'zeros']

MaskType = np.bool_

nomask = MaskType(0)

def doc_note(initialdoc, note):

"""

Adds a Notes section to an existing docstring.

"""

if initialdoc is None:

return

if note is None:

return initialdoc

newdoc = """

%s

Notes

-----

%s

"""

return newdoc % (initialdoc, note)

def get_object_signature(obj):

"""

Get the signature from obj

"""

try:

sig = formatargspec(*getargspec(obj))

except TypeError as errmsg:

sig = ''

# msg = "Unable to retrieve the signature of %s '%s'\n"\

# "(Initial error message: %s)"

# warnings.warn(msg % (type(obj),

# getattr(obj, '__name__', '???'),

# errmsg))

return sig

#####--------------------------------------------------------------------------

#---- --- Exceptions ---

#####--------------------------------------------------------------------------

class MAError(Exception):

"""Class for masked array related errors."""

pass

class MaskError(MAError):

"Class for mask related errors."

pass

#####--------------------------------------------------------------------------

#---- --- 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' : '???',

'U' : 'N/A',

'M8[D]' : np.datetime64('NaT', 'D'),

'M8[us]' : np.datetime64('NaT', 'us')

}

max_filler = ntypes._minvals

max_filler.update([(k, -np.inf) for k in [np.float32, np.float64]])

min_filler = ntypes._maxvals

min_filler.update([(k, +np.inf) for k in [np.float32, np.float64]])

if 'float128' in ntypes.typeDict:

max_filler.update([(np.float128, -np.inf)])

min_filler.update([(np.float128, +np.inf)])

def default_fill_value(obj):

"""

Return the default fill value for the argument object.

The default filling value depends on the datatype of the input

array or the type of the input scalar:

======== ========

datatype default

======== ========

bool True

int 999999

float 1.e20

complex 1.e20+0j

object '?'

string 'N/A'

======== ========

Parameters

----------

obj : ndarray, dtype or scalar

The array data-type or scalar for which the default fill value

is returned.

Returns

-------

fill_value : scalar

The default fill value.

Examples

--------

>>> np.ma.default_fill_value(1)

999999

>>> np.ma.default_fill_value(np.array([1.1, 2., np.pi]))

1e+20

>>> np.ma.default_fill_value(np.dtype(complex))

(1e+20+0j)

"""

if hasattr(obj, 'dtype'):

defval = _check_fill_value(None, obj.dtype)

elif isinstance(obj, np.dtype):

if obj.subdtype:

defval = default_filler.get(obj.subdtype[0].kind, '?')

elif obj.kind == 'M':

defval = default_filler.get(obj.str[1:], '?')

else:

defval = default_filler.get(obj.kind, '?')

elif isinstance(obj, float):

defval = default_filler['f']

elif isinstance(obj, int) or isinstance(obj, long):

defval = default_filler['i']

elif isinstance(obj, str):

defval = default_filler['S']

elif isinstance(obj, unicode):

defval = default_filler['U']

elif isinstance(obj, complex):

defval = default_filler['c']

else:

defval = default_filler['O']

return defval

def _recursive_extremum_fill_value(ndtype, extremum):

names = ndtype.names

if names:

deflist = []

for name in names:

fval = _recursive_extremum_fill_value(ndtype[name], extremum)

deflist.append(fval)

return tuple(deflist)

return extremum[ndtype]

def minimum_fill_value(obj):

"""

Return the maximum value that can be represented by the dtype of an object.

This function is useful for calculating a fill value suitable for

taking the minimum of an array with a given dtype.

Parameters

----------

obj : ndarray or dtype

An object that can be queried for it's numeric type.

Returns

-------

val : scalar

The maximum representable value.

Raises

------

TypeError

If `obj` isn't a suitable numeric type.

See Also

--------

maximum_fill_value : The inverse function.

set_fill_value : Set the filling value of a masked array.

MaskedArray.fill_value : Return current fill value.

Examples

--------

>>> import numpy.ma as ma

>>> a = np.int8()

>>> ma.minimum_fill_value(a)

127

>>> a = np.int32()

>>> ma.minimum_fill_value(a)

2147483647

An array of numeric data can also be passed.

>>> a = np.array([1, 2, 3], dtype=np.int8)

>>> ma.minimum_fill_value(a)

127

>>> a = np.array([1, 2, 3], dtype=np.float32)

>>> ma.minimum_fill_value(a)

inf

"""

errmsg = "Unsuitable type for calculating minimum."

if hasattr(obj, 'dtype'):

return _recursive_extremum_fill_value(obj.dtype, min_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, np.dtype):

return min_filler[obj]

else:

raise TypeError(errmsg)

def maximum_fill_value(obj):

"""

Return the minimum value that can be represented by the dtype of an object.

This function is useful for calculating a fill value suitable for

taking the maximum of an array with a given dtype.

Parameters

----------

obj : {ndarray, dtype}

An object that can be queried for it's numeric type.

Returns

-------

val : scalar

The minimum representable value.

Raises

------

TypeError

If `obj` isn't a suitable numeric type.

See Also

--------

minimum_fill_value : The inverse function.

set_fill_value : Set the filling value of a masked array.

MaskedArray.fill_value : Return current fill value.

Examples

--------

>>> import numpy.ma as ma

>>> a = np.int8()

>>> ma.maximum_fill_value(a)

-128

>>> a = np.int32()

>>> ma.maximum_fill_value(a)

-2147483648

An array of numeric data can also be passed.

>>> a = np.array([1, 2, 3], dtype=np.int8)

>>> ma.maximum_fill_value(a)

-128

>>> a = np.array([1, 2, 3], dtype=np.float32)

>>> ma.maximum_fill_value(a)

-inf

"""

errmsg = "Unsuitable type for calculating maximum."

if hasattr(obj, 'dtype'):

return _recursive_extremum_fill_value(obj.dtype, max_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, np.dtype):

return max_filler[obj]

else:

raise TypeError(errmsg)

def _recursive_set_default_fill_value(dtypedescr):

deflist = []

for currentdescr in dtypedescr:

currenttype = currentdescr[1]

if isinstance(currenttype, list):

deflist.append(tuple(_recursive_set_default_fill_value(currenttype)))

else:

deflist.append(default_fill_value(np.dtype(currenttype)))

return tuple(deflist)

def _recursive_set_fill_value(fillvalue, dtypedescr):

fillvalue = np.resize(fillvalue, len(dtypedescr))

output_value = []

for (fval, descr) in zip(fillvalue, dtypedescr):

cdtype = descr[1]

if isinstance(cdtype, list):

output_value.append(tuple(_recursive_set_fill_value(fval, cdtype)))

else:

output_value.append(np.array(fval, dtype=cdtype).item())

return tuple(output_value)

def _check_fill_value(fill_value, ndtype):

"""

Private function validating the given `fill_value` for the given dtype.

If fill_value is None, it is set to the default corresponding to the dtype

if this latter is standard (no fields). If the datatype is flexible (named

fields), fill_value is set to a tuple whose elements are the default fill

values corresponding to each field.

If fill_value is not None, its value is forced to the given dtype.

"""

ndtype = np.dtype(ndtype)

fields = ndtype.fields

if fill_value is None:

if fields:

descr = ndtype.descr

fill_value = np.array(_recursive_set_default_fill_value(descr),

dtype=ndtype,)

else:

fill_value = default_fill_value(ndtype)

elif fields:

fdtype = [(_[0], _[1]) for _ in ndtype.descr]

if isinstance(fill_value, (ndarray, np.void)):

try:

fill_value = np.array(fill_value, copy=False, dtype=fdtype)

except ValueError:

err_msg = "Unable to transform %s to dtype %s"

raise ValueError(err_msg % (fill_value, fdtype))

else:

descr = ndtype.descr

fill_value = np.asarray(fill_value, dtype=object)

fill_value = np.array(_recursive_set_fill_value(fill_value, descr),

dtype=ndtype)

else:

if isinstance(fill_value, basestring) and (ndtype.char not in 'SVU'):

fill_value = default_fill_value(ndtype)

else:

# In case we want to convert 1e+20 to int...

try:

fill_value = np.array(fill_value, copy=False, dtype=ndtype)#.item()

except OverflowError:

fill_value = default_fill_value(ndtype)

return np.array(fill_value)

def set_fill_value(a, fill_value):

"""

Set the filling value of a, if a is a masked array.

This function changes the fill value of the masked array `a` in place.

If `a` is not a masked array, the function returns silently, without

doing anything.

Parameters

----------

a : array_like

Input array.

fill_value : dtype

Filling value. A consistency test is performed to make sure

the value is compatible with the dtype of `a`.

Returns

-------

None

Nothing returned by this function.

See Also

--------

maximum_fill_value : Return the default fill value for a dtype.

MaskedArray.fill_value : Return current fill value.

MaskedArray.set_fill_value : Equivalent method.

Examples

--------

>>> import numpy.ma as ma

>>> a = np.arange(5)

>>> a

array([0, 1, 2, 3, 4])

>>> a = ma.masked_where(a < 3, a)

>>> a

masked_array(data = [-- -- -- 3 4],

mask = [ True True True False False],

fill_value=999999)

>>> ma.set_fill_value(a, -999)

>>> a

masked_array(data = [-- -- -- 3 4],

mask = [ True True True False False],

fill_value=-999)

Nothing happens if `a` is not a masked array.

>>> a = range(5)

>>> a

[0, 1, 2, 3, 4]

>>> ma.set_fill_value(a, 100)

>>> a

[0, 1, 2, 3, 4]

>>> a = np.arange(5)

>>> a

array([0, 1, 2, 3, 4])

>>> ma.set_fill_value(a, 100)

>>> a

array([0, 1, 2, 3, 4])

"""

if isinstance(a, MaskedArray):

a.set_fill_value(fill_value)

return

def get_fill_value(a):

"""

Return the filling value of a, if any. Otherwise, returns the

default filling value for that type.

"""

if isinstance(a, MaskedArray):

result = a.fill_value

else:

result = default_fill_value(a)

return result

def common_fill_value(a, b):

"""

Return the common filling value of two masked arrays, if any.

If ``a.fill_value == b.fill_value``, return the fill value,

otherwise return None.

Parameters

----------

a, b : MaskedArray

The masked arrays for which to compare fill values.

Returns

-------

fill_value : scalar or None

The common fill value, or None.

Examples

--------

>>> x = np.ma.array([0, 1.], fill_value=3)

>>> y = np.ma.array([0, 1.], fill_value=3)

>>> np.ma.common_fill_value(x, y)

3.0

"""

t1 = get_fill_value(a)

t2 = get_fill_value(b)

if t1 == t2:

return t1

return None

#####--------------------------------------------------------------------------

def filled(a, fill_value=None):

"""

Return input as an array with masked data replaced by a fill value.

If `a` is not a `MaskedArray`, `a` itself is returned.

If `a` is a `MaskedArray` and `fill_value` is None, `fill_value` is set to

``a.fill_value``.

Parameters

----------

a : MaskedArray or array_like

An input object.

fill_value : scalar, optional

Filling value. Default is None.

Returns

-------

a : ndarray

The filled array.

See Also

--------

compressed

Examples

--------

>>> x = np.ma.array(np.arange(9).reshape(3, 3), mask=[[1, 0, 0],

... [1, 0, 0],

... [0, 0, 0]])

>>> x.filled()

array([[999999, 1, 2],

[999999, 4, 5],

[ 6, 7, 8]])

"""

if hasattr(a, 'filled'):

return a.filled(fill_value)

elif isinstance(a, ndarray):

# Should we check for contiguity ? and a.flags['CONTIGUOUS']:

return a

elif isinstance(a, dict):

return np.array(a, 'O')

else:

return np.array(a)

#####--------------------------------------------------------------------------

def get_masked_subclass(*arrays):

"""

Return the youngest subclass of MaskedArray from a list of (masked) arrays.

In case of siblings, the first listed 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

# Don't return MaskedConstant as result: revert to MaskedArray

if rcls.__name__ == 'MaskedConstant':

return MaskedArray

return rcls

#####--------------------------------------------------------------------------

def getdata(a, subok=True):

"""

Return the data of a masked array as an ndarray.

Return the data of `a` (if any) as an ndarray if `a` is a ``MaskedArray``,

else return `a` as a ndarray or subclass (depending on `subok`) if not.

Parameters

----------

a : array_like

Input ``MaskedArray``, alternatively a ndarray or a subclass thereof.

subok : bool

Whether to force the output to be a `pure` ndarray (False) or to

return a subclass of ndarray if appropriate (True, default).

See Also

--------

getmask : Return the mask of a masked array, or nomask.

getmaskarray : Return the mask of a masked array, or full array of False.

Examples

--------

>>> import numpy.ma as ma

>>> a = ma.masked_equal([[1,2],[3,4]], 2)

>>> a

masked_array(data =

[[1 --]

[3 4]],

mask =

[[False True]

[False False]],

fill_value=999999)

>>> ma.getdata(a)

array([[1, 2],

[3, 4]])

Equivalently use the ``MaskedArray`` `data` attribute.

>>> a.data

array([[1, 2],

[3, 4]])

"""

try:

data = a._data

except AttributeError:

data = np.array(a, copy=False, subok=subok)

if not subok:

return data.view(ndarray)

return data

get_data = getdata

def fix_invalid(a, mask=nomask, copy=True, fill_value=None):

"""

Return input with invalid data masked and replaced by a fill value.

Invalid data means values of `nan`, `inf`, etc.

Parameters

----------

a : array_like

Input array, a (subclass of) ndarray.

copy : bool, optional

Whether to use a copy of `a` (True) or to fix `a` in place (False).

Default is True.

fill_value : scalar, optional

Value used for fixing invalid data. Default is None, in which case

the ``a.fill_value`` is used.

Returns

-------

b : MaskedArray

The input array with invalid entries fixed.

Notes

-----

A copy is performed by default.

Examples

--------

>>> x = np.ma.array([1., -1, np.nan, np.inf], mask=[1] + [0]*3)

>>> x

masked_array(data = [-- -1.0 nan inf],

mask = [ True False False False],

fill_value = 1e+20)

>>> np.ma.fix_invalid(x)

masked_array(data = [-- -1.0 -- --],

mask = [ True False True True],

fill_value = 1e+20)

>>> fixed = np.ma.fix_invalid(x)

>>> fixed.data

array([ 1.00000000e+00, -1.00000000e+00, 1.00000000e+20,

1.00000000e+20])

>>> x.data

array([ 1., -1., NaN, Inf])

"""

a = masked_array(a, copy=copy, mask=mask, subok=True)

#invalid = (numpy.isnan(a._data) | numpy.isinf(a._data))

invalid = np.logical_not(np.isfinite(a._data))

if not invalid.any():

return a

a._mask |= invalid

if fill_value is None:

fill_value = a.fill_value

a._data[invalid] = fill_value

return a

#####--------------------------------------------------------------------------

#---- --- Ufuncs ---

#####--------------------------------------------------------------------------

ufunc_domain = {}

ufunc_fills = {}

class _DomainCheckInterval:

"""

Define a valid interval, so that :

``domain_check_interval(a,b)(x) == True`` where

``x < a`` or ``x > b``.

"""

def __init__(self, a, b):

"domain_check_interval(a,b)(x) = true where x < a or y > b"

if (a > b):

(a, b) = (b, a)

self.a = a

self.b = b

def __call__ (self, x):

"Execute the call behavior."

return umath.logical_or(umath.greater (x, self.b),

umath.less(x, self.a))

class _DomainTan:

"""Define a valid interval for the `tan` function, so that:

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

"""

def __init__(self, eps):

"domain_tan(eps) = true where abs(cos(x)) < eps)"

self.eps = eps

def __call__ (self, x):

"Executes the call behavior."

return umath.less(umath.absolute(umath.cos(x)), self.eps)

class _DomainSafeDivide:

"""Define a domain for safe division."""

def __init__ (self, tolerance=None):

self.tolerance = tolerance

def __call__ (self, a, b):

# Delay the selection of the tolerance to here in order to reduce numpy

# import times. The calculation of these parameters is a substantial

# component of numpy's import time.

if self.tolerance is None:

self.tolerance = np.finfo(float).tiny

return umath.absolute(a) * self.tolerance >= umath.absolute(b)

class _DomainGreater:

"""DomainGreater(v)(x) is True where x <= v."""

def __init__(self, critical_value):

"DomainGreater(v)(x) = true where x <= v"

self.critical_value = critical_value

def __call__ (self, x):

"Executes the call behavior."

return umath.less_equal(x, self.critical_value)

class _DomainGreaterEqual:

"""DomainGreaterEqual(v)(x) is True where x < v."""

def __init__(self, critical_value):

"DomainGreaterEqual(v)(x) = true where x < v"

self.critical_value = critical_value

def __call__ (self, x):

"Executes the call behavior."

return umath.less(x, self.critical_value)

#..............................................................................

class _MaskedUnaryOperation:

"""

Defines masked version of unary operations, where invalid values are

pre-masked.

Parameters

----------

mufunc : callable

The function for which to define a masked version. Made available

as ``_MaskedUnaryOperation.f``.

fill : scalar, optional

Filling value, default is 0.

domain : class instance

Domain for the function. Should be one of the ``_Domain*``

classes. Default is None.

"""

def __init__ (self, mufunc, fill=0, domain=None):

""" _MaskedUnaryOperation(aufunc, fill=0, domain=None)

aufunc(fill) must be defined

self(x) returns aufunc(x)

with masked values where domain(x) is true or getmask(x) is true.

"""

self.f = mufunc

self.fill = fill

self.domain = domain

self.__doc__ = getattr(mufunc, "__doc__", str(mufunc))

self.__name__ = getattr(mufunc, "__name__", str(mufunc))

ufunc_domain[mufunc] = domain

ufunc_fills[mufunc] = fill

#

def __call__ (self, a, *args, **kwargs):

"Execute the call behavior."

d = getdata(a)

# Case 1.1. : Domained function

if self.domain is not None:

with np.errstate():

np.seterr(divide='ignore', invalid='ignore')

result = self.f(d, *args, **kwargs)

# Make a mask

m = ~umath.isfinite(result)

m |= self.domain(d)

m |= getmask(a)

# Case 1.2. : Function without a domain

else:

# Get the result and the mask

result = self.f(d, *args, **kwargs)

m = getmask(a)

# Case 2.1. : The result is scalarscalar

if not result.ndim:

if m:

return masked

return result

# Case 2.2. The result is an array

# We need to fill the invalid data back w/ the input

# Now, that's plain silly: in C, we would just skip the element and keep

# the original, but we do have to do it that way in Python

if m is not nomask:

# In case result has a lower dtype than the inputs (as in equal)

try:

np.copyto(result, d, where=m)

except TypeError:

pass

# Transform to

if isinstance(a, MaskedArray):

subtype = type(a)

else:

subtype = MaskedArray

result = result.view(subtype)

result._mask = m

result._update_from(a)

return result

#

def __str__ (self):

return "Masked version of %s. [Invalid values are masked]" % str(self.f)

class _MaskedBinaryOperation:

"""

Define masked version of binary operations, where invalid

values are pre-masked.

Parameters

----------

mbfunc : function

The function for which to define a masked version. Made available

as ``_MaskedBinaryOperation.f``.

domain : class instance

Default domain for the function. Should be one of the ``_Domain*``

classes. Default is None.

fillx : scalar, optional

Filling value for the first argument, default is 0.

filly : scalar, optional

Filling value for the second argument, default is 0.

"""

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."

# Get the data, as ndarray

(da, db) = (getdata(a, subok=False), getdata(b, subok=False))

# Get the mask

(ma, mb) = (getmask(a), getmask(b))

if ma is nomask:

if mb is nomask:

m = nomask

else:

m = umath.logical_or(getmaskarray(a), mb)

elif mb is nomask:

m = umath.logical_or(ma, getmaskarray(b))

else:

m = umath.logical_or(ma, mb)

# Get the result

with np.errstate():

np.seterr(divide='ignore', invalid='ignore')

result = self.f(da, db, *args, **kwargs)

# check it worked

if result is NotImplemented:

return NotImplemented

# Case 1. : scalar

if not result.ndim:

if m:

return masked

return result

# Case 2. : array

# Revert result to da where masked

if m.any():

np.copyto(result, 0, casting='unsafe', where=m)

# This only makes sense if the operation preserved the dtype

if result.dtype == da.dtype:

result += m * da

# Transforms to a (subclass of) MaskedArray

result = result.view(get_masked_subclass(a, b))

result._mask = m

# Update the optional info from the inputs

if isinstance(b, MaskedArray):

if isinstance(a, MaskedArray):

result._update_from(a)

else:

result._update_from(b)

elif isinstance(a, MaskedArray):

result._update_from(a)

return result

def reduce(self, target, axis=0, dtype=None):

"""Reduce `target` along the given `axis`."""

if isinstance(target, MaskedArray):

tclass = type(target)

else:

tclass = MaskedArray

m = getmask(target)

t = filled(target, self.filly)

if t.shape == ():

t = t.reshape(1)

if m is not nomask:

m = make_mask(m, copy=1)

m.shape = (1,)

if m is nomask:

return self.f.reduce(t, axis).view(tclass)

t = t.view(tclass)

t._mask = m

tr = self.f.reduce(getdata(t), axis, dtype=dtype or t.dtype)

mr = umath.logical_and.reduce(m, axis)

tr = tr.view(tclass)

if mr.ndim > 0:

tr._mask = mr

return tr

elif mr:

return masked

return tr

def outer (self, a, b):

"""Return the function applied to the outer product of a and b.

"""

ma = getmask(a)

mb = getmask(b)

if ma is nomask and mb is nomask:

m = nomask

else:

ma = getmaskarray(a)

mb = getmaskarray(b)

m = umath.logical_or.outer(ma, mb)