~jtaylor/ubuntu/precise/python-numpy/multiarch-fix-818867

« back to all changes in this revision

Viewing changes to numpy/oldnumeric/ma.py

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik, Riku Voipio, Tiziano Zito, Carlos Galisteo, Ondrej Certik
  • Date: 2008-07-08 15:08:16 UTC
  • mfrom: (0.1.21 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080708150816-ekf992jcp2k1eua3
Tags: 1:1.1.0-3
[ Riku Voipio ]
* debian/control: atlas is not available on armel, and after a quick look
  neither on alpha. I'd also suggest dropping
  libatlas-sse-dev|libatlas-sse2-dev|libatlas-3dnow-dev alternative combo
  away, these are potentially dangerous on buildd's. Ondrej: dropped.
  (Closes: #489568)

[ Tiziano Zito ]
* patch: build _dotblas.c when ATLAS is not installed, build-conflict with
  atlas, build-depend on blas+lapack only, as it used to be (Closes: #489726)

[ Carlos Galisteo ]
* debian/control
  - Added Homepage field.

[ Ondrej Certik ]
* Checked the package on i386 and amd64, both with and without atlas, all
  tests run and the numpy package is faster if atlas is around. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
# Incompatibility in that getmask and a.mask returns nomask
2
 
#  instead of None
3
 
 
4
 
from numpy.core.ma import *
5
 
import numpy.core.ma as nca
 
1
"""MA: a facility for dealing with missing observations
 
2
MA is generally used as a numpy.array look-alike.
 
3
by Paul F. Dubois.
 
4
 
 
5
Copyright 1999, 2000, 2001 Regents of the University of California.
 
6
Released for unlimited redistribution.
 
7
Adapted for numpy_core 2005 by Travis Oliphant and
 
8
(mainly) Paul Dubois.
 
9
 
 
10
"""
 
11
import types, sys
 
12
 
 
13
import numpy.core.umath as umath
 
14
import numpy.core.fromnumeric as fromnumeric
 
15
from numpy.core.numeric import newaxis, ndarray, inf
 
16
from numpy.core.fromnumeric import amax, amin
 
17
from numpy.core.numerictypes import bool_, typecodes
 
18
import numpy.core.numeric as numeric
 
19
import warnings
 
20
 
 
21
# Ufunc domain lookup for __array_wrap__
 
22
ufunc_domain = {}
 
23
# Ufunc fills lookup for __array__
 
24
ufunc_fills = {}
 
25
 
 
26
MaskType = bool_
 
27
nomask = MaskType(0)
 
28
divide_tolerance = 1.e-35
 
29
 
 
30
class MAError (Exception):
 
31
    def __init__ (self, args=None):
 
32
        "Create an exception"
 
33
 
 
34
        # The .args attribute must be a tuple.
 
35
        if not isinstance(args, tuple):
 
36
            args = (args,)
 
37
        self.args = args
 
38
    def __str__(self):
 
39
        "Calculate the string representation"
 
40
        return str(self.args[0])
 
41
    __repr__ = __str__
 
42
 
 
43
class _MaskedPrintOption:
 
44
    "One instance of this class, masked_print_option, is created."
 
45
    def __init__ (self, display):
 
46
        "Create the masked print option object."
 
47
        self.set_display(display)
 
48
        self._enabled = 1
 
49
 
 
50
    def display (self):
 
51
        "Show what prints for masked values."
 
52
        return self._display
 
53
 
 
54
    def set_display (self, s):
 
55
        "set_display(s) sets what prints for masked values."
 
56
        self._display = s
 
57
 
 
58
    def enabled (self):
 
59
        "Is the use of the display value enabled?"
 
60
        return self._enabled
 
61
 
 
62
    def enable(self, flag=1):
 
63
        "Set the enabling flag to flag."
 
64
        self._enabled = flag
 
65
 
 
66
    def __str__ (self):
 
67
        return str(self._display)
 
68
 
 
69
    __repr__ = __str__
 
70
 
 
71
#if you single index into a masked location you get this object.
 
72
masked_print_option = _MaskedPrintOption('--')
 
73
 
 
74
# Use single element arrays or scalars.
 
75
default_real_fill_value = 1.e20
 
76
default_complex_fill_value = 1.e20 + 0.0j
 
77
default_character_fill_value = '-'
 
78
default_integer_fill_value = 999999
 
79
default_object_fill_value = '?'
 
80
 
 
81
def default_fill_value (obj):
 
82
    "Function to calculate default fill value for an object."
 
83
    if isinstance(obj, types.FloatType):
 
84
        return default_real_fill_value
 
85
    elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType):
 
86
        return default_integer_fill_value
 
87
    elif isinstance(obj, types.StringType):
 
88
        return default_character_fill_value
 
89
    elif isinstance(obj, types.ComplexType):
 
90
        return default_complex_fill_value
 
91
    elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray):
 
92
        x = obj.dtype.char
 
93
        if x in typecodes['Float']:
 
94
            return default_real_fill_value
 
95
        if x in typecodes['Integer']:
 
96
            return default_integer_fill_value
 
97
        if x in typecodes['Complex']:
 
98
            return default_complex_fill_value
 
99
        if x in typecodes['Character']:
 
100
            return default_character_fill_value
 
101
        if x in typecodes['UnsignedInteger']:
 
102
            return umath.absolute(default_integer_fill_value)
 
103
        return default_object_fill_value
 
104
    else:
 
105
        return default_object_fill_value
 
106
 
 
107
def minimum_fill_value (obj):
 
108
    "Function to calculate default fill value suitable for taking minima."
 
109
    if isinstance(obj, types.FloatType):
 
110
        return numeric.inf
 
111
    elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType):
 
112
        return sys.maxint
 
113
    elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray):
 
114
        x = obj.dtype.char
 
115
        if x in typecodes['Float']:
 
116
            return numeric.inf
 
117
        if x in typecodes['Integer']:
 
118
            return sys.maxint
 
119
        if x in typecodes['UnsignedInteger']:
 
120
            return sys.maxint
 
121
    else:
 
122
        raise TypeError, 'Unsuitable type for calculating minimum.'
 
123
 
 
124
def maximum_fill_value (obj):
 
125
    "Function to calculate default fill value suitable for taking maxima."
 
126
    if isinstance(obj, types.FloatType):
 
127
        return -inf
 
128
    elif isinstance(obj, types.IntType) or isinstance(obj, types.LongType):
 
129
        return -sys.maxint
 
130
    elif isinstance(obj, MaskedArray) or isinstance(obj, ndarray):
 
131
        x = obj.dtype.char
 
132
        if x in typecodes['Float']:
 
133
            return -inf
 
134
        if x in typecodes['Integer']:
 
135
            return -sys.maxint
 
136
        if x in typecodes['UnsignedInteger']:
 
137
            return 0
 
138
    else:
 
139
        raise TypeError, 'Unsuitable type for calculating maximum.'
 
140
 
 
141
def set_fill_value (a, fill_value):
 
142
    "Set fill value of a if it is a masked array."
 
143
    if isMaskedArray(a):
 
144
        a.set_fill_value (fill_value)
 
145
 
 
146
def getmask (a):
 
147
    """Mask of values in a; could be nomask.
 
148
       Returns nomask if a is not a masked array.
 
149
       To get an array for sure use getmaskarray."""
 
150
    if isinstance(a, MaskedArray):
 
151
        return a.raw_mask()
 
152
    else:
 
153
        return nomask
 
154
 
 
155
def getmaskarray (a):
 
156
    """Mask of values in a; an array of zeros if mask is nomask
 
157
     or not a masked array, and is a byte-sized integer.
 
158
     Do not try to add up entries, for example.
 
159
    """
 
160
    m = getmask(a)
 
161
    if m is nomask:
 
162
        return make_mask_none(shape(a))
 
163
    else:
 
164
        return m
 
165
 
 
166
def is_mask (m):
 
167
    """Is m a legal mask? Does not check contents, only type.
 
168
    """
 
169
    try:
 
170
        return m.dtype.type is MaskType
 
171
    except AttributeError:
 
172
        return False
 
173
 
 
174
def make_mask (m, copy=0, flag=0):
 
175
    """make_mask(m, copy=0, flag=0)
 
176
       return m as a mask, creating a copy if necessary or requested.
 
177
       Can accept any sequence of integers or nomask. Does not check
 
178
       that contents must be 0s and 1s.
 
179
       if flag, return nomask if m contains no true elements.
 
180
    """
 
181
    if m is nomask:
 
182
        return nomask
 
183
    elif isinstance(m, ndarray):
 
184
        if m.dtype.type is MaskType:
 
185
            if copy:
 
186
                result = numeric.array(m, dtype=MaskType, copy=copy)
 
187
            else:
 
188
                result = m
 
189
        else:
 
190
            result = m.astype(MaskType)
 
191
    else:
 
192
        result = filled(m, True).astype(MaskType)
 
193
 
 
194
    if flag and not fromnumeric.sometrue(fromnumeric.ravel(result)):
 
195
        return nomask
 
196
    else:
 
197
        return result
 
198
 
 
199
def make_mask_none (s):
 
200
    "Return a mask of all zeros of shape s."
 
201
    result = numeric.zeros(s, dtype=MaskType)
 
202
    result.shape = s
 
203
    return result
 
204
 
 
205
def mask_or (m1, m2):
 
206
    """Logical or of the mask candidates m1 and m2, treating nomask as false.
 
207
       Result may equal m1 or m2 if the other is nomask.
 
208
     """
 
209
    if m1 is nomask: return make_mask(m2)
 
210
    if m2 is nomask: return make_mask(m1)
 
211
    if m1 is m2 and is_mask(m1): return m1
 
212
    return make_mask(umath.logical_or(m1, m2))
 
213
 
 
214
def filled (a, value = None):
 
215
    """a as a contiguous numeric array with any masked areas replaced by value
 
216
    if value is None or the special element "masked", get_fill_value(a)
 
217
    is used instead.
 
218
 
 
219
    If a is already a contiguous numeric array, a itself is returned.
 
220
 
 
221
    filled(a) can be used to be sure that the result is numeric when
 
222
    passing an object a to other software ignorant of MA, in particular to
 
223
    numeric itself.
 
224
    """
 
225
    if isinstance(a, MaskedArray):
 
226
        return a.filled(value)
 
227
    elif isinstance(a, ndarray) and a.flags['CONTIGUOUS']:
 
228
        return a
 
229
    elif isinstance(a, types.DictType):
 
230
        return numeric.array(a, 'O')
 
231
    else:
 
232
        return numeric.array(a)
 
233
 
 
234
def get_fill_value (a):
 
235
    """
 
236
    The fill value of a, if it has one; otherwise, the default fill value
 
237
    for that type.
 
238
    """
 
239
    if isMaskedArray(a):
 
240
        result = a.fill_value()
 
241
    else:
 
242
        result = default_fill_value(a)
 
243
    return result
 
244
 
 
245
def common_fill_value (a, b):
 
246
    "The common fill_value of a and b, if there is one, or None"
 
247
    t1 = get_fill_value(a)
 
248
    t2 = get_fill_value(b)
 
249
    if t1 == t2: return t1
 
250
    return None
 
251
 
 
252
# Domain functions return 1 where the argument(s) are not in the domain.
 
253
class domain_check_interval:
 
254
    "domain_check_interval(a,b)(x) = true where x < a or y > b"
 
255
    def __init__(self, y1, y2):
 
256
        "domain_check_interval(a,b)(x) = true where x < a or y > b"
 
257
        self.y1 = y1
 
258
        self.y2 = y2
 
259
 
 
260
    def __call__ (self, x):
 
261
        "Execute the call behavior."
 
262
        return umath.logical_or(umath.greater (x, self.y2),
 
263
                                   umath.less(x, self.y1)
 
264
                                  )
 
265
 
 
266
class domain_tan:
 
267
    "domain_tan(eps) = true where abs(cos(x)) < eps)"
 
268
    def __init__(self, eps):
 
269
        "domain_tan(eps) = true where abs(cos(x)) < eps)"
 
270
        self.eps = eps
 
271
 
 
272
    def __call__ (self, x):
 
273
        "Execute the call behavior."
 
274
        return umath.less(umath.absolute(umath.cos(x)), self.eps)
 
275
 
 
276
class domain_greater:
 
277
    "domain_greater(v)(x) = true where x <= v"
 
278
    def __init__(self, critical_value):
 
279
        "domain_greater(v)(x) = true where x <= v"
 
280
        self.critical_value = critical_value
 
281
 
 
282
    def __call__ (self, x):
 
283
        "Execute the call behavior."
 
284
        return umath.less_equal (x, self.critical_value)
 
285
 
 
286
class domain_greater_equal:
 
287
    "domain_greater_equal(v)(x) = true where x < v"
 
288
    def __init__(self, critical_value):
 
289
        "domain_greater_equal(v)(x) = true where x < v"
 
290
        self.critical_value = critical_value
 
291
 
 
292
    def __call__ (self, x):
 
293
        "Execute the call behavior."
 
294
        return umath.less (x, self.critical_value)
 
295
 
 
296
class masked_unary_operation:
 
297
    def __init__ (self, aufunc, fill=0, domain=None):
 
298
        """ masked_unary_operation(aufunc, fill=0, domain=None)
 
299
            aufunc(fill) must be defined
 
300
            self(x) returns aufunc(x)
 
301
            with masked values where domain(x) is true or getmask(x) is true.
 
302
        """
 
303
        self.f = aufunc
 
304
        self.fill = fill
 
305
        self.domain = domain
 
306
        self.__doc__ = getattr(aufunc, "__doc__", str(aufunc))
 
307
        self.__name__ = getattr(aufunc, "__name__", str(aufunc))
 
308
        ufunc_domain[aufunc] = domain
 
309
        ufunc_fills[aufunc] = fill,
 
310
 
 
311
    def __call__ (self, a, *args, **kwargs):
 
312
        "Execute the call behavior."
 
313
# numeric tries to return scalars rather than arrays when given scalars.
 
314
        m = getmask(a)
 
315
        d1 = filled(a, self.fill)
 
316
        if self.domain is not None:
 
317
            m = mask_or(m, self.domain(d1))
 
318
        result = self.f(d1, *args, **kwargs)
 
319
        return masked_array(result, m)
 
320
 
 
321
    def __str__ (self):
 
322
        return "Masked version of " + str(self.f)
 
323
 
 
324
 
 
325
class domain_safe_divide:
 
326
    def __init__ (self, tolerance=divide_tolerance):
 
327
        self.tolerance = tolerance
 
328
    def __call__ (self, a, b):
 
329
        return umath.absolute(a) * self.tolerance >= umath.absolute(b)
 
330
 
 
331
class domained_binary_operation:
 
332
    """Binary operations that have a domain, like divide. These are complicated
 
333
       so they are a separate class. They have no reduce, outer or accumulate.
 
334
    """
 
335
    def __init__ (self, abfunc, domain, fillx=0, filly=0):
 
336
        """abfunc(fillx, filly) must be defined.
 
337
           abfunc(x, filly) = x for all x to enable reduce.
 
338
        """
 
339
        self.f = abfunc
 
340
        self.domain = domain
 
341
        self.fillx = fillx
 
342
        self.filly = filly
 
343
        self.__doc__ = getattr(abfunc, "__doc__", str(abfunc))
 
344
        self.__name__ = getattr(abfunc, "__name__", str(abfunc))
 
345
        ufunc_domain[abfunc] = domain
 
346
        ufunc_fills[abfunc] = fillx, filly
 
347
 
 
348
    def __call__(self, a, b):
 
349
        "Execute the call behavior."
 
350
        ma = getmask(a)
 
351
        mb = getmask(b)
 
352
        d1 = filled(a, self.fillx)
 
353
        d2 = filled(b, self.filly)
 
354
        t = self.domain(d1, d2)
 
355
 
 
356
        if fromnumeric.sometrue(t, None):
 
357
            d2 = where(t, self.filly, d2)
 
358
            mb = mask_or(mb, t)
 
359
        m = mask_or(ma, mb)
 
360
        result =  self.f(d1, d2)
 
361
        return masked_array(result, m)
 
362
 
 
363
    def __str__ (self):
 
364
        return "Masked version of " + str(self.f)
 
365
 
 
366
class masked_binary_operation:
 
367
    def __init__ (self, abfunc, fillx=0, filly=0):
 
368
        """abfunc(fillx, filly) must be defined.
 
369
           abfunc(x, filly) = x for all x to enable reduce.
 
370
        """
 
371
        self.f = abfunc
 
372
        self.fillx = fillx
 
373
        self.filly = filly
 
374
        self.__doc__ = getattr(abfunc, "__doc__", str(abfunc))
 
375
        ufunc_domain[abfunc] = None
 
376
        ufunc_fills[abfunc] = fillx, filly
 
377
 
 
378
    def __call__ (self, a, b, *args, **kwargs):
 
379
        "Execute the call behavior."
 
380
        m = mask_or(getmask(a), getmask(b))
 
381
        d1 = filled(a, self.fillx)
 
382
        d2 = filled(b, self.filly)
 
383
        result = self.f(d1, d2, *args, **kwargs)
 
384
        if isinstance(result, ndarray) \
 
385
               and m.ndim != 0 \
 
386
               and m.shape != result.shape:
 
387
            m = mask_or(getmaskarray(a), getmaskarray(b))
 
388
        return masked_array(result, m)
 
389
 
 
390
    def reduce (self, target, axis=0, dtype=None):
 
391
        """Reduce target along the given axis with this function."""
 
392
        m = getmask(target)
 
393
        t = filled(target, self.filly)
 
394
        if t.shape == ():
 
395
            t = t.reshape(1)
 
396
            if m is not nomask:
 
397
                m = make_mask(m, copy=1)
 
398
                m.shape = (1,)
 
399
        if m is nomask:
 
400
            t = self.f.reduce(t, axis)
 
401
        else:
 
402
            t = masked_array (t, m)
 
403
            # XXX: "or t.dtype" below is a workaround for what appears
 
404
            # XXX: to be a bug in reduce.
 
405
            t = self.f.reduce(filled(t, self.filly), axis,
 
406
                              dtype=dtype or t.dtype)
 
407
            m = umath.logical_and.reduce(m, axis)
 
408
        if isinstance(t, ndarray):
 
409
            return masked_array(t, m, get_fill_value(target))
 
410
        elif m:
 
411
            return masked
 
412
        else:
 
413
            return t
 
414
 
 
415
    def outer (self, a, b):
 
416
        "Return the function applied to the outer product of a and b."
 
417
        ma = getmask(a)
 
418
        mb = getmask(b)
 
419
        if ma is nomask and mb is nomask:
 
420
            m = nomask
 
421
        else:
 
422
            ma = getmaskarray(a)
 
423
            mb = getmaskarray(b)
 
424
            m = logical_or.outer(ma, mb)
 
425
        d = self.f.outer(filled(a, self.fillx), filled(b, self.filly))
 
426
        return masked_array(d, m)
 
427
 
 
428
    def accumulate (self, target, axis=0):
 
429
        """Accumulate target along axis after filling with y fill value."""
 
430
        t = filled(target, self.filly)
 
431
        return masked_array (self.f.accumulate (t, axis))
 
432
    def __str__ (self):
 
433
        return "Masked version of " + str(self.f)
 
434
 
 
435
sqrt = masked_unary_operation(umath.sqrt, 0.0, domain_greater_equal(0.0))
 
436
log = masked_unary_operation(umath.log, 1.0, domain_greater(0.0))
 
437
log10 = masked_unary_operation(umath.log10, 1.0, domain_greater(0.0))
 
438
exp = masked_unary_operation(umath.exp)
 
439
conjugate = masked_unary_operation(umath.conjugate)
 
440
sin = masked_unary_operation(umath.sin)
 
441
cos = masked_unary_operation(umath.cos)
 
442
tan = masked_unary_operation(umath.tan, 0.0, domain_tan(1.e-35))
 
443
arcsin = masked_unary_operation(umath.arcsin, 0.0, domain_check_interval(-1.0, 1.0))
 
444
arccos = masked_unary_operation(umath.arccos, 0.0, domain_check_interval(-1.0, 1.0))
 
445
arctan = masked_unary_operation(umath.arctan)
 
446
# Missing from numeric
 
447
arcsinh = masked_unary_operation(umath.arcsinh)
 
448
arccosh = masked_unary_operation(umath.arccosh, 1.0, domain_greater_equal(1.0))
 
449
arctanh = masked_unary_operation(umath.arctanh, 0.0, domain_check_interval(-1.0+1e-15, 1.0-1e-15))
 
450
sinh = masked_unary_operation(umath.sinh)
 
451
cosh = masked_unary_operation(umath.cosh)
 
452
tanh = masked_unary_operation(umath.tanh)
 
453
absolute = masked_unary_operation(umath.absolute)
 
454
fabs = masked_unary_operation(umath.fabs)
 
455
negative = masked_unary_operation(umath.negative)
 
456
 
 
457
def nonzero(a):
 
458
    """returns the indices of the elements of a which are not zero
 
459
    and not masked
 
460
    """
 
461
    return numeric.asarray(filled(a, 0).nonzero())
 
462
 
 
463
around = masked_unary_operation(fromnumeric.round_)
 
464
floor = masked_unary_operation(umath.floor)
 
465
ceil = masked_unary_operation(umath.ceil)
 
466
logical_not = masked_unary_operation(umath.logical_not)
 
467
 
 
468
add = masked_binary_operation(umath.add)
 
469
subtract = masked_binary_operation(umath.subtract)
 
470
subtract.reduce = None
 
471
multiply = masked_binary_operation(umath.multiply, 1, 1)
 
472
divide = domained_binary_operation(umath.divide, domain_safe_divide(), 0, 1)
 
473
true_divide = domained_binary_operation(umath.true_divide, domain_safe_divide(), 0, 1)
 
474
floor_divide = domained_binary_operation(umath.floor_divide, domain_safe_divide(), 0, 1)
 
475
remainder = domained_binary_operation(umath.remainder, domain_safe_divide(), 0, 1)
 
476
fmod = domained_binary_operation(umath.fmod, domain_safe_divide(), 0, 1)
 
477
hypot = masked_binary_operation(umath.hypot)
 
478
arctan2 = masked_binary_operation(umath.arctan2, 0.0, 1.0)
 
479
arctan2.reduce = None
 
480
equal = masked_binary_operation(umath.equal)
 
481
equal.reduce = None
 
482
not_equal = masked_binary_operation(umath.not_equal)
 
483
not_equal.reduce = None
 
484
less_equal = masked_binary_operation(umath.less_equal)
 
485
less_equal.reduce = None
 
486
greater_equal = masked_binary_operation(umath.greater_equal)
 
487
greater_equal.reduce = None
 
488
less = masked_binary_operation(umath.less)
 
489
less.reduce = None
 
490
greater = masked_binary_operation(umath.greater)
 
491
greater.reduce = None
 
492
logical_and = masked_binary_operation(umath.logical_and)
 
493
alltrue = masked_binary_operation(umath.logical_and, 1, 1).reduce
 
494
logical_or = masked_binary_operation(umath.logical_or)
 
495
sometrue = logical_or.reduce
 
496
logical_xor = masked_binary_operation(umath.logical_xor)
 
497
bitwise_and = masked_binary_operation(umath.bitwise_and)
 
498
bitwise_or = masked_binary_operation(umath.bitwise_or)
 
499
bitwise_xor = masked_binary_operation(umath.bitwise_xor)
 
500
 
 
501
def rank (object):
 
502
    return fromnumeric.rank(filled(object))
 
503
 
 
504
def shape (object):
 
505
    return fromnumeric.shape(filled(object))
 
506
 
 
507
def size (object, axis=None):
 
508
    return fromnumeric.size(filled(object), axis)
 
509
 
 
510
class MaskedArray (object):
 
511
    """Arrays with possibly masked values.
 
512
       Masked values of 1 exclude the corresponding element from
 
513
       any computation.
 
514
 
 
515
       Construction:
 
516
           x = array(data, dtype=None, copy=True, order=False,
 
517
                     mask = nomask, fill_value=None)
 
518
 
 
519
       If copy=False, every effort is made not to copy the data:
 
520
           If data is a MaskedArray, and argument mask=nomask,
 
521
           then the candidate data is data.data and the
 
522
           mask used is data.mask. If data is a numeric array,
 
523
           it is used as the candidate raw data.
 
524
           If dtype is not None and
 
525
           is != data.dtype.char then a data copy is required.
 
526
           Otherwise, the candidate is used.
 
527
 
 
528
       If a data copy is required, raw data stored is the result of:
 
529
       numeric.array(data, dtype=dtype.char, copy=copy)
 
530
 
 
531
       If mask is nomask there are no masked values. Otherwise mask must
 
532
       be convertible to an array of booleans with the same shape as x.
 
533
 
 
534
       fill_value is used to fill in masked values when necessary,
 
535
       such as when printing and in method/function filled().
 
536
       The fill_value is not used for computation within this module.
 
537
    """
 
538
    __array_priority__ = 10.1
 
539
    def __init__(self, data, dtype=None, copy=True, order=False,
 
540
                 mask=nomask, fill_value=None):
 
541
        """array(data, dtype=None, copy=True, order=False, mask=nomask, fill_value=None)
 
542
           If data already a numeric array, its dtype becomes the default value of dtype.
 
543
        """
 
544
        if dtype is None:
 
545
            tc = None
 
546
        else:
 
547
            tc = numeric.dtype(dtype)
 
548
        need_data_copied = copy
 
549
        if isinstance(data, MaskedArray):
 
550
            c = data.data
 
551
            if tc is None:
 
552
                tc = c.dtype
 
553
            elif tc != c.dtype:
 
554
                need_data_copied = True
 
555
            if mask is nomask:
 
556
                mask = data.mask
 
557
            elif mask is not nomask: #attempting to change the mask
 
558
                need_data_copied = True
 
559
 
 
560
        elif isinstance(data, ndarray):
 
561
            c = data
 
562
            if tc is None:
 
563
                tc = c.dtype
 
564
            elif tc != c.dtype:
 
565
                need_data_copied = True
 
566
        else:
 
567
            need_data_copied = False #because I'll do it now
 
568
            c = numeric.array(data, dtype=tc, copy=True, order=order)
 
569
            tc = c.dtype
 
570
 
 
571
        if need_data_copied:
 
572
            if tc == c.dtype:
 
573
                self._data = numeric.array(c, dtype=tc, copy=True, order=order)
 
574
            else:
 
575
                self._data = c.astype(tc)
 
576
        else:
 
577
            self._data = c
 
578
 
 
579
        if mask is nomask:
 
580
            self._mask = nomask
 
581
            self._shared_mask = 0
 
582
        else:
 
583
            self._mask = make_mask (mask)
 
584
            if self._mask is nomask:
 
585
                self._shared_mask = 0
 
586
            else:
 
587
                self._shared_mask = (self._mask is mask)
 
588
                nm = size(self._mask)
 
589
                nd = size(self._data)
 
590
                if nm != nd:
 
591
                    if nm == 1:
 
592
                        self._mask = fromnumeric.resize(self._mask, self._data.shape)
 
593
                        self._shared_mask = 0
 
594
                    elif nd == 1:
 
595
                        self._data = fromnumeric.resize(self._data, self._mask.shape)
 
596
                        self._data.shape = self._mask.shape
 
597
                    else:
 
598
                        raise MAError, "Mask and data not compatible."
 
599
                elif nm == 1 and shape(self._mask) != shape(self._data):
 
600
                    self.unshare_mask()
 
601
                    self._mask.shape = self._data.shape
 
602
 
 
603
        self.set_fill_value(fill_value)
 
604
 
 
605
    def __array__ (self, t=None, context=None):
 
606
        "Special hook for numeric. Converts to numeric if possible."
 
607
        if self._mask is not nomask:
 
608
            if fromnumeric.ravel(self._mask).any():
 
609
                if context is None:
 
610
                    warnings.warn("Cannot automatically convert masked array to "\
 
611
                                  "numeric because data\n    is masked in one or "\
 
612
                                  "more locations.");
 
613
                    return self._data
 
614
                    #raise MAError, \
 
615
                    #      """Cannot automatically convert masked array to numeric because data
 
616
                    #      is masked in one or more locations.
 
617
                    #      """
 
618
                else:
 
619
                    func, args, i = context
 
620
                    fills = ufunc_fills.get(func)
 
621
                    if fills is None:
 
622
                        raise MAError, "%s not known to ma" % func
 
623
                    return self.filled(fills[i])
 
624
            else:  # Mask is all false
 
625
                   # Optimize to avoid future invocations of this section.
 
626
                self._mask = nomask
 
627
                self._shared_mask = 0
 
628
        if t:
 
629
            return self._data.astype(t)
 
630
        else:
 
631
            return self._data
 
632
 
 
633
    def __array_wrap__ (self, array, context=None):
 
634
        """Special hook for ufuncs.
 
635
 
 
636
        Wraps the numpy array and sets the mask according to
 
637
        context.
 
638
        """
 
639
        if context is None:
 
640
            return MaskedArray(array, copy=False, mask=nomask)
 
641
        func, args = context[:2]
 
642
        domain = ufunc_domain[func]
 
643
        m = reduce(mask_or, [getmask(a) for a in args])
 
644
        if domain is not None:
 
645
            m = mask_or(m, domain(*[getattr(a, '_data', a)
 
646
                                    for a in args]))
 
647
        if m is not nomask:
 
648
            try:
 
649
                shape = array.shape
 
650
            except AttributeError:
 
651
                pass
 
652
            else:
 
653
                if m.shape != shape:
 
654
                    m = reduce(mask_or, [getmaskarray(a) for a in args])
 
655
 
 
656
        return MaskedArray(array, copy=False, mask=m)
 
657
 
 
658
    def _get_shape(self):
 
659
        "Return the current shape."
 
660
        return self._data.shape
 
661
 
 
662
    def _set_shape (self, newshape):
 
663
        "Set the array's shape."
 
664
        self._data.shape = newshape
 
665
        if self._mask is not nomask:
 
666
            self._mask = self._mask.copy()
 
667
            self._mask.shape = newshape
 
668
 
 
669
    def _get_flat(self):
 
670
        """Calculate the flat value.
 
671
        """
 
672
        if self._mask is nomask:
 
673
            return masked_array(self._data.ravel(), mask=nomask,
 
674
                                fill_value = self.fill_value())
 
675
        else:
 
676
            return masked_array(self._data.ravel(),
 
677
                                mask=self._mask.ravel(),
 
678
                                fill_value = self.fill_value())
 
679
 
 
680
    def _set_flat (self, value):
 
681
        "x.flat = value"
 
682
        y = self.ravel()
 
683
        y[:] = value
 
684
 
 
685
    def _get_real(self):
 
686
        "Get the real part of a complex array."
 
687
        if self._mask is nomask:
 
688
            return masked_array(self._data.real, mask=nomask,
 
689
                            fill_value = self.fill_value())
 
690
        else:
 
691
            return masked_array(self._data.real, mask=self._mask,
 
692
                            fill_value = self.fill_value())
 
693
 
 
694
    def _set_real (self, value):
 
695
        "x.real = value"
 
696
        y = self.real
 
697
        y[...] = value
 
698
 
 
699
    def _get_imaginary(self):
 
700
        "Get the imaginary part of a complex array."
 
701
        if self._mask is nomask:
 
702
            return masked_array(self._data.imag, mask=nomask,
 
703
                            fill_value = self.fill_value())
 
704
        else:
 
705
            return masked_array(self._data.imag, mask=self._mask,
 
706
                            fill_value = self.fill_value())
 
707
 
 
708
    def _set_imaginary (self, value):
 
709
        "x.imaginary = value"
 
710
        y = self.imaginary
 
711
        y[...] = value
 
712
 
 
713
    def __str__(self):
 
714
        """Calculate the str representation, using masked for fill if
 
715
           it is enabled. Otherwise fill with fill value.
 
716
        """
 
717
        if masked_print_option.enabled():
 
718
            f = masked_print_option
 
719
            # XXX: Without the following special case masked
 
720
            # XXX: would print as "[--]", not "--". Can we avoid
 
721
            # XXX: checks for masked by choosing a different value
 
722
            # XXX: for the masked singleton? 2005-01-05 -- sasha
 
723
            if self is masked:
 
724
                return str(f)
 
725
            m = self._mask
 
726
            if m is not nomask and m.shape == () and m:
 
727
                return str(f)
 
728
            # convert to object array to make filled work
 
729
            self = self.astype(object)
 
730
        else:
 
731
            f = self.fill_value()
 
732
        res = self.filled(f)
 
733
        return str(res)
 
734
 
 
735
    def __repr__(self):
 
736
        """Calculate the repr representation, using masked for fill if
 
737
           it is enabled. Otherwise fill with fill value.
 
738
        """
 
739
        with_mask = """\
 
740
array(data =
 
741
 %(data)s,
 
742
      mask =
 
743
 %(mask)s,
 
744
      fill_value=%(fill)s)
 
745
"""
 
746
        with_mask1 = """\
 
747
array(data = %(data)s,
 
748
      mask = %(mask)s,
 
749
      fill_value=%(fill)s)
 
750
"""
 
751
        without_mask = """array(
 
752
 %(data)s)"""
 
753
        without_mask1 = """array(%(data)s)"""
 
754
 
 
755
        n = len(self.shape)
 
756
        if self._mask is nomask:
 
757
            if n <= 1:
 
758
                return without_mask1 % {'data':str(self.filled())}
 
759
            return without_mask % {'data':str(self.filled())}
 
760
        else:
 
761
            if n <= 1:
 
762
                return with_mask % {
 
763
                    'data': str(self.filled()),
 
764
                    'mask': str(self._mask),
 
765
                    'fill': str(self.fill_value())
 
766
                    }
 
767
            return with_mask % {
 
768
                'data': str(self.filled()),
 
769
                'mask': str(self._mask),
 
770
                'fill': str(self.fill_value())
 
771
                }
 
772
        without_mask1 = """array(%(data)s)"""
 
773
        if self._mask is nomask:
 
774
            return without_mask % {'data':str(self.filled())}
 
775
        else:
 
776
            return with_mask % {
 
777
                'data': str(self.filled()),
 
778
                'mask': str(self._mask),
 
779
                'fill': str(self.fill_value())
 
780
                }
 
781
 
 
782
    def __float__(self):
 
783
        "Convert self to float."
 
784
        self.unmask()
 
785
        if self._mask is not nomask:
 
786
            raise MAError, 'Cannot convert masked element to a Python float.'
 
787
        return float(self.data.item())
 
788
 
 
789
    def __int__(self):
 
790
        "Convert self to int."
 
791
        self.unmask()
 
792
        if self._mask is not nomask:
 
793
            raise MAError, 'Cannot convert masked element to a Python int.'
 
794
        return int(self.data.item())
 
795
 
 
796
    def __getitem__(self, i):
 
797
        "Get item described by i. Not a copy as in previous versions."
 
798
        self.unshare_mask()
 
799
        m = self._mask
 
800
        dout = self._data[i]
 
801
        if m is nomask:
 
802
            try:
 
803
                if dout.size == 1:
 
804
                    return dout
 
805
                else:
 
806
                    return masked_array(dout, fill_value=self._fill_value)
 
807
            except AttributeError:
 
808
                return dout
 
809
        mi = m[i]
 
810
        if mi.size == 1:
 
811
            if mi:
 
812
                return masked
 
813
            else:
 
814
                return dout
 
815
        else:
 
816
            return masked_array(dout, mi, fill_value=self._fill_value)
 
817
 
 
818
# --------
 
819
# setitem and setslice notes
 
820
# note that if value is masked, it means to mask those locations.
 
821
# setting a value changes the mask to match the value in those locations.
 
822
 
 
823
    def __setitem__(self, index, value):
 
824
        "Set item described by index. If value is masked, mask those locations."
 
825
        d = self._data
 
826
        if self is masked:
 
827
            raise MAError, 'Cannot alter masked elements.'
 
828
        if value is masked:
 
829
            if self._mask is nomask:
 
830
                self._mask = make_mask_none(d.shape)
 
831
                self._shared_mask = False
 
832
            else:
 
833
                self.unshare_mask()
 
834
            self._mask[index] = True
 
835
            return
 
836
        m = getmask(value)
 
837
        value = filled(value).astype(d.dtype)
 
838
        d[index] = value
 
839
        if m is nomask:
 
840
            if self._mask is not nomask:
 
841
                self.unshare_mask()
 
842
                self._mask[index] = False
 
843
        else:
 
844
            if self._mask is nomask:
 
845
                self._mask = make_mask_none(d.shape)
 
846
                self._shared_mask = True
 
847
            else:
 
848
                self.unshare_mask()
 
849
            self._mask[index] = m
 
850
 
 
851
    def __nonzero__(self):
 
852
        """returns true if any element is non-zero or masked
 
853
 
 
854
        """
 
855
        # XXX: This changes bool conversion logic from MA.
 
856
        # XXX: In MA bool(a) == len(a) != 0, but in numpy
 
857
        # XXX: scalars do not have len
 
858
        m = self._mask
 
859
        d = self._data
 
860
        return bool(m is not nomask and m.any()
 
861
                    or d is not nomask and d.any())
 
862
 
 
863
    def __len__ (self):
 
864
        """Return length of first dimension. This is weird but Python's
 
865
         slicing behavior depends on it."""
 
866
        return len(self._data)
 
867
 
 
868
    def __and__(self, other):
 
869
        "Return bitwise_and"
 
870
        return bitwise_and(self, other)
 
871
 
 
872
    def __or__(self, other):
 
873
        "Return bitwise_or"
 
874
        return bitwise_or(self, other)
 
875
 
 
876
    def __xor__(self, other):
 
877
        "Return bitwise_xor"
 
878
        return bitwise_xor(self, other)
 
879
 
 
880
    __rand__ = __and__
 
881
    __ror__ = __or__
 
882
    __rxor__ = __xor__
 
883
 
 
884
    def __abs__(self):
 
885
        "Return absolute(self)"
 
886
        return absolute(self)
 
887
 
 
888
    def __neg__(self):
 
889
        "Return negative(self)"
 
890
        return negative(self)
 
891
 
 
892
    def __pos__(self):
 
893
        "Return array(self)"
 
894
        return array(self)
 
895
 
 
896
    def __add__(self, other):
 
897
        "Return add(self, other)"
 
898
        return add(self, other)
 
899
 
 
900
    __radd__ = __add__
 
901
 
 
902
    def __mod__ (self, other):
 
903
        "Return remainder(self, other)"
 
904
        return remainder(self, other)
 
905
 
 
906
    def __rmod__ (self, other):
 
907
        "Return remainder(other, self)"
 
908
        return remainder(other, self)
 
909
 
 
910
    def __lshift__ (self, n):
 
911
        return left_shift(self, n)
 
912
 
 
913
    def __rshift__ (self, n):
 
914
        return right_shift(self, n)
 
915
 
 
916
    def __sub__(self, other):
 
917
        "Return subtract(self, other)"
 
918
        return subtract(self, other)
 
919
 
 
920
    def __rsub__(self, other):
 
921
        "Return subtract(other, self)"
 
922
        return subtract(other, self)
 
923
 
 
924
    def __mul__(self, other):
 
925
        "Return multiply(self, other)"
 
926
        return multiply(self, other)
 
927
 
 
928
    __rmul__ = __mul__
 
929
 
 
930
    def __div__(self, other):
 
931
        "Return divide(self, other)"
 
932
        return divide(self, other)
 
933
 
 
934
    def __rdiv__(self, other):
 
935
        "Return divide(other, self)"
 
936
        return divide(other, self)
 
937
 
 
938
    def __truediv__(self, other):
 
939
        "Return divide(self, other)"
 
940
        return true_divide(self, other)
 
941
 
 
942
    def __rtruediv__(self, other):
 
943
        "Return divide(other, self)"
 
944
        return true_divide(other, self)
 
945
 
 
946
    def __floordiv__(self, other):
 
947
        "Return divide(self, other)"
 
948
        return floor_divide(self, other)
 
949
 
 
950
    def __rfloordiv__(self, other):
 
951
        "Return divide(other, self)"
 
952
        return floor_divide(other, self)
 
953
 
 
954
    def __pow__(self, other, third=None):
 
955
        "Return power(self, other, third)"
 
956
        return power(self, other, third)
 
957
 
 
958
    def __sqrt__(self):
 
959
        "Return sqrt(self)"
 
960
        return sqrt(self)
 
961
 
 
962
    def __iadd__(self, other):
 
963
        "Add other to self in place."
 
964
        t = self._data.dtype.char
 
965
        f = filled(other, 0)
 
966
        t1 = f.dtype.char
 
967
        if t == t1:
 
968
            pass
 
969
        elif t in typecodes['Integer']:
 
970
            if t1 in typecodes['Integer']:
 
971
                f = f.astype(t)
 
972
            else:
 
973
                raise TypeError, 'Incorrect type for in-place operation.'
 
974
        elif t in typecodes['Float']:
 
975
            if t1 in typecodes['Integer']:
 
976
                f = f.astype(t)
 
977
            elif t1 in typecodes['Float']:
 
978
                f = f.astype(t)
 
979
            else:
 
980
                raise TypeError, 'Incorrect type for in-place operation.'
 
981
        elif t in typecodes['Complex']:
 
982
            if t1 in typecodes['Integer']:
 
983
                f = f.astype(t)
 
984
            elif t1 in typecodes['Float']:
 
985
                f = f.astype(t)
 
986
            elif t1 in typecodes['Complex']:
 
987
                f = f.astype(t)
 
988
            else:
 
989
                raise TypeError, 'Incorrect type for in-place operation.'
 
990
        else:
 
991
            raise TypeError, 'Incorrect type for in-place operation.'
 
992
 
 
993
        if self._mask is nomask:
 
994
            self._data += f
 
995
            m = getmask(other)
 
996
            self._mask = m
 
997
            self._shared_mask = m is not nomask
 
998
        else:
 
999
            result = add(self, masked_array(f, mask=getmask(other)))
 
1000
            self._data = result.data
 
1001
            self._mask = result.mask
 
1002
            self._shared_mask = 1
 
1003
        return self
 
1004
 
 
1005
    def __imul__(self, other):
 
1006
        "Add other to self in place."
 
1007
        t = self._data.dtype.char
 
1008
        f = filled(other, 0)
 
1009
        t1 = f.dtype.char
 
1010
        if t == t1:
 
1011
            pass
 
1012
        elif t in typecodes['Integer']:
 
1013
            if t1 in typecodes['Integer']:
 
1014
                f = f.astype(t)
 
1015
            else:
 
1016
                raise TypeError, 'Incorrect type for in-place operation.'
 
1017
        elif t in typecodes['Float']:
 
1018
            if t1 in typecodes['Integer']:
 
1019
                f = f.astype(t)
 
1020
            elif t1 in typecodes['Float']:
 
1021
                f = f.astype(t)
 
1022
            else:
 
1023
                raise TypeError, 'Incorrect type for in-place operation.'
 
1024
        elif t in typecodes['Complex']:
 
1025
            if t1 in typecodes['Integer']:
 
1026
                f = f.astype(t)
 
1027
            elif t1 in typecodes['Float']:
 
1028
                f = f.astype(t)
 
1029
            elif t1 in typecodes['Complex']:
 
1030
                f = f.astype(t)
 
1031
            else:
 
1032
                raise TypeError, 'Incorrect type for in-place operation.'
 
1033
        else:
 
1034
            raise TypeError, 'Incorrect type for in-place operation.'
 
1035
 
 
1036
        if self._mask is nomask:
 
1037
            self._data *= f
 
1038
            m = getmask(other)
 
1039
            self._mask = m
 
1040
            self._shared_mask = m is not nomask
 
1041
        else:
 
1042
            result = multiply(self, masked_array(f, mask=getmask(other)))
 
1043
            self._data = result.data
 
1044
            self._mask = result.mask
 
1045
            self._shared_mask = 1
 
1046
        return self
 
1047
 
 
1048
    def __isub__(self, other):
 
1049
        "Subtract other from self in place."
 
1050
        t = self._data.dtype.char
 
1051
        f = filled(other, 0)
 
1052
        t1 = f.dtype.char
 
1053
        if t == t1:
 
1054
            pass
 
1055
        elif t in typecodes['Integer']:
 
1056
            if t1 in typecodes['Integer']:
 
1057
                f = f.astype(t)
 
1058
            else:
 
1059
                raise TypeError, 'Incorrect type for in-place operation.'
 
1060
        elif t in typecodes['Float']:
 
1061
            if t1 in typecodes['Integer']:
 
1062
                f = f.astype(t)
 
1063
            elif t1 in typecodes['Float']:
 
1064
                f = f.astype(t)
 
1065
            else:
 
1066
                raise TypeError, 'Incorrect type for in-place operation.'
 
1067
        elif t in typecodes['Complex']:
 
1068
            if t1 in typecodes['Integer']:
 
1069
                f = f.astype(t)
 
1070
            elif t1 in typecodes['Float']:
 
1071
                f = f.astype(t)
 
1072
            elif t1 in typecodes['Complex']:
 
1073
                f = f.astype(t)
 
1074
            else:
 
1075
                raise TypeError, 'Incorrect type for in-place operation.'
 
1076
        else:
 
1077
            raise TypeError, 'Incorrect type for in-place operation.'
 
1078
 
 
1079
        if self._mask is nomask:
 
1080
            self._data -= f
 
1081
            m = getmask(other)
 
1082
            self._mask = m
 
1083
            self._shared_mask = m is not nomask
 
1084
        else:
 
1085
            result = subtract(self, masked_array(f, mask=getmask(other)))
 
1086
            self._data = result.data
 
1087
            self._mask = result.mask
 
1088
            self._shared_mask = 1
 
1089
        return self
 
1090
 
 
1091
 
 
1092
 
 
1093
    def __idiv__(self, other):
 
1094
        "Divide self by other in place."
 
1095
        t = self._data.dtype.char
 
1096
        f = filled(other, 0)
 
1097
        t1 = f.dtype.char
 
1098
        if t == t1:
 
1099
            pass
 
1100
        elif t in typecodes['Integer']:
 
1101
            if t1 in typecodes['Integer']:
 
1102
                f = f.astype(t)
 
1103
            else:
 
1104
                raise TypeError, 'Incorrect type for in-place operation.'
 
1105
        elif t in typecodes['Float']:
 
1106
            if t1 in typecodes['Integer']:
 
1107
                f = f.astype(t)
 
1108
            elif t1 in typecodes['Float']:
 
1109
                f = f.astype(t)
 
1110
            else:
 
1111
                raise TypeError, 'Incorrect type for in-place operation.'
 
1112
        elif t in typecodes['Complex']:
 
1113
            if t1 in typecodes['Integer']:
 
1114
                f = f.astype(t)
 
1115
            elif t1 in typecodes['Float']:
 
1116
                f = f.astype(t)
 
1117
            elif t1 in typecodes['Complex']:
 
1118
                f = f.astype(t)
 
1119
            else:
 
1120
                raise TypeError, 'Incorrect type for in-place operation.'
 
1121
        else:
 
1122
            raise TypeError, 'Incorrect type for in-place operation.'
 
1123
        mo = getmask(other)
 
1124
        result = divide(self, masked_array(f, mask=mo))
 
1125
        self._data = result.data
 
1126
        dm = result.raw_mask()
 
1127
        if dm is not self._mask:
 
1128
            self._mask = dm
 
1129
            self._shared_mask = 1
 
1130
        return self
 
1131
 
 
1132
    def __eq__(self, other):
 
1133
        return equal(self,other)
 
1134
 
 
1135
    def __ne__(self, other):
 
1136
        return not_equal(self,other)
 
1137
 
 
1138
    def __lt__(self, other):
 
1139
        return less(self,other)
 
1140
 
 
1141
    def __le__(self, other):
 
1142
        return less_equal(self,other)
 
1143
 
 
1144
    def __gt__(self, other):
 
1145
        return greater(self,other)
 
1146
 
 
1147
    def __ge__(self, other):
 
1148
        return greater_equal(self,other)
 
1149
 
 
1150
    def astype (self, tc):
 
1151
        "return self as array of given type."
 
1152
        d = self._data.astype(tc)
 
1153
        return array(d, mask=self._mask)
 
1154
 
 
1155
    def byte_swapped(self):
 
1156
        """Returns the raw data field, byte_swapped. Included for consistency
 
1157
         with numeric but doesn't make sense in this context.
 
1158
        """
 
1159
        return self._data.byte_swapped()
 
1160
 
 
1161
    def compressed (self):
 
1162
        "A 1-D array of all the non-masked data."
 
1163
        d = fromnumeric.ravel(self._data)
 
1164
        if self._mask is nomask:
 
1165
            return array(d)
 
1166
        else:
 
1167
            m = 1 - fromnumeric.ravel(self._mask)
 
1168
            c = fromnumeric.compress(m, d)
 
1169
            return array(c, copy=0)
 
1170
 
 
1171
    def count (self, axis = None):
 
1172
        "Count of the non-masked elements in a, or along a certain axis."
 
1173
        m = self._mask
 
1174
        s = self._data.shape
 
1175
        ls = len(s)
 
1176
        if m is nomask:
 
1177
            if ls == 0:
 
1178
                return 1
 
1179
            if ls == 1:
 
1180
                return s[0]
 
1181
            if axis is None:
 
1182
                return reduce(lambda x, y:x*y, s)
 
1183
            else:
 
1184
                n = s[axis]
 
1185
                t = list(s)
 
1186
                del t[axis]
 
1187
                return ones(t) * n
 
1188
        if axis is None:
 
1189
            w = fromnumeric.ravel(m).astype(int)
 
1190
            n1 = size(w)
 
1191
            if n1 == 1:
 
1192
                n2 = w[0]
 
1193
            else:
 
1194
                n2 = umath.add.reduce(w)
 
1195
            return n1 - n2
 
1196
        else:
 
1197
            n1 = size(m, axis)
 
1198
            n2 = sum(m.astype(int), axis)
 
1199
            return n1 - n2
 
1200
 
 
1201
    def dot (self, other):
 
1202
        "s.dot(other) = innerproduct(s, other)"
 
1203
        return innerproduct(self, other)
 
1204
 
 
1205
    def fill_value(self):
 
1206
        "Get the current fill value."
 
1207
        return self._fill_value
 
1208
 
 
1209
    def filled (self, fill_value=None):
 
1210
        """A numeric array with masked values filled. If fill_value is None,
 
1211
           use self.fill_value().
 
1212
 
 
1213
           If mask is nomask, copy data only if not contiguous.
 
1214
           Result is always a contiguous, numeric array.
 
1215
# Is contiguous really necessary now?
 
1216
        """
 
1217
        d = self._data
 
1218
        m = self._mask
 
1219
        if m is nomask:
 
1220
            if d.flags['CONTIGUOUS']:
 
1221
                return d
 
1222
            else:
 
1223
                return d.copy()
 
1224
        else:
 
1225
            if fill_value is None:
 
1226
                value = self._fill_value
 
1227
            else:
 
1228
                value = fill_value
 
1229
 
 
1230
            if self is masked:
 
1231
                result = numeric.array(value)
 
1232
            else:
 
1233
                try:
 
1234
                    result = numeric.array(d, dtype=d.dtype, copy=1)
 
1235
                    result[m] = value
 
1236
                except (TypeError, AttributeError):
 
1237
                    #ok, can't put that value in here
 
1238
                    value = numeric.array(value, dtype=object)
 
1239
                    d = d.astype(object)
 
1240
                    result = fromnumeric.choose(m, (d, value))
 
1241
            return result
 
1242
 
 
1243
    def ids (self):
 
1244
        """Return the ids of the data and mask areas"""
 
1245
        return (id(self._data), id(self._mask))
 
1246
 
 
1247
    def iscontiguous (self):
 
1248
        "Is the data contiguous?"
 
1249
        return self._data.flags['CONTIGUOUS']
 
1250
 
 
1251
    def itemsize(self):
 
1252
        "Item size of each data item."
 
1253
        return self._data.itemsize
 
1254
 
 
1255
 
 
1256
    def outer(self, other):
 
1257
        "s.outer(other) = outerproduct(s, other)"
 
1258
        return outerproduct(self, other)
 
1259
 
 
1260
    def put (self, values):
 
1261
        """Set the non-masked entries of self to filled(values).
 
1262
           No change to mask
 
1263
        """
 
1264
        iota = numeric.arange(self.size)
 
1265
        d = self._data
 
1266
        if self._mask is nomask:
 
1267
            ind = iota
 
1268
        else:
 
1269
            ind = fromnumeric.compress(1 - self._mask, iota)
 
1270
        d[ind] =  filled(values).astype(d.dtype)
 
1271
 
 
1272
    def putmask (self, values):
 
1273
        """Set the masked entries of self to filled(values).
 
1274
           Mask changed to nomask.
 
1275
        """
 
1276
        d = self._data
 
1277
        if self._mask is not nomask:
 
1278
            d[self._mask] = filled(values).astype(d.dtype)
 
1279
            self._shared_mask = 0
 
1280
            self._mask = nomask
 
1281
 
 
1282
    def ravel (self):
 
1283
        """Return a 1-D view of self."""
 
1284
        if self._mask is nomask:
 
1285
            return masked_array(self._data.ravel())
 
1286
        else:
 
1287
            return masked_array(self._data.ravel(), self._mask.ravel())
 
1288
 
 
1289
    def raw_data (self):
 
1290
        """ Obsolete; use data property instead.
 
1291
            The raw data; portions may be meaningless.
 
1292
            May be noncontiguous. Expert use only."""
 
1293
        return self._data
 
1294
    data = property(fget=raw_data,
 
1295
           doc="The data, but values at masked locations are meaningless.")
 
1296
 
 
1297
    def raw_mask (self):
 
1298
        """ Obsolete; use mask property instead.
 
1299
            May be noncontiguous. Expert use only.
 
1300
        """
 
1301
        return self._mask
 
1302
    mask = property(fget=raw_mask,
 
1303
           doc="The mask, may be nomask. Values where mask true are meaningless.")
 
1304
 
 
1305
    def reshape (self, *s):
 
1306
        """This array reshaped to shape s"""
 
1307
        d = self._data.reshape(*s)
 
1308
        if self._mask is nomask:
 
1309
            return masked_array(d)
 
1310
        else:
 
1311
            m = self._mask.reshape(*s)
 
1312
        return masked_array(d, m)
 
1313
 
 
1314
    def set_fill_value (self, v=None):
 
1315
        "Set the fill value to v. Omit v to restore default."
 
1316
        if v is None:
 
1317
            v = default_fill_value (self.raw_data())
 
1318
        self._fill_value = v
 
1319
 
 
1320
    def _get_ndim(self):
 
1321
        return self._data.ndim
 
1322
    ndim = property(_get_ndim, doc=numeric.ndarray.ndim.__doc__)
 
1323
 
 
1324
    def _get_size (self):
 
1325
        return self._data.size
 
1326
    size = property(fget=_get_size, doc="Number of elements in the array.")
 
1327
## CHECK THIS: signature of numeric.array.size?
 
1328
 
 
1329
    def _get_dtype(self):
 
1330
        return self._data.dtype
 
1331
    dtype = property(fget=_get_dtype, doc="type of the array elements.")
 
1332
 
 
1333
    def item(self, *args):
 
1334
        "Return Python scalar if possible"
 
1335
        if self._mask is not nomask:
 
1336
            m = self._mask.item(*args)
 
1337
            try:
 
1338
                if m[0]:
 
1339
                    return masked
 
1340
            except IndexError:
 
1341
                return masked
 
1342
        return self._data.item(*args)
 
1343
 
 
1344
    def itemset(self, *args):
 
1345
        "Set Python scalar into array"
 
1346
        item = args[-1]
 
1347
        args = args[:-1]
 
1348
        self[args] = item
 
1349
 
 
1350
    def tolist(self, fill_value=None):
 
1351
        "Convert to list"
 
1352
        return self.filled(fill_value).tolist()
 
1353
 
 
1354
    def tostring(self, fill_value=None):
 
1355
        "Convert to string"
 
1356
        return self.filled(fill_value).tostring()
 
1357
 
 
1358
    def unmask (self):
 
1359
        "Replace the mask by nomask if possible."
 
1360
        if self._mask is nomask: return
 
1361
        m = make_mask(self._mask, flag=1)
 
1362
        if m is nomask:
 
1363
            self._mask = nomask
 
1364
            self._shared_mask = 0
 
1365
 
 
1366
    def unshare_mask (self):
 
1367
        "If currently sharing mask, make a copy."
 
1368
        if self._shared_mask:
 
1369
            self._mask = make_mask (self._mask, copy=1, flag=0)
 
1370
            self._shared_mask = 0
 
1371
 
 
1372
    def _get_ctypes(self):
 
1373
        return self._data.ctypes
 
1374
 
 
1375
    def _get_T(self):
 
1376
        if (self.ndim < 2):
 
1377
            return self
 
1378
        return self.transpose()
 
1379
 
 
1380
    shape = property(_get_shape, _set_shape,
 
1381
           doc = 'tuple giving the shape of the array')
 
1382
 
 
1383
    flat = property(_get_flat, _set_flat,
 
1384
           doc = 'Access array in flat form.')
 
1385
 
 
1386
    real = property(_get_real, _set_real,
 
1387
           doc = 'Access the real part of the array')
 
1388
 
 
1389
    imaginary = property(_get_imaginary, _set_imaginary,
 
1390
           doc = 'Access the imaginary part of the array')
 
1391
 
 
1392
    imag = imaginary
 
1393
 
 
1394
    ctypes = property(_get_ctypes, None, doc="ctypes")
 
1395
 
 
1396
    T = property(_get_T, None, doc="get transpose")
 
1397
 
 
1398
#end class MaskedArray
 
1399
 
 
1400
array = MaskedArray
 
1401
 
 
1402
def isMaskedArray (x):
 
1403
    "Is x a masked array, that is, an instance of MaskedArray?"
 
1404
    return isinstance(x, MaskedArray)
 
1405
 
 
1406
isarray = isMaskedArray
 
1407
isMA = isMaskedArray  #backward compatibility
 
1408
 
 
1409
def allclose (a, b, fill_value=1, rtol=1.e-5, atol=1.e-8):
 
1410
    """ Returns true if all components of a and b are equal
 
1411
        subject to given tolerances.
 
1412
        If fill_value is 1, masked values considered equal.
 
1413
        If fill_value is 0, masked values considered unequal.
 
1414
        The relative error rtol should be positive and << 1.0
 
1415
        The absolute error atol comes into play for those elements
 
1416
        of b that are very small or zero; it says how small a must be also.
 
1417
    """
 
1418
    m = mask_or(getmask(a), getmask(b))
 
1419
    d1 = filled(a)
 
1420
    d2 = filled(b)
 
1421
    x = filled(array(d1, copy=0, mask=m), fill_value).astype(float)
 
1422
    y = filled(array(d2, copy=0, mask=m), 1).astype(float)
 
1423
    d = umath.less_equal(umath.absolute(x-y), atol + rtol * umath.absolute(y))
 
1424
    return fromnumeric.alltrue(fromnumeric.ravel(d))
 
1425
 
 
1426
def allequal (a, b, fill_value=1):
 
1427
    """
 
1428
        True if all entries of  a and b are equal, using
 
1429
        fill_value as a truth value where either or both are masked.
 
1430
    """
 
1431
    m = mask_or(getmask(a), getmask(b))
 
1432
    if m is nomask:
 
1433
        x = filled(a)
 
1434
        y = filled(b)
 
1435
        d = umath.equal(x, y)
 
1436
        return fromnumeric.alltrue(fromnumeric.ravel(d))
 
1437
    elif fill_value:
 
1438
        x = filled(a)
 
1439
        y = filled(b)
 
1440
        d = umath.equal(x, y)
 
1441
        dm = array(d, mask=m, copy=0)
 
1442
        return fromnumeric.alltrue(fromnumeric.ravel(filled(dm, 1)))
 
1443
    else:
 
1444
        return 0
 
1445
 
 
1446
def masked_values (data, value, rtol=1.e-5, atol=1.e-8, copy=1):
 
1447
    """
 
1448
       masked_values(data, value, rtol=1.e-5, atol=1.e-8)
 
1449
       Create a masked array; mask is nomask if possible.
 
1450
       If copy==0, and otherwise possible, result
 
1451
       may share data values with original array.
 
1452
       Let d = filled(data, value). Returns d
 
1453
       masked where abs(data-value)<= atol + rtol * abs(value)
 
1454
       if d is of a floating point type. Otherwise returns
 
1455
       masked_object(d, value, copy)
 
1456
    """
 
1457
    abs = umath.absolute
 
1458
    d = filled(data, value)
 
1459
    if issubclass(d.dtype.type, numeric.floating):
 
1460
        m = umath.less_equal(abs(d-value), atol+rtol*abs(value))
 
1461
        m = make_mask(m, flag=1)
 
1462
        return array(d, mask = m, copy=copy,
 
1463
                      fill_value=value)
 
1464
    else:
 
1465
        return masked_object(d, value, copy=copy)
 
1466
 
 
1467
def masked_object (data, value, copy=1):
 
1468
    "Create array masked where exactly data equal to value"
 
1469
    d = filled(data, value)
 
1470
    dm = make_mask(umath.equal(d, value), flag=1)
 
1471
    return array(d, mask=dm, copy=copy, fill_value=value)
 
1472
 
 
1473
def arange(start, stop=None, step=1, dtype=None):
 
1474
    """Just like range() except it returns a array whose type can be specified
 
1475
    by the keyword argument dtype.
 
1476
    """
 
1477
    return array(numeric.arange(start, stop, step, dtype))
 
1478
 
 
1479
arrayrange = arange
 
1480
 
 
1481
def fromstring (s, t):
 
1482
    "Construct a masked array from a string. Result will have no mask."
 
1483
    return masked_array(numeric.fromstring(s, t))
 
1484
 
 
1485
def left_shift (a, n):
 
1486
    "Left shift n bits"
 
1487
    m = getmask(a)
 
1488
    if m is nomask:
 
1489
        d = umath.left_shift(filled(a), n)
 
1490
        return masked_array(d)
 
1491
    else:
 
1492
        d = umath.left_shift(filled(a, 0), n)
 
1493
        return masked_array(d, m)
 
1494
 
 
1495
def right_shift (a, n):
 
1496
    "Right shift n bits"
 
1497
    m = getmask(a)
 
1498
    if m is nomask:
 
1499
        d = umath.right_shift(filled(a), n)
 
1500
        return masked_array(d)
 
1501
    else:
 
1502
        d = umath.right_shift(filled(a, 0), n)
 
1503
        return masked_array(d, m)
 
1504
 
 
1505
def resize (a, new_shape):
 
1506
    """resize(a, new_shape) returns a new array with the specified shape.
 
1507
    The original array's total size can be any size."""
 
1508
    m = getmask(a)
 
1509
    if m is not nomask:
 
1510
        m = fromnumeric.resize(m, new_shape)
 
1511
    result = array(fromnumeric.resize(filled(a), new_shape), mask=m)
 
1512
    result.set_fill_value(get_fill_value(a))
 
1513
    return result
 
1514
 
 
1515
def new_repeat(a, repeats, axis=None):
 
1516
    """repeat elements of a repeats times along axis
 
1517
       repeats is a sequence of length a.shape[axis]
 
1518
       telling how many times to repeat each element.
 
1519
    """
 
1520
    af = filled(a)
 
1521
    if isinstance(repeats, types.IntType):
 
1522
        if axis is None:
 
1523
            num = af.size
 
1524
        else:
 
1525
            num = af.shape[axis]
 
1526
        repeats = tuple([repeats]*num)
 
1527
 
 
1528
    m = getmask(a)
 
1529
    if m is not nomask:
 
1530
        m = fromnumeric.repeat(m, repeats, axis)
 
1531
    d = fromnumeric.repeat(af, repeats, axis)
 
1532
    result = masked_array(d, m)
 
1533
    result.set_fill_value(get_fill_value(a))
 
1534
    return result
 
1535
 
 
1536
 
 
1537
 
 
1538
def identity(n):
 
1539
    """identity(n) returns the identity matrix of shape n x n.
 
1540
    """
 
1541
    return array(numeric.identity(n))
 
1542
 
 
1543
def indices (dimensions, dtype=None):
 
1544
    """indices(dimensions,dtype=None) returns an array representing a grid
 
1545
    of indices with row-only, and column-only variation.
 
1546
    """
 
1547
    return array(numeric.indices(dimensions, dtype))
 
1548
 
 
1549
def zeros (shape, dtype=float):
 
1550
    """zeros(n, dtype=float) =
 
1551
     an array of all zeros of the given length or shape."""
 
1552
    return array(numeric.zeros(shape, dtype))
 
1553
 
 
1554
def ones (shape, dtype=float):
 
1555
    """ones(n, dtype=float) =
 
1556
     an array of all ones of the given length or shape."""
 
1557
    return array(numeric.ones(shape, dtype))
 
1558
 
 
1559
def count (a, axis = None):
 
1560
    "Count of the non-masked elements in a, or along a certain axis."
 
1561
    a = masked_array(a)
 
1562
    return a.count(axis)
 
1563
 
 
1564
def power (a, b, third=None):
 
1565
    "a**b"
 
1566
    if third is not None:
 
1567
        raise MAError, "3-argument power not supported."
 
1568
    ma = getmask(a)
 
1569
    mb = getmask(b)
 
1570
    m = mask_or(ma, mb)
 
1571
    fa = filled(a, 1)
 
1572
    fb = filled(b, 1)
 
1573
    if fb.dtype.char in typecodes["Integer"]:
 
1574
        return masked_array(umath.power(fa, fb), m)
 
1575
    md = make_mask(umath.less(fa, 0), flag=1)
 
1576
    m = mask_or(m, md)
 
1577
    if m is nomask:
 
1578
        return masked_array(umath.power(fa, fb))
 
1579
    else:
 
1580
        fa = numeric.where(m, 1, fa)
 
1581
        return masked_array(umath.power(fa, fb), m)
 
1582
 
 
1583
def masked_array (a, mask=nomask, fill_value=None):
 
1584
    """masked_array(a, mask=nomask) =
 
1585
       array(a, mask=mask, copy=0, fill_value=fill_value)
 
1586
    """
 
1587
    return array(a, mask=mask, copy=0, fill_value=fill_value)
 
1588
 
 
1589
def sum (target, axis=None, dtype=None):
 
1590
    if axis is None:
 
1591
        target = ravel(target)
 
1592
        axis = 0
 
1593
    return add.reduce(target, axis, dtype)
 
1594
 
 
1595
def product (target, axis=None, dtype=None):
 
1596
    if axis is None:
 
1597
        target = ravel(target)
 
1598
        axis = 0
 
1599
    return multiply.reduce(target, axis, dtype)
 
1600
 
 
1601
def new_average (a, axis=None, weights=None, returned = 0):
 
1602
    """average(a, axis=None, weights=None)
 
1603
       Computes average along indicated axis.
 
1604
       If axis is None, average over the entire array
 
1605
       Inputs can be integer or floating types; result is of type float.
 
1606
 
 
1607
       If weights are given, result is sum(a*weights,axis=0)/(sum(weights,axis=0)*1.0)
 
1608
       weights must have a's shape or be the 1-d with length the size
 
1609
       of a in the given axis.
 
1610
 
 
1611
       If returned, return a tuple: the result and the sum of the weights
 
1612
       or count of values. Results will have the same shape.
 
1613
 
 
1614
       masked values in the weights will be set to 0.0
 
1615
    """
 
1616
    a = masked_array(a)
 
1617
    mask = a.mask
 
1618
    ash = a.shape
 
1619
    if ash == ():
 
1620
        ash = (1,)
 
1621
    if axis is None:
 
1622
        if mask is nomask:
 
1623
            if weights is None:
 
1624
                n = add.reduce(a.raw_data().ravel())
 
1625
                d = reduce(lambda x, y: x * y, ash, 1.0)
 
1626
            else:
 
1627
                w = filled(weights, 0.0).ravel()
 
1628
                n = umath.add.reduce(a.raw_data().ravel() * w)
 
1629
                d = umath.add.reduce(w)
 
1630
                del w
 
1631
        else:
 
1632
            if weights is None:
 
1633
                n = add.reduce(a.ravel())
 
1634
                w = fromnumeric.choose(mask, (1.0, 0.0)).ravel()
 
1635
                d = umath.add.reduce(w)
 
1636
                del w
 
1637
            else:
 
1638
                w = array(filled(weights, 0.0), float, mask=mask).ravel()
 
1639
                n = add.reduce(a.ravel() * w)
 
1640
                d = add.reduce(w)
 
1641
                del w
 
1642
    else:
 
1643
        if mask is nomask:
 
1644
            if weights is None:
 
1645
                d = ash[axis] * 1.0
 
1646
                n = umath.add.reduce(a.raw_data(), axis)
 
1647
            else:
 
1648
                w = filled(weights, 0.0)
 
1649
                wsh = w.shape
 
1650
                if wsh == ():
 
1651
                    wsh = (1,)
 
1652
                if wsh == ash:
 
1653
                    w = numeric.array(w, float, copy=0)
 
1654
                    n = add.reduce(a*w, axis)
 
1655
                    d = add.reduce(w, axis)
 
1656
                    del w
 
1657
                elif wsh == (ash[axis],):
 
1658
                    r = [newaxis]*len(ash)
 
1659
                    r[axis] = slice(None, None, 1)
 
1660
                    w = eval ("w["+ repr(tuple(r)) + "] * ones(ash, float)")
 
1661
                    n = add.reduce(a*w, axis)
 
1662
                    d = add.reduce(w, axis)
 
1663
                    del w, r
 
1664
                else:
 
1665
                    raise ValueError, 'average: weights wrong shape.'
 
1666
        else:
 
1667
            if weights is None:
 
1668
                n = add.reduce(a, axis)
 
1669
                w = numeric.choose(mask, (1.0, 0.0))
 
1670
                d = umath.add.reduce(w, axis)
 
1671
                del w
 
1672
            else:
 
1673
                w = filled(weights, 0.0)
 
1674
                wsh = w.shape
 
1675
                if wsh == ():
 
1676
                    wsh = (1,)
 
1677
                if wsh == ash:
 
1678
                    w = array(w, float, mask=mask, copy=0)
 
1679
                    n = add.reduce(a*w, axis)
 
1680
                    d = add.reduce(w, axis)
 
1681
                elif wsh == (ash[axis],):
 
1682
                    r = [newaxis]*len(ash)
 
1683
                    r[axis] = slice(None, None, 1)
 
1684
                    w = eval ("w["+ repr(tuple(r)) + "] * masked_array(ones(ash, float), mask)")
 
1685
                    n = add.reduce(a*w, axis)
 
1686
                    d = add.reduce(w, axis)
 
1687
                else:
 
1688
                    raise ValueError, 'average: weights wrong shape.'
 
1689
                del w
 
1690
    #print n, d, repr(mask), repr(weights)
 
1691
    if n is masked or d is masked: return masked
 
1692
    result = divide (n, d)
 
1693
    del n
 
1694
 
 
1695
    if isinstance(result, MaskedArray):
 
1696
        result.unmask()
 
1697
        if returned:
 
1698
            if not isinstance(d, MaskedArray):
 
1699
                d = masked_array(d)
 
1700
            if not d.shape == result.shape:
 
1701
                d = ones(result.shape, float) * d
 
1702
            d.unmask()
 
1703
    if returned:
 
1704
        return result, d
 
1705
    else:
 
1706
        return result
 
1707
 
 
1708
def where (condition, x, y):
 
1709
    """where(condition, x, y) is x where condition is nonzero, y otherwise.
 
1710
       condition must be convertible to an integer array.
 
1711
       Answer is always the shape of condition.
 
1712
       The type depends on x and y. It is integer if both x and y are
 
1713
       the value masked.
 
1714
    """
 
1715
    fc = filled(not_equal(condition, 0), 0)
 
1716
    xv = filled(x)
 
1717
    xm = getmask(x)
 
1718
    yv = filled(y)
 
1719
    ym = getmask(y)
 
1720
    d = numeric.choose(fc, (yv, xv))
 
1721
    md = numeric.choose(fc, (ym, xm))
 
1722
    m = getmask(condition)
 
1723
    m = make_mask(mask_or(m, md), copy=0, flag=1)
 
1724
    return masked_array(d, m)
 
1725
 
 
1726
def choose (indices, t, out=None, mode='raise'):
 
1727
    "Returns array shaped like indices with elements chosen from t"
 
1728
    def fmask (x):
 
1729
        if x is masked: return 1
 
1730
        return filled(x)
 
1731
    def nmask (x):
 
1732
        if x is masked: return 1
 
1733
        m = getmask(x)
 
1734
        if m is nomask: return 0
 
1735
        return m
 
1736
    c = filled(indices, 0)
 
1737
    masks = [nmask(x) for x in t]
 
1738
    a = [fmask(x) for x in t]
 
1739
    d = numeric.choose(c, a)
 
1740
    m = numeric.choose(c, masks)
 
1741
    m = make_mask(mask_or(m, getmask(indices)), copy=0, flag=1)
 
1742
    return masked_array(d, m)
 
1743
 
 
1744
def masked_where(condition, x, copy=1):
 
1745
    """Return x as an array masked where condition is true.
 
1746
       Also masked where x or condition masked.
 
1747
    """
 
1748
    cm = filled(condition,1)
 
1749
    m = mask_or(getmask(x), cm)
 
1750
    return array(filled(x), copy=copy, mask=m)
 
1751
 
 
1752
def masked_greater(x, value, copy=1):
 
1753
    "masked_greater(x, value) = x masked where x > value"
 
1754
    return masked_where(greater(x, value), x, copy)
 
1755
 
 
1756
def masked_greater_equal(x, value, copy=1):
 
1757
    "masked_greater_equal(x, value) = x masked where x >= value"
 
1758
    return masked_where(greater_equal(x, value), x, copy)
 
1759
 
 
1760
def masked_less(x, value, copy=1):
 
1761
    "masked_less(x, value) = x masked where x < value"
 
1762
    return masked_where(less(x, value), x, copy)
 
1763
 
 
1764
def masked_less_equal(x, value, copy=1):
 
1765
    "masked_less_equal(x, value) = x masked where x <= value"
 
1766
    return masked_where(less_equal(x, value), x, copy)
 
1767
 
 
1768
def masked_not_equal(x, value, copy=1):
 
1769
    "masked_not_equal(x, value) = x masked where x != value"
 
1770
    d = filled(x, 0)
 
1771
    c = umath.not_equal(d, value)
 
1772
    m = mask_or(c, getmask(x))
 
1773
    return array(d, mask=m, copy=copy)
 
1774
 
 
1775
def masked_equal(x, value, copy=1):
 
1776
    """masked_equal(x, value) = x masked where x == value
 
1777
       For floating point consider masked_values(x, value) instead.
 
1778
    """
 
1779
    d = filled(x, 0)
 
1780
    c = umath.equal(d, value)
 
1781
    m = mask_or(c, getmask(x))
 
1782
    return array(d, mask=m, copy=copy)
 
1783
 
 
1784
def masked_inside(x, v1, v2, copy=1):
 
1785
    """x with mask of all values of x that are inside [v1,v2]
 
1786
       v1 and v2 can be given in either order.
 
1787
    """
 
1788
    if v2 < v1:
 
1789
        t = v2
 
1790
        v2 = v1
 
1791
        v1 = t
 
1792
    d = filled(x, 0)
 
1793
    c = umath.logical_and(umath.less_equal(d, v2), umath.greater_equal(d, v1))
 
1794
    m = mask_or(c, getmask(x))
 
1795
    return array(d, mask = m, copy=copy)
 
1796
 
 
1797
def masked_outside(x, v1, v2, copy=1):
 
1798
    """x with mask of all values of x that are outside [v1,v2]
 
1799
       v1 and v2 can be given in either order.
 
1800
    """
 
1801
    if v2 < v1:
 
1802
        t = v2
 
1803
        v2 = v1
 
1804
        v1 = t
 
1805
    d = filled(x, 0)
 
1806
    c = umath.logical_or(umath.less(d, v1), umath.greater(d, v2))
 
1807
    m = mask_or(c, getmask(x))
 
1808
    return array(d, mask = m, copy=copy)
 
1809
 
 
1810
def reshape (a, *newshape):
 
1811
    "Copy of a with a new shape."
 
1812
    m = getmask(a)
 
1813
    d = filled(a).reshape(*newshape)
 
1814
    if m is nomask:
 
1815
        return masked_array(d)
 
1816
    else:
 
1817
        return masked_array(d, mask=numeric.reshape(m, *newshape))
 
1818
 
 
1819
def ravel (a):
 
1820
    "a as one-dimensional, may share data and mask"
 
1821
    m = getmask(a)
 
1822
    d = fromnumeric.ravel(filled(a))
 
1823
    if m is nomask:
 
1824
        return masked_array(d)
 
1825
    else:
 
1826
        return masked_array(d, mask=numeric.ravel(m))
 
1827
 
 
1828
def concatenate (arrays, axis=0):
 
1829
    "Concatenate the arrays along the given axis"
 
1830
    d = []
 
1831
    for x in arrays:
 
1832
        d.append(filled(x))
 
1833
    d = numeric.concatenate(d, axis)
 
1834
    for x in arrays:
 
1835
        if getmask(x) is not nomask: break
 
1836
    else:
 
1837
        return masked_array(d)
 
1838
    dm = []
 
1839
    for x in arrays:
 
1840
        dm.append(getmaskarray(x))
 
1841
    dm = numeric.concatenate(dm, axis)
 
1842
    return masked_array(d, mask=dm)
 
1843
 
 
1844
def swapaxes (a, axis1, axis2):
 
1845
    m = getmask(a)
 
1846
    d = masked_array(a).data
 
1847
    if m is nomask:
 
1848
        return masked_array(data=numeric.swapaxes(d, axis1, axis2))
 
1849
    else:
 
1850
        return masked_array(data=numeric.swapaxes(d, axis1, axis2),
 
1851
                            mask=numeric.swapaxes(m, axis1, axis2),)
 
1852
 
 
1853
 
 
1854
def new_take (a, indices, axis=None, out=None, mode='raise'):
 
1855
    "returns selection of items from a."
 
1856
    m = getmask(a)
 
1857
    # d = masked_array(a).raw_data()
 
1858
    d = masked_array(a).data
 
1859
    if m is nomask:
 
1860
        return masked_array(numeric.take(d, indices, axis))
 
1861
    else:
 
1862
        return masked_array(numeric.take(d, indices, axis),
 
1863
                     mask = numeric.take(m, indices, axis))
 
1864
 
 
1865
def transpose(a, axes=None):
 
1866
    "reorder dimensions per tuple axes"
 
1867
    m = getmask(a)
 
1868
    d = filled(a)
 
1869
    if m is nomask:
 
1870
        return masked_array(numeric.transpose(d, axes))
 
1871
    else:
 
1872
        return masked_array(numeric.transpose(d, axes),
 
1873
                     mask = numeric.transpose(m, axes))
 
1874
 
 
1875
 
 
1876
def put(a, indices, values, mode='raise'):
 
1877
    """sets storage-indexed locations to corresponding values.
 
1878
 
 
1879
    Values and indices are filled if necessary.
 
1880
 
 
1881
    """
 
1882
    d = a.raw_data()
 
1883
    ind = filled(indices)
 
1884
    v = filled(values)
 
1885
    numeric.put (d, ind, v)
 
1886
    m = getmask(a)
 
1887
    if m is not nomask:
 
1888
        a.unshare_mask()
 
1889
        numeric.put(a.raw_mask(), ind, 0)
 
1890
 
 
1891
def putmask(a, mask, values):
 
1892
    "putmask(a, mask, values) sets a where mask is true."
 
1893
    if mask is nomask:
 
1894
        return
 
1895
    numeric.putmask(a.raw_data(), mask, values)
 
1896
    m = getmask(a)
 
1897
    if m is nomask: return
 
1898
    a.unshare_mask()
 
1899
    numeric.putmask(a.raw_mask(), mask, 0)
 
1900
 
 
1901
def inner(a, b):
 
1902
    """inner(a,b) returns the dot product of two arrays, which has
 
1903
    shape a.shape[:-1] + b.shape[:-1] with elements computed by summing the
 
1904
    product of the elements from the last dimensions of a and b.
 
1905
    Masked elements are replace by zeros.
 
1906
    """
 
1907
    fa = filled(a, 0)
 
1908
    fb = filled(b, 0)
 
1909
    if len(fa.shape) == 0: fa.shape = (1,)
 
1910
    if len(fb.shape) == 0: fb.shape = (1,)
 
1911
    return masked_array(numeric.inner(fa, fb))
 
1912
 
 
1913
innerproduct = inner
 
1914
 
 
1915
def outer(a, b):
 
1916
    """outer(a,b) = {a[i]*b[j]}, has shape (len(a),len(b))"""
 
1917
    fa = filled(a, 0).ravel()
 
1918
    fb = filled(b, 0).ravel()
 
1919
    d = numeric.outer(fa, fb)
 
1920
    ma = getmask(a)
 
1921
    mb = getmask(b)
 
1922
    if ma is nomask and mb is nomask:
 
1923
        return masked_array(d)
 
1924
    ma = getmaskarray(a)
 
1925
    mb = getmaskarray(b)
 
1926
    m = make_mask(1-numeric.outer(1-ma, 1-mb), copy=0)
 
1927
    return masked_array(d, m)
 
1928
 
 
1929
outerproduct = outer
 
1930
 
 
1931
def dot(a, b):
 
1932
    """dot(a,b) returns matrix-multiplication between a and b.  The product-sum
 
1933
    is over the last dimension of a and the second-to-last dimension of b.
 
1934
    Masked values are replaced by zeros. See also innerproduct.
 
1935
    """
 
1936
    return innerproduct(filled(a, 0), numeric.swapaxes(filled(b, 0), -1, -2))
 
1937
 
 
1938
def compress(condition, x, dimension=-1, out=None):
 
1939
    """Select those parts of x for which condition is true.
 
1940
       Masked values in condition are considered false.
 
1941
    """
 
1942
    c = filled(condition, 0)
 
1943
    m = getmask(x)
 
1944
    if m is not nomask:
 
1945
        m = numeric.compress(c, m, dimension)
 
1946
    d = numeric.compress(c, filled(x), dimension)
 
1947
    return masked_array(d, m)
 
1948
 
 
1949
class _minimum_operation:
 
1950
    "Object to calculate minima"
 
1951
    def __init__ (self):
 
1952
        """minimum(a, b) or minimum(a)
 
1953
           In one argument case returns the scalar minimum.
 
1954
        """
 
1955
        pass
 
1956
 
 
1957
    def __call__ (self, a, b=None):
 
1958
        "Execute the call behavior."
 
1959
        if b is None:
 
1960
            m = getmask(a)
 
1961
            if m is nomask:
 
1962
                d = amin(filled(a).ravel())
 
1963
                return d
 
1964
            ac = a.compressed()
 
1965
            if len(ac) == 0:
 
1966
                return masked
 
1967
            else:
 
1968
                return amin(ac.raw_data())
 
1969
        else:
 
1970
            return where(less(a, b), a, b)
 
1971
 
 
1972
    def reduce (self, target, axis=0):
 
1973
        """Reduce target along the given axis."""
 
1974
        m = getmask(target)
 
1975
        if m is nomask:
 
1976
            t = filled(target)
 
1977
            return masked_array (umath.minimum.reduce (t, axis))
 
1978
        else:
 
1979
            t = umath.minimum.reduce(filled(target, minimum_fill_value(target)), axis)
 
1980
            m = umath.logical_and.reduce(m, axis)
 
1981
            return masked_array(t, m, get_fill_value(target))
 
1982
 
 
1983
    def outer (self, a, b):
 
1984
        "Return the function applied to the outer product of a and b."
 
1985
        ma = getmask(a)
 
1986
        mb = getmask(b)
 
1987
        if ma is nomask and mb is nomask:
 
1988
            m = nomask
 
1989
        else:
 
1990
            ma = getmaskarray(a)
 
1991
            mb = getmaskarray(b)
 
1992
            m = logical_or.outer(ma, mb)
 
1993
        d = umath.minimum.outer(filled(a), filled(b))
 
1994
        return masked_array(d, m)
 
1995
 
 
1996
minimum = _minimum_operation ()
 
1997
 
 
1998
class _maximum_operation:
 
1999
    "Object to calculate maxima"
 
2000
    def __init__ (self):
 
2001
        """maximum(a, b) or maximum(a)
 
2002
           In one argument case returns the scalar maximum.
 
2003
        """
 
2004
        pass
 
2005
 
 
2006
    def __call__ (self, a, b=None):
 
2007
        "Execute the call behavior."
 
2008
        if b is None:
 
2009
            m = getmask(a)
 
2010
            if m is nomask:
 
2011
                d = amax(filled(a).ravel())
 
2012
                return d
 
2013
            ac = a.compressed()
 
2014
            if len(ac) == 0:
 
2015
                return masked
 
2016
            else:
 
2017
                return amax(ac.raw_data())
 
2018
        else:
 
2019
            return where(greater(a, b), a, b)
 
2020
 
 
2021
    def reduce (self, target, axis=0):
 
2022
        """Reduce target along the given axis."""
 
2023
        m = getmask(target)
 
2024
        if m is nomask:
 
2025
            t = filled(target)
 
2026
            return masked_array (umath.maximum.reduce (t, axis))
 
2027
        else:
 
2028
            t = umath.maximum.reduce(filled(target, maximum_fill_value(target)), axis)
 
2029
            m = umath.logical_and.reduce(m, axis)
 
2030
            return masked_array(t, m, get_fill_value(target))
 
2031
 
 
2032
    def outer (self, a, b):
 
2033
        "Return the function applied to the outer product of a and b."
 
2034
        ma = getmask(a)
 
2035
        mb = getmask(b)
 
2036
        if ma is nomask and mb is nomask:
 
2037
            m = nomask
 
2038
        else:
 
2039
            ma = getmaskarray(a)
 
2040
            mb = getmaskarray(b)
 
2041
            m = logical_or.outer(ma, mb)
 
2042
        d = umath.maximum.outer(filled(a), filled(b))
 
2043
        return masked_array(d, m)
 
2044
 
 
2045
maximum = _maximum_operation ()
 
2046
 
 
2047
def sort (x, axis = -1, fill_value=None):
 
2048
    """If x does not have a mask, return a masked array formed from the
 
2049
       result of numeric.sort(x, axis).
 
2050
       Otherwise, fill x with fill_value. Sort it.
 
2051
       Set a mask where the result is equal to fill_value.
 
2052
       Note that this may have unintended consequences if the data contains the
 
2053
       fill value at a non-masked site.
 
2054
 
 
2055
       If fill_value is not given the default fill value for x's type will be
 
2056
       used.
 
2057
    """
 
2058
    if fill_value is None:
 
2059
        fill_value = default_fill_value (x)
 
2060
    d = filled(x, fill_value)
 
2061
    s = fromnumeric.sort(d, axis)
 
2062
    if getmask(x) is nomask:
 
2063
        return masked_array(s)
 
2064
    return masked_values(s, fill_value, copy=0)
 
2065
 
 
2066
def diagonal(a, k = 0, axis1=0, axis2=1):
 
2067
    """diagonal(a,k=0,axis1=0, axis2=1) = the k'th diagonal of a"""
 
2068
    d = fromnumeric.diagonal(filled(a), k, axis1, axis2)
 
2069
    m = getmask(a)
 
2070
    if m is nomask:
 
2071
        return masked_array(d, m)
 
2072
    else:
 
2073
        return masked_array(d, fromnumeric.diagonal(m, k, axis1, axis2))
 
2074
 
 
2075
def trace (a, offset=0, axis1=0, axis2=1, dtype=None, out=None):
 
2076
    """trace(a,offset=0, axis1=0, axis2=1) returns the sum along diagonals
 
2077
    (defined by the last two dimenions) of the array.
 
2078
    """
 
2079
    return diagonal(a, offset, axis1, axis2).sum(dtype=dtype)
 
2080
 
 
2081
def argsort (x, axis = -1, out=None, fill_value=None):
 
2082
    """Treating masked values as if they have the value fill_value,
 
2083
       return sort indices for sorting along given axis.
 
2084
       if fill_value is None, use get_fill_value(x)
 
2085
       Returns a numpy array.
 
2086
    """
 
2087
    d = filled(x, fill_value)
 
2088
    return fromnumeric.argsort(d, axis)
 
2089
 
 
2090
def argmin (x, axis = -1, out=None, fill_value=None):
 
2091
    """Treating masked values as if they have the value fill_value,
 
2092
       return indices for minimum values along given axis.
 
2093
       if fill_value is None, use get_fill_value(x).
 
2094
       Returns a numpy array if x has more than one dimension.
 
2095
       Otherwise, returns a scalar index.
 
2096
    """
 
2097
    d = filled(x, fill_value)
 
2098
    return fromnumeric.argmin(d, axis)
 
2099
 
 
2100
def argmax (x, axis = -1, out=None, fill_value=None):
 
2101
    """Treating masked values as if they have the value fill_value,
 
2102
       return sort indices for maximum along given axis.
 
2103
       if fill_value is None, use -get_fill_value(x) if it exists.
 
2104
       Returns a numpy array if x has more than one dimension.
 
2105
       Otherwise, returns a scalar index.
 
2106
    """
 
2107
    if fill_value is None:
 
2108
        fill_value = default_fill_value (x)
 
2109
        try:
 
2110
            fill_value = - fill_value
 
2111
        except:
 
2112
            pass
 
2113
    d = filled(x, fill_value)
 
2114
    return fromnumeric.argmax(d, axis)
 
2115
 
 
2116
def fromfunction (f, s):
 
2117
    """apply f to s to create array as in umath."""
 
2118
    return masked_array(numeric.fromfunction(f, s))
 
2119
 
 
2120
def asarray(data, dtype=None):
 
2121
    """asarray(data, dtype) = array(data, dtype, copy=0)
 
2122
    """
 
2123
    if isinstance(data, MaskedArray) and \
 
2124
        (dtype is None or dtype == data.dtype):
 
2125
        return data
 
2126
    return array(data, dtype=dtype, copy=0)
 
2127
 
 
2128
# Add methods to support ndarray interface
 
2129
# XXX: I is better to to change the masked_*_operation adaptors
 
2130
# XXX: to wrap ndarray methods directly to create ma.array methods.
 
2131
from types import MethodType
 
2132
def _m(f):
 
2133
    return MethodType(f, None, array)
 
2134
def not_implemented(*args, **kwds):
 
2135
    raise NotImplementedError, "not yet implemented for numpy.ma arrays"
 
2136
array.all = _m(alltrue)
 
2137
array.any = _m(sometrue)
 
2138
array.argmax = _m(argmax)
 
2139
array.argmin = _m(argmin)
 
2140
array.argsort = _m(argsort)
 
2141
array.base = property(_m(not_implemented))
 
2142
array.byteswap = _m(not_implemented)
 
2143
 
 
2144
def _choose(self, *args, **kwds):
 
2145
    return choose(self, args)
 
2146
array.choose = _m(_choose)
 
2147
del _choose
 
2148
 
 
2149
def _clip(self,a_min,a_max,out=None):
 
2150
    return MaskedArray(data = self.data.clip(asarray(a_min).data,
 
2151
                                             asarray(a_max).data),
 
2152
                       mask = mask_or(self.mask,
 
2153
                                      mask_or(getmask(a_min),getmask(a_max))))
 
2154
array.clip = _m(_clip)
 
2155
 
 
2156
def _compress(self, cond, axis=None, out=None):
 
2157
    return compress(cond, self, axis)
 
2158
array.compress = _m(_compress)
 
2159
del _compress
 
2160
 
 
2161
array.conj = array.conjugate = _m(conjugate)
 
2162
array.copy = _m(not_implemented)
 
2163
 
 
2164
def _cumprod(self, axis=None, dtype=None, out=None):
 
2165
    m = self.mask
 
2166
    if m is not nomask:
 
2167
        m = umath.logical_or.accumulate(self.mask, axis)
 
2168
    return MaskedArray(data = self.filled(1).cumprod(axis, dtype), mask=m)
 
2169
array.cumprod = _m(_cumprod)
 
2170
 
 
2171
def _cumsum(self, axis=None, dtype=None, out=None):
 
2172
    m = self.mask
 
2173
    if m is not nomask:
 
2174
        m = umath.logical_or.accumulate(self.mask, axis)
 
2175
    return MaskedArray(data=self.filled(0).cumsum(axis, dtype), mask=m)
 
2176
array.cumsum = _m(_cumsum)
 
2177
 
 
2178
array.diagonal = _m(diagonal)
 
2179
array.dump = _m(not_implemented)
 
2180
array.dumps = _m(not_implemented)
 
2181
array.fill = _m(not_implemented)
 
2182
array.flags = property(_m(not_implemented))
 
2183
array.flatten = _m(ravel)
 
2184
array.getfield = _m(not_implemented)
 
2185
 
 
2186
def _max(a, axis=None, out=None):
 
2187
    if out is not None:
 
2188
        raise TypeError("Output arrays Unsupported for masked arrays")
 
2189
    if axis is None:
 
2190
        return maximum(a)
 
2191
    else:
 
2192
        return maximum.reduce(a, axis)
 
2193
array.max = _m(_max)
 
2194
del _max
 
2195
def _min(a, axis=None, out=None):
 
2196
    if out is not None:
 
2197
        raise TypeError("Output arrays Unsupported for masked arrays")
 
2198
    if axis is None:
 
2199
        return minimum(a)
 
2200
    else:
 
2201
        return minimum.reduce(a, axis)
 
2202
array.min = _m(_min)
 
2203
del _min
 
2204
array.mean = _m(new_average)
 
2205
array.nbytes = property(_m(not_implemented))
 
2206
array.newbyteorder = _m(not_implemented)
 
2207
array.nonzero = _m(nonzero)
 
2208
array.prod = _m(product)
 
2209
 
 
2210
def _ptp(a,axis=None,out=None):
 
2211
    return a.max(axis,out)-a.min(axis)
 
2212
array.ptp = _m(_ptp)
 
2213
array.repeat = _m(new_repeat)
 
2214
array.resize = _m(resize)
 
2215
array.searchsorted = _m(not_implemented)
 
2216
array.setfield = _m(not_implemented)
 
2217
array.setflags = _m(not_implemented)
 
2218
array.sort = _m(not_implemented)  # NB: ndarray.sort is inplace
 
2219
 
 
2220
def _squeeze(self):
 
2221
    try:
 
2222
        result = MaskedArray(data = self.data.squeeze(),
 
2223
                             mask = self.mask.squeeze())
 
2224
    except AttributeError:
 
2225
        result = _wrapit(self, 'squeeze')
 
2226
    return result
 
2227
array.squeeze = _m(_squeeze)
 
2228
 
 
2229
array.strides = property(_m(not_implemented))
 
2230
array.sum = _m(sum)
 
2231
def _swapaxes(self,axis1,axis2):
 
2232
    return MaskedArray(data = self.data.swapaxes(axis1, axis2),
 
2233
                       mask = self.mask.swapaxes(axis1, axis2))
 
2234
array.swapaxes = _m(_swapaxes)
 
2235
array.take = _m(new_take)
 
2236
array.tofile = _m(not_implemented)
 
2237
array.trace = _m(trace)
 
2238
array.transpose = _m(transpose)
 
2239
 
 
2240
def _var(self,axis=None,dtype=None, out=None):
 
2241
    if axis is None:
 
2242
        return numeric.asarray(self.compressed()).var()
 
2243
    a = self.swapaxes(axis,0)
 
2244
    a = a - a.mean(axis=0)
 
2245
    a *= a
 
2246
    a /= a.count(axis=0)
 
2247
    return a.swapaxes(0,axis).sum(axis)
 
2248
def _std(self,axis=None, dtype=None, out=None):
 
2249
    return (self.var(axis,dtype))**0.5
 
2250
array.var = _m(_var)
 
2251
array.std = _m(_std)
 
2252
 
 
2253
array.view =  _m(not_implemented)
 
2254
array.round = _m(around)
 
2255
del _m, MethodType, not_implemented
 
2256
 
 
2257
 
 
2258
masked = MaskedArray(0, int, mask=1)
6
2259
 
7
2260
def repeat(a, repeats, axis=0):
8
 
    return nca.repeat(a, repeats, axis)
 
2261
    return new_repeat(a, repeats, axis)
9
2262
 
10
2263
def average(a, axis=0, weights=None, returned=0):
11
 
    return nca.average(a, axis, weights, returned)
 
2264
    return new_average(a, axis, weights, returned)
12
2265
 
13
2266
def take(a, indices, axis=0):
14
 
    return nca.average(a, indices, axis=0)
15
 
 
 
2267
    return new_take(a, indices, axis)