source: MGET/Branches/Jason/PythonPackage/src/GeoEco/AggregatedModules/numpy-1.6.2/core/numeric.py @ 1025

Revision 1025, 68.1 KB checked in by jjr8, 8 years ago (diff)
  • Incremented build number
  • Added support for Python 2.7 (this has not been tested yet!)
  • Upgraded MGET's internal copy of GDAL to 1.9.1
  • Added support for numpy 1.6.2
  • Upgraded MGET's internal copy of numpy to 1.6.2 on Python 2.5, 2.6, 2.7
  • Upgraded MGET's internal copy of lxml to 2.2.8
Line 
1__all__ = ['newaxis', 'ndarray', 'flatiter', 'nditer', 'nested_iters', 'ufunc',
2           'arange', 'array', 'zeros', 'count_nonzero', 'empty', 'broadcast',
3           'dtype', 'fromstring', 'fromfile', 'frombuffer',
4           'int_asbuffer', 'where', 'argwhere',
5           'concatenate', 'fastCopyAndTranspose', 'lexsort', 'set_numeric_ops',
6           'can_cast', 'promote_types', 'min_scalar_type', 'result_type',
7           'asarray', 'asanyarray', 'ascontiguousarray', 'asfortranarray',
8           'isfortran', 'empty_like', 'zeros_like',
9           'correlate', 'convolve', 'inner', 'dot', 'einsum', 'outer', 'vdot',
10           'alterdot', 'restoredot', 'roll', 'rollaxis', 'cross', 'tensordot',
11           'array2string', 'get_printoptions', 'set_printoptions',
12           'array_repr', 'array_str', 'set_string_function',
13           'little_endian', 'require',
14           'fromiter', 'array_equal', 'array_equiv',
15           'indices', 'fromfunction',
16           'load', 'loads', 'isscalar', 'binary_repr', 'base_repr',
17           'ones', 'identity', 'allclose', 'compare_chararrays', 'putmask',
18           'seterr', 'geterr', 'setbufsize', 'getbufsize',
19           'seterrcall', 'geterrcall', 'errstate', 'flatnonzero',
20           'Inf', 'inf', 'infty', 'Infinity',
21           'nan', 'NaN', 'False_', 'True_', 'bitwise_not',
22           'CLIP', 'RAISE', 'WRAP', 'MAXDIMS', 'BUFSIZE', 'ALLOW_THREADS',
23           'ComplexWarning']
24
25import sys
26import warnings
27import multiarray
28import umath
29from umath import *
30import numerictypes
31from numerictypes import *
32
33if sys.version_info[0] < 3:
34    __all__.extend(['getbuffer', 'newbuffer'])
35
36class ComplexWarning(RuntimeWarning):
37    """
38    The warning raised when casting a complex dtype to a real dtype.
39
40    As implemented, casting a complex number to a real discards its imaginary
41    part, but this behavior may not be what the user actually wants.
42
43    """
44    pass
45
46bitwise_not = invert
47
48CLIP = multiarray.CLIP
49WRAP = multiarray.WRAP
50RAISE = multiarray.RAISE
51MAXDIMS = multiarray.MAXDIMS
52ALLOW_THREADS = multiarray.ALLOW_THREADS
53BUFSIZE = multiarray.BUFSIZE
54
55ndarray = multiarray.ndarray
56flatiter = multiarray.flatiter
57nditer = multiarray.nditer
58nested_iters = multiarray.nested_iters
59broadcast = multiarray.broadcast
60dtype = multiarray.dtype
61ufunc = type(sin)
62
63
64# originally from Fernando Perez's IPython
65def zeros_like(a, dtype=None, order='K', subok=True):
66    """
67    Return an array of zeros with the same shape and type as a given array.
68
69    With default parameters, is equivalent to ``a.copy().fill(0)``.
70
71    Parameters
72    ----------
73    a : array_like
74        The shape and data-type of `a` define these same attributes of
75        the returned array.
76    dtype : data-type, optional
77        Overrides the data type of the result.
78    order : {'C', 'F', 'A', or 'K'}, optional
79        Overrides the memory layout of the result. 'C' means C-order,
80        'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous,
81        'C' otherwise. 'K' means match the layout of `a` as closely
82        as possible.
83
84    Returns
85    -------
86    out : ndarray
87        Array of zeros with the same shape and type as `a`.
88
89    See Also
90    --------
91    ones_like : Return an array of ones with shape and type of input.
92    empty_like : Return an empty array with shape and type of input.
93    zeros : Return a new array setting values to zero.
94    ones : Return a new array setting values to one.
95    empty : Return a new uninitialized array.
96
97    Examples
98    --------
99    >>> x = np.arange(6)
100    >>> x = x.reshape((2, 3))
101    >>> x
102    array([[0, 1, 2],
103           [3, 4, 5]])
104    >>> np.zeros_like(x)
105    array([[0, 0, 0],
106           [0, 0, 0]])
107
108    >>> y = np.arange(3, dtype=np.float)
109    >>> y
110    array([ 0.,  1.,  2.])
111    >>> np.zeros_like(y)
112    array([ 0.,  0.,  0.])
113
114    """
115    res = empty_like(a, dtype=dtype, order=order, subok=subok)
116    res.fill(0)
117    return res
118
119# end Fernando's utilities
120
121
122def extend_all(module):
123    adict = {}
124    for a in __all__:
125        adict[a] = 1
126    try:
127        mall = getattr(module, '__all__')
128    except AttributeError:
129        mall = [k for k in module.__dict__.keys() if not k.startswith('_')]
130    for a in mall:
131        if a not in adict:
132            __all__.append(a)
133
134extend_all(umath)
135extend_all(numerictypes)
136
137newaxis = None
138
139
140arange = multiarray.arange
141array = multiarray.array
142zeros = multiarray.zeros
143count_nonzero = multiarray.count_nonzero
144empty = multiarray.empty
145empty_like = multiarray.empty_like
146fromstring = multiarray.fromstring
147fromiter = multiarray.fromiter
148fromfile = multiarray.fromfile
149frombuffer = multiarray.frombuffer
150if sys.version_info[0] < 3:
151    newbuffer = multiarray.newbuffer
152    getbuffer = multiarray.getbuffer
153int_asbuffer = multiarray.int_asbuffer
154where = multiarray.where
155concatenate = multiarray.concatenate
156fastCopyAndTranspose = multiarray._fastCopyAndTranspose
157set_numeric_ops = multiarray.set_numeric_ops
158can_cast = multiarray.can_cast
159promote_types = multiarray.promote_types
160min_scalar_type = multiarray.min_scalar_type
161result_type = multiarray.result_type
162lexsort = multiarray.lexsort
163compare_chararrays = multiarray.compare_chararrays
164putmask = multiarray.putmask
165einsum = multiarray.einsum
166
167def asarray(a, dtype=None, order=None):
168    """
169    Convert the input to an array.
170
171    Parameters
172    ----------
173    a : array_like
174        Input data, in any form that can be converted to an array.  This
175        includes lists, lists of tuples, tuples, tuples of tuples, tuples
176        of lists and ndarrays.
177    dtype : data-type, optional
178        By default, the data-type is inferred from the input data.
179    order : {'C', 'F'}, optional
180        Whether to use row-major ('C') or column-major ('F' for FORTRAN)
181        memory representation.  Defaults to 'C'.
182
183    Returns
184    -------
185    out : ndarray
186        Array interpretation of `a`.  No copy is performed if the input
187        is already an ndarray.  If `a` is a subclass of ndarray, a base
188        class ndarray is returned.
189
190    See Also
191    --------
192    asanyarray : Similar function which passes through subclasses.
193    ascontiguousarray : Convert input to a contiguous array.
194    asfarray : Convert input to a floating point ndarray.
195    asfortranarray : Convert input to an ndarray with column-major
196                     memory order.
197    asarray_chkfinite : Similar function which checks input for NaNs and Infs.
198    fromiter : Create an array from an iterator.
199    fromfunction : Construct an array by executing a function on grid
200                   positions.
201
202    Examples
203    --------
204    Convert a list into an array:
205
206    >>> a = [1, 2]
207    >>> np.asarray(a)
208    array([1, 2])
209
210    Existing arrays are not copied:
211
212    >>> a = np.array([1, 2])
213    >>> np.asarray(a) is a
214    True
215
216    If `dtype` is set, array is copied only if dtype does not match:
217
218    >>> a = np.array([1, 2], dtype=np.float32)
219    >>> np.asarray(a, dtype=np.float32) is a
220    True
221    >>> np.asarray(a, dtype=np.float64) is a
222    False
223
224    Contrary to `asanyarray`, ndarray subclasses are not passed through:
225
226    >>> issubclass(np.matrix, np.ndarray)
227    True
228    >>> a = np.matrix([[1, 2]])
229    >>> np.asarray(a) is a
230    False
231    >>> np.asanyarray(a) is a
232    True
233
234    """
235    return array(a, dtype, copy=False, order=order)
236
237def asanyarray(a, dtype=None, order=None):
238    """
239    Convert the input to an ndarray, but pass ndarray subclasses through.
240
241    Parameters
242    ----------
243    a : array_like
244        Input data, in any form that can be converted to an array.  This
245        includes scalars, lists, lists of tuples, tuples, tuples of tuples,
246        tuples of lists, and ndarrays.
247    dtype : data-type, optional
248        By default, the data-type is inferred from the input data.
249    order : {'C', 'F'}, optional
250        Whether to use row-major ('C') or column-major ('F') memory
251        representation.  Defaults to 'C'.
252
253    Returns
254    -------
255    out : ndarray or an ndarray subclass
256        Array interpretation of `a`.  If `a` is an ndarray or a subclass
257        of ndarray, it is returned as-is and no copy is performed.
258
259    See Also
260    --------
261    asarray : Similar function which always returns ndarrays.
262    ascontiguousarray : Convert input to a contiguous array.
263    asfarray : Convert input to a floating point ndarray.
264    asfortranarray : Convert input to an ndarray with column-major
265                     memory order.
266    asarray_chkfinite : Similar function which checks input for NaNs and
267                        Infs.
268    fromiter : Create an array from an iterator.
269    fromfunction : Construct an array by executing a function on grid
270                   positions.
271
272    Examples
273    --------
274    Convert a list into an array:
275
276    >>> a = [1, 2]
277    >>> np.asanyarray(a)
278    array([1, 2])
279
280    Instances of `ndarray` subclasses are passed through as-is:
281
282    >>> a = np.matrix([1, 2])
283    >>> np.asanyarray(a) is a
284    True
285
286    """
287    return array(a, dtype, copy=False, order=order, subok=True)
288
289def ascontiguousarray(a, dtype=None):
290    """
291    Return a contiguous array in memory (C order).
292
293    Parameters
294    ----------
295    a : array_like
296        Input array.
297    dtype : str or dtype object, optional
298        Data-type of returned array.
299
300    Returns
301    -------
302    out : ndarray
303        Contiguous array of same shape and content as `a`, with type `dtype`
304        if specified.
305
306    See Also
307    --------
308    asfortranarray : Convert input to an ndarray with column-major
309                     memory order.
310    require : Return an ndarray that satisfies requirements.
311    ndarray.flags : Information about the memory layout of the array.
312
313    Examples
314    --------
315    >>> x = np.arange(6).reshape(2,3)
316    >>> np.ascontiguousarray(x, dtype=np.float32)
317    array([[ 0.,  1.,  2.],
318           [ 3.,  4.,  5.]], dtype=float32)
319    >>> x.flags['C_CONTIGUOUS']
320    True
321
322    """
323    return array(a, dtype, copy=False, order='C', ndmin=1)
324
325def asfortranarray(a, dtype=None):
326    """
327    Return an array laid out in Fortran order in memory.
328
329    Parameters
330    ----------
331    a : array_like
332        Input array.
333    dtype : str or dtype object, optional
334        By default, the data-type is inferred from the input data.
335
336    Returns
337    -------
338    out : ndarray
339        The input `a` in Fortran, or column-major, order.
340
341    See Also
342    --------
343    ascontiguousarray : Convert input to a contiguous (C order) array.
344    asanyarray : Convert input to an ndarray with either row or
345        column-major memory order.
346    require : Return an ndarray that satisfies requirements.
347    ndarray.flags : Information about the memory layout of the array.
348
349    Examples
350    --------
351    >>> x = np.arange(6).reshape(2,3)
352    >>> y = np.asfortranarray(x)
353    >>> x.flags['F_CONTIGUOUS']
354    False
355    >>> y.flags['F_CONTIGUOUS']
356    True
357
358    """
359    return array(a, dtype, copy=False, order='F', ndmin=1)
360
361def require(a, dtype=None, requirements=None):
362    """
363    Return an ndarray of the provided type that satisfies requirements.
364
365    This function is useful to be sure that an array with the correct flags
366    is returned for passing to compiled code (perhaps through ctypes).
367
368    Parameters
369    ----------
370    a : array_like
371       The object to be converted to a type-and-requirement-satisfying array.
372    dtype : data-type
373       The required data-type, the default data-type is float64).
374    requirements : str or list of str
375       The requirements list can be any of the following
376
377       * 'F_CONTIGUOUS' ('F') - ensure a Fortran-contiguous array
378       * 'C_CONTIGUOUS' ('C') - ensure a C-contiguous array
379       * 'ALIGNED' ('A')      - ensure a data-type aligned array
380       * 'WRITEABLE' ('W')    - ensure a writable array
381       * 'OWNDATA' ('O')      - ensure an array that owns its own data
382
383    See Also
384    --------
385    asarray : Convert input to an ndarray.
386    asanyarray : Convert to an ndarray, but pass through ndarray subclasses.
387    ascontiguousarray : Convert input to a contiguous array.
388    asfortranarray : Convert input to an ndarray with column-major
389                     memory order.
390    ndarray.flags : Information about the memory layout of the array.
391
392    Notes
393    -----
394    The returned array will be guaranteed to have the listed requirements
395    by making a copy if needed.
396
397    Examples
398    --------
399    >>> x = np.arange(6).reshape(2,3)
400    >>> x.flags
401      C_CONTIGUOUS : True
402      F_CONTIGUOUS : False
403      OWNDATA : False
404      WRITEABLE : True
405      ALIGNED : True
406      UPDATEIFCOPY : False
407
408    >>> y = np.require(x, dtype=np.float32, requirements=['A', 'O', 'W', 'F'])
409    >>> y.flags
410      C_CONTIGUOUS : False
411      F_CONTIGUOUS : True
412      OWNDATA : True
413      WRITEABLE : True
414      ALIGNED : True
415      UPDATEIFCOPY : False
416
417    """
418    if requirements is None:
419        requirements = []
420    else:
421        requirements = [x.upper() for x in requirements]
422
423    if not requirements:
424        return asanyarray(a, dtype=dtype)
425
426    if 'ENSUREARRAY' in requirements or 'E' in requirements:
427        subok = False
428    else:
429        subok = True
430
431    arr = array(a, dtype=dtype, copy=False, subok=subok)
432
433    copychar = 'A'
434    if 'FORTRAN' in requirements or \
435       'F_CONTIGUOUS' in requirements or \
436       'F' in requirements:
437        copychar = 'F'
438    elif 'CONTIGUOUS' in requirements or \
439         'C_CONTIGUOUS' in requirements or \
440         'C' in requirements:
441        copychar = 'C'
442
443    for prop in requirements:
444        if not arr.flags[prop]:
445            arr = arr.copy(copychar)
446            break
447    return arr
448
449def isfortran(a):
450    """
451    Returns True if array is arranged in Fortran-order in memory
452    and dimension > 1.
453
454    Parameters
455    ----------
456    a : ndarray
457        Input array.
458
459
460    Examples
461    --------
462
463    np.array allows to specify whether the array is written in C-contiguous
464    order (last index varies the fastest), or FORTRAN-contiguous order in
465    memory (first index varies the fastest).
466
467    >>> a = np.array([[1, 2, 3], [4, 5, 6]], order='C')
468    >>> a
469    array([[1, 2, 3],
470           [4, 5, 6]])
471    >>> np.isfortran(a)
472    False
473
474    >>> b = np.array([[1, 2, 3], [4, 5, 6]], order='FORTRAN')
475    >>> b
476    array([[1, 2, 3],
477           [4, 5, 6]])
478    >>> np.isfortran(b)
479    True
480
481
482    The transpose of a C-ordered array is a FORTRAN-ordered array.
483
484    >>> a = np.array([[1, 2, 3], [4, 5, 6]], order='C')
485    >>> a
486    array([[1, 2, 3],
487           [4, 5, 6]])
488    >>> np.isfortran(a)
489    False
490    >>> b = a.T
491    >>> b
492    array([[1, 4],
493           [2, 5],
494           [3, 6]])
495    >>> np.isfortran(b)
496    True
497
498    1-D arrays always evaluate as False.
499
500    >>> np.isfortran(np.array([1, 2], order='FORTRAN'))
501    False
502
503    """
504    return a.flags.fnc
505
506def argwhere(a):
507    """
508    Find the indices of array elements that are non-zero, grouped by element.
509
510    Parameters
511    ----------
512    a : array_like
513        Input data.
514
515    Returns
516    -------
517    index_array : ndarray
518        Indices of elements that are non-zero. Indices are grouped by element.
519
520    See Also
521    --------
522    where, nonzero
523
524    Notes
525    -----
526    ``np.argwhere(a)`` is the same as ``np.transpose(np.nonzero(a))``.
527
528    The output of ``argwhere`` is not suitable for indexing arrays.
529    For this purpose use ``where(a)`` instead.
530
531    Examples
532    --------
533    >>> x = np.arange(6).reshape(2,3)
534    >>> x
535    array([[0, 1, 2],
536           [3, 4, 5]])
537    >>> np.argwhere(x>1)
538    array([[0, 2],
539           [1, 0],
540           [1, 1],
541           [1, 2]])
542
543    """
544    return transpose(asanyarray(a).nonzero())
545
546def flatnonzero(a):
547    """
548    Return indices that are non-zero in the flattened version of a.
549
550    This is equivalent to a.ravel().nonzero()[0].
551
552    Parameters
553    ----------
554    a : ndarray
555        Input array.
556
557    Returns
558    -------
559    res : ndarray
560        Output array, containing the indices of the elements of `a.ravel()`
561        that are non-zero.
562
563    See Also
564    --------
565    nonzero : Return the indices of the non-zero elements of the input array.
566    ravel : Return a 1-D array containing the elements of the input array.
567
568    Examples
569    --------
570    >>> x = np.arange(-2, 3)
571    >>> x
572    array([-2, -1,  0,  1,  2])
573    >>> np.flatnonzero(x)
574    array([0, 1, 3, 4])
575
576    Use the indices of the non-zero elements as an index array to extract
577    these elements:
578
579    >>> x.ravel()[np.flatnonzero(x)]
580    array([-2, -1,  1,  2])
581
582    """
583    return a.ravel().nonzero()[0]
584
585_mode_from_name_dict = {'v': 0,
586                        's' : 1,
587                        'f' : 2}
588
589def _mode_from_name(mode):
590    if isinstance(mode, type("")):
591        return _mode_from_name_dict[mode.lower()[0]]
592    return mode
593
594def correlate(a, v, mode='valid', old_behavior=False):
595    """
596    Cross-correlation of two 1-dimensional sequences.
597
598    This function computes the correlation as generally defined in signal
599    processing texts::
600
601        z[k] = sum_n a[n] * conj(v[n+k])
602
603    with a and v sequences being zero-padded where necessary and conj being
604    the conjugate.
605
606    Parameters
607    ----------
608    a, v : array_like
609        Input sequences.
610    mode : {'valid', 'same', 'full'}, optional
611        Refer to the `convolve` docstring.  Note that the default
612        is `valid`, unlike `convolve`, which uses `full`.
613    old_behavior : bool
614        If True, uses the old behavior from Numeric, (correlate(a,v) == correlate(v,
615        a), and the conjugate is not taken for complex arrays). If False, uses
616        the conventional signal processing definition (see note).
617
618    See Also
619    --------
620    convolve : Discrete, linear convolution of two one-dimensional sequences.
621
622    Examples
623    --------
624    >>> np.correlate([1, 2, 3], [0, 1, 0.5])
625    array([ 3.5])
626    >>> np.correlate([1, 2, 3], [0, 1, 0.5], "same")
627    array([ 2. ,  3.5,  3. ])
628    >>> np.correlate([1, 2, 3], [0, 1, 0.5], "full")
629    array([ 0.5,  2. ,  3.5,  3. ,  0. ])
630
631    """
632    mode = _mode_from_name(mode)
633# the old behavior should be made available under a different name, see thread
634# http://thread.gmane.org/gmane.comp.python.numeric.general/12609/focus=12630
635    if old_behavior:
636        warnings.warn("""
637The old behavior of correlate was deprecated for 1.4.0, and will be completely removed
638for NumPy 2.0.
639
640The new behavior fits the conventional definition of correlation: inputs are
641never swapped, and the second argument is conjugated for complex arrays.""",
642            DeprecationWarning)
643        return multiarray.correlate(a,v,mode)
644    else:
645        return multiarray.correlate2(a,v,mode)
646
647def convolve(a,v,mode='full'):
648    """
649    Returns the discrete, linear convolution of two one-dimensional sequences.
650
651    The convolution operator is often seen in signal processing, where it
652    models the effect of a linear time-invariant system on a signal [1]_.  In
653    probability theory, the sum of two independent random variables is
654    distributed according to the convolution of their individual
655    distributions.
656
657    Parameters
658    ----------
659    a : (N,) array_like
660        First one-dimensional input array.
661    v : (M,) array_like
662        Second one-dimensional input array.
663    mode : {'full', 'valid', 'same'}, optional
664        'full':
665          By default, mode is 'full'.  This returns the convolution
666          at each point of overlap, with an output shape of (N+M-1,). At
667          the end-points of the convolution, the signals do not overlap
668          completely, and boundary effects may be seen.
669
670        'same':
671          Mode `same` returns output of length ``max(M, N)``.  Boundary
672          effects are still visible.
673
674        'valid':
675          Mode `valid` returns output of length
676          ``max(M, N) - min(M, N) + 1``.  The convolution product is only given
677          for points where the signals overlap completely.  Values outside
678          the signal boundary have no effect.
679
680    Returns
681    -------
682    out : ndarray
683        Discrete, linear convolution of `a` and `v`.
684
685    See Also
686    --------
687    scipy.signal.fftconvolve : Convolve two arrays using the Fast Fourier
688                               Transform.
689    scipy.linalg.toeplitz : Used to construct the convolution operator.
690
691    Notes
692    -----
693    The discrete convolution operation is defined as
694
695    .. math:: (f * g)[n] = \\sum_{m = -\\infty}^{\\infty} f[m] g[n - m]
696
697    It can be shown that a convolution :math:`x(t) * y(t)` in time/space
698    is equivalent to the multiplication :math:`X(f) Y(f)` in the Fourier
699    domain, after appropriate padding (padding is necessary to prevent
700    circular convolution).  Since multiplication is more efficient (faster)
701    than convolution, the function `scipy.signal.fftconvolve` exploits the
702    FFT to calculate the convolution of large data-sets.
703
704    References
705    ----------
706    .. [1] Wikipedia, "Convolution", http://en.wikipedia.org/wiki/Convolution.
707
708    Examples
709    --------
710    Note how the convolution operator flips the second array
711    before "sliding" the two across one another:
712
713    >>> np.convolve([1, 2, 3], [0, 1, 0.5])
714    array([ 0. ,  1. ,  2.5,  4. ,  1.5])
715
716    Only return the middle values of the convolution.
717    Contains boundary effects, where zeros are taken
718    into account:
719
720    >>> np.convolve([1,2,3],[0,1,0.5], 'same')
721    array([ 1. ,  2.5,  4. ])
722
723    The two arrays are of the same length, so there
724    is only one position where they completely overlap:
725
726    >>> np.convolve([1,2,3],[0,1,0.5], 'valid')
727    array([ 2.5])
728
729    """
730    a,v = array(a, ndmin=1),array(v, ndmin=1)
731    if (len(v) > len(a)):
732        a, v = v, a
733    if len(a) == 0 :
734        raise ValueError('a cannot be empty')
735    if len(v) == 0 :
736        raise ValueError('v cannot be empty')
737    mode = _mode_from_name(mode)
738    return multiarray.correlate(a, v[::-1], mode)
739
740def outer(a,b):
741    """
742    Compute the outer product of two vectors.
743
744    Given two vectors, ``a = [a0, a1, ..., aM]`` and
745    ``b = [b0, b1, ..., bN]``,
746    the outer product [1]_ is::
747
748      [[a0*b0  a0*b1 ... a0*bN ]
749       [a1*b0    .
750       [ ...          .
751       [aM*b0            aM*bN ]]
752
753    Parameters
754    ----------
755    a, b : array_like, shape (M,), (N,)
756        First and second input vectors.  Inputs are flattened if they
757        are not already 1-dimensional.
758
759    Returns
760    -------
761    out : ndarray, shape (M, N)
762        ``out[i, j] = a[i] * b[j]``
763
764    See also
765    --------
766    inner, einsum
767
768    References
769    ----------
770    .. [1] : G. H. Golub and C. F. van Loan, *Matrix Computations*, 3rd
771             ed., Baltimore, MD, Johns Hopkins University Press, 1996,
772             pg. 8.
773
774    Examples
775    --------
776    Make a (*very* coarse) grid for computing a Mandelbrot set:
777
778    >>> rl = np.outer(np.ones((5,)), np.linspace(-2, 2, 5))
779    >>> rl
780    array([[-2., -1.,  0.,  1.,  2.],
781           [-2., -1.,  0.,  1.,  2.],
782           [-2., -1.,  0.,  1.,  2.],
783           [-2., -1.,  0.,  1.,  2.],
784           [-2., -1.,  0.,  1.,  2.]])
785    >>> im = np.outer(1j*np.linspace(2, -2, 5), np.ones((5,)))
786    >>> im
787    array([[ 0.+2.j,  0.+2.j,  0.+2.j,  0.+2.j,  0.+2.j],
788           [ 0.+1.j,  0.+1.j,  0.+1.j,  0.+1.j,  0.+1.j],
789           [ 0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j],
790           [ 0.-1.j,  0.-1.j,  0.-1.j,  0.-1.j,  0.-1.j],
791           [ 0.-2.j,  0.-2.j,  0.-2.j,  0.-2.j,  0.-2.j]])
792    >>> grid = rl + im
793    >>> grid
794    array([[-2.+2.j, -1.+2.j,  0.+2.j,  1.+2.j,  2.+2.j],
795           [-2.+1.j, -1.+1.j,  0.+1.j,  1.+1.j,  2.+1.j],
796           [-2.+0.j, -1.+0.j,  0.+0.j,  1.+0.j,  2.+0.j],
797           [-2.-1.j, -1.-1.j,  0.-1.j,  1.-1.j,  2.-1.j],
798           [-2.-2.j, -1.-2.j,  0.-2.j,  1.-2.j,  2.-2.j]])
799
800    An example using a "vector" of letters:
801
802    >>> x = np.array(['a', 'b', 'c'], dtype=object)
803    >>> np.outer(x, [1, 2, 3])
804    array([[a, aa, aaa],
805           [b, bb, bbb],
806           [c, cc, ccc]], dtype=object)
807
808    """
809    a = asarray(a)
810    b = asarray(b)
811    return a.ravel()[:,newaxis]*b.ravel()[newaxis,:]
812
813# try to import blas optimized dot if available
814try:
815    # importing this changes the dot function for basic 4 types
816    # to blas-optimized versions.
817    from _dotblas import dot, vdot, inner, alterdot, restoredot
818except ImportError:
819    # docstrings are in add_newdocs.py
820    inner = multiarray.inner
821    dot = multiarray.dot
822    def vdot(a, b):
823        return dot(asarray(a).ravel().conj(), asarray(b).ravel())
824    def alterdot():
825        pass
826    def restoredot():
827        pass
828
829def tensordot(a, b, axes=2):
830    """
831    Compute tensor dot product along specified axes for arrays >= 1-D.
832
833    Given two tensors (arrays of dimension greater than or equal to one),
834    ``a`` and ``b``, and an array_like object containing two array_like
835    objects, ``(a_axes, b_axes)``, sum the products of ``a``'s and ``b``'s
836    elements (components) over the axes specified by ``a_axes`` and
837    ``b_axes``. The third argument can be a single non-negative
838    integer_like scalar, ``N``; if it is such, then the last ``N``
839    dimensions of ``a`` and the first ``N`` dimensions of ``b`` are summed
840    over.
841
842    Parameters
843    ----------
844    a, b : array_like, len(shape) >= 1
845        Tensors to "dot".
846
847    axes : variable type
848
849    * integer_like scalar
850      Number of axes to sum over (applies to both arrays); or
851
852    * array_like, shape = (2,), both elements array_like
853      Axes to be summed over, first sequence applying to ``a``, second
854      to ``b``.
855
856    See Also
857    --------
858    dot, einsum
859
860    Notes
861    -----
862    When there is more than one axis to sum over - and they are not the last
863    (first) axes of ``a`` (``b``) - the argument ``axes`` should consist of
864    two sequences of the same length, with the first axis to sum over given
865    first in both sequences, the second axis second, and so forth.
866
867    Examples
868    --------
869    A "traditional" example:
870
871    >>> a = np.arange(60.).reshape(3,4,5)
872    >>> b = np.arange(24.).reshape(4,3,2)
873    >>> c = np.tensordot(a,b, axes=([1,0],[0,1]))
874    >>> c.shape
875    (5, 2)
876    >>> c
877    array([[ 4400.,  4730.],
878           [ 4532.,  4874.],
879           [ 4664.,  5018.],
880           [ 4796.,  5162.],
881           [ 4928.,  5306.]])
882    >>> # A slower but equivalent way of computing the same...
883    >>> d = np.zeros((5,2))
884    >>> for i in range(5):
885    ...   for j in range(2):
886    ...     for k in range(3):
887    ...       for n in range(4):
888    ...         d[i,j] += a[k,n,i] * b[n,k,j]
889    >>> c == d
890    array([[ True,  True],
891           [ True,  True],
892           [ True,  True],
893           [ True,  True],
894           [ True,  True]], dtype=bool)
895
896    An extended example taking advantage of the overloading of + and \\*:
897
898    >>> a = np.array(range(1, 9))
899    >>> a.shape = (2, 2, 2)
900    >>> A = np.array(('a', 'b', 'c', 'd'), dtype=object)
901    >>> A.shape = (2, 2)
902    >>> a; A
903    array([[[1, 2],
904            [3, 4]],
905           [[5, 6],
906            [7, 8]]])
907    array([[a, b],
908           [c, d]], dtype=object)
909
910    >>> np.tensordot(a, A) # third argument default is 2
911    array([abbcccdddd, aaaaabbbbbbcccccccdddddddd], dtype=object)
912
913    >>> np.tensordot(a, A, 1)
914    array([[[acc, bdd],
915            [aaacccc, bbbdddd]],
916           [[aaaaacccccc, bbbbbdddddd],
917            [aaaaaaacccccccc, bbbbbbbdddddddd]]], dtype=object)
918
919    >>> np.tensordot(a, A, 0) # "Left for reader" (result too long to incl.)
920    array([[[[[a, b],
921              [c, d]],
922              ...
923
924    >>> np.tensordot(a, A, (0, 1))
925    array([[[abbbbb, cddddd],
926            [aabbbbbb, ccdddddd]],
927           [[aaabbbbbbb, cccddddddd],
928            [aaaabbbbbbbb, ccccdddddddd]]], dtype=object)
929
930    >>> np.tensordot(a, A, (2, 1))
931    array([[[abb, cdd],
932            [aaabbbb, cccdddd]],
933           [[aaaaabbbbbb, cccccdddddd],
934            [aaaaaaabbbbbbbb, cccccccdddddddd]]], dtype=object)
935
936    >>> np.tensordot(a, A, ((0, 1), (0, 1)))
937    array([abbbcccccddddddd, aabbbbccccccdddddddd], dtype=object)
938
939    >>> np.tensordot(a, A, ((2, 1), (1, 0)))
940    array([acccbbdddd, aaaaacccccccbbbbbbdddddddd], dtype=object)
941
942    """
943    try:
944        iter(axes)
945    except:
946        axes_a = range(-axes,0)
947        axes_b = range(0,axes)
948    else:
949        axes_a, axes_b = axes
950    try:
951        na = len(axes_a)
952        axes_a = list(axes_a)
953    except TypeError:
954        axes_a = [axes_a]
955        na = 1
956    try:
957        nb = len(axes_b)
958        axes_b = list(axes_b)
959    except TypeError:
960        axes_b = [axes_b]
961        nb = 1
962
963    a, b = asarray(a), asarray(b)
964    as_ = a.shape
965    nda = len(a.shape)
966    bs = b.shape
967    ndb = len(b.shape)
968    equal = True
969    if (na != nb): equal = False
970    else:
971        for k in xrange(na):
972            if as_[axes_a[k]] != bs[axes_b[k]]:
973                equal = False
974                break
975            if axes_a[k] < 0:
976                axes_a[k] += nda
977            if axes_b[k] < 0:
978                axes_b[k] += ndb
979    if not equal:
980        raise ValueError, "shape-mismatch for sum"
981
982    # Move the axes to sum over to the end of "a"
983    # and to the front of "b"
984    notin = [k for k in range(nda) if k not in axes_a]
985    newaxes_a = notin + axes_a
986    N2 = 1
987    for axis in axes_a:
988        N2 *= as_[axis]
989    newshape_a = (-1, N2)
990    olda = [as_[axis] for axis in notin]
991
992    notin = [k for k in range(ndb) if k not in axes_b]
993    newaxes_b = axes_b + notin
994    N2 = 1
995    for axis in axes_b:
996        N2 *= bs[axis]
997    newshape_b = (N2, -1)
998    oldb = [bs[axis] for axis in notin]
999
1000    at = a.transpose(newaxes_a).reshape(newshape_a)
1001    bt = b.transpose(newaxes_b).reshape(newshape_b)
1002    res = dot(at, bt)
1003    return res.reshape(olda + oldb)
1004
1005def roll(a, shift, axis=None):
1006    """
1007    Roll array elements along a given axis.
1008
1009    Elements that roll beyond the last position are re-introduced at
1010    the first.
1011
1012    Parameters
1013    ----------
1014    a : array_like
1015        Input array.
1016    shift : int
1017        The number of places by which elements are shifted.
1018    axis : int, optional
1019        The axis along which elements are shifted.  By default, the array
1020        is flattened before shifting, after which the original
1021        shape is restored.
1022
1023    Returns
1024    -------
1025    res : ndarray
1026        Output array, with the same shape as `a`.
1027
1028    See Also
1029    --------
1030    rollaxis : Roll the specified axis backwards, until it lies in a
1031               given position.
1032
1033    Examples
1034    --------
1035    >>> x = np.arange(10)
1036    >>> np.roll(x, 2)
1037    array([8, 9, 0, 1, 2, 3, 4, 5, 6, 7])
1038
1039    >>> x2 = np.reshape(x, (2,5))
1040    >>> x2
1041    array([[0, 1, 2, 3, 4],
1042           [5, 6, 7, 8, 9]])
1043    >>> np.roll(x2, 1)
1044    array([[9, 0, 1, 2, 3],
1045           [4, 5, 6, 7, 8]])
1046    >>> np.roll(x2, 1, axis=0)
1047    array([[5, 6, 7, 8, 9],
1048           [0, 1, 2, 3, 4]])
1049    >>> np.roll(x2, 1, axis=1)
1050    array([[4, 0, 1, 2, 3],
1051           [9, 5, 6, 7, 8]])
1052
1053    """
1054    a = asanyarray(a)
1055    if axis is None:
1056        n = a.size
1057        reshape = True
1058    else:
1059        n = a.shape[axis]
1060        reshape = False
1061    shift %= n
1062    indexes = concatenate((arange(n-shift,n),arange(n-shift)))
1063    res = a.take(indexes, axis)
1064    if reshape:
1065        return res.reshape(a.shape)
1066    else:
1067        return res
1068
1069def rollaxis(a, axis, start=0):
1070    """
1071    Roll the specified axis backwards, until it lies in a given position.
1072
1073    Parameters
1074    ----------
1075    a : ndarray
1076        Input array.
1077    axis : int
1078        The axis to roll backwards.  The positions of the other axes do not
1079        change relative to one another.
1080    start : int, optional
1081        The axis is rolled until it lies before this position.  The default,
1082        0, results in a "complete" roll.
1083
1084    Returns
1085    -------
1086    res : ndarray
1087        Output array.
1088
1089    See Also
1090    --------
1091    roll : Roll the elements of an array by a number of positions along a
1092        given axis.
1093
1094    Examples
1095    --------
1096    >>> a = np.ones((3,4,5,6))
1097    >>> np.rollaxis(a, 3, 1).shape
1098    (3, 6, 4, 5)
1099    >>> np.rollaxis(a, 2).shape
1100    (5, 3, 4, 6)
1101    >>> np.rollaxis(a, 1, 4).shape
1102    (3, 5, 6, 4)
1103
1104    """
1105    n = a.ndim
1106    if axis < 0:
1107        axis += n
1108    if start < 0:
1109        start += n
1110    msg = 'rollaxis: %s (%d) must be >=0 and < %d'
1111    if not (0 <= axis < n):
1112        raise ValueError, msg % ('axis', axis, n)
1113    if not (0 <= start < n+1):
1114        raise ValueError, msg % ('start', start, n+1)
1115    if (axis < start): # it's been removed
1116        start -= 1
1117    if axis==start:
1118        return a
1119    axes = range(0,n)
1120    axes.remove(axis)
1121    axes.insert(start, axis)
1122    return a.transpose(axes)
1123
1124# fix hack in scipy which imports this function
1125def _move_axis_to_0(a, axis):
1126    return rollaxis(a, axis, 0)
1127
1128def cross(a, b, axisa=-1, axisb=-1, axisc=-1, axis=None):
1129    """
1130    Return the cross product of two (arrays of) vectors.
1131
1132    The cross product of `a` and `b` in :math:`R^3` is a vector perpendicular
1133    to both `a` and `b`.  If `a` and `b` are arrays of vectors, the vectors
1134    are defined by the last axis of `a` and `b` by default, and these axes
1135    can have dimensions 2 or 3.  Where the dimension of either `a` or `b` is
1136    2, the third component of the input vector is assumed to be zero and the
1137    cross product calculated accordingly.  In cases where both input vectors
1138    have dimension 2, the z-component of the cross product is returned.
1139
1140    Parameters
1141    ----------
1142    a : array_like
1143        Components of the first vector(s).
1144    b : array_like
1145        Components of the second vector(s).
1146    axisa : int, optional
1147        Axis of `a` that defines the vector(s).  By default, the last axis.
1148    axisb : int, optional
1149        Axis of `b` that defines the vector(s).  By default, the last axis.
1150    axisc : int, optional
1151        Axis of `c` containing the cross product vector(s).  By default, the
1152        last axis.
1153    axis : int, optional
1154        If defined, the axis of `a`, `b` and `c` that defines the vector(s)
1155        and cross product(s).  Overrides `axisa`, `axisb` and `axisc`.
1156
1157    Returns
1158    -------
1159    c : ndarray
1160        Vector cross product(s).
1161
1162    Raises
1163    ------
1164    ValueError
1165        When the dimension of the vector(s) in `a` and/or `b` does not
1166        equal 2 or 3.
1167
1168    See Also
1169    --------
1170    inner : Inner product
1171    outer : Outer product.
1172    ix_ : Construct index arrays.
1173
1174    Examples
1175    --------
1176    Vector cross-product.
1177
1178    >>> x = [1, 2, 3]
1179    >>> y = [4, 5, 6]
1180    >>> np.cross(x, y)
1181    array([-3,  6, -3])
1182
1183    One vector with dimension 2.
1184
1185    >>> x = [1, 2]
1186    >>> y = [4, 5, 6]
1187    >>> np.cross(x, y)
1188    array([12, -6, -3])
1189
1190    Equivalently:
1191
1192    >>> x = [1, 2, 0]
1193    >>> y = [4, 5, 6]
1194    >>> np.cross(x, y)
1195    array([12, -6, -3])
1196
1197    Both vectors with dimension 2.
1198
1199    >>> x = [1,2]
1200    >>> y = [4,5]
1201    >>> np.cross(x, y)
1202    -3
1203
1204    Multiple vector cross-products. Note that the direction of the cross
1205    product vector is defined by the `right-hand rule`.
1206
1207    >>> x = np.array([[1,2,3], [4,5,6]])
1208    >>> y = np.array([[4,5,6], [1,2,3]])
1209    >>> np.cross(x, y)
1210    array([[-3,  6, -3],
1211           [ 3, -6,  3]])
1212
1213    The orientation of `c` can be changed using the `axisc` keyword.
1214
1215    >>> np.cross(x, y, axisc=0)
1216    array([[-3,  3],
1217           [ 6, -6],
1218           [-3,  3]])
1219
1220    Change the vector definition of `x` and `y` using `axisa` and `axisb`.
1221
1222    >>> x = np.array([[1,2,3], [4,5,6], [7, 8, 9]])
1223    >>> y = np.array([[7, 8, 9], [4,5,6], [1,2,3]])
1224    >>> np.cross(x, y)
1225    array([[ -6,  12,  -6],
1226           [  0,   0,   0],
1227           [  6, -12,   6]])
1228    >>> np.cross(x, y, axisa=0, axisb=0)
1229    array([[-24,  48, -24],
1230           [-30,  60, -30],
1231           [-36,  72, -36]])
1232
1233    """
1234    if axis is not None:
1235        axisa,axisb,axisc=(axis,)*3
1236    a = asarray(a).swapaxes(axisa, 0)
1237    b = asarray(b).swapaxes(axisb, 0)
1238    msg = "incompatible dimensions for cross product\n"\
1239          "(dimension must be 2 or 3)"
1240    if (a.shape[0] not in [2,3]) or (b.shape[0] not in [2,3]):
1241        raise ValueError(msg)
1242    if a.shape[0] == 2:
1243        if (b.shape[0] == 2):
1244            cp = a[0]*b[1] - a[1]*b[0]
1245            if cp.ndim == 0:
1246                return cp
1247            else:
1248                return cp.swapaxes(0, axisc)
1249        else:
1250            x = a[1]*b[2]
1251            y = -a[0]*b[2]
1252            z = a[0]*b[1] - a[1]*b[0]
1253    elif a.shape[0] == 3:
1254        if (b.shape[0] == 3):
1255            x = a[1]*b[2] - a[2]*b[1]
1256            y = a[2]*b[0] - a[0]*b[2]
1257            z = a[0]*b[1] - a[1]*b[0]
1258        else:
1259            x = -a[2]*b[1]
1260            y = a[2]*b[0]
1261            z = a[0]*b[1] - a[1]*b[0]
1262    cp = array([x,y,z])
1263    if cp.ndim == 1:
1264        return cp
1265    else:
1266        return cp.swapaxes(0,axisc)
1267
1268
1269#Use numarray's printing function
1270from arrayprint import array2string, get_printoptions, set_printoptions
1271
1272_typelessdata = [int_, float_, complex_]
1273if issubclass(intc, int):
1274    _typelessdata.append(intc)
1275
1276if issubclass(longlong, int):
1277    _typelessdata.append(longlong)
1278
1279def array_repr(arr, max_line_width=None, precision=None, suppress_small=None):
1280    """
1281    Return the string representation of an array.
1282
1283    Parameters
1284    ----------
1285    arr : ndarray
1286        Input array.
1287    max_line_width : int, optional
1288        The maximum number of columns the string should span. Newline
1289        characters split the string appropriately after array elements.
1290    precision : int, optional
1291        Floating point precision. Default is the current printing precision
1292        (usually 8), which can be altered using `set_printoptions`.
1293    suppress_small : bool, optional
1294        Represent very small numbers as zero, default is False. Very small
1295        is defined by `precision`, if the precision is 8 then
1296        numbers smaller than 5e-9 are represented as zero.
1297
1298    Returns
1299    -------
1300    string : str
1301      The string representation of an array.
1302
1303    See Also
1304    --------
1305    array_str, array2string, set_printoptions
1306
1307    Examples
1308    --------
1309    >>> np.array_repr(np.array([1,2]))
1310    'array([1, 2])'
1311    >>> np.array_repr(np.ma.array([0.]))
1312    'MaskedArray([ 0.])'
1313    >>> np.array_repr(np.array([], np.int32))
1314    'array([], dtype=int32)'
1315
1316    >>> x = np.array([1e-6, 4e-7, 2, 3])
1317    >>> np.array_repr(x, precision=6, suppress_small=True)
1318    'array([ 0.000001,  0.      ,  2.      ,  3.      ])'
1319
1320    """
1321    if arr.size > 0 or arr.shape==(0,):
1322        lst = array2string(arr, max_line_width, precision, suppress_small,
1323                           ', ', "array(")
1324    else: # show zero-length shape unless it is (0,)
1325        lst = "[], shape=%s" % (repr(arr.shape),)
1326    typeless = arr.dtype.type in _typelessdata
1327
1328    if arr.__class__ is not ndarray:
1329        cName= arr.__class__.__name__
1330    else:
1331        cName = "array"
1332    if typeless and arr.size:
1333        return cName + "(%s)" % lst
1334    else:
1335        typename=arr.dtype.name
1336        lf = ''
1337        if issubclass(arr.dtype.type, flexible):
1338            if arr.dtype.names:
1339                typename = "%s" % str(arr.dtype)
1340            else:
1341                typename = "'%s'" % str(arr.dtype)
1342            lf = '\n'+' '*len("array(")
1343        return cName + "(%s, %sdtype=%s)" % (lst, lf, typename)
1344
1345def array_str(a, max_line_width=None, precision=None, suppress_small=None):
1346    """
1347    Return a string representation of the data in an array.
1348
1349    The data in the array is returned as a single string.  This function is
1350    similar to `array_repr`, the difference being that `array_repr` also
1351    returns information on the kind of array and its data type.
1352
1353    Parameters
1354    ----------
1355    a : ndarray
1356        Input array.
1357    max_line_width : int, optional
1358        Inserts newlines if text is longer than `max_line_width`.  The
1359        default is, indirectly, 75.
1360    precision : int, optional
1361        Floating point precision.  Default is the current printing precision
1362        (usually 8), which can be altered using `set_printoptions`.
1363    suppress_small : bool, optional
1364        Represent numbers "very close" to zero as zero; default is False.
1365        Very close is defined by precision: if the precision is 8, e.g.,
1366        numbers smaller (in absolute value) than 5e-9 are represented as
1367        zero.
1368
1369    See Also
1370    --------
1371    array2string, array_repr, set_printoptions
1372
1373    Examples
1374    --------
1375    >>> np.array_str(np.arange(3))
1376    '[0 1 2]'
1377
1378    """
1379    return array2string(a, max_line_width, precision, suppress_small, ' ', "", str)
1380
1381def set_string_function(f, repr=True):
1382    """
1383    Set a Python function to be used when pretty printing arrays.
1384
1385    Parameters
1386    ----------
1387    f : function or None
1388        Function to be used to pretty print arrays. The function should expect
1389        a single array argument and return a string of the representation of
1390        the array. If None, the function is reset to the default NumPy function
1391        to print arrays.
1392    repr : bool, optional
1393        If True (default), the function for pretty printing (``__repr__``)
1394        is set, if False the function that returns the default string
1395        representation (``__str__``) is set.
1396
1397    See Also
1398    --------
1399    set_printoptions, get_printoptions
1400
1401    Examples
1402    --------
1403    >>> def pprint(arr):
1404    ...     return 'HA! - What are you going to do now?'
1405    ...
1406    >>> np.set_string_function(pprint)
1407    >>> a = np.arange(10)
1408    >>> a
1409    HA! - What are you going to do now?
1410    >>> print a
1411    [0 1 2 3 4 5 6 7 8 9]
1412
1413    We can reset the function to the default:
1414
1415    >>> np.set_string_function(None)
1416    >>> a
1417    array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
1418
1419    `repr` affects either pretty printing or normal string representation.
1420    Note that ``__repr__`` is still affected by setting ``__str__``
1421    because the width of each array element in the returned string becomes
1422    equal to the length of the result of ``__str__()``.
1423
1424    >>> x = np.arange(4)
1425    >>> np.set_string_function(lambda x:'random', repr=False)
1426    >>> x.__str__()
1427    'random'
1428    >>> x.__repr__()
1429    'array([     0,      1,      2,      3])'
1430
1431    """
1432    if f is None:
1433        if repr:
1434            return multiarray.set_string_function(array_repr, 1)
1435        else:
1436            return multiarray.set_string_function(array_str, 0)
1437    else:
1438        return multiarray.set_string_function(f, repr)
1439
1440set_string_function(array_str, 0)
1441set_string_function(array_repr, 1)
1442
1443little_endian = (sys.byteorder == 'little')
1444
1445
1446def indices(dimensions, dtype=int):
1447    """
1448    Return an array representing the indices of a grid.
1449
1450    Compute an array where the subarrays contain index values 0,1,...
1451    varying only along the corresponding axis.
1452
1453    Parameters
1454    ----------
1455    dimensions : sequence of ints
1456        The shape of the grid.
1457    dtype : dtype, optional
1458        Data type of the result.
1459
1460    Returns
1461    -------
1462    grid : ndarray
1463        The array of grid indices,
1464        ``grid.shape = (len(dimensions),) + tuple(dimensions)``.
1465
1466    See Also
1467    --------
1468    mgrid, meshgrid
1469
1470    Notes
1471    -----
1472    The output shape is obtained by prepending the number of dimensions
1473    in front of the tuple of dimensions, i.e. if `dimensions` is a tuple
1474    ``(r0, ..., rN-1)`` of length ``N``, the output shape is
1475    ``(N,r0,...,rN-1)``.
1476
1477    The subarrays ``grid[k]`` contains the N-D array of indices along the
1478    ``k-th`` axis. Explicitly::
1479
1480        grid[k,i0,i1,...,iN-1] = ik
1481
1482    Examples
1483    --------
1484    >>> grid = np.indices((2, 3))
1485    >>> grid.shape
1486    (2, 2, 3)
1487    >>> grid[0]        # row indices
1488    array([[0, 0, 0],
1489           [1, 1, 1]])
1490    >>> grid[1]        # column indices
1491    array([[0, 1, 2],
1492           [0, 1, 2]])
1493
1494    The indices can be used as an index into an array.
1495
1496    >>> x = np.arange(20).reshape(5, 4)
1497    >>> row, col = np.indices((2, 3))
1498    >>> x[row, col]
1499    array([[0, 1, 2],
1500           [4, 5, 6]])
1501
1502    Note that it would be more straightforward in the above example to
1503    extract the required elements directly with ``x[:2, :3]``.
1504
1505    """
1506    dimensions = tuple(dimensions)
1507    N = len(dimensions)
1508    if N == 0:
1509        return array([],dtype=dtype)
1510    res = empty((N,)+dimensions, dtype=dtype)
1511    for i, dim in enumerate(dimensions):
1512        tmp = arange(dim,dtype=dtype)
1513        tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
1514        newdim = dimensions[:i] + (1,)+ dimensions[i+1:]
1515        val = zeros(newdim, dtype)
1516        add(tmp, val, res[i])
1517    return res
1518
1519def fromfunction(function, shape, **kwargs):
1520    """
1521    Construct an array by executing a function over each coordinate.
1522
1523    The resulting array therefore has a value ``fn(x, y, z)`` at
1524    coordinate ``(x, y, z)``.
1525
1526    Parameters
1527    ----------
1528    function : callable
1529        The function is called with N parameters, each of which
1530        represents the coordinates of the array varying along a
1531        specific axis.  For example, if `shape` were ``(2, 2)``, then
1532        the parameters would be two arrays, ``[[0, 0], [1, 1]]`` and
1533        ``[[0, 1], [0, 1]]``.  `function` must be capable of operating on
1534        arrays, and should return a scalar value.
1535    shape : (N,) tuple of ints
1536        Shape of the output array, which also determines the shape of
1537        the coordinate arrays passed to `function`.
1538    dtype : data-type, optional
1539        Data-type of the coordinate arrays passed to `function`.
1540        By default, `dtype` is float.
1541
1542    Returns
1543    -------
1544    out : any
1545        The result of the call to `function` is passed back directly.
1546        Therefore the type and shape of `out` is completely determined by
1547        `function`.
1548
1549    See Also
1550    --------
1551    indices, meshgrid
1552
1553    Notes
1554    -----
1555    Keywords other than `shape` and `dtype` are passed to `function`.
1556
1557    Examples
1558    --------
1559    >>> np.fromfunction(lambda i, j: i == j, (3, 3), dtype=int)
1560    array([[ True, False, False],
1561           [False,  True, False],
1562           [False, False,  True]], dtype=bool)
1563
1564    >>> np.fromfunction(lambda i, j: i + j, (3, 3), dtype=int)
1565    array([[0, 1, 2],
1566           [1, 2, 3],
1567           [2, 3, 4]])
1568
1569    """
1570    dtype = kwargs.pop('dtype', float)
1571    args = indices(shape, dtype=dtype)
1572    return function(*args,**kwargs)
1573
1574def isscalar(num):
1575    """
1576    Returns True if the type of `num` is a scalar type.
1577
1578    Parameters
1579    ----------
1580    num : any
1581        Input argument, can be of any type and shape.
1582
1583    Returns
1584    -------
1585    val : bool
1586        True if `num` is a scalar type, False if it is not.
1587
1588    Examples
1589    --------
1590    >>> np.isscalar(3.1)
1591    True
1592    >>> np.isscalar([3.1])
1593    False
1594    >>> np.isscalar(False)
1595    True
1596
1597    """
1598    if isinstance(num, generic):
1599        return True
1600    else:
1601        return type(num) in ScalarType
1602
1603_lkup = {
1604    '0':'0000',
1605    '1':'0001',
1606    '2':'0010',
1607    '3':'0011',
1608    '4':'0100',
1609    '5':'0101',
1610    '6':'0110',
1611    '7':'0111',
1612    '8':'1000',
1613    '9':'1001',
1614    'a':'1010',
1615    'b':'1011',
1616    'c':'1100',
1617    'd':'1101',
1618    'e':'1110',
1619    'f':'1111',
1620    'A':'1010',
1621    'B':'1011',
1622    'C':'1100',
1623    'D':'1101',
1624    'E':'1110',
1625    'F':'1111',
1626    'L':''}
1627
1628def binary_repr(num, width=None):
1629    """
1630    Return the binary representation of the input number as a string.
1631
1632    For negative numbers, if width is not given, a minus sign is added to the
1633    front. If width is given, the two's complement of the number is
1634    returned, with respect to that width.
1635
1636    In a two's-complement system negative numbers are represented by the two's
1637    complement of the absolute value. This is the most common method of
1638    representing signed integers on computers [1]_. A N-bit two's-complement
1639    system can represent every integer in the range
1640    :math:`-2^{N-1}` to :math:`+2^{N-1}-1`.
1641
1642    Parameters
1643    ----------
1644    num : int
1645        Only an integer decimal number can be used.
1646    width : int, optional
1647        The length of the returned string if `num` is positive, the length of
1648        the two's complement if `num` is negative.
1649
1650    Returns
1651    -------
1652    bin : str
1653        Binary representation of `num` or two's complement of `num`.
1654
1655    See Also
1656    --------
1657    base_repr: Return a string representation of a number in the given base
1658               system.
1659
1660    Notes
1661    -----
1662    `binary_repr` is equivalent to using `base_repr` with base 2, but about 25x
1663    faster.
1664
1665    References
1666    ----------
1667    .. [1] Wikipedia, "Two's complement",
1668        http://en.wikipedia.org/wiki/Two's_complement
1669
1670    Examples
1671    --------
1672    >>> np.binary_repr(3)
1673    '11'
1674    >>> np.binary_repr(-3)
1675    '-11'
1676    >>> np.binary_repr(3, width=4)
1677    '0011'
1678
1679    The two's complement is returned when the input number is negative and
1680    width is specified:
1681
1682    >>> np.binary_repr(-3, width=4)
1683    '1101'
1684
1685    """
1686    sign = ''
1687    if num < 0:
1688        if width is None:
1689            sign = '-'
1690            num = -num
1691        else:
1692            # replace num with its 2-complement
1693            num = 2**width + num
1694    elif num == 0:
1695        return '0'*(width or 1)
1696    ostr = hex(num)
1697    bin = ''.join([_lkup[ch] for ch in ostr[2:]])
1698    bin = bin.lstrip('0')
1699    if width is not None:
1700        bin = bin.zfill(width)
1701    return sign + bin
1702
1703def base_repr(number, base=2, padding=0):
1704    """
1705    Return a string representation of a number in the given base system.
1706
1707    Parameters
1708    ----------
1709    number : int
1710        The value to convert. Only positive values are handled.
1711    base : int, optional
1712        Convert `number` to the `base` number system. The valid range is 2-36,
1713        the default value is 2.
1714    padding : int, optional
1715        Number of zeros padded on the left. Default is 0 (no padding).
1716
1717    Returns
1718    -------
1719    out : str
1720        String representation of `number` in `base` system.
1721
1722    See Also
1723    --------
1724    binary_repr : Faster version of `base_repr` for base 2.
1725
1726    Examples
1727    --------
1728    >>> np.base_repr(5)
1729    '101'
1730    >>> np.base_repr(6, 5)
1731    '11'
1732    >>> np.base_repr(7, base=5, padding=3)
1733    '00012'
1734
1735    >>> np.base_repr(10, base=16)
1736    'A'
1737    >>> np.base_repr(32, base=16)
1738    '20'
1739
1740    """
1741    digits = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'
1742    if base > len(digits):
1743        raise ValueError("Bases greater than 36 not handled in base_repr.")
1744
1745    num = abs(number)
1746    res = []
1747    while num:
1748        res.append(digits[num % base])
1749        num //= base
1750    if padding:
1751        res.append('0' * padding)
1752    if number < 0:
1753        res.append('-')
1754    return ''.join(reversed(res or '0'))
1755
1756from cPickle import load, loads
1757_cload = load
1758_file = open
1759
1760def load(file):
1761    """
1762    Wrapper around cPickle.load which accepts either a file-like object or
1763    a filename.
1764
1765    Note that the NumPy binary format is not based on pickle/cPickle anymore.
1766    For details on the preferred way of loading and saving files, see `load`
1767    and `save`.
1768
1769    See Also
1770    --------
1771    load, save
1772
1773    """
1774    if isinstance(file, type("")):
1775        file = _file(file,"rb")
1776    return _cload(file)
1777
1778# These are all essentially abbreviations
1779# These might wind up in a special abbreviations module
1780
1781def _maketup(descr, val):
1782    dt = dtype(descr)
1783    # Place val in all scalar tuples:
1784    fields = dt.fields
1785    if fields is None:
1786        return val
1787    else:
1788        res = [_maketup(fields[name][0],val) for name in dt.names]
1789        return tuple(res)
1790
1791def ones(shape, dtype=None, order='C'):
1792    """
1793    Return a new array of given shape and type, filled with ones.
1794
1795    Please refer to the documentation for `zeros` for further details.
1796
1797    See Also
1798    --------
1799    zeros, ones_like
1800
1801    Examples
1802    --------
1803    >>> np.ones(5)
1804    array([ 1.,  1.,  1.,  1.,  1.])
1805
1806    >>> np.ones((5,), dtype=np.int)
1807    array([1, 1, 1, 1, 1])
1808
1809    >>> np.ones((2, 1))
1810    array([[ 1.],
1811           [ 1.]])
1812
1813    >>> s = (2,2)
1814    >>> np.ones(s)
1815    array([[ 1.,  1.],
1816           [ 1.,  1.]])
1817
1818    """
1819    a = empty(shape, dtype, order)
1820    try:
1821        a.fill(1)
1822        # Above is faster now after addition of fast loops.
1823        #a = zeros(shape, dtype, order)
1824        #a+=1
1825    except TypeError:
1826        obj = _maketup(dtype, 1)
1827        a.fill(obj)
1828    return a
1829
1830def identity(n, dtype=None):
1831    """
1832    Return the identity array.
1833
1834    The identity array is a square array with ones on
1835    the main diagonal.
1836
1837    Parameters
1838    ----------
1839    n : int
1840        Number of rows (and columns) in `n` x `n` output.
1841    dtype : data-type, optional
1842        Data-type of the output.  Defaults to ``float``.
1843
1844    Returns
1845    -------
1846    out : ndarray
1847        `n` x `n` array with its main diagonal set to one,
1848        and all other elements 0.
1849
1850    Examples
1851    --------
1852    >>> np.identity(3)
1853    array([[ 1.,  0.,  0.],
1854           [ 0.,  1.,  0.],
1855           [ 0.,  0.,  1.]])
1856
1857    """
1858    a = zeros((n,n), dtype=dtype)
1859    a.flat[::n+1] = 1
1860    return a
1861
1862def allclose(a, b, rtol=1.e-5, atol=1.e-8):
1863    """
1864    Returns True if two arrays are element-wise equal within a tolerance.
1865
1866    The tolerance values are positive, typically very small numbers.  The
1867    relative difference (`rtol` * abs(`b`)) and the absolute difference
1868    `atol` are added together to compare against the absolute difference
1869    between `a` and `b`.
1870
1871    If either array contains one or more NaNs, False is returned.
1872    Infs are treated as equal if they are in the same place and of the same
1873    sign in both arrays.
1874
1875    Parameters
1876    ----------
1877    a, b : array_like
1878        Input arrays to compare.
1879    rtol : float
1880        The relative tolerance parameter (see Notes).
1881    atol : float
1882        The absolute tolerance parameter (see Notes).
1883
1884    Returns
1885    -------
1886    y : bool
1887        Returns True if the two arrays are equal within the given
1888        tolerance; False otherwise.
1889
1890    See Also
1891    --------
1892    all, any, alltrue, sometrue
1893
1894    Notes
1895    -----
1896    If the following equation is element-wise True, then allclose returns
1897    True.
1898
1899     absolute(`a` - `b`) <= (`atol` + `rtol` * absolute(`b`))
1900
1901    The above equation is not symmetric in `a` and `b`, so that
1902    `allclose(a, b)` might be different from `allclose(b, a)` in
1903    some rare cases.
1904
1905    Examples
1906    --------
1907    >>> np.allclose([1e10,1e-7], [1.00001e10,1e-8])
1908    False
1909    >>> np.allclose([1e10,1e-8], [1.00001e10,1e-9])
1910    True
1911    >>> np.allclose([1e10,1e-8], [1.0001e10,1e-9])
1912    False
1913    >>> np.allclose([1.0, np.nan], [1.0, np.nan])
1914    False
1915
1916    """
1917    x = array(a, copy=False, ndmin=1)
1918    y = array(b, copy=False, ndmin=1)
1919
1920    if any(isnan(x)) or any(isnan(y)):
1921        return False
1922
1923    xinf = isinf(x)
1924    yinf = isinf(y)
1925    if any(xinf) or any(yinf):
1926        # Check that x and y have inf's only in the same positions
1927        if not all(xinf == yinf):
1928            return False
1929        # Check that sign of inf's in x and y is the same
1930        if not all(x[xinf] == y[xinf]):
1931            return False
1932
1933        x = x[~xinf]
1934        y = y[~xinf]
1935
1936    return all(less_equal(abs(x-y), atol + rtol * abs(y)))
1937
1938def array_equal(a1, a2):
1939    """
1940    True if two arrays have the same shape and elements, False otherwise.
1941
1942    Parameters
1943    ----------
1944    a1, a2 : array_like
1945        Input arrays.
1946
1947    Returns
1948    -------
1949    b : bool
1950        Returns True if the arrays are equal.
1951
1952    See Also
1953    --------
1954    allclose: Returns True if two arrays are element-wise equal within a
1955              tolerance.
1956    array_equiv: Returns True if input arrays are shape consistent and all
1957                 elements equal.
1958
1959    Examples
1960    --------
1961    >>> np.array_equal([1, 2], [1, 2])
1962    True
1963    >>> np.array_equal(np.array([1, 2]), np.array([1, 2]))
1964    True
1965    >>> np.array_equal([1, 2], [1, 2, 3])
1966    False
1967    >>> np.array_equal([1, 2], [1, 4])
1968    False
1969
1970    """
1971    try:
1972        a1, a2 = asarray(a1), asarray(a2)
1973    except:
1974        return False
1975    if a1.shape != a2.shape:
1976        return False
1977    return bool(logical_and.reduce(equal(a1,a2).ravel()))
1978
1979def array_equiv(a1, a2):
1980    """
1981    Returns True if input arrays are shape consistent and all elements equal.
1982
1983    Shape consistent means they are either the same shape, or one input array
1984    can be broadcasted to create the same shape as the other one.
1985
1986    Parameters
1987    ----------
1988    a1, a2 : array_like
1989        Input arrays.
1990
1991    Returns
1992    -------
1993    out : bool
1994        True if equivalent, False otherwise.
1995
1996    Examples
1997    --------
1998    >>> np.array_equiv([1, 2], [1, 2])
1999    True
2000    >>> np.array_equiv([1, 2], [1, 3])
2001    False
2002
2003    Showing the shape equivalence:
2004
2005    >>> np.array_equiv([1, 2], [[1, 2], [1, 2]])
2006    True
2007    >>> np.array_equiv([1, 2], [[1, 2, 1, 2], [1, 2, 1, 2]])
2008    False
2009
2010    >>> np.array_equiv([1, 2], [[1, 2], [1, 3]])
2011    False
2012
2013    """
2014    try:
2015        a1, a2 = asarray(a1), asarray(a2)
2016    except:
2017        return False
2018    try:
2019        return bool(logical_and.reduce(equal(a1,a2).ravel()))
2020    except ValueError:
2021        return False
2022
2023
2024_errdict = {"ignore":ERR_IGNORE,
2025            "warn":ERR_WARN,
2026            "raise":ERR_RAISE,
2027            "call":ERR_CALL,
2028            "print":ERR_PRINT,
2029            "log":ERR_LOG}
2030
2031_errdict_rev = {}
2032for key in _errdict.keys():
2033    _errdict_rev[_errdict[key]] = key
2034del key
2035
2036def seterr(all=None, divide=None, over=None, under=None, invalid=None):
2037    """
2038    Set how floating-point errors are handled.
2039
2040    Note that operations on integer scalar types (such as `int16`) are
2041    handled like floating point, and are affected by these settings.
2042
2043    Parameters
2044    ----------
2045    all : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional
2046        Set treatment for all types of floating-point errors at once:
2047
2048        - ignore: Take no action when the exception occurs.
2049        - warn: Print a `RuntimeWarning` (via the Python `warnings` module).
2050        - raise: Raise a `FloatingPointError`.
2051        - call: Call a function specified using the `seterrcall` function.
2052        - print: Print a warning directly to ``stdout``.
2053        - log: Record error in a Log object specified by `seterrcall`.
2054
2055        The default is not to change the current behavior.
2056    divide : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional
2057        Treatment for division by zero.
2058    over : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional
2059        Treatment for floating-point overflow.
2060    under : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional
2061        Treatment for floating-point underflow.
2062    invalid : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional
2063        Treatment for invalid floating-point operation.
2064
2065    Returns
2066    -------
2067    old_settings : dict
2068        Dictionary containing the old settings.
2069
2070    See also
2071    --------
2072    seterrcall : Set a callback function for the 'call' mode.
2073    geterr, geterrcall
2074
2075    Notes
2076    -----
2077    The floating-point exceptions are defined in the IEEE 754 standard [1]:
2078
2079    - Division by zero: infinite result obtained from finite numbers.
2080    - Overflow: result too large to be expressed.
2081    - Underflow: result so close to zero that some precision
2082      was lost.
2083    - Invalid operation: result is not an expressible number, typically
2084      indicates that a NaN was produced.
2085
2086    .. [1] http://en.wikipedia.org/wiki/IEEE_754
2087
2088    Examples
2089    --------
2090    >>> old_settings = np.seterr(all='ignore')  #seterr to known value
2091    >>> np.seterr(over='raise')
2092    {'over': 'ignore', 'divide': 'ignore', 'invalid': 'ignore',
2093     'under': 'ignore'}
2094    >>> np.seterr(all='ignore')  # reset to default
2095    {'over': 'raise', 'divide': 'ignore', 'invalid': 'ignore', 'under': 'ignore'}
2096
2097    >>> np.int16(32000) * np.int16(3)
2098    30464
2099    >>> old_settings = np.seterr(all='warn', over='raise')
2100    >>> np.int16(32000) * np.int16(3)
2101    Traceback (most recent call last):
2102      File "<stdin>", line 1, in <module>
2103    FloatingPointError: overflow encountered in short_scalars
2104
2105    >>> old_settings = np.seterr(all='print')
2106    >>> np.geterr()
2107    {'over': 'print', 'divide': 'print', 'invalid': 'print', 'under': 'print'}
2108    >>> np.int16(32000) * np.int16(3)
2109    Warning: overflow encountered in short_scalars
2110    30464
2111
2112    """
2113
2114    pyvals = umath.geterrobj()
2115    old = geterr()
2116
2117    if divide is None: divide = all or old['divide']
2118    if over is None: over = all or old['over']
2119    if under is None: under = all or old['under']
2120    if invalid is None: invalid = all or old['invalid']
2121
2122    maskvalue = ((_errdict[divide] << SHIFT_DIVIDEBYZERO) +
2123                 (_errdict[over] << SHIFT_OVERFLOW ) +
2124                 (_errdict[under] << SHIFT_UNDERFLOW) +
2125                 (_errdict[invalid] << SHIFT_INVALID))
2126
2127    pyvals[1] = maskvalue
2128    umath.seterrobj(pyvals)
2129    return old
2130
2131
2132def geterr():
2133    """
2134    Get the current way of handling floating-point errors.
2135
2136    Returns
2137    -------
2138    res : dict
2139        A dictionary with keys "divide", "over", "under", and "invalid",
2140        whose values are from the strings "ignore", "print", "log", "warn",
2141        "raise", and "call". The keys represent possible floating-point
2142        exceptions, and the values define how these exceptions are handled.
2143
2144    See Also
2145    --------
2146    geterrcall, seterr, seterrcall
2147
2148    Notes
2149    -----
2150    For complete documentation of the types of floating-point exceptions and
2151    treatment options, see `seterr`.
2152
2153    Examples
2154    --------
2155    >>> np.geterr()
2156    {'over': 'warn', 'divide': 'warn', 'invalid': 'warn',
2157    'under': 'ignore'}
2158    >>> np.arange(3.) / np.arange(3.)
2159    array([ NaN,   1.,   1.])
2160
2161    >>> oldsettings = np.seterr(all='warn', over='raise')
2162    >>> np.geterr()
2163    {'over': 'raise', 'divide': 'warn', 'invalid': 'warn', 'under': 'warn'}
2164    >>> np.arange(3.) / np.arange(3.)
2165    __main__:1: RuntimeWarning: invalid value encountered in divide
2166    array([ NaN,   1.,   1.])
2167
2168    """
2169    maskvalue = umath.geterrobj()[1]
2170    mask = 7
2171    res = {}
2172    val = (maskvalue >> SHIFT_DIVIDEBYZERO) & mask
2173    res['divide'] = _errdict_rev[val]
2174    val = (maskvalue >> SHIFT_OVERFLOW) & mask
2175    res['over'] = _errdict_rev[val]
2176    val = (maskvalue >> SHIFT_UNDERFLOW) & mask
2177    res['under'] = _errdict_rev[val]
2178    val = (maskvalue >> SHIFT_INVALID) & mask
2179    res['invalid'] = _errdict_rev[val]
2180    return res
2181
2182def setbufsize(size):
2183    """
2184    Set the size of the buffer used in ufuncs.
2185
2186    Parameters
2187    ----------
2188    size : int
2189        Size of buffer.
2190
2191    """
2192    if size > 10e6:
2193        raise ValueError, "Buffer size, %s, is too big." % size
2194    if size < 5:
2195        raise ValueError, "Buffer size, %s, is too small." %size
2196    if size % 16 != 0:
2197        raise ValueError, "Buffer size, %s, is not a multiple of 16." %size
2198
2199    pyvals = umath.geterrobj()
2200    old = getbufsize()
2201    pyvals[0] = size
2202    umath.seterrobj(pyvals)
2203    return old
2204
2205def getbufsize():
2206    """Return the size of the buffer used in ufuncs.
2207    """
2208    return umath.geterrobj()[0]
2209
2210def seterrcall(func):
2211    """
2212    Set the floating-point error callback function or log object.
2213
2214    There are two ways to capture floating-point error messages.  The first
2215    is to set the error-handler to 'call', using `seterr`.  Then, set
2216    the function to call using this function.
2217
2218    The second is to set the error-handler to 'log', using `seterr`.
2219    Floating-point errors then trigger a call to the 'write' method of
2220    the provided object.
2221
2222    Parameters
2223    ----------
2224    func : callable f(err, flag) or object with write method
2225        Function to call upon floating-point errors ('call'-mode) or
2226        object whose 'write' method is used to log such message ('log'-mode).
2227
2228        The call function takes two arguments. The first is the
2229        type of error (one of "divide", "over", "under", or "invalid"),
2230        and the second is the status flag.  The flag is a byte, whose
2231        least-significant bits indicate the status::
2232
2233          [0 0 0 0 invalid over under invalid]
2234
2235        In other words, ``flags = divide + 2*over + 4*under + 8*invalid``.
2236
2237        If an object is provided, its write method should take one argument,
2238        a string.
2239
2240    Returns
2241    -------
2242    h : callable, log instance or None
2243        The old error handler.
2244
2245    See Also
2246    --------
2247    seterr, geterr, geterrcall
2248
2249    Examples
2250    --------
2251    Callback upon error:
2252
2253    >>> def err_handler(type, flag):
2254    ...     print "Floating point error (%s), with flag %s" % (type, flag)
2255    ...
2256
2257    >>> saved_handler = np.seterrcall(err_handler)
2258    >>> save_err = np.seterr(all='call')
2259
2260    >>> np.array([1, 2, 3]) / 0.0
2261    Floating point error (divide by zero), with flag 1
2262    array([ Inf,  Inf,  Inf])
2263
2264    >>> np.seterrcall(saved_handler)
2265    <function err_handler at 0x...>
2266    >>> np.seterr(**save_err)
2267    {'over': 'call', 'divide': 'call', 'invalid': 'call', 'under': 'call'}
2268
2269    Log error message:
2270
2271    >>> class Log(object):
2272    ...     def write(self, msg):
2273    ...         print "LOG: %s" % msg
2274    ...
2275
2276    >>> log = Log()
2277    >>> saved_handler = np.seterrcall(log)
2278    >>> save_err = np.seterr(all='log')
2279
2280    >>> np.array([1, 2, 3]) / 0.0
2281    LOG: Warning: divide by zero encountered in divide
2282    <BLANKLINE>
2283    array([ Inf,  Inf,  Inf])
2284
2285    >>> np.seterrcall(saved_handler)
2286    <__main__.Log object at 0x...>
2287    >>> np.seterr(**save_err)
2288    {'over': 'log', 'divide': 'log', 'invalid': 'log', 'under': 'log'}
2289
2290    """
2291    if func is not None and not callable(func):
2292        if not hasattr(func, 'write') or not callable(func.write):
2293            raise ValueError, "Only callable can be used as callback"
2294    pyvals = umath.geterrobj()
2295    old = geterrcall()
2296    pyvals[2] = func
2297    umath.seterrobj(pyvals)
2298    return old
2299
2300def geterrcall():
2301    """
2302    Return the current callback function used on floating-point errors.
2303
2304    When the error handling for a floating-point error (one of "divide",
2305    "over", "under", or "invalid") is set to 'call' or 'log', the function
2306    that is called or the log instance that is written to is returned by
2307    `geterrcall`. This function or log instance has been set with
2308    `seterrcall`.
2309
2310    Returns
2311    -------
2312    errobj : callable, log instance or None
2313        The current error handler. If no handler was set through `seterrcall`,
2314        ``None`` is returned.
2315
2316    See Also
2317    --------
2318    seterrcall, seterr, geterr
2319
2320    Notes
2321    -----
2322    For complete documentation of the types of floating-point exceptions and
2323    treatment options, see `seterr`.
2324
2325    Examples
2326    --------
2327    >>> np.geterrcall()  # we did not yet set a handler, returns None
2328
2329    >>> oldsettings = np.seterr(all='call')
2330    >>> def err_handler(type, flag):
2331    ...     print "Floating point error (%s), with flag %s" % (type, flag)
2332    >>> oldhandler = np.seterrcall(err_handler)
2333    >>> np.array([1, 2, 3]) / 0.0
2334    Floating point error (divide by zero), with flag 1
2335    array([ Inf,  Inf,  Inf])
2336
2337    >>> cur_handler = np.geterrcall()
2338    >>> cur_handler is err_handler
2339    True
2340
2341    """
2342    return umath.geterrobj()[2]
2343
2344class _unspecified(object):
2345    pass
2346_Unspecified = _unspecified()
2347
2348class errstate(object):
2349    """
2350    errstate(**kwargs)
2351
2352    Context manager for floating-point error handling.
2353
2354    Using an instance of `errstate` as a context manager allows statements in
2355    that context to execute with a known error handling behavior. Upon entering
2356    the context the error handling is set with `seterr` and `seterrcall`, and
2357    upon exiting it is reset to what it was before.
2358
2359    Parameters
2360    ----------
2361    kwargs : {divide, over, under, invalid}
2362        Keyword arguments. The valid keywords are the possible floating-point
2363        exceptions. Each keyword should have a string value that defines the
2364        treatment for the particular error. Possible values are
2365        {'ignore', 'warn', 'raise', 'call', 'print', 'log'}.
2366
2367    See Also
2368    --------
2369    seterr, geterr, seterrcall, geterrcall
2370
2371    Notes
2372    -----
2373    The ``with`` statement was introduced in Python 2.5, and can only be used
2374    there by importing it: ``from __future__ import with_statement``. In
2375    earlier Python versions the ``with`` statement is not available.
2376
2377    For complete documentation of the types of floating-point exceptions and
2378    treatment options, see `seterr`.
2379
2380    Examples
2381    --------
2382    >>> from __future__ import with_statement  # use 'with' in Python 2.5
2383    >>> olderr = np.seterr(all='ignore')  # Set error handling to known state.
2384
2385    >>> np.arange(3) / 0.
2386    array([ NaN,  Inf,  Inf])
2387    >>> with np.errstate(divide='warn'):
2388    ...     np.arange(3) / 0.
2389    ...
2390    __main__:2: RuntimeWarning: divide by zero encountered in divide
2391    array([ NaN,  Inf,  Inf])
2392
2393    >>> np.sqrt(-1)
2394    nan
2395    >>> with np.errstate(invalid='raise'):
2396    ...     np.sqrt(-1)
2397    Traceback (most recent call last):
2398      File "<stdin>", line 2, in <module>
2399    FloatingPointError: invalid value encountered in sqrt
2400
2401    Outside the context the error handling behavior has not changed:
2402
2403    >>> np.geterr()
2404    {'over': 'warn', 'divide': 'warn', 'invalid': 'warn',
2405    'under': 'ignore'}
2406
2407    """
2408    # Note that we don't want to run the above doctests because they will fail
2409    # without a from __future__ import with_statement
2410    def __init__(self, **kwargs):
2411        self.call = kwargs.pop('call',_Unspecified)
2412        self.kwargs = kwargs
2413    def __enter__(self):
2414        self.oldstate = seterr(**self.kwargs)
2415        if self.call is not _Unspecified:
2416            self.oldcall = seterrcall(self.call)
2417    def __exit__(self, *exc_info):
2418        seterr(**self.oldstate)
2419        if self.call is not _Unspecified:
2420            seterrcall(self.oldcall)
2421
2422def _setdef():
2423    defval = [UFUNC_BUFSIZE_DEFAULT, ERR_DEFAULT2, None]
2424    umath.seterrobj(defval)
2425
2426# set the default values
2427_setdef()
2428
2429Inf = inf = infty = Infinity = PINF
2430nan = NaN = NAN
2431False_ = bool_(False)
2432True_ = bool_(True)
2433
2434import fromnumeric
2435from fromnumeric import *
2436extend_all(fromnumeric)
Note: See TracBrowser for help on using the repository browser.