diff options
Diffstat (limited to 'numpy')
71 files changed, 1806 insertions, 752 deletions
diff --git a/numpy/_import_tools.py b/numpy/_import_tools.py index 526217359..9b1942e8d 100644 --- a/numpy/_import_tools.py +++ b/numpy/_import_tools.py @@ -2,6 +2,7 @@ from __future__ import division, absolute_import, print_function import os import sys +import warnings __all__ = ['PackageLoader'] @@ -162,7 +163,10 @@ class PackageLoader(object): postpone= : bool when True, don't load packages [default: False] - """ + """ + warnings.warn('pkgload and PackageLoader are obsolete ' + 'and will be removed in a future version of numpy', + DeprecationWarning) frame = self.parent_frame self.info_modules = {} if options.get('force', False): diff --git a/numpy/core/bscript b/numpy/core/bscript index 682f35e8e..5df5a3f8a 100644 --- a/numpy/core/bscript +++ b/numpy/core/bscript @@ -60,27 +60,13 @@ def define_no_smp(): #-------------------------------- # Checking SMP and thread options #-------------------------------- - # Python 2.3 causes a segfault when - # trying to re-acquire the thread-state - # which is done in error-handling - # ufunc code. NPY_ALLOW_C_API and friends - # cause the segfault. So, we disable threading - # for now. - if sys.version[:5] < '2.4.2': - nosmp = 1 - else: - # Perhaps a fancier check is in order here. - # so that threads are only enabled if there - # are actually multiple CPUS? -- but - # threaded code can be nice even on a single - # CPU so that long-calculating code doesn't - # block. - try: - nosmp = os.environ['NPY_NOSMP'] - nosmp = 1 - except KeyError: - nosmp = 0 - return nosmp == 1 + # Perhaps a fancier check is in order here. + # so that threads are only enabled if there + # are actually multiple CPUS? -- but + # threaded code can be nice even on a single + # CPU so that long-calculating code doesn't + # block. + return 'NPY_NOSMP' in os.environ def write_numpy_config(conf): subst_dict = {} diff --git a/numpy/core/code_generators/ufunc_docstrings.py b/numpy/core/code_generators/ufunc_docstrings.py index 67c01cc67..1abbb666b 100644 --- a/numpy/core/code_generators/ufunc_docstrings.py +++ b/numpy/core/code_generators/ufunc_docstrings.py @@ -932,8 +932,8 @@ add_newdoc('numpy.core.umath', 'divide', Returns ------- y : {ndarray, scalar} - The quotient `x1/x2`, element-wise. Returns a scalar if - both `x1` and `x2` are scalars. + The quotient ``x1/x2``, element-wise. Returns a scalar if + both ``x1`` and ``x2`` are scalars. See Also -------- @@ -942,13 +942,13 @@ add_newdoc('numpy.core.umath', 'divide', Notes ----- - Equivalent to `x1` / `x2` in terms of array-broadcasting. + Equivalent to ``x1`` / ``x2`` in terms of array-broadcasting. - Behavior on division by zero can be changed using `seterr`. + Behavior on division by zero can be changed using ``seterr``. - When both `x1` and `x2` are of an integer type, `divide` will return - integers and throw away the fractional part. Moreover, division by zero - always yields zero in integer arithmetic. + In Python 2, when both ``x1`` and ``x2`` are of an integer type, + ``divide`` will behave like ``floor_divide``. In Python 3, it behaves + like ``true_divide``. Examples -------- @@ -961,20 +961,20 @@ add_newdoc('numpy.core.umath', 'divide', [ Inf, 4. , 2.5], [ Inf, 7. , 4. ]]) - Note the behavior with integer types: + Note the behavior with integer types (Python 2 only): >>> np.divide(2, 4) 0 >>> np.divide(2, 4.) 0.5 - Division by zero always yields zero in integer arithmetic, and does not - raise an exception or a warning: + Division by zero always yields zero in integer arithmetic (again, + Python 2 only), and does not raise an exception or a warning: >>> np.divide(np.array([0, 1], dtype=int), np.array([0, 0], dtype=int)) array([0, 0]) - Division by zero can, however, be caught using `seterr`: + Division by zero can, however, be caught using ``seterr``: >>> old_err_state = np.seterr(divide='raise') >>> np.divide(1, 0) diff --git a/numpy/core/fromnumeric.py b/numpy/core/fromnumeric.py index 1cb1ed8bc..84a10bf04 100644 --- a/numpy/core/fromnumeric.py +++ b/numpy/core/fromnumeric.py @@ -513,6 +513,12 @@ def transpose(a, axes=None): See Also -------- rollaxis + argsort + + Notes + ----- + Use `transpose(a, argsort(axes))` to invert the transposition of tensors + when using the `axes` keyword argument. Examples -------- @@ -1082,6 +1088,9 @@ def resize(a, new_shape): Examples -------- >>> a=np.array([[0,1],[2,3]]) + >>> np.resize(a,(2,3)) + array([[0, 1, 2], + [3, 0, 1]]) >>> np.resize(a,(1,4)) array([[0, 1, 2, 3]]) >>> np.resize(a,(2,4)) @@ -2095,8 +2104,14 @@ def amax(a, axis=None, out=None, keepdims=False): ---------- a : array_like Input data. - axis : int, optional - Axis along which to operate. By default, flattened input is used. + axis : None or int or tuple of ints, optional + Axis or axes along which to operate. By default, flattened input is + used. + + .. versionadded: 1.7.0 + + If this is a tuple of ints, the maximum is selected over multiple axes, + instead of a single axis or all the axes as before. out : ndarray, optional Alternative output array in which to place the result. Must be of the same shape and buffer length as the expected output. @@ -2179,8 +2194,14 @@ def amin(a, axis=None, out=None, keepdims=False): ---------- a : array_like Input data. - axis : int, optional - Axis along which to operate. By default, flattened input is used. + axis : None or int or tuple of ints, optional + Axis or axes along which to operate. By default, flattened input is + used. + + .. versionadded: 1.7.0 + + If this is a tuple of ints, the minimum is selected over multiple axes, + instead of a single axis or all the axes as before. out : ndarray, optional Alternative output array in which to place the result. Must be of the same shape and buffer length as the expected output. @@ -2693,9 +2714,14 @@ def mean(a, axis=None, dtype=None, out=None, keepdims=False): a : array_like Array containing numbers whose mean is desired. If `a` is not an array, a conversion is attempted. - axis : int, optional - Axis along which the means are computed. The default is to compute - the mean of the flattened array. + axis : None or int or tuple of ints, optional + Axis or axes along which the means are computed. The default is to + compute the mean of the flattened array. + + .. versionadded: 1.7.0 + + If this is a tuple of ints, a mean is performed over multiple axes, + instead of a single axis or all the axes as before. dtype : data-type, optional Type to use in computing the mean. For integer inputs, the default is `float64`; for floating point inputs, it is the same as the @@ -2778,9 +2804,14 @@ def std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False): ---------- a : array_like Calculate the standard deviation of these values. - axis : int, optional - Axis along which the standard deviation is computed. The default is - to compute the standard deviation of the flattened array. + axis : None or int or tuple of ints, optional + Axis or axes along which the standard deviation is computed. The + default is to compute the standard deviation of the flattened array. + + .. versionadded: 1.7.0 + + If this is a tuple of ints, a standard deviation is performed over + multiple axes, instead of a single axis or all the axes as before. dtype : dtype, optional Type to use in computing the standard deviation. For arrays of integer type the default is float64, for arrays of float types it is @@ -2881,9 +2912,14 @@ def var(a, axis=None, dtype=None, out=None, ddof=0, a : array_like Array containing numbers whose variance is desired. If `a` is not an array, a conversion is attempted. - axis : int, optional - Axis along which the variance is computed. The default is to compute - the variance of the flattened array. + axis : None or int or tuple of ints, optional + Axis or axes along which the variance is computed. The default is to + compute the variance of the flattened array. + + .. versionadded: 1.7.0 + + If this is a tuple of ints, a variance is performed over multiple axes, + instead of a single axis or all the axes as before. dtype : data-type, optional Type to use in computing the variance. For arrays of integer type the default is `float32`; for arrays of float types it is the same as diff --git a/numpy/core/include/numpy/ndarrayobject.h b/numpy/core/include/numpy/ndarrayobject.h index b8c7c3a2d..fbaaeacea 100644 --- a/numpy/core/include/numpy/ndarrayobject.h +++ b/numpy/core/include/numpy/ndarrayobject.h @@ -14,6 +14,7 @@ extern "C" CONFUSE_EMACS everything when you're typing */ #endif +#include <Python.h> #include "ndarraytypes.h" /* Includes the "function" C-API -- these are all stored in a @@ -50,14 +51,26 @@ extern "C" CONFUSE_EMACS #define PyArray_CheckScalar(m) (PyArray_IsScalar(m, Generic) || \ PyArray_IsZeroDim(m)) - +#if PY_MAJOR_VERSION >= 3 +#define PyArray_IsPythonNumber(obj) \ + (PyFloat_Check(obj) || PyComplex_Check(obj) || \ + PyLong_Check(obj) || PyBool_Check(obj)) +#define PyArray_IsIntegerScalar(obj) (PyLong_Check(obj) \ + || PyArray_IsScalar((obj), Integer)) +#define PyArray_IsPythonScalar(obj) \ + (PyArray_IsPythonNumber(obj) || PyBytes_Check(obj) || \ + PyUnicode_Check(obj)) +#else #define PyArray_IsPythonNumber(obj) \ (PyInt_Check(obj) || PyFloat_Check(obj) || PyComplex_Check(obj) || \ PyLong_Check(obj) || PyBool_Check(obj)) - +#define PyArray_IsIntegerScalar(obj) (PyInt_Check(obj) \ + || PyLong_Check(obj) \ + || PyArray_IsScalar((obj), Integer)) #define PyArray_IsPythonScalar(obj) \ (PyArray_IsPythonNumber(obj) || PyString_Check(obj) || \ PyUnicode_Check(obj)) +#endif #define PyArray_IsAnyScalar(obj) \ (PyArray_IsScalar(obj, Generic) || PyArray_IsPythonScalar(obj)) @@ -65,10 +78,6 @@ extern "C" CONFUSE_EMACS #define PyArray_CheckAnyScalar(obj) (PyArray_IsPythonScalar(obj) || \ PyArray_CheckScalar(obj)) -#define PyArray_IsIntegerScalar(obj) (PyInt_Check(obj) \ - || PyLong_Check(obj) \ - || PyArray_IsScalar((obj), Integer)) - #define PyArray_GETCONTIGUOUS(m) (PyArray_ISCONTIGUOUS(m) ? \ Py_INCREF(m), (m) : \ diff --git a/numpy/core/include/numpy/npy_common.h b/numpy/core/include/numpy/npy_common.h index 5cba8c9d2..92b03d20c 100644 --- a/numpy/core/include/numpy/npy_common.h +++ b/numpy/core/include/numpy/npy_common.h @@ -133,6 +133,7 @@ extern long long __cdecl _ftelli64(FILE *); #else #define npy_ftell ftell #endif + #include <sys/types.h> #define npy_lseek lseek #define npy_off_t off_t @@ -264,18 +265,9 @@ typedef unsigned PY_LONG_LONG npy_ulonglong; # ifdef _MSC_VER # define NPY_LONGLONG_FMT "I64d" # define NPY_ULONGLONG_FMT "I64u" -# elif defined(__APPLE__) || defined(__FreeBSD__) -/* "%Ld" only parses 4 bytes -- "L" is floating modifier on MacOS X/BSD */ +# else # define NPY_LONGLONG_FMT "lld" # define NPY_ULONGLONG_FMT "llu" -/* - another possible variant -- *quad_t works on *BSD, but is deprecated: - #define LONGLONG_FMT "qd" - #define ULONGLONG_FMT "qu" -*/ -# else -# define NPY_LONGLONG_FMT "Ld" -# define NPY_ULONGLONG_FMT "Lu" # endif # ifdef _MSC_VER # define NPY_LONGLONG_SUFFIX(x) (x##i64) diff --git a/numpy/core/include/numpy/npy_cpu.h b/numpy/core/include/numpy/npy_cpu.h index 24d4ce1fc..60abae4e0 100644 --- a/numpy/core/include/numpy/npy_cpu.h +++ b/numpy/core/include/numpy/npy_cpu.h @@ -20,6 +20,7 @@ #define _NPY_CPUARCH_H_ #include "numpyconfig.h" +#include <string.h> /* for memcpy */ #if defined( __i386__ ) || defined(i386) || defined(_M_IX86) /* @@ -80,38 +81,7 @@ information about your platform (OS, CPU and compiler) #endif -/* - This "white-lists" the architectures that we know don't require - pointer alignment. We white-list, since the memcpy version will - work everywhere, whereas assignment will only work where pointer - dereferencing doesn't require alignment. - - TODO: There may be more architectures we can white list. -*/ -#if defined(NPY_CPU_X86) || defined(NPY_CPU_AMD64) - #define NPY_COPY_PYOBJECT_PTR(dst, src) (*((PyObject **)(dst)) = *((PyObject **)(src))) -#else - #if NPY_SIZEOF_PY_INTPTR_T == 4 - #define NPY_COPY_PYOBJECT_PTR(dst, src) \ - ((char*)(dst))[0] = ((char*)(src))[0]; \ - ((char*)(dst))[1] = ((char*)(src))[1]; \ - ((char*)(dst))[2] = ((char*)(src))[2]; \ - ((char*)(dst))[3] = ((char*)(src))[3]; - #elif NPY_SIZEOF_PY_INTPTR_T == 8 - #define NPY_COPY_PYOBJECT_PTR(dst, src) \ - ((char*)(dst))[0] = ((char*)(src))[0]; \ - ((char*)(dst))[1] = ((char*)(src))[1]; \ - ((char*)(dst))[2] = ((char*)(src))[2]; \ - ((char*)(dst))[3] = ((char*)(src))[3]; \ - ((char*)(dst))[4] = ((char*)(src))[4]; \ - ((char*)(dst))[5] = ((char*)(src))[5]; \ - ((char*)(dst))[6] = ((char*)(src))[6]; \ - ((char*)(dst))[7] = ((char*)(src))[7]; - #else - #error Unknown architecture, please report this to numpy maintainers with \ - information about your platform (OS, CPU and compiler) - #endif -#endif +#define NPY_COPY_PYOBJECT_PTR(dst, src) memcpy(dst, src, sizeof(PyObject *)) #if (defined(NPY_CPU_X86) || defined(NPY_CPU_AMD64)) #define NPY_CPU_HAVE_UNALIGNED_ACCESS 1 diff --git a/numpy/core/memmap.py b/numpy/core/memmap.py index b1c96ee29..4b10f361c 100644 --- a/numpy/core/memmap.py +++ b/numpy/core/memmap.py @@ -111,6 +111,11 @@ class memmap(ndarray): certain size depending on the platform. This size is always < 2GB even on 64-bit systems. + When a memmap causes a file to be created or extended beyond its + current size in the filesystem, the contents of the new part are + unspecified. On systems with POSIX filesystem semantics, the extended + part will be filled with zero bytes. + Examples -------- >>> data = np.arange(12, dtype='float32') diff --git a/numpy/core/numerictypes.py b/numpy/core/numerictypes.py index 1545bc734..0c03cce89 100644 --- a/numpy/core/numerictypes.py +++ b/numpy/core/numerictypes.py @@ -670,7 +670,7 @@ def issubclass_(arg1, arg2): Determine if a class is a subclass of a second class. `issubclass_` is equivalent to the Python built-in ``issubclass``, - except that it returns False instead of raising a TypeError is one + except that it returns False instead of raising a TypeError if one of the arguments is not a class. Parameters diff --git a/numpy/core/records.py b/numpy/core/records.py index d0f82a25c..bf4d835ea 100644 --- a/numpy/core/records.py +++ b/numpy/core/records.py @@ -71,7 +71,6 @@ _byteorderconv = {'b':'>', # are equally allowed numfmt = nt.typeDict -_typestr = nt._typestr def find_duplicate(list): """Find duplication in a list, return a list of duplicated elements""" @@ -268,7 +267,7 @@ class record(nt.void): """Pretty-print all fields.""" # pretty-print all fields names = self.dtype.names - maxlen = max([len(name) for name in names]) + maxlen = max(len(name) for name in names) rows = [] fmt = '%% %ds: %%s' % maxlen for name in names: @@ -527,15 +526,12 @@ def fromarrays(arrayList, dtype=None, shape=None, formats=None, if formats is None and dtype is None: # go through each object in the list to see if it is an ndarray # and determine the formats. - formats = '' + formats = [] for obj in arrayList: if not isinstance(obj, ndarray): raise ValueError("item in the array list must be an ndarray.") - formats += _typestr[obj.dtype.type] - if issubclass(obj.dtype.type, nt.flexible): - formats += repr(obj.itemsize) - formats += ',' - formats = formats[:-1] + formats.append(obj.dtype.str) + formats = ','.join(formats) if dtype is not None: descr = sb.dtype(dtype) diff --git a/numpy/core/setup.py b/numpy/core/setup.py index 9cb9d7361..a51eb690b 100644 --- a/numpy/core/setup.py +++ b/numpy/core/setup.py @@ -78,27 +78,13 @@ def is_npy_no_signal(): def is_npy_no_smp(): """Return True if the NPY_NO_SMP symbol must be defined in public header (when SMP support cannot be reliably enabled).""" - # Python 2.3 causes a segfault when - # trying to re-acquire the thread-state - # which is done in error-handling - # ufunc code. NPY_ALLOW_C_API and friends - # cause the segfault. So, we disable threading - # for now. - if sys.version[:5] < '2.4.2': - nosmp = 1 - else: - # Perhaps a fancier check is in order here. - # so that threads are only enabled if there - # are actually multiple CPUS? -- but - # threaded code can be nice even on a single - # CPU so that long-calculating code doesn't - # block. - try: - nosmp = os.environ['NPY_NOSMP'] - nosmp = 1 - except KeyError: - nosmp = 0 - return nosmp == 1 + # Perhaps a fancier check is in order here. + # so that threads are only enabled if there + # are actually multiple CPUS? -- but + # threaded code can be nice even on a single + # CPU so that long-calculating code doesn't + # block. + return 'NPY_NOSMP' in os.environ def win32_checks(deflist): from numpy.distutils.misc_util import get_build_architecture @@ -279,11 +265,11 @@ def check_types(config_cmd, ext, build_dir): expected['long'] = [8, 4] expected['float'] = [4] expected['double'] = [8] - expected['long double'] = [8, 12, 16] - expected['Py_intptr_t'] = [4, 8] + expected['long double'] = [16, 12, 8] + expected['Py_intptr_t'] = [8, 4] expected['PY_LONG_LONG'] = [8] expected['long long'] = [8] - expected['off_t'] = [4, 8] + expected['off_t'] = [8, 4] # Check we have the python header (-dev* packages on Linux) result = config_cmd.check_header('Python.h') @@ -323,7 +309,8 @@ def check_types(config_cmd, ext, build_dir): # definition is binary compatible with C99 complex type (check done at # build time in npy_common.h) complex_def = "struct {%s __x; %s __y;}" % (type, type) - res = config_cmd.check_type_size(complex_def, expected=2*expected[type]) + res = config_cmd.check_type_size(complex_def, + expected=[2 * x for x in expected[type]]) if res >= 0: public_defines.append(('NPY_SIZEOF_COMPLEX_%s' % sym2def(type), '%d' % res)) else: diff --git a/numpy/core/setup_common.py b/numpy/core/setup_common.py index 3dc1cecf5..e51797c03 100644 --- a/numpy/core/setup_common.py +++ b/numpy/core/setup_common.py @@ -172,10 +172,20 @@ def check_long_double_representation(cmd): body = LONG_DOUBLE_REPRESENTATION_SRC % {'type': 'long double'} # We need to use _compile because we need the object filename - src, object = cmd._compile(body, None, None, 'c') + src, obj = cmd._compile(body, None, None, 'c') try: - type = long_double_representation(pyod(object)) - return type + ltype = long_double_representation(pyod(obj)) + return ltype + except ValueError: + # try linking to support CC="gcc -flto" or icc -ipo + # struct needs to be volatile so it isn't optimized away + body = body.replace('struct', 'volatile struct') + body += "int main(void) { return 0; }\n" + src, obj = cmd._compile(body, None, None, 'c') + cmd.temp_files.append("_configtest") + cmd.compiler.link_executable([obj], "_configtest") + ltype = long_double_representation(pyod("_configtest")) + return ltype finally: cmd._clean() diff --git a/numpy/core/src/multiarray/arraytypes.c.src b/numpy/core/src/multiarray/arraytypes.c.src index 1b9bc9a8c..89b5404b4 100644 --- a/numpy/core/src/multiarray/arraytypes.c.src +++ b/numpy/core/src/multiarray/arraytypes.c.src @@ -4259,8 +4259,12 @@ datetime_dtype_metadata_clone(NpyAuxData *data) /* * Initializes the c_metadata field for the _builtin_descrs DATETIME * and TIMEDELTA. + * + * must not be static, gcc 4.1.2 on redhat 5 then miscompiles this function + * see gh-5163 + * */ -static int +NPY_NO_EXPORT int initialize_builtin_datetime_metadata(void) { PyArray_DatetimeDTypeMetaData *data1, *data2; diff --git a/numpy/core/src/multiarray/buffer.c b/numpy/core/src/multiarray/buffer.c index 6ba98d1cd..7f7607e1f 100644 --- a/numpy/core/src/multiarray/buffer.c +++ b/numpy/core/src/multiarray/buffer.c @@ -629,6 +629,8 @@ array_getbuffer(PyObject *obj, Py_buffer *view, int flags) { PyArrayObject *self; _buffer_info_t *info = NULL; + int i; + Py_ssize_t sd; self = (PyArrayObject*)obj; @@ -703,6 +705,31 @@ array_getbuffer(PyObject *obj, Py_buffer *view, int flags) } if ((flags & PyBUF_STRIDES) == PyBUF_STRIDES) { view->strides = info->strides; + +#ifdef NPY_RELAXED_STRIDES_CHECKING + /* + * If NPY_RELAXED_STRIDES_CHECKING is on, the array may be + * contiguous, but it won't look that way to Python when it + * tries to determine contiguity by looking at the strides + * (since one of the elements may be -1). In that case, just + * regenerate strides from shape. + */ + if (PyArray_CHKFLAGS(self, NPY_ARRAY_C_CONTIGUOUS) && + !((flags & PyBUF_F_CONTIGUOUS) == PyBUF_F_CONTIGUOUS)) { + sd = view->itemsize; + for (i = view->ndim-1; i >= 0; --i) { + view->strides[i] = sd; + sd *= view->shape[i]; + } + } + else if (PyArray_CHKFLAGS(self, NPY_ARRAY_F_CONTIGUOUS)) { + sd = view->itemsize; + for (i = 0; i < view->ndim; ++i) { + view->strides[i] = sd; + sd *= view->shape[i]; + } + } +#endif } else { view->strides = NULL; @@ -922,7 +949,7 @@ _descriptor_from_pep3118_format_fast(char *s, PyObject **result) *result = (PyObject*)PyArray_DescrNewByteorder(descr, byte_order); Py_DECREF(descr); } - + return 1; } diff --git a/numpy/core/src/multiarray/datetime.c b/numpy/core/src/multiarray/datetime.c index a5be9f670..a870650fc 100644 --- a/numpy/core/src/multiarray/datetime.c +++ b/numpy/core/src/multiarray/datetime.c @@ -3047,7 +3047,7 @@ cast_timedelta_to_timedelta(PyArray_DatetimeMetaData *src_meta, * Returns true if the object is something that is best considered * a Datetime, false otherwise. */ -static npy_bool +static NPY_GCC_NONNULL(1) npy_bool is_any_numpy_datetime(PyObject *obj) { return (PyArray_IsScalar(obj, Datetime) || @@ -3296,7 +3296,8 @@ datetime_arange(PyObject *start, PyObject *stop, PyObject *step, } } else { - if (is_any_numpy_datetime(start) || is_any_numpy_datetime(stop)) { + if ((start && is_any_numpy_datetime(start)) || + is_any_numpy_datetime(stop)) { type_nums[0] = NPY_DATETIME; } else { diff --git a/numpy/core/src/multiarray/item_selection.c b/numpy/core/src/multiarray/item_selection.c index b2bf17f4c..cd0ae1680 100644 --- a/numpy/core/src/multiarray/item_selection.c +++ b/numpy/core/src/multiarray/item_selection.c @@ -265,6 +265,11 @@ PyArray_PutTo(PyArrayObject *self, PyObject* values0, PyObject *indices0, "put: first argument must be an array"); return NULL; } + + if (PyArray_FailUnlessWriteable(self, "put: output array") < 0) { + return NULL; + } + if (!PyArray_ISCONTIGUOUS(self)) { PyArrayObject *obj; int flags = NPY_ARRAY_CARRAY | NPY_ARRAY_UPDATEIFCOPY; diff --git a/numpy/core/src/multiarray/multiarray_tests.c.src b/numpy/core/src/multiarray/multiarray_tests.c.src index a22319cfe..77e699562 100644 --- a/numpy/core/src/multiarray/multiarray_tests.c.src +++ b/numpy/core/src/multiarray/multiarray_tests.c.src @@ -2,6 +2,22 @@ #include <Python.h> #include "numpy/arrayobject.h" +/* test PyArray_IsPythonScalar, before including private py3 compat header */ +static PyObject * +IsPythonScalar(PyObject * dummy, PyObject *args) +{ + PyObject *arg = NULL; + if (!PyArg_ParseTuple(args, "O", &arg)) { + return NULL; + } + if (PyArray_IsPythonScalar(arg)) { + Py_RETURN_TRUE; + } + else { + Py_RETURN_FALSE; + } +} + #include "npy_pycompat.h" /* @@ -860,6 +876,9 @@ test_nditer_too_large(PyObject *NPY_UNUSED(self), PyObject *args) { static PyMethodDef Multiarray_TestsMethods[] = { + {"IsPythonScalar", + IsPythonScalar, + METH_VARARGS, NULL}, {"test_neighborhood_iterator", test_neighborhood_iterator, METH_VARARGS, NULL}, diff --git a/numpy/core/src/multiarray/multiarraymodule.c b/numpy/core/src/multiarray/multiarraymodule.c index 51c594781..1bf0613e4 100644 --- a/numpy/core/src/multiarray/multiarraymodule.c +++ b/numpy/core/src/multiarray/multiarraymodule.c @@ -253,7 +253,7 @@ PyArray_As1D(PyObject **op, char **ptr, int *d1, int typecode) { npy_intp newd1; PyArray_Descr *descr; - char msg[] = "PyArray_As1D: use PyArray_AsCArray."; + static const char msg[] = "PyArray_As1D: use PyArray_AsCArray."; if (DEPRECATE(msg) < 0) { return -1; @@ -274,7 +274,7 @@ PyArray_As2D(PyObject **op, char ***ptr, int *d1, int *d2, int typecode) { npy_intp newdims[2]; PyArray_Descr *descr; - char msg[] = "PyArray_As1D: use PyArray_AsCArray."; + static const char msg[] = "PyArray_As1D: use PyArray_AsCArray."; if (DEPRECATE(msg) < 0) { return -1; @@ -346,8 +346,9 @@ PyArray_ConcatenateArrays(int narrays, PyArrayObject **arrays, int axis) } if (ndim == 1 && axis != 0) { - char msg[] = "axis != 0 for ndim == 1; this will raise an error in " - "future versions of numpy"; + static const char msg[] = "axis != 0 for ndim == 1; " + "this will raise an error in " + "future versions of numpy"; if (DEPRECATE(msg) < 0) { return NULL; } diff --git a/numpy/core/src/npymath/npy_math.c.src b/numpy/core/src/npymath/npy_math.c.src index ba5c58fb4..b7f28bb39 100644 --- a/numpy/core/src/npymath/npy_math.c.src +++ b/numpy/core/src/npymath/npy_math.c.src @@ -491,7 +491,11 @@ double npy_log2(double x) #ifndef HAVE_CBRT@C@ @type@ npy_cbrt@c@(@type@ x) { - if (x < 0) { + /* don't set invalid flag */ + if (npy_isnan(x)) { + return NPY_NAN; + } + else if (x < 0) { return -npy_pow@c@(-x, 1. / 3.); } else { diff --git a/numpy/core/src/npymath/npy_math_complex.c.src b/numpy/core/src/npymath/npy_math_complex.c.src index 920f107b8..9cbfd32ae 100644 --- a/numpy/core/src/npymath/npy_math_complex.c.src +++ b/numpy/core/src/npymath/npy_math_complex.c.src @@ -247,7 +247,8 @@ #ifdef HAVE_@KIND@@C@ @type@ npy_@kind@@c@(@ctype@ z) { - __@ctype@_to_c99_cast z1 = {z}; + __@ctype@_to_c99_cast z1; + z1.npy_z = z; return @kind@@c@(z1.c99_z); } #endif @@ -260,8 +261,9 @@ #ifdef HAVE_@KIND@@C@ @ctype@ npy_@kind@@c@(@ctype@ z) { - __@ctype@_to_c99_cast z1 = {z}; + __@ctype@_to_c99_cast z1; __@ctype@_to_c99_cast ret; + z1.npy_z = z; ret.c99_z = @kind@@c@(z1.c99_z); return ret.npy_z; } @@ -275,9 +277,11 @@ #ifdef HAVE_@KIND@@C@ @ctype@ npy_@kind@@c@(@ctype@ x, @ctype@ y) { - __@ctype@_to_c99_cast xcast = {x}; - __@ctype@_to_c99_cast ycast = {y}; + __@ctype@_to_c99_cast xcast; + __@ctype@_to_c99_cast ycast; __@ctype@_to_c99_cast ret; + xcast.npy_z = x; + ycast.npy_z = y; ret.c99_z = @kind@@c@(xcast.c99_z, ycast.c99_z); return ret.npy_z; } diff --git a/numpy/core/src/npymath/npy_math_private.h b/numpy/core/src/npymath/npy_math_private.h index b3b1690be..284d203bf 100644 --- a/numpy/core/src/npymath/npy_math_private.h +++ b/numpy/core/src/npymath/npy_math_private.h @@ -485,6 +485,24 @@ do { \ * support is available */ #ifdef NPY_USE_C99_COMPLEX + +/* Microsoft C defines _MSC_VER */ +#ifdef _MSC_VER +typedef union { + npy_cdouble npy_z; + _Dcomplex c99_z; +} __npy_cdouble_to_c99_cast; + +typedef union { + npy_cfloat npy_z; + _Fcomplex c99_z; +} __npy_cfloat_to_c99_cast; + +typedef union { + npy_clongdouble npy_z; + _Lcomplex c99_z; +} __npy_clongdouble_to_c99_cast; +#else /* !_MSC_VER */ typedef union { npy_cdouble npy_z; complex double c99_z; @@ -499,7 +517,9 @@ typedef union { npy_clongdouble npy_z; complex long double c99_z; } __npy_clongdouble_to_c99_cast; -#else +#endif /* !_MSC_VER */ + +#else /* !NPY_USE_C99_COMPLEX */ typedef union { npy_cdouble npy_z; npy_cdouble c99_z; @@ -514,6 +534,6 @@ typedef union { npy_clongdouble npy_z; npy_clongdouble c99_z; } __npy_clongdouble_to_c99_cast; -#endif +#endif /* !NPY_USE_C99_COMPLEX */ #endif /* !_NPY_MATH_PRIVATE_H_ */ diff --git a/numpy/core/src/private/npy_config.h b/numpy/core/src/private/npy_config.h index 71d448ee9..882913e2f 100644 --- a/numpy/core/src/private/npy_config.h +++ b/numpy/core/src/private/npy_config.h @@ -21,16 +21,6 @@ */ #define NPY_MAX_COPY_ALIGNMENT 16 -/* Safe to use ldexp and frexp for long double for MSVC builds */ -#if (NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE) || defined(_MSC_VER) - #ifdef HAVE_LDEXP - #define HAVE_LDEXPL 1 - #endif - #ifdef HAVE_FREXP - #define HAVE_FREXPL 1 - #endif -#endif - /* Disable broken Sun Workshop Pro math functions */ #ifdef __SUNPRO_C #undef HAVE_ATAN2 diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src index 14076f3fa..a69fc8147 100644 --- a/numpy/core/src/umath/loops.c.src +++ b/numpy/core/src/umath/loops.c.src @@ -47,20 +47,23 @@ ***************************************************************************** */ +/* unary loop input and output contiguous */ +#define IS_UNARY_CONT(tin, tout) (steps[0] == sizeof(tin) && \ + steps[1] == sizeof(tout)) #define IS_BINARY_REDUCE ((args[0] == args[2])\ && (steps[0] == steps[2])\ && (steps[0] == 0)) -/* binary loop input and output continous */ +/* binary loop input and output contiguous */ #define IS_BINARY_CONT(tin, tout) (steps[0] == sizeof(tin) && \ steps[1] == sizeof(tin) && \ steps[2] == sizeof(tout)) -/* binary loop input and output continous with first scalar */ +/* binary loop input and output contiguous with first scalar */ #define IS_BINARY_CONT_S1(tin, tout) (steps[0] == 0 && \ steps[1] == sizeof(tin) && \ steps[2] == sizeof(tout)) -/* binary loop input and output continous with second scalar */ +/* binary loop input and output contiguous with second scalar */ #define IS_BINARY_CONT_S2(tin, tout) (steps[0] == sizeof(tin) && \ steps[1] == 0 && \ steps[2] == sizeof(tout)) @@ -79,6 +82,33 @@ npy_intp i;\ for(i = 0; i < n; i++, ip1 += is1, op1 += os1) +/* + * loop with contiguous specialization + * op should be the code working on `tin in` and + * storing the result in `tout * out` + * combine with NPY_GCC_OPT_3 to allow autovectorization + * should only be used where its worthwhile to avoid code bloat + */ +#define UNARY_LOOP_FAST(tin, tout, op) \ + do { \ + /* condition allows compiler to optimize the generic macro */ \ + if (IS_UNARY_CONT(tin, tout)) { \ + UNARY_LOOP { \ + const tin in = *(tin *)ip1; \ + tout * out = (tout *)op1; \ + op; \ + } \ + } \ + else { \ + UNARY_LOOP { \ + const tin in = *(tin *)ip1; \ + tout * out = (tout *)op1; \ + op; \ + } \ + } \ + } \ + while (0) + #define UNARY_LOOP_TWO_OUT\ char *ip1 = args[0], *op1 = args[1], *op2 = args[2];\ npy_intp is1 = steps[0], os1 = steps[1], os2 = steps[2];\ @@ -93,6 +123,51 @@ npy_intp i;\ for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1) +/* + * loop with contiguous specialization + * op should be the code working on `tin in1`, `tin in2` and + * storing the result in `tout * out` + * combine with NPY_GCC_OPT_3 to allow autovectorization + * should only be used where its worthwhile to avoid code bloat + */ +#define BINARY_LOOP_FAST(tin, tout, op) \ + do { \ + /* condition allows compiler to optimize the generic macro */ \ + if (IS_BINARY_CONT(tin, tout)) { \ + BINARY_LOOP { \ + const tin in1 = *(tin *)ip1; \ + const tin in2 = *(tin *)ip2; \ + tout * out = (tout *)op1; \ + op; \ + } \ + } \ + else if (IS_BINARY_CONT_S1(tin, tout)) { \ + const tin in1 = *(tin *)args[0]; \ + BINARY_LOOP { \ + const tin in2 = *(tin *)ip2; \ + tout * out = (tout *)op1; \ + op; \ + } \ + } \ + else if (IS_BINARY_CONT_S2(tin, tout)) { \ + const tin in2 = *(tin *)args[1]; \ + BINARY_LOOP { \ + const tin in1 = *(tin *)ip1; \ + tout * out = (tout *)op1; \ + op; \ + } \ + } \ + else { \ + BINARY_LOOP { \ + const tin in1 = *(tin *)ip1; \ + const tin in2 = *(tin *)ip2; \ + tout * out = (tout *)op1; \ + op; \ + } \ + } \ + } \ + while (0) + #define BINARY_REDUCE_LOOP_INNER\ char *ip2 = args[1]; \ npy_intp is2 = steps[1]; \ @@ -703,58 +778,40 @@ NPY_NO_EXPORT void } } -NPY_NO_EXPORT void +NPY_NO_EXPORT NPY_GCC_OPT_3 void @TYPE@_square(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(data)) { - UNARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - *((@type@ *)op1) = in1*in1; - } + UNARY_LOOP_FAST(@type@, @type@, *out = in * in); } NPY_NO_EXPORT void @TYPE@_reciprocal(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(data)) { - UNARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - *((@type@ *)op1) = (@type@)(1.0/in1); - } + UNARY_LOOP_FAST(@type@, @type@, *out = 1.0 / in); } NPY_NO_EXPORT void @TYPE@_conjugate(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { - UNARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - *((@type@ *)op1) = in1; - } + UNARY_LOOP_FAST(@type@, @type@, *out = in); } -NPY_NO_EXPORT void +NPY_NO_EXPORT NPY_GCC_OPT_3 void @TYPE@_negative(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { - UNARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - *((@type@ *)op1) = (@type@)(-(@type@)in1); - } + UNARY_LOOP_FAST(@type@, @type@, *out = -in); } -NPY_NO_EXPORT void +NPY_NO_EXPORT NPY_GCC_OPT_3 void @TYPE@_logical_not(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { - UNARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - *((npy_bool *)op1) = !in1; - } + UNARY_LOOP_FAST(@type@, npy_bool, *out = !in); } -NPY_NO_EXPORT void +NPY_NO_EXPORT NPY_GCC_OPT_3 void @TYPE@_invert(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { - UNARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - *((@type@ *)op1) = ~in1; - } + UNARY_LOOP_FAST(@type@, @type@, *out = ~in); } /**begin repeat1 @@ -764,7 +821,7 @@ NPY_NO_EXPORT void * #OP = +, -,*, &, |, ^, <<, >># */ -NPY_NO_EXPORT void +NPY_NO_EXPORT NPY_GCC_OPT_3 void @TYPE@_@kind@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { if(IS_BINARY_REDUCE) { @@ -774,11 +831,7 @@ NPY_NO_EXPORT void *((@type@ *)iop1) = io1; } else { - BINARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - const @type@ in2 = *(@type@ *)ip2; - *((@type@ *)op1) = in1 @OP@ in2; - } + BINARY_LOOP_FAST(@type@, @type@, *out = in1 @OP@ in2); } } @@ -797,39 +850,7 @@ NPY_NO_EXPORT NPY_GCC_OPT_3 void * gcc vectorization of this is not good (PR60575) but manual integer * vectorization is too tedious to be worthwhile */ - if (IS_BINARY_CONT(@type@, npy_bool)) { - npy_intp i, n = dimensions[0]; - @type@ * a = (@type@ *)args[0], * b = (@type@ *)args[1]; - npy_bool * o = (npy_bool *)args[2]; - for (i = 0; i < n; i++) { - o[i] = a[i] @OP@ b[i]; - } - } - else if (IS_BINARY_CONT_S1(@type@, npy_bool)) { - npy_intp i, n = dimensions[0]; - @type@ a = *(@type@ *)args[0]; - @type@ * b = (@type@ *)args[1]; - npy_bool * o = (npy_bool *)args[2]; - for (i = 0; i < n; i++) { - o[i] = a @OP@ b[i]; - } - } - else if (IS_BINARY_CONT_S2(@type@, npy_bool)) { - npy_intp i, n = dimensions[0]; - @type@ * a = (@type@ *)args[0]; - @type@ b = *(@type@*)args[1]; - npy_bool * o = (npy_bool *)args[2]; - for (i = 0; i < n; i++) { - o[i] = a[i] @OP@ b; - } - } - else { - BINARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - const @type@ in2 = *(@type@ *)ip2; - *((npy_bool *)op1) = in1 @OP@ in2; - } - } + BINARY_LOOP_FAST(@type@, npy_bool, *out = in1 @OP@ in2); } /**end repeat1**/ @@ -914,22 +935,16 @@ NPY_NO_EXPORT void * #type = npy_byte, npy_short, npy_int, npy_long, npy_longlong# */ -NPY_NO_EXPORT void +NPY_NO_EXPORT NPY_GCC_OPT_3 void @TYPE@_absolute(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { - UNARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - *((@type@ *)op1) = (in1 >= 0) ? in1 : -in1; - } + UNARY_LOOP_FAST(@type@, @type@, *out = (in >= 0) ? in : -in); } -NPY_NO_EXPORT void +NPY_NO_EXPORT NPY_GCC_OPT_3 void @TYPE@_sign(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { - UNARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - *((@type@ *)op1) = in1 > 0 ? 1 : (in1 < 0 ? -1 : 0); - } + UNARY_LOOP_FAST(@type@, @type@, *out = in > 0 ? 1 : (in < 0 ? -1 : 0)); } NPY_NO_EXPORT void @@ -997,13 +1012,10 @@ NPY_NO_EXPORT void } } -NPY_NO_EXPORT void +NPY_NO_EXPORT NPY_GCC_OPT_3 void @TYPE@_sign(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { - UNARY_LOOP { - const @type@ in1 = *(@type@ *)ip1; - *((@type@ *)op1) = in1 > 0 ? 1 : 0; - } + UNARY_LOOP_FAST(@type@, @type@, *out = in > 0 ? 1 : 0); } NPY_NO_EXPORT void diff --git a/numpy/core/src/umath/ufunc_object.c b/numpy/core/src/umath/ufunc_object.c index b12fbf57f..dc5065f14 100644 --- a/numpy/core/src/umath/ufunc_object.c +++ b/numpy/core/src/umath/ufunc_object.c @@ -1901,60 +1901,67 @@ PyUFunc_GeneralizedFunction(PyUFuncObject *ufunc, } /* - * Validate the core dimensions of all the operands, - * and collect all of the labeled core dimension sizes - * into the array 'core_dim_sizes'. Initialize them to - * 1, for example in the case where the operand broadcasts - * to a core dimension, it won't be visited. + * Validate the core dimensions of all the operands, and collect all of + * the labelled core dimensions into 'core_dim_sizes'. + * + * The behavior has been changed in Numpy 1.10.0, and the following + * requirements must be fulfilled or an error will be raised: + * * Arguments, both input and output, must have at least as many + * dimensions as the corresponding number of core dimensions. In + * previous versions, 1's were prepended to the shape as needed. + * * Core dimensions with same labels must have exactly matching sizes. + * In previous versions, core dimensions of size 1 would broadcast + * against other core dimensions with the same label. + * * All core dimensions must have their size specified by a passed in + * input or output argument. In previous versions, core dimensions in + * an output argument that were not specified in an input argument, + * and whose size could not be inferred from a passed in output + * argument, would have their size set to 1. */ for (i = 0; i < ufunc->core_num_dim_ix; ++i) { - core_dim_sizes[i] = 1; + core_dim_sizes[i] = -1; } for (i = 0; i < nop; ++i) { if (op[i] != NULL) { int dim_offset = ufunc->core_offsets[i]; int num_dims = ufunc->core_num_dims[i]; int core_start_dim = PyArray_NDIM(op[i]) - num_dims; - /* Make sure any output operand has enough dimensions */ - if (i >= nin && core_start_dim < 0) { + + /* Check if operands have enough dimensions */ + if (core_start_dim < 0) { PyErr_Format(PyExc_ValueError, - "%s: Output operand %d does not have enough dimensions " - "(has %d, gufunc core with signature %s " - "requires %d)", - ufunc_name, i - nin, PyArray_NDIM(op[i]), + "%s: %s operand %d does not have enough " + "dimensions (has %d, gufunc core with " + "signature %s requires %d)", + ufunc_name, i < nin ? "Input" : "Output", + i < nin ? i : i - nin, PyArray_NDIM(op[i]), ufunc->core_signature, num_dims); retval = -1; goto fail; } /* - * Make sure each core dimension matches all other core - * dimensions with the same label - * - * NOTE: For input operands, core_start_dim may be negative. - * In that case, the operand is being broadcast onto - * core dimensions. For example, a scalar will broadcast - * to fit any core signature. + * Make sure every core dimension exactly matches all other core + * dimensions with the same label. */ - if (core_start_dim >= 0) { - idim = 0; - } else { - idim = -core_start_dim; - } - for (; idim < num_dims; ++idim) { - int core_dim_index = ufunc->core_dim_ixs[dim_offset + idim]; + for (idim = 0; idim < num_dims; ++idim) { + int core_dim_index = ufunc->core_dim_ixs[dim_offset+idim]; npy_intp op_dim_size = - PyArray_SHAPE(op[i])[core_start_dim + idim]; - if (core_dim_sizes[core_dim_index] == 1) { + PyArray_DIM(op[i], core_start_dim+idim); + + if (core_dim_sizes[core_dim_index] == -1) { core_dim_sizes[core_dim_index] = op_dim_size; - } else if ((i >= nin || op_dim_size != 1) && - core_dim_sizes[core_dim_index] != op_dim_size) { + } + else if (op_dim_size != core_dim_sizes[core_dim_index]) { PyErr_Format(PyExc_ValueError, - "%s: Operand %d has a mismatch in its core " - "dimension %d, with gufunc signature %s " - "(size %zd is different from %zd)", - ufunc_name, i, idim, ufunc->core_signature, - op_dim_size, core_dim_sizes[core_dim_index]); + "%s: %s operand %d has a mismatch in its " + "core dimension %d, with gufunc " + "signature %s (size %zd is different " + "from %zd)", + ufunc_name, i < nin ? "Input" : "Output", + i < nin ? i : i - nin, idim, + ufunc->core_signature, op_dim_size, + core_dim_sizes[core_dim_index]); retval = -1; goto fail; } @@ -1962,6 +1969,44 @@ PyUFunc_GeneralizedFunction(PyUFuncObject *ufunc, } } + /* + * Make sure no core dimension is unspecified. + */ + for (i = 0; i < ufunc->core_num_dim_ix; ++i) { + if (core_dim_sizes[i] == -1) { + break; + } + } + if (i != ufunc->core_num_dim_ix) { + /* + * There is at least one core dimension missing, find in which + * operand it comes up first (it has to be an output operand). + */ + const int missing_core_dim = i; + int out_op; + for (out_op = nin; out_op < nop; ++out_op) { + int first_idx = ufunc->core_offsets[out_op]; + int last_idx = first_idx + ufunc->core_num_dims[out_op]; + for (i = first_idx; i < last_idx; ++i) { + if (ufunc->core_dim_ixs[i] == missing_core_dim) { + break; + } + } + if (i < last_idx) { + /* Change index offsets for error message */ + out_op -= nin; + i -= first_idx; + break; + } + } + PyErr_Format(PyExc_ValueError, + "%s: Output operand %d has core dimension %d " + "unspecified, with gufunc signature %s", + ufunc_name, out_op, i, ufunc->core_signature); + retval = -1; + goto fail; + } + /* Fill in the initial part of 'iter_shape' */ for (idim = 0; idim < broadcast_ndim; ++idim) { iter_shape[idim] = -1; @@ -3926,18 +3971,19 @@ _find_array_wrap(PyObject *args, PyObject *kwds, PyObject *with_wrap[NPY_MAXARGS], *wraps[NPY_MAXARGS]; PyObject *obj, *wrap = NULL; - /* If a 'subok' parameter is passed and isn't True, don't wrap */ + /* + * If a 'subok' parameter is passed and isn't True, don't wrap but put None + * into slots with out arguments which means return the out argument + */ if (kwds != NULL && (obj = PyDict_GetItem(kwds, npy_um_str_subok)) != NULL) { if (obj != Py_True) { - for (i = 0; i < nout; i++) { - output_wrap[i] = NULL; - } - return; + /* skip search for wrap members */ + goto handle_out; } } - nargs = PyTuple_GET_SIZE(args); + for (i = 0; i < nin; i++) { obj = PyTuple_GET_ITEM(args, i); if (PyArray_CheckExact(obj) || PyArray_IsAnyScalar(obj)) { @@ -3995,6 +4041,8 @@ _find_array_wrap(PyObject *args, PyObject *kwds, * exact ndarray so that no PyArray_Return is * done in that case. */ +handle_out: + nargs = PyTuple_GET_SIZE(args); for (i = 0; i < nout; i++) { int j = nin + i; int incref = 1; diff --git a/numpy/core/src/umath/umath_tests.c.src b/numpy/core/src/umath/umath_tests.c.src index 46d06288d..33d9846bd 100644 --- a/numpy/core/src/umath/umath_tests.c.src +++ b/numpy/core/src/umath/umath_tests.c.src @@ -38,6 +38,9 @@ INIT_OUTER_LOOP_3 \ npy_intp s3 = *steps++; +#define BEGIN_OUTER_LOOP_2 \ + for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1) { + #define BEGIN_OUTER_LOOP_3 \ for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2) { @@ -58,13 +61,14 @@ char *inner1d_signature = "(i),(i)->()"; /**begin repeat #TYPE=LONG,DOUBLE# - #typ=npy_long, npy_double# + #typ=npy_long,npy_double# */ /* * This implements the function * out[n] = sum_i { in1[n, i] * in2[n, i] }. */ + static void @TYPE@_inner1d(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { @@ -91,7 +95,7 @@ char *innerwt_signature = "(i),(i),(i)->()"; /**begin repeat #TYPE=LONG,DOUBLE# - #typ=npy_long, npy_double# + #typ=npy_long,npy_double# */ @@ -127,7 +131,7 @@ char *matrix_multiply_signature = "(m,n),(n,p)->(m,p)"; /**begin repeat #TYPE=FLOAT,DOUBLE,LONG# - #typ=npy_float, npy_double,npy_long# + #typ=npy_float,npy_double,npy_long# */ /* @@ -135,7 +139,6 @@ char *matrix_multiply_signature = "(m,n),(n,p)->(m,p)"; * out[k, m, p] = sum_n { in1[k, m, n] * in2[k, n, p] }. */ - static void @TYPE@_matrix_multiply(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func)) { @@ -177,27 +180,66 @@ static void /**end repeat**/ -/* The following lines were generated using a slightly modified - version of code_generators/generate_umath.py and adding these - lines to defdict: - -defdict = { -'inner1d' : - Ufunc(2, 1, None_, - r'''inner on the last dimension and broadcast on the rest \n" - " \"(i),(i)->()\" \n''', - TD('ld'), - ), -'innerwt' : - Ufunc(3, 1, None_, - r'''inner1d with a weight argument \n" - " \"(i),(i),(i)->()\" \n''', - TD('ld'), - ), -} +char *euclidean_pdist_signature = "(n,d)->(p)"; +/**begin repeat + + #TYPE=FLOAT,DOUBLE# + #typ=npy_float,npy_double# + #sqrt_func=sqrtf,sqrt# */ +/* + * This implements the function + * out[j*(2*n-3-j)+k-1] = sum_d { (in1[j, d] - in1[k, d])^2 } + * with 0 < k < j < n, i.e. computes all unique pairwise euclidean distances. + */ + +static void +@TYPE@_euclidean_pdist(char **args, npy_intp *dimensions, npy_intp *steps, + void *NPY_UNUSED(func)) +{ + INIT_OUTER_LOOP_2 + npy_intp len_n = *dimensions++; + npy_intp len_d = *dimensions++; + npy_intp len_p = *dimensions; + npy_intp stride_n = *steps++; + npy_intp stride_d = *steps++; + npy_intp stride_p = *steps; + + assert(len_n * (len_n - 1) / 2 == len_p); + + BEGIN_OUTER_LOOP_2 + const char *data_this = (const char *)args[0]; + char *data_out = args[1]; + npy_intp n; + for (n = 0; n < len_n; ++n) { + const char *data_that = data_this + stride_n; + npy_intp nn; + for (nn = n + 1; nn < len_n; ++nn) { + const char *ptr_this = data_this; + const char *ptr_that = data_that; + @typ@ out = 0; + npy_intp d; + for (d = 0; d < len_d; ++d) { + const @typ@ delta = *(const @typ@ *)ptr_this - + *(const @typ@ *)ptr_that; + out += delta * delta; + ptr_this += stride_d; + ptr_that += stride_d; + } + *(@typ@ *)data_out = @sqrt_func@(out); + data_that += stride_n; + data_out += stride_p; + } + data_this += stride_n; + } + END_OUTER_LOOP +} + +/**end repeat**/ + + static PyUFuncGenericFunction inner1d_functions[] = { LONG_inner1d, DOUBLE_inner1d }; static void * inner1d_data[] = { (void *)NULL, (void *)NULL }; static char inner1d_signatures[] = { NPY_LONG, NPY_LONG, NPY_LONG, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE }; @@ -208,39 +250,49 @@ static PyUFuncGenericFunction matrix_multiply_functions[] = { LONG_matrix_multip static void *matrix_multiply_data[] = { (void *)NULL, (void *)NULL, (void *)NULL }; static char matrix_multiply_signatures[] = { NPY_LONG, NPY_LONG, NPY_LONG, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE }; +static PyUFuncGenericFunction euclidean_pdist_functions[] = + { FLOAT_euclidean_pdist, DOUBLE_euclidean_pdist }; +static void *eucldiean_pdist_data[] = { (void *)NULL, (void *)NULL }; +static char euclidean_pdist_signatures[] = { NPY_FLOAT, NPY_FLOAT, + NPY_DOUBLE, NPY_DOUBLE }; + + static void addUfuncs(PyObject *dictionary) { PyObject *f; - f = PyUFunc_FromFuncAndDataAndSignature(inner1d_functions, inner1d_data, inner1d_signatures, 2, - 2, 1, PyUFunc_None, "inner1d", - "inner on the last dimension and broadcast on the rest \n"\ - " \"(i),(i)->()\" \n", - 0, inner1d_signature); + f = PyUFunc_FromFuncAndDataAndSignature(inner1d_functions, inner1d_data, + inner1d_signatures, 2, 2, 1, PyUFunc_None, "inner1d", + "inner on the last dimension and broadcast on the rest \n" + " \"(i),(i)->()\" \n", + 0, inner1d_signature); PyDict_SetItemString(dictionary, "inner1d", f); Py_DECREF(f); - f = PyUFunc_FromFuncAndDataAndSignature(innerwt_functions, innerwt_data, innerwt_signatures, 2, - 3, 1, PyUFunc_None, "innerwt", - "inner1d with a weight argument \n"\ - " \"(i),(i),(i)->()\" \n", - 0, innerwt_signature); + f = PyUFunc_FromFuncAndDataAndSignature(innerwt_functions, innerwt_data, + innerwt_signatures, 2, 3, 1, PyUFunc_None, "innerwt", + "inner1d with a weight argument \n" + " \"(i),(i),(i)->()\" \n", + 0, innerwt_signature); PyDict_SetItemString(dictionary, "innerwt", f); Py_DECREF(f); f = PyUFunc_FromFuncAndDataAndSignature(matrix_multiply_functions, - matrix_multiply_data, matrix_multiply_signatures, - 3, 2, 1, PyUFunc_None, "matrix_multiply", - "matrix multiplication on last two dimensions \n"\ - " \"(m,n),(n,p)->(m,p)\" \n", - 0, matrix_multiply_signature); + matrix_multiply_data, matrix_multiply_signatures, + 3, 2, 1, PyUFunc_None, "matrix_multiply", + "matrix multiplication on last two dimensions \n" + " \"(m,n),(n,p)->(m,p)\" \n", + 0, matrix_multiply_signature); PyDict_SetItemString(dictionary, "matrix_multiply", f); Py_DECREF(f); + f = PyUFunc_FromFuncAndDataAndSignature(euclidean_pdist_functions, + eucldiean_pdist_data, euclidean_pdist_signatures, + 2, 1, 1, PyUFunc_None, "euclidean_pdist", + "pairwise euclidean distance on last two dimensions \n" + " \"(n,d)->(p)\" \n", + 0, euclidean_pdist_signature); + PyDict_SetItemString(dictionary, "euclidean_pdist", f); + Py_DECREF(f); } -/* - End of auto-generated code. -*/ - - static PyObject * UMath_Tests_test_signature(PyObject *NPY_UNUSED(dummy), PyObject *args) diff --git a/numpy/core/tests/test_datetime.py b/numpy/core/tests/test_datetime.py index bf0ba6807..4e432f885 100644 --- a/numpy/core/tests/test_datetime.py +++ b/numpy/core/tests/test_datetime.py @@ -1412,6 +1412,11 @@ class TestDateTime(TestCase): np.datetime64('2012-02-03T14Z', 's'), np.timedelta64(5, 'Y')) + def test_datetime_arange_no_dtype(self): + d = np.array('2010-01-04', dtype="M8[D]") + assert_equal(np.arange(d, d + 1), d) + assert_raises(ValueError, np.arange, d) + def test_timedelta_arange(self): a = np.arange(3, 10, dtype='m8') assert_equal(a.dtype, np.dtype('m8')) @@ -1430,6 +1435,11 @@ class TestDateTime(TestCase): assert_raises(TypeError, np.arange, np.timedelta64(0, 'Y'), np.timedelta64(5, 'D')) + def test_timedelta_arange_no_dtype(self): + d = np.array(5, dtype="m8[D]") + assert_equal(np.arange(d, d + 1), d) + assert_raises(ValueError, np.arange, d) + def test_datetime_maximum_reduce(self): a = np.array(['2010-01-02', '1999-03-14', '1833-03'], dtype='M8[D]') assert_equal(np.maximum.reduce(a).dtype, np.dtype('M8[D]')) diff --git a/numpy/core/tests/test_dtype.py b/numpy/core/tests/test_dtype.py index 18660351c..2621c8696 100644 --- a/numpy/core/tests/test_dtype.py +++ b/numpy/core/tests/test_dtype.py @@ -508,13 +508,8 @@ class TestDtypeAttributeDeletion(object): "isbuiltin", "isnative", "isalignedstruct", "fields", "metadata", "hasobject"] - if sys.version[:3] == '2.4': - error = TypeError - else: - error = AttributeError - for s in attr: - assert_raises(error, delattr, dt, s) + assert_raises(AttributeError, delattr, dt, s) def test_dtype_writable_attributes_deletion(self): diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py index f5a8bc4d2..1e47a2297 100644 --- a/numpy/core/tests/test_multiarray.py +++ b/numpy/core/tests/test_multiarray.py @@ -1757,6 +1757,40 @@ class TestMethods(TestCase): a.diagonal() assert_(sys.getrefcount(a) < 50) + def test_put(self): + icodes = np.typecodes['AllInteger'] + fcodes = np.typecodes['AllFloat'] + for dt in icodes + fcodes + 'O': + tgt = np.array([0, 1, 0, 3, 0, 5], dtype=dt) + + # test 1-d + a = np.zeros(6, dtype=dt) + a.put([1, 3, 5], [1, 3, 5]) + assert_equal(a, tgt) + + # test 2-d + a = np.zeros((2, 3), dtype=dt) + a.put([1, 3, 5], [1, 3, 5]) + assert_equal(a, tgt.reshape(2, 3)) + + for dt in '?': + tgt = np.array([False, True, False, True, False, True], dtype=dt) + + # test 1-d + a = np.zeros(6, dtype=dt) + a.put([1, 3, 5], [True]*3) + assert_equal(a, tgt) + + # test 2-d + a = np.zeros((2, 3), dtype=dt) + a.put([1, 3, 5], [True]*3) + assert_equal(a, tgt.reshape(2, 3)) + + # check must be writeable + a = np.zeros(6) + a.flags.writeable = False + assert_raises(ValueError, a.put, [1, 3, 5], [1, 3, 5]) + def test_ravel(self): a = np.array([[0, 1], [2, 3]]) assert_equal(a.ravel(), [0, 1, 2, 3]) @@ -1887,7 +1921,7 @@ class TestMethods(TestCase): a = np.array([1-1j, 1, 2.0, 'f'], object) assert_raises(AttributeError, lambda: a.conj()) - assert_raises(AttributeError, lambda: a.conjugate()) + assert_raises(AttributeError, lambda: a.conjugate()) class TestBinop(object): @@ -2128,6 +2162,16 @@ class TestBinop(object): assert_(isinstance(obj2, SomeClass2)) +class TestCAPI(TestCase): + def test_IsPythonScalar(self): + from numpy.core.multiarray_tests import IsPythonScalar + assert_(IsPythonScalar(b'foobar')) + assert_(IsPythonScalar(1)) + assert_(IsPythonScalar(2**80)) + assert_(IsPythonScalar(2.)) + assert_(IsPythonScalar("a")) + + class TestSubscripting(TestCase): def test_test_zero_rank(self): x = array([1, 2, 3]) @@ -4282,6 +4326,31 @@ class TestNewBufferProtocol(object): x3 = np.arange(dt3.itemsize, dtype=np.int8).view(dt3) self._check_roundtrip(x3) + def test_relaxed_strides(self): + # Test that relaxed strides are converted to non-relaxed + c = np.ones((1, 10, 10), dtype='i8') + + # Check for NPY_RELAXED_STRIDES_CHECKING: + if np.ones((10, 1), order="C").flags.f_contiguous: + c.strides = (-1, 80, 8) + + assert memoryview(c).strides == (800, 80, 8) + + # Writing C-contiguous data to a BytesIO buffer should work + fd = io.BytesIO() + fd.write(c.data) + + fortran = c.T + assert memoryview(fortran).strides == (8, 80, 800) + + arr = np.ones((1, 10)) + if arr.flags.f_contiguous: + shape, strides = get_buffer_info(arr, ['F_CONTIGUOUS']) + assert_(strides[0] == 8) + arr = np.ones((10, 1), order='F') + shape, strides = get_buffer_info(arr, ['C_CONTIGUOUS']) + assert_(strides[-1] == 8) + class TestArrayAttributeDeletion(object): diff --git a/numpy/core/tests/test_nditer.py b/numpy/core/tests/test_nditer.py index 0055c038b..dbf2cfaae 100644 --- a/numpy/core/tests/test_nditer.py +++ b/numpy/core/tests/test_nditer.py @@ -2481,13 +2481,8 @@ def test_iter_non_writable_attribute_deletion(): "iterationneedsapi", "has_multi_index", "has_index", "dtypes", "ndim", "nop", "itersize", "finished"] - if sys.version[:3] == '2.4': - error = TypeError - else: - error = AttributeError - for s in attr: - assert_raises(error, delattr, it, s) + assert_raises(AttributeError, delattr, it, s) def test_iter_writable_attribute_deletion(): diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py index 46e864495..ea145ef81 100644 --- a/numpy/core/tests/test_numeric.py +++ b/numpy/core/tests/test_numeric.py @@ -1697,6 +1697,20 @@ class TestStdVar(TestCase): assert_almost_equal(std(self.A, ddof=2)**2, self.real_var*len(self.A)/float(len(self.A)-2)) + def test_out_scalar(self): + d = np.arange(10) + out = np.array(0.) + r = np.std(d, out=out) + assert_(r is out) + assert_array_equal(r, out) + r = np.var(d, out=out) + assert_(r is out) + assert_array_equal(r, out) + r = np.mean(d, out=out) + assert_(r is out) + assert_array_equal(r, out) + + class TestStdVarComplex(TestCase): def test_basic(self): A = array([1, 1.j, -1, -1.j]) diff --git a/numpy/core/tests/test_records.py b/numpy/core/tests/test_records.py index 8c9ce5c70..355e5480a 100644 --- a/numpy/core/tests/test_records.py +++ b/numpy/core/tests/test_records.py @@ -1,5 +1,6 @@ from __future__ import division, absolute_import, print_function +import sys from os import path import numpy as np from numpy.testing import * @@ -15,6 +16,14 @@ class TestFromrecords(TestCase): r = np.rec.fromrecords([[456, 'dbe', 1.2], [2, 'de', 1.3]], names='col1,col2,col3') assert_equal(r[0].item(), (456, 'dbe', 1.2)) + assert_equal(r['col1'].dtype.kind, 'i') + if sys.version_info[0] >= 3: + assert_equal(r['col2'].dtype.kind, 'U') + assert_equal(r['col2'].dtype.itemsize, 12) + else: + assert_equal(r['col2'].dtype.kind, 'S') + assert_equal(r['col2'].dtype.itemsize, 3) + assert_equal(r['col3'].dtype.kind, 'f') def test_method_array(self): r = np.rec.array(asbytes('abcdefg') * 100, formats='i2,a3,i4', shape=3, byteorder='big') diff --git a/numpy/core/tests/test_ufunc.py b/numpy/core/tests/test_ufunc.py index eacc266be..a285d5334 100644 --- a/numpy/core/tests/test_ufunc.py +++ b/numpy/core/tests/test_ufunc.py @@ -388,21 +388,18 @@ class TestUfunc(TestCase): msg = "extend & broadcast loop dimensions" b = np.arange(4).reshape((2, 2)) assert_array_equal(umt.inner1d(a, b), np.sum(a*b, axis=-1), err_msg=msg) - msg = "broadcast in core dimensions" + # Broadcast in core dimensions should fail a = np.arange(8).reshape((4, 2)) b = np.arange(4).reshape((4, 1)) - assert_array_equal(umt.inner1d(a, b), np.sum(a*b, axis=-1), err_msg=msg) - msg = "extend & broadcast core and loop dimensions" + assert_raises(ValueError, umt.inner1d, a, b) + # Extend core dimensions should fail a = np.arange(8).reshape((4, 2)) b = np.array(7) - assert_array_equal(umt.inner1d(a, b), np.sum(a*b, axis=-1), err_msg=msg) - msg = "broadcast should fail" + assert_raises(ValueError, umt.inner1d, a, b) + # Broadcast should fail a = np.arange(2).reshape((2, 1, 1)) b = np.arange(3).reshape((3, 1, 1)) - try: - ret = umt.inner1d(a, b) - assert_equal(ret, None, err_msg=msg) - except ValueError: None + assert_raises(ValueError, umt.inner1d, a, b) def test_type_cast(self): msg = "type cast" @@ -542,8 +539,8 @@ class TestUfunc(TestCase): a2 = d2.transpose(p2)[s2] ref = ref and a1.base != None ref = ref and a2.base != None - if broadcastable(a1.shape[-1], a2.shape[-2]) and \ - broadcastable(a1.shape[0], a2.shape[0]): + if (a1.shape[-1] == a2.shape[-2] and + broadcastable(a1.shape[0], a2.shape[0])): assert_array_almost_equal( umt.matrix_multiply(a1, a2), np.sum(a2[..., np.newaxis].swapaxes(-3, -1) * @@ -553,6 +550,16 @@ class TestUfunc(TestCase): assert_equal(ref, True, err_msg="reference check") + def test_euclidean_pdist(self): + a = np.arange(12, dtype=np.float).reshape(4, 3) + out = np.empty((a.shape[0] * (a.shape[0] - 1) // 2,), dtype=a.dtype) + umt.euclidean_pdist(a, out) + b = np.sqrt(np.sum((a[:, None] - a)**2, axis=-1)) + b = b[~np.tri(a.shape[0], dtype=bool)] + assert_almost_equal(out, b) + # An output array is required to determine p with signature (n,d)->(p) + assert_raises(ValueError, umt.euclidean_pdist, a) + def test_object_logical(self): a = np.array([3, None, True, False, "test", ""], dtype=object) assert_equal(np.logical_or(a, None), diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py index f8d6520d7..c71b7b658 100644 --- a/numpy/core/tests/test_umath.py +++ b/numpy/core/tests/test_umath.py @@ -36,6 +36,36 @@ class TestConstants(TestCase): def test_euler_gamma(self): assert_allclose(ncu.euler_gamma, 0.5772156649015329, 1e-15) +class TestOut(TestCase): + def test_out_subok(self): + for b in (True, False): + aout = np.array(0.5) + + r = np.add(aout, 2, out=aout) + assert_(r is aout) + assert_array_equal(r, aout) + + r = np.add(aout, 2, out=aout, subok=b) + assert_(r is aout) + assert_array_equal(r, aout) + + r = np.add(aout, 2, aout, subok=False) + assert_(r is aout) + assert_array_equal(r, aout) + + d = np.ones(5) + o1 = np.zeros(5) + o2 = np.zeros(5, dtype=np.int32) + r1, r2 = np.frexp(d, o1, o2, subok=b) + assert_(r1 is o1) + assert_array_equal(r1, o1) + assert_(r2 is o2) + assert_array_equal(r2, o2) + + r1, r2 = np.frexp(d, out=o1, subok=b) + assert_(r1 is o1) + assert_array_equal(r1, o1) + class TestDivision(TestCase): def test_division_int(self): @@ -759,6 +789,17 @@ class TestBool(TestCase): assert_equal(np.bitwise_xor(arg1, arg2), out) +class TestInt(TestCase): + def test_logical_not(self): + x = np.ones(10, dtype=np.int16) + o = np.ones(10 * 2, dtype=np.bool) + tgt = o.copy() + tgt[::2] = False + os = o[::2] + assert_array_equal(np.logical_not(x, out=os), False) + assert_array_equal(o, tgt) + + class TestFloatingPoint(TestCase): def test_floating_point(self): assert_equal(ncu.FLOATING_POINT_SUPPORT, 1) diff --git a/numpy/distutils/ccompiler.py b/numpy/distutils/ccompiler.py index 8484685c0..d7fe702a6 100644 --- a/numpy/distutils/ccompiler.py +++ b/numpy/distutils/ccompiler.py @@ -16,7 +16,7 @@ from distutils.version import LooseVersion from numpy.distutils import log from numpy.distutils.exec_command import exec_command from numpy.distutils.misc_util import cyg2win32, is_sequence, mingw32, \ - quote_args + quote_args, get_num_build_jobs from numpy.distutils.compat import get_exception @@ -165,9 +165,10 @@ def CCompiler_compile(self, sources, output_dir=None, macros=None, return [] # FIXME:RELATIVE_IMPORT if sys.version_info[0] < 3: - from .fcompiler import FCompiler + from .fcompiler import FCompiler, is_f_file, has_f90_header else: - from numpy.distutils.fcompiler import FCompiler + from numpy.distutils.fcompiler import (FCompiler, is_f_file, + has_f90_header) if isinstance(self, FCompiler): display = [] for fc in ['f77', 'f90', 'fix']: @@ -189,20 +190,45 @@ def CCompiler_compile(self, sources, output_dir=None, macros=None, display += "\nextra options: '%s'" % (' '.join(extra_postargs)) log.info(display) - # build any sources in same order as they were originally specified - # especially important for fortran .f90 files using modules + def single_compile(args): + obj, (src, ext) = args + self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts) + if isinstance(self, FCompiler): objects_to_build = list(build.keys()) + f77_objects, other_objects = [], [] for obj in objects: if obj in objects_to_build: src, ext = build[obj] if self.compiler_type=='absoft': obj = cyg2win32(obj) src = cyg2win32(src) - self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts) + if is_f_file(src) and not has_f90_header(src): + f77_objects.append((obj, (src, ext))) + else: + other_objects.append((obj, (src, ext))) + + # f77 objects can be built in parallel + build_items = f77_objects + # build f90 modules serial, module files are generated during + # compilation and may be used by files later in the list so the + # ordering is important + for o in other_objects: + single_compile(o) + else: + build_items = build.items() + + jobs = get_num_build_jobs() + if len(build) > 1 and jobs > 1: + # build parallel + import multiprocessing.pool + pool = multiprocessing.pool.ThreadPool(jobs) + pool.map(single_compile, build_items) + pool.close() else: - for obj, (src, ext) in build.items(): - self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts) + # build serial + for o in build_items: + single_compile(o) # Return *all* object filenames, not just the ones we just built. return objects diff --git a/numpy/distutils/command/build.py b/numpy/distutils/command/build.py index b6912be15..f7249ae81 100644 --- a/numpy/distutils/command/build.py +++ b/numpy/distutils/command/build.py @@ -16,6 +16,8 @@ class build(old_build): user_options = old_build.user_options + [ ('fcompiler=', None, "specify the Fortran compiler type"), + ('jobs=', 'j', + "number of parallel jobs"), ] help_options = old_build.help_options + [ @@ -26,8 +28,14 @@ class build(old_build): def initialize_options(self): old_build.initialize_options(self) self.fcompiler = None + self.jobs = None def finalize_options(self): + if self.jobs: + try: + self.jobs = int(self.jobs) + except ValueError: + raise ValueError("--jobs/-j argument must be an integer") build_scripts = self.build_scripts old_build.finalize_options(self) plat_specifier = ".%s-%s" % (get_platform(), sys.version[0:3]) diff --git a/numpy/distutils/command/build_clib.py b/numpy/distutils/command/build_clib.py index 84ca87250..6e65a3bfb 100644 --- a/numpy/distutils/command/build_clib.py +++ b/numpy/distutils/command/build_clib.py @@ -30,6 +30,8 @@ class build_clib(old_build_clib): ('fcompiler=', None, "specify the Fortran compiler type"), ('inplace', 'i', 'Build in-place'), + ('jobs=', 'j', + "number of parallel jobs"), ] boolean_options = old_build_clib.boolean_options + ['inplace'] @@ -38,7 +40,16 @@ class build_clib(old_build_clib): old_build_clib.initialize_options(self) self.fcompiler = None self.inplace = 0 - return + self.jobs = None + + def finalize_options(self): + if self.jobs: + try: + self.jobs = int(self.jobs) + except ValueError: + raise ValueError("--jobs/-j argument must be an integer") + old_build_clib.finalize_options(self) + self.set_undefined_options('build', ('jobs', 'jobs')) def have_f_sources(self): for (lib_name, build_info) in self.libraries: diff --git a/numpy/distutils/command/build_ext.py b/numpy/distutils/command/build_ext.py index b48e4227a..59c453607 100644 --- a/numpy/distutils/command/build_ext.py +++ b/numpy/distutils/command/build_ext.py @@ -34,6 +34,8 @@ class build_ext (old_build_ext): user_options = old_build_ext.user_options + [ ('fcompiler=', None, "specify the Fortran compiler type"), + ('jobs=', 'j', + "number of parallel jobs"), ] help_options = old_build_ext.help_options + [ @@ -44,12 +46,19 @@ class build_ext (old_build_ext): def initialize_options(self): old_build_ext.initialize_options(self) self.fcompiler = None + self.jobs = None def finalize_options(self): + if self.jobs: + try: + self.jobs = int(self.jobs) + except ValueError: + raise ValueError("--jobs/-j argument must be an integer") incl_dirs = self.include_dirs old_build_ext.finalize_options(self) if incl_dirs is not None: self.include_dirs.extend(self.distribution.include_dirs or []) + self.set_undefined_options('build', ('jobs', 'jobs')) def run(self): if not self.extensions: @@ -407,11 +416,6 @@ class build_ext (old_build_ext): if ext.language=='c++' and cxx_compiler is not None: linker = cxx_compiler.link_shared_object - if sys.version[:3]>='2.3': - kws = {'target_lang':ext.language} - else: - kws = {} - linker(objects, ext_filename, libraries=libraries, library_dirs=library_dirs, @@ -419,7 +423,8 @@ class build_ext (old_build_ext): extra_postargs=extra_args, export_symbols=self.get_export_symbols(ext), debug=self.debug, - build_temp=self.build_temp,**kws) + build_temp=self.build_temp, + target_lang=ext.language) def _add_dummy_mingwex_sym(self, c_sources): build_src = self.get_finalized_command("build_src").build_src diff --git a/numpy/distutils/misc_util.py b/numpy/distutils/misc_util.py index c146178f0..19103f2c1 100644 --- a/numpy/distutils/misc_util.py +++ b/numpy/distutils/misc_util.py @@ -13,6 +13,10 @@ import shutil import distutils from distutils.errors import DistutilsError +try: + from threading import local as tlocal +except ImportError: + from dummy_threading import local as tlocal try: set @@ -31,7 +35,8 @@ __all__ = ['Configuration', 'get_numpy_include_dirs', 'default_config_dict', 'get_script_files', 'get_lib_source_files', 'get_data_files', 'dot_join', 'get_frame', 'minrelpath', 'njoin', 'is_sequence', 'is_string', 'as_list', 'gpaths', 'get_language', - 'quote_args', 'get_build_architecture', 'get_info', 'get_pkg_info'] + 'quote_args', 'get_build_architecture', 'get_info', 'get_pkg_info', + 'get_num_build_jobs'] class InstallableLib(object): """ @@ -60,6 +65,36 @@ class InstallableLib(object): self.build_info = build_info self.target_dir = target_dir + +def get_num_build_jobs(): + """ + Get number of parallel build jobs set by the --jobs command line argument + of setup.py + If the command did not receive a setting the environment variable + NPY_NUM_BUILD_JOBS checked and if that is unset it returns 1. + + Returns + ------- + out : int + number of parallel jobs that can be run + + """ + from numpy.distutils.core import get_distribution + envjobs = int(os.environ.get("NPY_NUM_BUILD_JOBS", 1)) + dist = get_distribution() + # may be None during configuration + if dist is None: + return envjobs + + # any of these three may have the job set, take the largest + cmdattr = (getattr(dist.get_command_obj('build'), 'jobs', None), + getattr(dist.get_command_obj('build_ext'), 'jobs', None), + getattr(dist.get_command_obj('build_clib'), 'jobs', None)) + if all(x is None for x in cmdattr): + return envjobs + else: + return max(x for x in cmdattr if x is not None) + def quote_args(args): # don't used _nt_quote_args as it does not check if # args items already have quotes or not. @@ -249,9 +284,9 @@ def gpaths(paths, local_path='', include_non_existing=True): return _fix_paths(paths, local_path, include_non_existing) -_temporary_directory = None def clean_up_temporary_directory(): - global _temporary_directory + tdata = tlocal() + _temporary_directory = getattr(tdata, 'tempdir', None) if not _temporary_directory: return try: @@ -261,13 +296,13 @@ def clean_up_temporary_directory(): _temporary_directory = None def make_temp_file(suffix='', prefix='', text=True): - global _temporary_directory - if not _temporary_directory: - _temporary_directory = tempfile.mkdtemp() + tdata = tlocal() + if not hasattr(tdata, 'tempdir'): + tdata.tempdir = tempfile.mkdtemp() atexit.register(clean_up_temporary_directory) fid, name = tempfile.mkstemp(suffix=suffix, prefix=prefix, - dir=_temporary_directory, + dir=tdata.tempdir, text=text) fo = os.fdopen(fid, 'w') return fo, name diff --git a/numpy/f2py/cb_rules.py b/numpy/f2py/cb_rules.py index f3bf848a7..5cf6b3010 100644 --- a/numpy/f2py/cb_rules.py +++ b/numpy/f2py/cb_rules.py @@ -357,11 +357,11 @@ cb_arg_rules=[ 'pyobjfrom': [{debugcapi:'\tfprintf(stderr,"debug-capi:cb:#varname#\\n");'}, {isintent_c: """\ \tif (#name#_nofargs>capi_i) { -\t\tPyArrayObject *tmp_arr = (PyArrayObject *)PyArray_New(&PyArray_Type,#rank#,#varname_i#_Dims,#atype#,NULL,(char*)#varname_i#,0,NPY_CARRAY,NULL); /*XXX: Hmm, what will destroy this array??? */ +\t\tPyArrayObject *tmp_arr = (PyArrayObject *)PyArray_New(&PyArray_Type,#rank#,#varname_i#_Dims,#atype#,NULL,(char*)#varname_i#,0,NPY_ARRAY_CARRAY,NULL); /*XXX: Hmm, what will destroy this array??? */ """, l_not(isintent_c): """\ \tif (#name#_nofargs>capi_i) { -\t\tPyArrayObject *tmp_arr = (PyArrayObject *)PyArray_New(&PyArray_Type,#rank#,#varname_i#_Dims,#atype#,NULL,(char*)#varname_i#,0,NPY_FARRAY,NULL); /*XXX: Hmm, what will destroy this array??? */ +\t\tPyArrayObject *tmp_arr = (PyArrayObject *)PyArray_New(&PyArray_Type,#rank#,#varname_i#_Dims,#atype#,NULL,(char*)#varname_i#,0,NPY_ARRAY_FARRAY,NULL); /*XXX: Hmm, what will destroy this array??? */ """, }, """ @@ -384,7 +384,7 @@ cb_arg_rules=[ \t\t\tfprintf(stderr,\"rv_cb_arr is NULL\\n\"); \t\t\tgoto capi_fail; \t\t} -\t\tMEMCOPY(#varname_i#,rv_cb_arr->data,PyArray_NBYTES(rv_cb_arr)); +\t\tMEMCOPY(#varname_i#,PyArray_DATA(rv_cb_arr),PyArray_NBYTES(rv_cb_arr)); \t\tif (capi_tmp != (PyObject *)rv_cb_arr) { \t\t\tPy_DECREF(rv_cb_arr); \t\t} diff --git a/numpy/f2py/cfuncs.py b/numpy/f2py/cfuncs.py index 7fb630697..01e189dc1 100644 --- a/numpy/f2py/cfuncs.py +++ b/numpy/f2py/cfuncs.py @@ -223,7 +223,7 @@ cppmacros['SWAP']="""\ \ta = b;\\ \tb = c;} """ -#cppmacros['ISCONTIGUOUS']='#define ISCONTIGUOUS(m) ((m)->flags & NPY_CONTIGUOUS)' +#cppmacros['ISCONTIGUOUS']='#define ISCONTIGUOUS(m) (PyArray_FLAGS(m) & NPY_ARRAY_C_CONTIGUOUS)' cppmacros['PRINTPYOBJERR']="""\ #define PRINTPYOBJERR(obj)\\ \tfprintf(stderr,\"#modulename#.error is related to \");\\ @@ -248,8 +248,8 @@ needs['len..']=['f2py_size'] cppmacros['len..']="""\ #define rank(var) var ## _Rank #define shape(var,dim) var ## _Dims[dim] -#define old_rank(var) (((PyArrayObject *)(capi_ ## var ## _tmp))->nd) -#define old_shape(var,dim) (((PyArrayObject *)(capi_ ## var ## _tmp))->dimensions[dim]) +#define old_rank(var) (PyArray_NDIM((PyArrayObject *)(capi_ ## var ## _tmp))) +#define old_shape(var,dim) PyArray_DIM(((PyArrayObject *)(capi_ ## var ## _tmp)),dim) #define fshape(var,dim) shape(var,rank(var)-dim-1) #define len(var) shape(var,0) #define flen(var) fshape(var,0) @@ -316,35 +316,35 @@ cppmacros['pyobj_from_string1size']='#define pyobj_from_string1size(v,len) (PyUS needs['TRYPYARRAYTEMPLATE']=['PRINTPYOBJERR'] cppmacros['TRYPYARRAYTEMPLATE']="""\ /* New SciPy */ -#define TRYPYARRAYTEMPLATECHAR case NPY_STRING: *(char *)(arr->data)=*v; break; -#define TRYPYARRAYTEMPLATELONG case NPY_LONG: *(long *)(arr->data)=*v; break; -#define TRYPYARRAYTEMPLATEOBJECT case NPY_OBJECT: (arr->descr->f->setitem)(pyobj_from_ ## ctype ## 1(*v),arr->data); break; +#define TRYPYARRAYTEMPLATECHAR case NPY_STRING: *(char *)(PyArray_DATA(arr))=*v; break; +#define TRYPYARRAYTEMPLATELONG case NPY_LONG: *(long *)(PyArray_DATA(arr))=*v; break; +#define TRYPYARRAYTEMPLATEOBJECT case NPY_OBJECT: (PyArray_DESCR(arr)->f->setitem)(pyobj_from_ ## ctype ## 1(*v),PyArray_DATA(arr)); break; #define TRYPYARRAYTEMPLATE(ctype,typecode) \\ PyArrayObject *arr = NULL;\\ if (!obj) return -2;\\ if (!PyArray_Check(obj)) return -1;\\ if (!(arr=(PyArrayObject *)obj)) {fprintf(stderr,\"TRYPYARRAYTEMPLATE:\");PRINTPYOBJERR(obj);return 0;}\\ - if (arr->descr->type==typecode) {*(ctype *)(arr->data)=*v; return 1;}\\ - switch (arr->descr->type_num) {\\ - case NPY_DOUBLE: *(double *)(arr->data)=*v; break;\\ - case NPY_INT: *(int *)(arr->data)=*v; break;\\ - case NPY_LONG: *(long *)(arr->data)=*v; break;\\ - case NPY_FLOAT: *(float *)(arr->data)=*v; break;\\ - case NPY_CDOUBLE: *(double *)(arr->data)=*v; break;\\ - case NPY_CFLOAT: *(float *)(arr->data)=*v; break;\\ - case NPY_BOOL: *(npy_bool *)(arr->data)=(*v!=0); break;\\ - case NPY_UBYTE: *(unsigned char *)(arr->data)=*v; break;\\ - case NPY_BYTE: *(signed char *)(arr->data)=*v; break;\\ - case NPY_SHORT: *(short *)(arr->data)=*v; break;\\ - case NPY_USHORT: *(npy_ushort *)(arr->data)=*v; break;\\ - case NPY_UINT: *(npy_uint *)(arr->data)=*v; break;\\ - case NPY_ULONG: *(npy_ulong *)(arr->data)=*v; break;\\ - case NPY_LONGLONG: *(npy_longlong *)(arr->data)=*v; break;\\ - case NPY_ULONGLONG: *(npy_ulonglong *)(arr->data)=*v; break;\\ - case NPY_LONGDOUBLE: *(npy_longdouble *)(arr->data)=*v; break;\\ - case NPY_CLONGDOUBLE: *(npy_longdouble *)(arr->data)=*v; break;\\ - case NPY_OBJECT: (arr->descr->f->setitem)(pyobj_from_ ## ctype ## 1(*v),arr->data, arr); break;\\ + if (PyArray_DESCR(arr)->type==typecode) {*(ctype *)(PyArray_DATA(arr))=*v; return 1;}\\ + switch (PyArray_TYPE(arr)) {\\ + case NPY_DOUBLE: *(double *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_INT: *(int *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_LONG: *(long *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_FLOAT: *(float *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_CDOUBLE: *(double *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_CFLOAT: *(float *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_BOOL: *(npy_bool *)(PyArray_DATA(arr))=(*v!=0); break;\\ + case NPY_UBYTE: *(unsigned char *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_BYTE: *(signed char *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_SHORT: *(short *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_USHORT: *(npy_ushort *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_UINT: *(npy_uint *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_ULONG: *(npy_ulong *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_LONGLONG: *(npy_longlong *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_ULONGLONG: *(npy_ulonglong *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_LONGDOUBLE: *(npy_longdouble *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_CLONGDOUBLE: *(npy_longdouble *)(PyArray_DATA(arr))=*v; break;\\ + case NPY_OBJECT: (PyArray_DESCR(arr)->f->setitem)(pyobj_from_ ## ctype ## 1(*v),PyArray_DATA(arr), arr); break;\\ default: return -2;\\ };\\ return 1 @@ -352,36 +352,36 @@ cppmacros['TRYPYARRAYTEMPLATE']="""\ needs['TRYCOMPLEXPYARRAYTEMPLATE']=['PRINTPYOBJERR'] cppmacros['TRYCOMPLEXPYARRAYTEMPLATE']="""\ -#define TRYCOMPLEXPYARRAYTEMPLATEOBJECT case NPY_OBJECT: (arr->descr->f->setitem)(pyobj_from_complex_ ## ctype ## 1((*v)),arr->data, arr); break; +#define TRYCOMPLEXPYARRAYTEMPLATEOBJECT case NPY_OBJECT: (PyArray_DESCR(arr)->f->setitem)(pyobj_from_complex_ ## ctype ## 1((*v)),PyArray_DATA(arr), arr); break; #define TRYCOMPLEXPYARRAYTEMPLATE(ctype,typecode)\\ PyArrayObject *arr = NULL;\\ if (!obj) return -2;\\ if (!PyArray_Check(obj)) return -1;\\ if (!(arr=(PyArrayObject *)obj)) {fprintf(stderr,\"TRYCOMPLEXPYARRAYTEMPLATE:\");PRINTPYOBJERR(obj);return 0;}\\ - if (arr->descr->type==typecode) {\\ - *(ctype *)(arr->data)=(*v).r;\\ - *(ctype *)(arr->data+sizeof(ctype))=(*v).i;\\ + if (PyArray_DESCR(arr)->type==typecode) {\\ + *(ctype *)(PyArray_DATA(arr))=(*v).r;\\ + *(ctype *)(PyArray_DATA(arr)+sizeof(ctype))=(*v).i;\\ return 1;\\ }\\ - switch (arr->descr->type_num) {\\ - case NPY_CDOUBLE: *(double *)(arr->data)=(*v).r;*(double *)(arr->data+sizeof(double))=(*v).i;break;\\ - case NPY_CFLOAT: *(float *)(arr->data)=(*v).r;*(float *)(arr->data+sizeof(float))=(*v).i;break;\\ - case NPY_DOUBLE: *(double *)(arr->data)=(*v).r; break;\\ - case NPY_LONG: *(long *)(arr->data)=(*v).r; break;\\ - case NPY_FLOAT: *(float *)(arr->data)=(*v).r; break;\\ - case NPY_INT: *(int *)(arr->data)=(*v).r; break;\\ - case NPY_SHORT: *(short *)(arr->data)=(*v).r; break;\\ - case NPY_UBYTE: *(unsigned char *)(arr->data)=(*v).r; break;\\ - case NPY_BYTE: *(signed char *)(arr->data)=(*v).r; break;\\ - case NPY_BOOL: *(npy_bool *)(arr->data)=((*v).r!=0 && (*v).i!=0); break;\\ - case NPY_USHORT: *(npy_ushort *)(arr->data)=(*v).r; break;\\ - case NPY_UINT: *(npy_uint *)(arr->data)=(*v).r; break;\\ - case NPY_ULONG: *(npy_ulong *)(arr->data)=(*v).r; break;\\ - case NPY_LONGLONG: *(npy_longlong *)(arr->data)=(*v).r; break;\\ - case NPY_ULONGLONG: *(npy_ulonglong *)(arr->data)=(*v).r; break;\\ - case NPY_LONGDOUBLE: *(npy_longdouble *)(arr->data)=(*v).r; break;\\ - case NPY_CLONGDOUBLE: *(npy_longdouble *)(arr->data)=(*v).r;*(npy_longdouble *)(arr->data+sizeof(npy_longdouble))=(*v).i;break;\\ - case NPY_OBJECT: (arr->descr->f->setitem)(pyobj_from_complex_ ## ctype ## 1((*v)),arr->data, arr); break;\\ + switch (PyArray_TYPE(arr)) {\\ + case NPY_CDOUBLE: *(double *)(PyArray_DATA(arr))=(*v).r;*(double *)(PyArray_DATA(arr)+sizeof(double))=(*v).i;break;\\ + case NPY_CFLOAT: *(float *)(PyArray_DATA(arr))=(*v).r;*(float *)(PyArray_DATA(arr)+sizeof(float))=(*v).i;break;\\ + case NPY_DOUBLE: *(double *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_LONG: *(long *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_FLOAT: *(float *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_INT: *(int *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_SHORT: *(short *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_UBYTE: *(unsigned char *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_BYTE: *(signed char *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_BOOL: *(npy_bool *)(PyArray_DATA(arr))=((*v).r!=0 && (*v).i!=0); break;\\ + case NPY_USHORT: *(npy_ushort *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_UINT: *(npy_uint *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_ULONG: *(npy_ulong *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_LONGLONG: *(npy_longlong *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_ULONGLONG: *(npy_ulonglong *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_LONGDOUBLE: *(npy_longdouble *)(PyArray_DATA(arr))=(*v).r; break;\\ + case NPY_CLONGDOUBLE: *(npy_longdouble *)(PyArray_DATA(arr))=(*v).r;*(npy_longdouble *)(PyArray_DATA(arr)+sizeof(npy_longdouble))=(*v).i;break;\\ + case NPY_OBJECT: (PyArray_DESCR(arr)->f->setitem)(pyobj_from_complex_ ## ctype ## 1((*v)),PyArray_DATA(arr), arr); break;\\ default: return -2;\\ };\\ return -1; @@ -391,11 +391,11 @@ cppmacros['TRYCOMPLEXPYARRAYTEMPLATE']="""\ ## \tif (PyArray_Check(obj)) arr = (PyArrayObject *)obj;\\ ## \telse arr = (PyArrayObject *)PyArray_ContiguousFromObject(obj,typenum,0,0);\\ ## \tif (arr) {\\ -## \t\tif (arr->descr->type_num==NPY_OBJECT) {\\ -## \t\t\tif (!ctype ## _from_pyobj(v,(arr->descr->getitem)(arr->data),\"\"))\\ +## \t\tif (PyArray_TYPE(arr)==NPY_OBJECT) {\\ +## \t\t\tif (!ctype ## _from_pyobj(v,(PyArray_DESCR(arr)->getitem)(PyArray_DATA(arr)),\"\"))\\ ## \t\t\tgoto capi_fail;\\ ## \t\t} else {\\ -## \t\t\t(arr->descr->cast[typenum])(arr->data,1,(char*)v,1,1);\\ +## \t\t\t(PyArray_DESCR(arr)->cast[typenum])(PyArray_DATA(arr),1,(char*)v,1,1);\\ ## \t\t}\\ ## \t\tif ((PyObject *)arr != obj) { Py_DECREF(arr); }\\ ## \t\treturn 1;\\ @@ -407,11 +407,11 @@ cppmacros['TRYCOMPLEXPYARRAYTEMPLATE']="""\ ## \tif (PyArray_Check(obj)) arr = (PyArrayObject *)obj;\\ ## \telse arr = (PyArrayObject *)PyArray_ContiguousFromObject(obj,typenum,0,0);\\ ## \tif (arr) {\\ -## \t\tif (arr->descr->type_num==NPY_OBJECT) {\\ -## \t\t\tif (!ctype ## _from_pyobj(v,(arr->descr->getitem)(arr->data),\"\"))\\ +## \t\tif (PyArray_TYPE(arr)==NPY_OBJECT) {\\ +## \t\t\tif (!ctype ## _from_pyobj(v,(PyArray_DESCR(arr)->getitem)(PyArray_DATA(arr)),\"\"))\\ ## \t\t\tgoto capi_fail;\\ ## \t\t} else {\\ -## \t\t\t(arr->descr->cast[typenum])((void *)(arr->data),1,(void *)(v),1,1);\\ +## \t\t\t(PyArray_DESCR(arr)->cast[typenum])((void *)(PyArray_DATA(arr)),1,(void *)(v),1,1);\\ ## \t\t}\\ ## \t\tif ((PyObject *)arr != obj) { Py_DECREF(arr); }\\ ## \t\treturn 1;\\ @@ -536,15 +536,15 @@ cppmacros['OLDPYNUM']="""\ cfuncs['calcarrindex']="""\ static int calcarrindex(int *i,PyArrayObject *arr) { \tint k,ii = i[0]; -\tfor (k=1; k < arr->nd; k++) -\t\tii += (ii*(arr->dimensions[k] - 1)+i[k]); /* assuming contiguous arr */ +\tfor (k=1; k < PyArray_NDIM(arr); k++) +\t\tii += (ii*(PyArray_DIM(arr,k) - 1)+i[k]); /* assuming contiguous arr */ \treturn ii; }""" cfuncs['calcarrindextr']="""\ static int calcarrindextr(int *i,PyArrayObject *arr) { -\tint k,ii = i[arr->nd-1]; -\tfor (k=1; k < arr->nd; k++) -\t\tii += (ii*(arr->dimensions[arr->nd-k-1] - 1)+i[arr->nd-k-1]); /* assuming contiguous arr */ +\tint k,ii = i[PyArray_NDIM(arr)-1]; +\tfor (k=1; k < PyArray_NDIM(arr); k++) +\t\tii += (ii*(PyArray_DIM(arr,PyArray_NDIM(arr)-k-1) - 1)+i[PyArray_NDIM(arr)-k-1]); /* assuming contiguous arr */ \treturn ii; }""" cfuncs['forcomb']="""\ @@ -592,7 +592,7 @@ cfuncs['try_pyarr_from_string']="""\ static int try_pyarr_from_string(PyObject *obj,const string str) { \tPyArrayObject *arr = NULL; \tif (PyArray_Check(obj) && (!((arr = (PyArrayObject *)obj) == NULL))) -\t\t{ STRINGCOPYN(arr->data,str,PyArray_NBYTES(arr)); } +\t\t{ STRINGCOPYN(PyArray_DATA(arr),str,PyArray_NBYTES(arr)); } \treturn 1; capi_fail: \tPRINTPYOBJERR(obj); @@ -623,9 +623,9 @@ fprintf(stderr,\"string_from_pyobj(str='%s',len=%d,inistr='%s',obj=%p)\\n\",(cha \t\t\tgoto capi_fail; \t\t} \t\tif (*len == -1) -\t\t\t*len = (arr->descr->elsize)*PyArray_SIZE(arr); +\t\t\t*len = (PyArray_ITEMSIZE(arr))*PyArray_SIZE(arr); \t\tSTRINGMALLOC(*str,*len); -\t\tSTRINGCOPYN(*str,arr->data,*len+1); +\t\tSTRINGCOPYN(*str,PyArray_DATA(arr),*len+1); \t\treturn 1; \t} \tif (PyString_Check(obj)) { diff --git a/numpy/f2py/crackfortran.py b/numpy/f2py/crackfortran.py index 893081126..0fde37bcf 100755 --- a/numpy/f2py/crackfortran.py +++ b/numpy/f2py/crackfortran.py @@ -1274,7 +1274,7 @@ def markinnerspaces(line): cb='' for c in line: if cb=='\\' and c in ['\\', '\'', '"']: - l=l+c; + l=l+c cb=c continue if f==0 and c in ['\'', '"']: cc=c; cc1={'\'':'"','"':'\''}[c] @@ -2198,8 +2198,10 @@ def analyzevars(block): if 'intent' not in vars[n]: vars[n]['intent']=[] for c in [x.strip() for x in markoutercomma(intent).split('@,@')]: - if not c in vars[n]['intent']: - vars[n]['intent'].append(c) + # Remove spaces so that 'in out' becomes 'inout' + tmp = c.replace(' ', '') + if tmp not in vars[n]['intent']: + vars[n]['intent'].append(tmp) intent=None if note: note=note.replace('\\n\\n', '\n\n') @@ -2220,7 +2222,7 @@ def analyzevars(block): if 'check' not in vars[n]: vars[n]['check']=[] for c in [x.strip() for x in markoutercomma(check).split('@,@')]: - if not c in vars[n]['check']: + if c not in vars[n]['check']: vars[n]['check'].append(c) check=None if dim and 'dimension' not in vars[n]: diff --git a/numpy/f2py/rules.py b/numpy/f2py/rules.py index 4c186712c..b46776777 100644 --- a/numpy/f2py/rules.py +++ b/numpy/f2py/rules.py @@ -109,6 +109,8 @@ module_rules={ * $D"""+"""ate:$ * Do not edit this file directly unless you know what you are doing!!! */ + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #ifdef __cplusplus extern \"C\" { #endif @@ -232,6 +234,7 @@ PyMODINIT_FUNC init#modulename#(void) { #ifdef __cplusplus } #endif +#undef NPY_NO_DEPRECATED_API """, 'separatorsfor':{'latexdoc':'\n\n', 'restdoc':'\n\n'}, @@ -1027,9 +1030,9 @@ if (#varname#_capi==Py_None) { # \tif ((#varname#_capi != Py_None) && PyArray_Check(#varname#_capi) \\ # \t\t&& (#varname#_capi != (PyObject *)capi_#varname#_tmp)) { -# \t\tif (((PyArrayObject *)#varname#_capi)->nd != capi_#varname#_tmp->nd) { -# \t\t\tif (#varname#_capi != capi_#varname#_tmp->base) -# \t\t\t\tcopy_ND_array((PyArrayObject *)capi_#varname#_tmp->base,(PyArrayObject *)#varname#_capi); +# \t\tif (PyArray_NDIM((PyArrayObject *)#varname#_capi) != PyArray_NDIM(capi_#varname#_tmp)) { +# \t\t\tif (#varname#_capi != PyArray_BASE(capi_#varname#_tmp)) +# \t\t\t\tcopy_ND_array(PyArray_BASE((PyArrayObject *)capi_#varname#_tmp),(PyArrayObject *)#varname#_capi); # \t\t} else # \t\t\tcopy_ND_array(capi_#varname#_tmp,(PyArrayObject *)#varname#_capi); # \t} @@ -1047,7 +1050,7 @@ if (#varname#_capi==Py_None) { \t\tif (!PyErr_Occurred()) \t\t\tPyErr_SetString(#modulename#_error,\"failed in converting #nth# `#varname#\' of #pyname# to C/Fortran array\" ); \t} else { -\t\t#varname# = (#ctype# *)(capi_#varname#_tmp->data); +\t\t#varname# = (#ctype# *)(PyArray_DATA(capi_#varname#_tmp)); """, {hasinitvalue:[ {isintent_nothide:'\tif (#varname#_capi == Py_None) {'}, @@ -1056,7 +1059,7 @@ if (#varname#_capi==Py_None) { """\ \t\tint *_i,capi_i=0; \t\tCFUNCSMESS(\"#name#: Initializing #varname#=#init#\\n\"); -\t\tif (initforcomb(capi_#varname#_tmp->dimensions,capi_#varname#_tmp->nd,1)) { +\t\tif (initforcomb(PyArray_DIMS(capi_#varname#_tmp),PyArray_NDIM(capi_#varname#_tmp),1)) { \t\t\twhile ((_i = nextforcomb())) \t\t\t\t#varname#[capi_i++] = #init#; /* fortran way */ \t\t} else { diff --git a/numpy/f2py/setup.py b/numpy/f2py/setup.py index 3c6140b53..a27a001a9 100644 --- a/numpy/f2py/setup.py +++ b/numpy/f2py/setup.py @@ -88,26 +88,24 @@ main() if __name__ == "__main__": config = configuration(top_path='') - version = config.get_version() print('F2PY Version', version) config = config.todict() - if sys.version[:3]>='2.3': - config['download_url'] = "http://cens.ioc.ee/projects/f2py2e/2.x"\ - "/F2PY-2-latest.tar.gz" - config['classifiers'] = [ - 'Development Status :: 5 - Production/Stable', - 'Intended Audience :: Developers', - 'Intended Audience :: Science/Research', - 'License :: OSI Approved :: NumPy License', - 'Natural Language :: English', - 'Operating System :: OS Independent', - 'Programming Language :: C', - 'Programming Language :: Fortran', - 'Programming Language :: Python', - 'Topic :: Scientific/Engineering', - 'Topic :: Software Development :: Code Generators', - ] + config['download_url'] = "http://cens.ioc.ee/projects/f2py2e/2.x"\ + "/F2PY-2-latest.tar.gz" + config['classifiers'] = [ + 'Development Status :: 5 - Production/Stable', + 'Intended Audience :: Developers', + 'Intended Audience :: Science/Research', + 'License :: OSI Approved :: NumPy License', + 'Natural Language :: English', + 'Operating System :: OS Independent', + 'Programming Language :: C', + 'Programming Language :: Fortran', + 'Programming Language :: Python', + 'Topic :: Scientific/Engineering', + 'Topic :: Software Development :: Code Generators', + ] setup(version=version, description = "F2PY - Fortran to Python Interface Generaton", author = "Pearu Peterson", diff --git a/numpy/f2py/src/fortranobject.c b/numpy/f2py/src/fortranobject.c index 001a4c7de..e88e83f0e 100644 --- a/numpy/f2py/src/fortranobject.c +++ b/numpy/f2py/src/fortranobject.c @@ -1,4 +1,5 @@ #define FORTRANOBJECT_C +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include "fortranobject.h" #ifdef __cplusplus @@ -58,11 +59,11 @@ PyFortranObject_New(FortranDataDef* defs, f2py_void_func init) { int n = fp->defs[i].rank-1; v = PyArray_New(&PyArray_Type, n, fp->defs[i].dims.d, NPY_STRING, NULL, fp->defs[i].data, fp->defs[i].dims.d[n], - NPY_FARRAY, NULL); + NPY_ARRAY_FARRAY, NULL); } else { v = PyArray_New(&PyArray_Type, fp->defs[i].rank, fp->defs[i].dims.d, - fp->defs[i].type, NULL, fp->defs[i].data, 0, NPY_FARRAY, + fp->defs[i].type, NULL, fp->defs[i].data, 0, NPY_ARRAY_FARRAY, NULL); } if (v==NULL) return NULL; @@ -273,7 +274,7 @@ fortran_getattr(PyFortranObject *fp, char *name) { k = fp->defs[i].rank; if (fp->defs[i].data !=NULL) { /* array is allocated */ PyObject *v = PyArray_New(&PyArray_Type, k, fp->defs[i].dims.d, - fp->defs[i].type, NULL, fp->defs[i].data, 0, NPY_FARRAY, + fp->defs[i].type, NULL, fp->defs[i].data, 0, NPY_ARRAY_FARRAY, NULL); if (v==NULL) return NULL; /* Py_INCREF(v); */ @@ -345,7 +346,7 @@ fortran_setattr(PyFortranObject *fp, char *name, PyObject *v) { for(k=0;k<fp->defs[i].rank;k++) dims[k]=-1; if ((arr = array_from_pyobj(fp->defs[i].type,dims,fp->defs[i].rank,F2PY_INTENT_IN,v))==NULL) return -1; - (*(fp->defs[i].func))(&fp->defs[i].rank,arr->dimensions,set_data,&flag); + (*(fp->defs[i].func))(&fp->defs[i].rank,PyArray_DIMS(arr),set_data,&flag); } else { /* deallocate */ for(k=0;k<fp->defs[i].rank;k++) dims[k]=0; (*(fp->defs[i].func))(&fp->defs[i].rank,dims,set_data,&flag); @@ -357,11 +358,11 @@ fortran_setattr(PyFortranObject *fp, char *name, PyObject *v) { return -1; } if (fp->defs[i].data!=NULL) { /* copy Python object to Fortran array */ - npy_intp s = PyArray_MultiplyList(fp->defs[i].dims.d,arr->nd); + npy_intp s = PyArray_MultiplyList(fp->defs[i].dims.d,PyArray_NDIM(arr)); if (s==-1) - s = PyArray_MultiplyList(arr->dimensions,arr->nd); + s = PyArray_MultiplyList(PyArray_DIMS(arr),PyArray_NDIM(arr)); if (s<0 || - (memcpy(fp->defs[i].data,arr->data,s*PyArray_ITEMSIZE(arr)))==NULL) { + (memcpy(fp->defs[i].data,PyArray_DATA(arr),s*PyArray_ITEMSIZE(arr)))==NULL) { if ((PyObject*)arr!=v) { Py_DECREF(arr); } @@ -616,8 +617,9 @@ void dump_dims(int rank, npy_intp* dims) { } printf("]\n"); } -void dump_attrs(const PyArrayObject* arr) { - int rank = arr->nd; +void dump_attrs(const PyArrayObject* obj) { + const PyArrayObject_fields *arr = (const PyArrayObject_fields*) obj; + int rank = PyArray_NDIM(arr); npy_intp size = PyArray_Size((PyObject *)arr); printf("\trank = %d, flags = %d, size = %" NPY_INTP_FMT "\n", rank,arr->flags,size); @@ -630,7 +632,9 @@ void dump_attrs(const PyArrayObject* arr) { #define SWAPTYPE(a,b,t) {t c; c = (a); (a) = (b); (b) = c; } -static int swap_arrays(PyArrayObject* arr1, PyArrayObject* arr2) { +static int swap_arrays(PyArrayObject* obj1, PyArrayObject* obj2) { + PyArrayObject_fields *arr1 = (PyArrayObject_fields*) obj1, + *arr2 = (PyArrayObject_fields*) obj2; SWAPTYPE(arr1->data,arr2->data,char*); SWAPTYPE(arr1->nd,arr2->nd,int); SWAPTYPE(arr1->dimensions,arr2->dimensions,npy_intp*); @@ -707,17 +711,17 @@ PyArrayObject* array_from_pyobj(const int type_num, if (intent & F2PY_INTENT_CACHE) { /* intent(cache) */ - if (PyArray_ISONESEGMENT(obj) - && PyArray_ITEMSIZE((PyArrayObject *)obj)>=elsize) { - if (check_and_fix_dimensions((PyArrayObject *)obj,rank,dims)) { + if (PyArray_ISONESEGMENT(arr) + && PyArray_ITEMSIZE(arr)>=elsize) { + if (check_and_fix_dimensions(arr,rank,dims)) { return NULL; /*XXX: set exception */ } if (intent & F2PY_INTENT_OUT) - Py_INCREF(obj); - return (PyArrayObject *)obj; + Py_INCREF(arr); + return arr; } strcpy(mess, "failed to initialize intent(cache) array"); - if (!PyArray_ISONESEGMENT(obj)) + if (!PyArray_ISONESEGMENT(arr)) strcat(mess, " -- input must be in one segment"); if (PyArray_ITEMSIZE(arr)<elsize) sprintf(mess+strlen(mess)," -- expected at least elsize=%d but got %d", @@ -766,7 +770,7 @@ PyArrayObject* array_from_pyobj(const int type_num, ); if (!(ARRAY_ISCOMPATIBLE(arr,type_num))) sprintf(mess+strlen(mess)," -- input '%c' not compatible to '%c'", - arr->descr->type,typechar); + PyArray_DESCR(arr)->type,typechar); if (!(F2PY_CHECK_ALIGNMENT(arr, intent))) sprintf(mess+strlen(mess)," -- input not %d-aligned", F2PY_GET_ALIGNMENT(intent)); PyErr_SetString(PyExc_ValueError,mess); @@ -777,7 +781,7 @@ PyArrayObject* array_from_pyobj(const int type_num, { PyArrayObject *retarr = (PyArrayObject *) \ - PyArray_New(&PyArray_Type, arr->nd, arr->dimensions, type_num, + PyArray_New(&PyArray_Type, PyArray_NDIM(arr), PyArray_DIMS(arr), type_num, NULL,NULL,0, !(intent&F2PY_INTENT_C), NULL); @@ -814,8 +818,8 @@ PyArrayObject* array_from_pyobj(const int type_num, F2PY_REPORT_ON_ARRAY_COPY_FROMANY; arr = (PyArrayObject *) \ PyArray_FromAny(obj,PyArray_DescrFromType(type_num), 0,0, - ((intent & F2PY_INTENT_C)?NPY_CARRAY:NPY_FARRAY) \ - | NPY_FORCECAST, NULL); + ((intent & F2PY_INTENT_C)?NPY_ARRAY_CARRAY:NPY_ARRAY_FARRAY) \ + | NPY_ARRAY_FORCECAST, NULL); if (arr==NULL) return NULL; if (check_and_fix_dimensions(arr,rank,dims)) @@ -836,20 +840,20 @@ int check_and_fix_dimensions(const PyArrayObject* arr,const int rank,npy_intp *d the dimensions from arr. It also checks that non-blank dims will match with the corresponding values in arr dimensions. */ - const npy_intp arr_size = (arr->nd)?PyArray_Size((PyObject *)arr):1; + const npy_intp arr_size = (PyArray_NDIM(arr))?PyArray_Size((PyObject *)arr):1; #ifdef DEBUG_COPY_ND_ARRAY dump_attrs(arr); printf("check_and_fix_dimensions:init: dims="); dump_dims(rank,dims); #endif - if (rank > arr->nd) { /* [1,2] -> [[1],[2]]; 1 -> [[1]] */ + if (rank > PyArray_NDIM(arr)) { /* [1,2] -> [[1],[2]]; 1 -> [[1]] */ npy_intp new_size = 1; int free_axe = -1; int i; npy_intp d; /* Fill dims where -1 or 0; check dimensions; calc new_size; */ - for(i=0;i<arr->nd;++i) { - d = arr->dimensions[i]; + for(i=0;i<PyArray_NDIM(arr);++i) { + d = PyArray_DIM(arr,i); if (dims[i] >= 0) { if (d>1 && dims[i]!=d) { fprintf(stderr,"%d-th dimension must be fixed to %" NPY_INTP_FMT @@ -863,7 +867,7 @@ int check_and_fix_dimensions(const PyArrayObject* arr,const int rank,npy_intp *d } new_size *= dims[i]; } - for(i=arr->nd;i<rank;++i) + for(i=PyArray_NDIM(arr);i<rank;++i) if (dims[i]>1) { fprintf(stderr,"%d-th dimension must be %" NPY_INTP_FMT " but got 0 (not defined).\n", @@ -883,12 +887,12 @@ int check_and_fix_dimensions(const PyArrayObject* arr,const int rank,npy_intp *d " indices)\n", new_size,arr_size); return 1; } - } else if (rank==arr->nd) { + } else if (rank==PyArray_NDIM(arr)) { npy_intp new_size = 1; int i; npy_intp d; for (i=0; i<rank; ++i) { - d = arr->dimensions[i]; + d = PyArray_DIM(arr,i); if (dims[i]>=0) { if (d > 1 && d!=dims[i]) { fprintf(stderr,"%d-th dimension must be fixed to %" NPY_INTP_FMT @@ -910,19 +914,19 @@ int check_and_fix_dimensions(const PyArrayObject* arr,const int rank,npy_intp *d npy_intp d; int effrank; npy_intp size; - for (i=0,effrank=0;i<arr->nd;++i) - if (arr->dimensions[i]>1) ++effrank; + for (i=0,effrank=0;i<PyArray_NDIM(arr);++i) + if (PyArray_DIM(arr,i)>1) ++effrank; if (dims[rank-1]>=0) if (effrank>rank) { fprintf(stderr,"too many axes: %d (effrank=%d), expected rank=%d\n", - arr->nd,effrank,rank); + PyArray_NDIM(arr),effrank,rank); return 1; } for (i=0,j=0;i<rank;++i) { - while (j<arr->nd && arr->dimensions[j]<2) ++j; - if (j>=arr->nd) d = 1; - else d = arr->dimensions[j++]; + while (j<PyArray_NDIM(arr) && PyArray_DIM(arr,j)<2) ++j; + if (j>=PyArray_NDIM(arr)) d = 1; + else d = PyArray_DIM(arr,j++); if (dims[i]>=0) { if (d>1 && d!=dims[i]) { fprintf(stderr,"%d-th dimension must be fixed to %" NPY_INTP_FMT @@ -935,20 +939,20 @@ int check_and_fix_dimensions(const PyArrayObject* arr,const int rank,npy_intp *d dims[i] = d; } - for (i=rank;i<arr->nd;++i) { /* [[1,2],[3,4]] -> [1,2,3,4] */ - while (j<arr->nd && arr->dimensions[j]<2) ++j; - if (j>=arr->nd) d = 1; - else d = arr->dimensions[j++]; + for (i=rank;i<PyArray_NDIM(arr);++i) { /* [[1,2],[3,4]] -> [1,2,3,4] */ + while (j<PyArray_NDIM(arr) && PyArray_DIM(arr,j)<2) ++j; + if (j>=PyArray_NDIM(arr)) d = 1; + else d = PyArray_DIM(arr,j++); dims[rank-1] *= d; } for (i=0,size=1;i<rank;++i) size *= dims[i]; if (size != arr_size) { fprintf(stderr,"unexpected array size: size=%" NPY_INTP_FMT ", arr_size=%" NPY_INTP_FMT ", rank=%d, effrank=%d, arr.nd=%d, dims=[", - size,arr_size,rank,effrank,arr->nd); + size,arr_size,rank,effrank,PyArray_NDIM(arr)); for (i=0;i<rank;++i) fprintf(stderr," %" NPY_INTP_FMT,dims[i]); fprintf(stderr," ], arr.dims=["); - for (i=0;i<arr->nd;++i) fprintf(stderr," %" NPY_INTP_FMT,arr->dimensions[i]); + for (i=0;i<PyArray_NDIM(arr);++i) fprintf(stderr," %" NPY_INTP_FMT,PyArray_DIM(arr,i)); fprintf(stderr," ]\n"); return 1; } @@ -1029,4 +1033,6 @@ F2PyCapsule_Check(PyObject *ptr) #ifdef __cplusplus } #endif + +#undef NPY_NO_DEPRECATED_API /************************* EOF fortranobject.c *******************************/ diff --git a/numpy/f2py/src/fortranobject.h b/numpy/f2py/src/fortranobject.h index 689f78c92..c9b54e259 100644 --- a/numpy/f2py/src/fortranobject.h +++ b/numpy/f2py/src/fortranobject.h @@ -119,7 +119,7 @@ int F2PyCapsule_Check(PyObject *ptr); #endif -#define ISCONTIGUOUS(m) ((m)->flags & NPY_CONTIGUOUS) +#define ISCONTIGUOUS(m) (PyArray_FLAGS(m) & NPY_ARRAY_C_CONTIGUOUS) #define F2PY_INTENT_IN 1 #define F2PY_INTENT_INOUT 2 #define F2PY_INTENT_OUT 4 diff --git a/numpy/f2py/tests/src/array_from_pyobj/wrapmodule.c b/numpy/f2py/tests/src/array_from_pyobj/wrapmodule.c index 9a21f2420..2da6a2c5d 100644 --- a/numpy/f2py/tests/src/array_from_pyobj/wrapmodule.c +++ b/numpy/f2py/tests/src/array_from_pyobj/wrapmodule.c @@ -90,22 +90,22 @@ static PyObject *f2py_rout_wrap_attrs(PyObject *capi_self, &PyArray_Type,&arr_capi)) return NULL; arr = (PyArrayObject *)arr_capi; - sprintf(s,"%p",arr->data); - dimensions = PyTuple_New(arr->nd); - strides = PyTuple_New(arr->nd); - for (i=0;i<arr->nd;++i) { - PyTuple_SetItem(dimensions,i,PyInt_FromLong(arr->dimensions[i])); - PyTuple_SetItem(strides,i,PyInt_FromLong(arr->strides[i])); + sprintf(s,"%p",PyArray_DATA(arr)); + dimensions = PyTuple_New(PyArray_NDIM(arr)); + strides = PyTuple_New(PyArray_NDIM(arr)); + for (i=0;i<PyArray_NDIM(arr);++i) { + PyTuple_SetItem(dimensions,i,PyInt_FromLong(PyArray_DIM(arr,i))); + PyTuple_SetItem(strides,i,PyInt_FromLong(PyArray_STRIDE(arr,i))); } - return Py_BuildValue("siOOO(cciii)ii",s,arr->nd, + return Py_BuildValue("siOOO(cciii)ii",s,PyArray_NDIM(arr), dimensions,strides, - (arr->base==NULL?Py_None:arr->base), - arr->descr->kind, - arr->descr->type, - arr->descr->type_num, - arr->descr->elsize, - arr->descr->alignment, - arr->flags, + (PyArray_BASE(arr)==NULL?Py_None:PyArray_BASE(arr)), + PyArray_DESCR(arr)->kind, + PyArray_DESCR(arr)->type, + PyArray_TYPE(arr), + PyArray_ITEMSIZE(arr), + PyArray_DESCR(arr)->alignment, + PyArray_FLAGS(arr), PyArray_ITEMSIZE(arr)); } @@ -190,24 +190,24 @@ PyMODINIT_FUNC inittest_array_from_pyobj_ext(void) { PyDict_SetItemString(d, "NPY_NOTYPE", PyInt_FromLong(NPY_NOTYPE)); PyDict_SetItemString(d, "NPY_USERDEF", PyInt_FromLong(NPY_USERDEF)); - PyDict_SetItemString(d, "CONTIGUOUS", PyInt_FromLong(NPY_CONTIGUOUS)); - PyDict_SetItemString(d, "FORTRAN", PyInt_FromLong(NPY_FORTRAN)); - PyDict_SetItemString(d, "OWNDATA", PyInt_FromLong(NPY_OWNDATA)); - PyDict_SetItemString(d, "FORCECAST", PyInt_FromLong(NPY_FORCECAST)); - PyDict_SetItemString(d, "ENSURECOPY", PyInt_FromLong(NPY_ENSURECOPY)); - PyDict_SetItemString(d, "ENSUREARRAY", PyInt_FromLong(NPY_ENSUREARRAY)); - PyDict_SetItemString(d, "ALIGNED", PyInt_FromLong(NPY_ALIGNED)); - PyDict_SetItemString(d, "WRITEABLE", PyInt_FromLong(NPY_WRITEABLE)); - PyDict_SetItemString(d, "UPDATEIFCOPY", PyInt_FromLong(NPY_UPDATEIFCOPY)); + PyDict_SetItemString(d, "CONTIGUOUS", PyInt_FromLong(NPY_ARRAY_C_CONTIGUOUS)); + PyDict_SetItemString(d, "FORTRAN", PyInt_FromLong(NPY_ARRAY_F_CONTIGUOUS)); + PyDict_SetItemString(d, "OWNDATA", PyInt_FromLong(NPY_ARRAY_OWNDATA)); + PyDict_SetItemString(d, "FORCECAST", PyInt_FromLong(NPY_ARRAY_FORCECAST)); + PyDict_SetItemString(d, "ENSURECOPY", PyInt_FromLong(NPY_ARRAY_ENSURECOPY)); + PyDict_SetItemString(d, "ENSUREARRAY", PyInt_FromLong(NPY_ARRAY_ENSUREARRAY)); + PyDict_SetItemString(d, "ALIGNED", PyInt_FromLong(NPY_ARRAY_ALIGNED)); + PyDict_SetItemString(d, "WRITEABLE", PyInt_FromLong(NPY_ARRAY_WRITEABLE)); + PyDict_SetItemString(d, "UPDATEIFCOPY", PyInt_FromLong(NPY_ARRAY_UPDATEIFCOPY)); - PyDict_SetItemString(d, "BEHAVED", PyInt_FromLong(NPY_BEHAVED)); - PyDict_SetItemString(d, "BEHAVED_NS", PyInt_FromLong(NPY_BEHAVED_NS)); - PyDict_SetItemString(d, "CARRAY", PyInt_FromLong(NPY_CARRAY)); - PyDict_SetItemString(d, "FARRAY", PyInt_FromLong(NPY_FARRAY)); - PyDict_SetItemString(d, "CARRAY_RO", PyInt_FromLong(NPY_CARRAY_RO)); - PyDict_SetItemString(d, "FARRAY_RO", PyInt_FromLong(NPY_FARRAY_RO)); - PyDict_SetItemString(d, "DEFAULT", PyInt_FromLong(NPY_DEFAULT)); - PyDict_SetItemString(d, "UPDATE_ALL", PyInt_FromLong(NPY_UPDATE_ALL)); + PyDict_SetItemString(d, "BEHAVED", PyInt_FromLong(NPY_ARRAY_BEHAVED)); + PyDict_SetItemString(d, "BEHAVED_NS", PyInt_FromLong(NPY_ARRAY_BEHAVED_NS)); + PyDict_SetItemString(d, "CARRAY", PyInt_FromLong(NPY_ARRAY_CARRAY)); + PyDict_SetItemString(d, "FARRAY", PyInt_FromLong(NPY_ARRAY_FARRAY)); + PyDict_SetItemString(d, "CARRAY_RO", PyInt_FromLong(NPY_ARRAY_CARRAY_RO)); + PyDict_SetItemString(d, "FARRAY_RO", PyInt_FromLong(NPY_ARRAY_FARRAY_RO)); + PyDict_SetItemString(d, "DEFAULT", PyInt_FromLong(NPY_ARRAY_DEFAULT)); + PyDict_SetItemString(d, "UPDATE_ALL", PyInt_FromLong(NPY_ARRAY_UPDATE_ALL)); if (PyErr_Occurred()) Py_FatalError("can't initialize module wrap"); diff --git a/numpy/f2py/tests/src/regression/inout.f90 b/numpy/f2py/tests/src/regression/inout.f90 new file mode 100644 index 000000000..80cdad90c --- /dev/null +++ b/numpy/f2py/tests/src/regression/inout.f90 @@ -0,0 +1,9 @@ +! Check that intent(in out) translates as intent(inout). +! The separation seems to be a common usage. + subroutine foo(x) + implicit none + real(4), intent(in out) :: x + dimension x(3) + x(1) = x(1) + x(2) + x(3) + return + end diff --git a/numpy/f2py/tests/test_regression.py b/numpy/f2py/tests/test_regression.py new file mode 100644 index 000000000..9bd3f3fe3 --- /dev/null +++ b/numpy/f2py/tests/test_regression.py @@ -0,0 +1,32 @@ +from __future__ import division, absolute_import, print_function + +import os +import math + +import numpy as np +from numpy.testing import dec, assert_raises, assert_equal + +import util + +def _path(*a): + return os.path.join(*((os.path.dirname(__file__),) + a)) + +class TestIntentInOut(util.F2PyTest): + # Check that intent(in out) translates as intent(inout) + sources = [_path('src', 'regression', 'inout.f90')] + + @dec.slow + def test_inout(self): + # non-contiguous should raise error + x = np.arange(6, dtype=np.float32)[::2] + assert_raises(ValueError, self.module.foo, x) + + # check values with contiguous array + x = np.arange(3, dtype=np.float32) + self.module.foo(x) + assert_equal(x, [3, 1, 2]) + + +if __name__ == "__main__": + import nose + nose.runmodule() diff --git a/numpy/fft/fftpack_litemodule.c b/numpy/fft/fftpack_litemodule.c index 95da3194f..9753047dc 100644 --- a/numpy/fft/fftpack_litemodule.c +++ b/numpy/fft/fftpack_litemodule.c @@ -6,11 +6,9 @@ static PyObject *ErrorObject; -/* ----------------------------------------------------- */ +static const char fftpack_cfftf__doc__[] = ""; -static char fftpack_cfftf__doc__[] = ""; - -PyObject * +static PyObject * fftpack_cfftf(PyObject *NPY_UNUSED(self), PyObject *args) { PyObject *op1, *op2; @@ -61,9 +59,9 @@ fail: return NULL; } -static char fftpack_cfftb__doc__[] = ""; +static const char fftpack_cfftb__doc__[] = ""; -PyObject * +static PyObject * fftpack_cfftb(PyObject *NPY_UNUSED(self), PyObject *args) { PyObject *op1, *op2; @@ -114,7 +112,7 @@ fail: return NULL; } -static char fftpack_cffti__doc__[] =""; +static const char fftpack_cffti__doc__[] = ""; static PyObject * fftpack_cffti(PyObject *NPY_UNUSED(self), PyObject *args) @@ -143,9 +141,9 @@ fftpack_cffti(PyObject *NPY_UNUSED(self), PyObject *args) return (PyObject *)op; } -static char fftpack_rfftf__doc__[] =""; +static const char fftpack_rfftf__doc__[] = ""; -PyObject * +static PyObject * fftpack_rfftf(PyObject *NPY_UNUSED(self), PyObject *args) { PyObject *op1, *op2; @@ -187,7 +185,6 @@ fftpack_rfftf(PyObject *NPY_UNUSED(self), PyObject *args) rptr = (double *)PyArray_DATA(ret); dptr = (double *)PyArray_DATA(data); - Py_BEGIN_ALLOW_THREADS; NPY_SIGINT_ON; for (i = 0; i < nrepeats; i++) { @@ -211,10 +208,9 @@ fail: return NULL; } -static char fftpack_rfftb__doc__[] =""; - +static const char fftpack_rfftb__doc__[] = ""; -PyObject * +static PyObject * fftpack_rfftb(PyObject *NPY_UNUSED(self), PyObject *args) { PyObject *op1, *op2; @@ -274,8 +270,7 @@ fail: return NULL; } - -static char fftpack_rffti__doc__[] =""; +static const char fftpack_rffti__doc__[] = ""; static PyObject * fftpack_rffti(PyObject *NPY_UNUSED(self), PyObject *args) @@ -316,11 +311,6 @@ static struct PyMethodDef fftpack_methods[] = { {NULL, NULL, 0, NULL} /* sentinel */ }; - -/* Initialization function for the module (*must* be called initfftpack) */ - -static char fftpack_module_documentation[] = "" ; - #if PY_MAJOR_VERSION >= 3 static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT, @@ -349,6 +339,8 @@ initfftpack_lite(void) #if PY_MAJOR_VERSION >= 3 m = PyModule_Create(&moduledef); #else + static const char fftpack_module_documentation[] = ""; + m = Py_InitModule4("fftpack_lite", fftpack_methods, fftpack_module_documentation, (PyObject*)NULL,PYTHON_API_VERSION); diff --git a/numpy/lib/arraysetops.py b/numpy/lib/arraysetops.py index 5c3b504de..d3b6119f4 100644 --- a/numpy/lib/arraysetops.py +++ b/numpy/lib/arraysetops.py @@ -241,6 +241,11 @@ def intersect1d(ar1, ar2, assume_unique=False): >>> np.intersect1d([1, 3, 4, 3], [3, 1, 2, 1]) array([1, 3]) + To intersect more than two arrays, use functools.reduce: + + >>> from functools import reduce + >>> reduce(np.intersect1d, ([1, 3, 4, 3], [3, 1, 2, 1], [6, 3, 4, 2])) + array([3]) """ if not assume_unique: # Might be faster than unique( intersect1d( ar1, ar2 ) )? @@ -421,6 +426,11 @@ def union1d(ar1, ar2): >>> np.union1d([-1, 0, 1], [-2, 0, 2]) array([-2, -1, 0, 1, 2]) + To find the union of more than two arrays, use functools.reduce: + + >>> from functools import reduce + >>> reduce(np.union1d, ([1, 3, 4, 3], [3, 1, 2, 1], [6, 3, 4, 2])) + array([1, 2, 3, 4, 6]) """ return unique(np.concatenate((ar1, ar2))) diff --git a/numpy/lib/format.py b/numpy/lib/format.py index 98743b6ad..b93f86ca3 100644 --- a/numpy/lib/format.py +++ b/numpy/lib/format.py @@ -141,7 +141,7 @@ import sys import io import warnings from numpy.lib.utils import safe_eval -from numpy.compat import asbytes, isfileobj, long, basestring +from numpy.compat import asbytes, asstr, isfileobj, long, basestring if sys.version_info[0] >= 3: import pickle @@ -298,7 +298,8 @@ def _write_array_header(fp, d, version=None): # can take advantage of our premature optimization. current_header_len = MAGIC_LEN + 2 + len(header) + 1 # 1 for the newline topad = 16 - (current_header_len % 16) - header = asbytes(header + ' '*topad + '\n') + header = header + ' '*topad + '\n' + header = asbytes(_filter_header(header)) if len(header) >= (256*256) and version == (1, 0): raise ValueError("header does not fit inside %s bytes required by the" @@ -410,6 +411,45 @@ def read_array_header_2_0(fp): """ _read_array_header(fp, version=(2, 0)) + +def _filter_header(s): + """Clean up 'L' in npz header ints. + + Cleans up the 'L' in strings representing integers. Needed to allow npz + headers produced in Python2 to be read in Python3. + + Parameters + ---------- + s : byte string + Npy file header. + + Returns + ------- + header : str + Cleaned up header. + + """ + import tokenize + if sys.version_info[0] >= 3: + from io import StringIO + else: + from StringIO import StringIO + + tokens = [] + last_token_was_number = False + for token in tokenize.generate_tokens(StringIO(asstr(s)).read): + token_type = token[0] + token_string = token[1] + if (last_token_was_number and + token_type == tokenize.NAME and + token_string == "L"): + continue + else: + tokens.append(token) + last_token_was_number = (token_type == tokenize.NUMBER) + return tokenize.untokenize(tokens) + + def _read_array_header(fp, version): """ see read_array_header_1_0 @@ -434,6 +474,7 @@ def _read_array_header(fp, version): # "shape" : tuple of int # "fortran_order" : bool # "descr" : dtype.descr + header = _filter_header(header) try: d = safe_eval(header) except SyntaxError as e: diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 8b384bfaa..36ce94bad 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -515,8 +515,7 @@ def average(a, axis=None, weights=None, returned=False): scl = avg.dtype.type(a.size/avg.size) else: a = a + 0.0 - wgt = np.array(weights, dtype=a.dtype, copy=0) - + wgt = np.asarray(weights) # Sanity checks if a.shape != wgt.shape: if axis is None: @@ -533,7 +532,7 @@ def average(a, axis=None, weights=None, returned=False): # setup wgt to broadcast along axis wgt = np.array(wgt, copy=0, ndmin=a.ndim).swapaxes(-1, axis) - scl = wgt.sum(axis=axis) + scl = wgt.sum(axis=axis, dtype=np.result_type(a.dtype, wgt.dtype)) if (scl == 0.0).any(): raise ZeroDivisionError( "Weights sum to zero, can't be normalized") @@ -883,28 +882,33 @@ def copy(a, order='K'): # Basic operations -def gradient(f, *varargs): +def gradient(f, *varargs, **kwargs): """ Return the gradient of an N-dimensional array. - + The gradient is computed using second order accurate central differences - in the interior and second order accurate one-sides (forward or backwards) - differences at the boundaries. The returned gradient hence has the same - shape as the input array. + in the interior and either first differences or second order accurate + one-sides (forward or backwards) differences at the boundaries. The + returned gradient hence has the same shape as the input array. Parameters ---------- f : array_like - An N-dimensional array containing samples of a scalar function. - `*varargs` : scalars - 0, 1, or N scalars specifying the sample distances in each direction, - that is: `dx`, `dy`, `dz`, ... The default distance is 1. + An N-dimensional array containing samples of a scalar function. + varargs : list of scalar, optional + N scalars specifying the sample distances for each dimension, + i.e. `dx`, `dy`, `dz`, ... Default distance: 1. + edge_order : {1, 2}, optional + Gradient is calculated using N\ :sup:`th` order accurate differences + at the boundaries. Default: 1. + + .. versionadded:: 1.9.1 Returns ------- gradient : ndarray - N arrays of the same shape as `f` giving the derivative of `f` with - respect to each dimension. + N arrays of the same shape as `f` giving the derivative of `f` with + respect to each dimension. Examples -------- @@ -916,15 +920,14 @@ def gradient(f, *varargs): >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float)) [array([[ 2., 2., -1.], - [ 2., 2., -1.]]), - array([[ 1. , 2.5, 4. ], - [ 1. , 1. , 1. ]])] + [ 2., 2., -1.]]), array([[ 1. , 2.5, 4. ], + [ 1. , 1. , 1. ]])] - >>> x = np.array([0,1,2,3,4]) - >>> dx = gradient(x) + >>> x = np.array([0, 1, 2, 3, 4]) + >>> dx = np.gradient(x) >>> y = x**2 - >>> gradient(y,dx) - array([0., 2., 4., 6., 8.]) + >>> np.gradient(y, dx, edge_order=2) + array([-0., 2., 4., 6., 8.]) """ f = np.asanyarray(f) N = len(f.shape) # number of dimensions @@ -939,6 +942,13 @@ def gradient(f, *varargs): raise SyntaxError( "invalid number of arguments") + edge_order = kwargs.pop('edge_order', 1) + if kwargs: + raise TypeError('"{}" are not valid keyword arguments.'.format( + '", "'.join(kwargs.keys()))) + if edge_order > 2: + raise ValueError("'edge_order' greater than 2 not supported") + # use central differences on interior and one-sided differences on the # endpoints. This preserves second order-accuracy over the full domain. @@ -978,7 +988,7 @@ def gradient(f, *varargs): "at least two elements are required.") # Numerical differentiation: 1st order edges, 2nd order interior - if y.shape[axis] == 2: + if y.shape[axis] == 2 or edge_order == 1: # Use first order differences for time data out = np.empty_like(y, dtype=otype) @@ -1026,7 +1036,8 @@ def gradient(f, *varargs): out[slice1] = (3.0*y[slice2] - 4.0*y[slice3] + y[slice4])/2.0 # divide by step size - outvals.append(out / dx[axis]) + out /= dx[axis] + outvals.append(out) # reset the slice object in this dimension to ":" slice1[axis] = slice(None) @@ -1102,7 +1113,7 @@ def diff(a, n=1, axis=-1): return a[slice1]-a[slice2] -def interp(x, xp, fp, left=None, right=None): +def interp(x, xp, fp, left=None, right=None, period=None): """ One-dimensional linear interpolation. @@ -1115,7 +1126,9 @@ def interp(x, xp, fp, left=None, right=None): The x-coordinates of the interpolated values. xp : 1-D sequence of floats - The x-coordinates of the data points, must be increasing. + The x-coordinates of the data points, must be increasing if argument + `period` is not specified. Otherwise, `xp` is internally sorted after + normalizing the periodic boundaries with ``xp = xp % period``. fp : 1-D sequence of floats The y-coordinates of the data points, same length as `xp`. @@ -1126,6 +1139,12 @@ def interp(x, xp, fp, left=None, right=None): right : float, optional Value to return for `x > xp[-1]`, default is `fp[-1]`. + period : None or float, optional + .. versionadded:: 1.10.0 + A period for the x-coordinates. This parameter allows the proper + interpolation of angular x-coordinates. Parameters `left` and `right` + are ignored if `period` is specified. + Returns ------- y : {float, ndarray} @@ -1135,6 +1154,8 @@ def interp(x, xp, fp, left=None, right=None): ------ ValueError If `xp` and `fp` have different length + If `xp` or `fp` are not 1-D sequences + If `period == 0` Notes ----- @@ -1144,7 +1165,6 @@ def interp(x, xp, fp, left=None, right=None): np.all(np.diff(xp) > 0) - Examples -------- >>> xp = [1, 2, 3] @@ -1170,13 +1190,51 @@ def interp(x, xp, fp, left=None, right=None): [<matplotlib.lines.Line2D object at 0x...>] >>> plt.show() + Interpolation with periodic x-coordinates: + + >>> x = [-180, -170, -185, 185, -10, -5, 0, 365] + >>> xp = [190, -190, 350, -350] + >>> fp = [5, 10, 3, 4] + >>> np.interp(x, xp, fp, period=360) + array([7.5, 5., 8.75, 6.25, 3., 3.25, 3.5, 3.75]) + """ - if isinstance(x, (float, int, number)): - return compiled_interp([x], xp, fp, left, right).item() - elif isinstance(x, np.ndarray) and x.ndim == 0: - return compiled_interp([x], xp, fp, left, right).item() + if period is None: + if isinstance(x, (float, int, number)): + return compiled_interp([x], xp, fp, left, right).item() + elif isinstance(x, np.ndarray) and x.ndim == 0: + return compiled_interp([x], xp, fp, left, right).item() + else: + return compiled_interp(x, xp, fp, left, right) else: - return compiled_interp(x, xp, fp, left, right) + if period == 0: + raise ValueError("period must be a non-zero value") + period = abs(period) + left = None + right = None + return_array = True + if isinstance(x, (float, int, number)): + return_array = False + x = [x] + x = np.asarray(x, dtype=np.float64) + xp = np.asarray(xp, dtype=np.float64) + fp = np.asarray(fp, dtype=np.float64) + if xp.ndim != 1 or fp.ndim != 1: + raise ValueError("Data points must be 1-D sequences") + if xp.shape[0] != fp.shape[0]: + raise ValueError("fp and xp are not of the same length") + # normalizing periodic boundaries + x = x % period + xp = xp % period + asort_xp = np.argsort(xp) + xp = xp[asort_xp] + fp = fp[asort_xp] + xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period)) + fp = np.concatenate((fp[-1:], fp, fp[0:1])) + if return_array: + return compiled_interp(x, xp, fp, left, right) + else: + return compiled_interp(x, xp, fp, left, right).item() def angle(z, deg=0): diff --git a/numpy/lib/npyio.py b/numpy/lib/npyio.py index 641203f34..cbef1a6e2 100644 --- a/numpy/lib/npyio.py +++ b/numpy/lib/npyio.py @@ -122,6 +122,14 @@ class BagObj(object): return object.__getattribute__(self, '_obj')[key] except KeyError: raise AttributeError(key) + + def __dir__(self): + """ + Enables dir(bagobj) to list the files in an NpzFile. + + This also enables tab-completion in an interpreter or IPython. + """ + return object.__getattribute__(self, '_obj').keys() def zipfile_factory(*args, **kwargs): diff --git a/numpy/lib/tests/data/python3.npy b/numpy/lib/tests/data/python3.npy Binary files differnew file mode 100644 index 000000000..7c6997dd6 --- /dev/null +++ b/numpy/lib/tests/data/python3.npy diff --git a/numpy/lib/tests/data/win64python2.npy b/numpy/lib/tests/data/win64python2.npy Binary files differnew file mode 100644 index 000000000..d9bc36af7 --- /dev/null +++ b/numpy/lib/tests/data/win64python2.npy diff --git a/numpy/lib/tests/test_format.py b/numpy/lib/tests/test_format.py index c09386789..ee77386bc 100644 --- a/numpy/lib/tests/test_format.py +++ b/numpy/lib/tests/test_format.py @@ -524,6 +524,16 @@ def test_compressed_roundtrip(): assert_array_equal(arr, arr1) +def test_python2_python3_interoperability(): + if sys.version_info[0] >= 3: + fname = 'win64python2.npy' + else: + fname = 'python3.npy' + path = os.path.join(os.path.dirname(__file__), 'data', fname) + data = np.load(path) + assert_array_equal(data, np.ones(2)) + + def test_version_2_0(): f = BytesIO() # requires more than 2 byte for header diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index c800f8347..80faf85a6 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -124,6 +124,11 @@ class TestAverage(TestCase): assert_array_equal(average(y1, weights=w2, axis=1), desired) assert_equal(average(y1, weights=w2), 5.) + y3 = rand(5).astype(np.float32) + w3 = rand(5).astype(np.float64) + + assert_(np.average(y3, weights=w3).dtype == np.result_type(y3, w3)) + def test_returned(self): y = np.array([[1, 2, 3], [4, 5, 6]]) @@ -526,8 +531,18 @@ class TestGradient(TestCase): def test_masked(self): # Make sure that gradient supports subclasses like masked arrays - x = np.ma.array([[1, 1], [3, 4]]) - assert_equal(type(gradient(x)[0]), type(x)) + x = np.ma.array([[1, 1], [3, 4]], + mask=[[False, False], [False, False]]) + out = gradient(x)[0] + assert_equal(type(out), type(x)) + # And make sure that the output and input don't have aliased mask + # arrays + assert_(x.mask is not out.mask) + # Also check that edge_order=2 doesn't alter the original mask + x2 = np.ma.arange(5) + x2[2] = np.ma.masked + np.gradient(x2, edge_order=2) + assert_array_equal(x2.mask, [False, False, True, False, False]) def test_datetime64(self): # Make sure gradient() can handle special types like datetime64 @@ -536,7 +551,7 @@ class TestGradient(TestCase): '1910-10-12', '1910-12-12', '1912-12-12'], dtype='datetime64[D]') dx = np.array( - [-7, -3, 0, 31, 61, 396, 1066], + [-5, -3, 0, 31, 61, 396, 731], dtype='timedelta64[D]') assert_array_equal(gradient(x), dx) assert_(dx.dtype == np.dtype('timedelta64[D]')) @@ -547,7 +562,7 @@ class TestGradient(TestCase): [-5, -3, 10, 12, 61, 321, 300], dtype='timedelta64[D]') dx = np.array( - [-3, 7, 7, 25, 154, 119, -161], + [2, 7, 7, 25, 154, 119, -21], dtype='timedelta64[D]') assert_array_equal(gradient(x), dx) assert_(dx.dtype == np.dtype('timedelta64[D]')) @@ -561,7 +576,7 @@ class TestGradient(TestCase): dx = x[1] - x[0] y = 2 * x ** 3 + 4 * x ** 2 + 2 * x analytical = 6 * x ** 2 + 8 * x + 2 - num_error = np.abs((np.gradient(y, dx) / analytical) - 1) + num_error = np.abs((np.gradient(y, dx, edge_order=2) / analytical) - 1) assert_(np.all(num_error < 0.03) == True) @@ -1604,6 +1619,9 @@ class TestInterp(TestCase): def test_exceptions(self): assert_raises(ValueError, interp, 0, [], []) assert_raises(ValueError, interp, 0, [0], [1, 2]) + assert_raises(ValueError, interp, 0, [0, 1], [1, 2], period=0) + assert_raises(ValueError, interp, 0, [], [], period=360) + assert_raises(ValueError, interp, 0, [0], [1, 2], period=360) def test_basic(self): x = np.linspace(0, 1, 5) @@ -1644,6 +1662,16 @@ class TestInterp(TestCase): fp = np.sin(xp) assert_almost_equal(np.interp(np.pi, xp, fp), 0.0) + def test_period(self): + x = [-180, -170, -185, 185, -10, -5, 0, 365] + xp = [190, -190, 350, -350] + fp = [5, 10, 3, 4] + y = [7.5, 5., 8.75, 6.25, 3., 3.25, 3.5, 3.75] + assert_almost_equal(np.interp(x, xp, fp, period=360), y) + x = np.array(x, order='F').reshape(2, -1) + y = np.array(y, order='C').reshape(2, -1) + assert_almost_equal(np.interp(x, xp, fp, period=360), y) + def compare_results(res, desired): for i in range(len(desired)): diff --git a/numpy/lib/tests/test_io.py b/numpy/lib/tests/test_io.py index 4038d6a7f..68b2018cd 100644 --- a/numpy/lib/tests/test_io.py +++ b/numpy/lib/tests/test_io.py @@ -216,6 +216,17 @@ class TestSavezLoad(RoundtripTest, TestCase): l = np.load(c) assert_equal(a, l['file_a']) assert_equal(b, l['file_b']) + + def test_BagObj(self): + a = np.array([[1, 2], [3, 4]], float) + b = np.array([[1 + 2j, 2 + 7j], [3 - 6j, 4 + 12j]], complex) + c = BytesIO() + np.savez(c, file_a=a, file_b=b) + c.seek(0) + l = np.load(c) + assert_equal(sorted(dir(l.f)), ['file_a','file_b']) + assert_equal(a, l.f.file_a) + assert_equal(b, l.f.file_b) def test_savez_filename_clashes(self): # Test that issue #852 is fixed diff --git a/numpy/lib/tests/test_nanfunctions.py b/numpy/lib/tests/test_nanfunctions.py index 3da6b5149..35ae86c20 100644 --- a/numpy/lib/tests/test_nanfunctions.py +++ b/numpy/lib/tests/test_nanfunctions.py @@ -645,6 +645,22 @@ class TestNanFunctions_Median(TestCase): assert_raises(IndexError, np.nanmedian, d, axis=(0, 4)) assert_raises(ValueError, np.nanmedian, d, axis=(1, 1)) + def test_float_special(self): + with warnings.catch_warnings(record=True): + warnings.simplefilter('ignore', RuntimeWarning) + a = np.array([[np.inf, np.nan], [np.nan, np.nan]]) + assert_equal(np.nanmedian(a, axis=0), [np.inf, np.nan]) + assert_equal(np.nanmedian(a, axis=1), [np.inf, np.nan]) + assert_equal(np.nanmedian(a), np.inf) + + # minimum fill value check + a = np.array([[np.nan, np.nan, np.inf], [np.nan, np.nan, np.inf]]) + assert_equal(np.nanmedian(a, axis=1), np.inf) + + # no mask path + a = np.array([[np.inf, np.inf], [np.inf, np.inf]]) + assert_equal(np.nanmedian(a, axis=1), np.inf) + class TestNanFunctions_Percentile(TestCase): diff --git a/numpy/lib/utils.py b/numpy/lib/utils.py index df0052493..519d0e9b9 100644 --- a/numpy/lib/utils.py +++ b/numpy/lib/utils.py @@ -4,6 +4,7 @@ import os import sys import types import re +import warnings from numpy.core.numerictypes import issubclass_, issubsctype, issubdtype from numpy.core import ndarray, ufunc, asarray @@ -1002,111 +1003,70 @@ class SafeEval(object): This includes strings with lists, dicts and tuples using the abstract syntax tree created by ``compiler.parse``. - For an example of usage, see `safe_eval`. + .. deprecated:: 1.10.0 See Also -------- safe_eval """ + def __init__(self): + warnings.warn("SafeEval is deprecated in 1.10 and will be removed.", + DeprecationWarning) - if sys.version_info[0] < 3: - def visit(self, node, **kw): - cls = node.__class__ - meth = getattr(self, 'visit'+cls.__name__, self.default) - return meth(node, **kw) + def visit(self, node): + cls = node.__class__ + meth = getattr(self, 'visit' + cls.__name__, self.default) + return meth(node) - def default(self, node, **kw): - raise SyntaxError("Unsupported source construct: %s" - % node.__class__) + def default(self, node): + raise SyntaxError("Unsupported source construct: %s" + % node.__class__) - def visitExpression(self, node, **kw): - for child in node.getChildNodes(): - return self.visit(child, **kw) + def visitExpression(self, node): + return self.visit(node.body) - def visitConst(self, node, **kw): - return node.value + def visitNum(self, node): + return node.n - def visitDict(self, node,**kw): - return dict( - [(self.visit(k), self.visit(v)) for k, v in node.items] - ) - - def visitTuple(self, node, **kw): - return tuple([self.visit(i) for i in node.nodes]) - - def visitList(self, node, **kw): - return [self.visit(i) for i in node.nodes] - - def visitUnaryAdd(self, node, **kw): - return +self.visit(node.getChildNodes()[0]) - - def visitUnarySub(self, node, **kw): - return -self.visit(node.getChildNodes()[0]) - - def visitName(self, node, **kw): - if node.name == 'False': - return False - elif node.name == 'True': - return True - elif node.name == 'None': - return None - else: - raise SyntaxError("Unknown name: %s" % node.name) - else: - - def visit(self, node): - cls = node.__class__ - meth = getattr(self, 'visit' + cls.__name__, self.default) - return meth(node) - - def default(self, node): - raise SyntaxError("Unsupported source construct: %s" - % node.__class__) - - def visitExpression(self, node): - return self.visit(node.body) - - def visitNum(self, node): - return node.n + def visitStr(self, node): + return node.s - def visitStr(self, node): - return node.s + def visitBytes(self, node): + return node.s - def visitBytes(self, node): - return node.s + def visitDict(self, node,**kw): + return dict([(self.visit(k), self.visit(v)) + for k, v in zip(node.keys, node.values)]) - def visitDict(self, node,**kw): - return dict([(self.visit(k), self.visit(v)) - for k, v in zip(node.keys, node.values)]) + def visitTuple(self, node): + return tuple([self.visit(i) for i in node.elts]) - def visitTuple(self, node): - return tuple([self.visit(i) for i in node.elts]) + def visitList(self, node): + return [self.visit(i) for i in node.elts] - def visitList(self, node): - return [self.visit(i) for i in node.elts] + def visitUnaryOp(self, node): + import ast + if isinstance(node.op, ast.UAdd): + return +self.visit(node.operand) + elif isinstance(node.op, ast.USub): + return -self.visit(node.operand) + else: + raise SyntaxError("Unknown unary op: %r" % node.op) + + def visitName(self, node): + if node.id == 'False': + return False + elif node.id == 'True': + return True + elif node.id == 'None': + return None + else: + raise SyntaxError("Unknown name: %s" % node.id) - def visitUnaryOp(self, node): - import ast - if isinstance(node.op, ast.UAdd): - return +self.visit(node.operand) - elif isinstance(node.op, ast.USub): - return -self.visit(node.operand) - else: - raise SyntaxError("Unknown unary op: %r" % node.op) - - def visitName(self, node): - if node.id == 'False': - return False - elif node.id == 'True': - return True - elif node.id == 'None': - return None - else: - raise SyntaxError("Unknown name: %s" % node.id) + def visitNameConstant(self, node): + return node.value - def visitNameConstant(self, node): - return node.value def safe_eval(source): """ @@ -1151,26 +1111,8 @@ def safe_eval(source): SyntaxError: Unsupported source construct: compiler.ast.CallFunc """ - # Local imports to speed up numpy's import time. - import warnings - - with warnings.catch_warnings(): - # compiler package is deprecated for 3.x, which is already solved - # here - warnings.simplefilter('ignore', DeprecationWarning) - try: - import compiler - except ImportError: - import ast as compiler - - walker = SafeEval() - try: - ast = compiler.parse(source, mode="eval") - except SyntaxError: - raise - try: - return walker.visit(ast) - except SyntaxError: - raise + # Local import to speed up numpy's import time. + import ast + return ast.literal_eval(source) #----------------------------------------------------------------------------- diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 6b2299fe7..43b2dc4dc 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -14,7 +14,7 @@ from __future__ import division, absolute_import, print_function __all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv', 'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'slogdet', 'det', 'svd', 'eig', 'eigh', 'lstsq', 'norm', 'qr', 'cond', 'matrix_rank', - 'LinAlgError'] + 'LinAlgError', 'multi_dot'] import warnings @@ -23,7 +23,7 @@ from numpy.core import ( csingle, cdouble, inexact, complexfloating, newaxis, ravel, all, Inf, dot, add, multiply, sqrt, maximum, fastCopyAndTranspose, sum, isfinite, size, finfo, errstate, geterrobj, longdouble, rollaxis, amin, amax, product, abs, - broadcast + broadcast, atleast_2d, intp, asanyarray ) from numpy.lib import triu, asfarray from numpy.linalg import lapack_lite, _umath_linalg @@ -1921,7 +1921,7 @@ def _multi_svd_norm(x, row_axis, col_axis, op): return result -def norm(x, ord=None, axis=None): +def norm(x, ord=None, axis=None, keepdims=False): """ Matrix or vector norm. @@ -1942,6 +1942,11 @@ def norm(x, ord=None, axis=None): axes that hold 2-D matrices, and the matrix norms of these matrices are computed. If `axis` is None then either a vector norm (when `x` is 1-D) or a matrix norm (when `x` is 2-D) is returned. + keepdims : bool, optional + .. versionadded:: 1.10.0 + If this is set to True, the axes which are normed over are left in the + result as dimensions with size one. With this option the result will + broadcast correctly against the original `x`. Returns ------- @@ -2053,35 +2058,43 @@ def norm(x, ord=None, axis=None): # Check the default case first and handle it immediately. if ord is None and axis is None: + ndim = x.ndim x = x.ravel(order='K') if isComplexType(x.dtype.type): sqnorm = dot(x.real, x.real) + dot(x.imag, x.imag) else: sqnorm = dot(x, x) - return sqrt(sqnorm) + ret = sqrt(sqnorm) + if keepdims: + ret = ret.reshape(ndim*[1]) + return ret # Normalize the `axis` argument to a tuple. nd = x.ndim if axis is None: axis = tuple(range(nd)) elif not isinstance(axis, tuple): + try: + axis = int(axis) + except: + raise TypeError("'axis' must be None, an integer or a tuple of integers") axis = (axis,) if len(axis) == 1: if ord == Inf: - return abs(x).max(axis=axis) + return abs(x).max(axis=axis, keepdims=keepdims) elif ord == -Inf: - return abs(x).min(axis=axis) + return abs(x).min(axis=axis, keepdims=keepdims) elif ord == 0: # Zero norm - return (x != 0).sum(axis=axis) + return (x != 0).sum(axis=axis, keepdims=keepdims) elif ord == 1: # special case for speedup - return add.reduce(abs(x), axis=axis) + return add.reduce(abs(x), axis=axis, keepdims=keepdims) elif ord is None or ord == 2: # special case for speedup s = (x.conj() * x).real - return sqrt(add.reduce(s, axis=axis)) + return sqrt(add.reduce(s, axis=axis, keepdims=keepdims)) else: try: ord + 1 @@ -2100,7 +2113,7 @@ def norm(x, ord=None, axis=None): # if the type changed, we can safely overwrite absx abs(absx, out=absx) absx **= ord - return add.reduce(absx, axis=axis) ** (1.0 / ord) + return add.reduce(absx, axis=axis, keepdims=keepdims) ** (1.0 / ord) elif len(axis) == 2: row_axis, col_axis = axis if not (-nd <= row_axis < nd and -nd <= col_axis < nd): @@ -2109,28 +2122,217 @@ def norm(x, ord=None, axis=None): if row_axis % nd == col_axis % nd: raise ValueError('Duplicate axes given.') if ord == 2: - return _multi_svd_norm(x, row_axis, col_axis, amax) + ret = _multi_svd_norm(x, row_axis, col_axis, amax) elif ord == -2: - return _multi_svd_norm(x, row_axis, col_axis, amin) + ret = _multi_svd_norm(x, row_axis, col_axis, amin) elif ord == 1: if col_axis > row_axis: col_axis -= 1 - return add.reduce(abs(x), axis=row_axis).max(axis=col_axis) + ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis) elif ord == Inf: if row_axis > col_axis: row_axis -= 1 - return add.reduce(abs(x), axis=col_axis).max(axis=row_axis) + ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis) elif ord == -1: if col_axis > row_axis: col_axis -= 1 - return add.reduce(abs(x), axis=row_axis).min(axis=col_axis) + ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis) elif ord == -Inf: if row_axis > col_axis: row_axis -= 1 - return add.reduce(abs(x), axis=col_axis).min(axis=row_axis) + ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis) elif ord in [None, 'fro', 'f']: - return sqrt(add.reduce((x.conj() * x).real, axis=axis)) + ret = sqrt(add.reduce((x.conj() * x).real, axis=axis)) else: raise ValueError("Invalid norm order for matrices.") + if keepdims: + ret_shape = list(x.shape) + ret_shape[axis[0]] = 1 + ret_shape[axis[1]] = 1 + ret = ret.reshape(ret_shape) + return ret else: raise ValueError("Improper number of dimensions to norm.") + + +# multi_dot + +def multi_dot(arrays): + """ + Compute the dot product of two or more arrays in a single function call, + while automatically selecting the fastest evaluation order. + + `multi_dot` chains `numpy.dot` and uses an optimal parenthesizations + of the matrices [1]_ [2]_. Depending on the shape of the matrices + this can speed up the multiplication a lot. + + The first and last argument can be 1-D and are treated respectively as + row and column vector. The other arguments must be 2-D. + + Think of `multi_dot` as:: + + def multi_dot(arrays): return functools.reduce(np.dot, arrays) + + + Parameters + ---------- + arrays : sequence of array_like + First and last argument can be 1-D and are treated respectively as + row and column vector, the other arguments must be 2-D. + + Returns + ------- + output : ndarray + Returns the dot product of the supplied arrays. + + See Also + -------- + dot : dot multiplication with two arguments. + + References + ---------- + + .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378 + .. [2] http://en.wikipedia.org/wiki/Matrix_chain_multiplication + + Examples + -------- + `multi_dot` allows you to write:: + + >>> import numpy as np + >>> # Prepare some data + >>> A = np.random.random(10000, 100) + >>> B = np.random.random(100, 1000) + >>> C = np.random.random(1000, 5) + >>> D = np.random.random(5, 333) + >>> # the actual dot multiplication + >>> multi_dot([A, B, C, D]) + + instead of:: + + >>> np.dot(np.dot(np.dot(A, B), C), D) + >>> # or + >>> A.dot(B).dot(C).dot(D) + + + Example: multiplication costs of different parenthesizations + ------------------------------------------------------------ + + The cost for a matrix multiplication can be calculated with the + following function:: + + def cost(A, B): return A.shape[0] * A.shape[1] * B.shape[1] + + Let's assume we have three matrices + :math:`A_{10x100}, B_{100x5}, C_{5x50}$`. + + The costs for the two different parenthesizations are as follows:: + + cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500 + cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000 + + """ + n = len(arrays) + # optimization only makes sense for len(arrays) > 2 + if n < 2: + raise ValueError("Expecting at least two arrays.") + elif n == 2: + return dot(arrays[0], arrays[1]) + + arrays = [asanyarray(a) for a in arrays] + + # save original ndim to reshape the result array into the proper form later + ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim + # Explicitly convert vectors to 2D arrays to keep the logic of the internal + # _multi_dot_* functions as simple as possible. + if arrays[0].ndim == 1: + arrays[0] = atleast_2d(arrays[0]) + if arrays[-1].ndim == 1: + arrays[-1] = atleast_2d(arrays[-1]).T + _assertRank2(*arrays) + + # _multi_dot_three is much faster than _multi_dot_matrix_chain_order + if n == 3: + result = _multi_dot_three(arrays[0], arrays[1], arrays[2]) + else: + order = _multi_dot_matrix_chain_order(arrays) + result = _multi_dot(arrays, order, 0, n - 1) + + # return proper shape + if ndim_first == 1 and ndim_last == 1: + return result[0, 0] # scalar + elif ndim_first == 1 or ndim_last == 1: + return result.ravel() # 1-D + else: + return result + + +def _multi_dot_three(A, B, C): + """ + Find best ordering for three arrays and do the multiplication. + + Doing in manually instead of using dynamic programing is + approximately 15 times faster. + + """ + # cost1 = cost((AB)C) + cost1 = (A.shape[0] * A.shape[1] * B.shape[1] + # (AB) + A.shape[0] * B.shape[1] * C.shape[1]) # (--)C + # cost2 = cost((AB)C) + cost2 = (B.shape[0] * B.shape[1] * C.shape[1] + # (BC) + A.shape[0] * A.shape[1] * C.shape[1]) # A(--) + + if cost1 < cost2: + return dot(dot(A, B), C) + else: + return dot(A, dot(B, C)) + + +def _multi_dot_matrix_chain_order(arrays, return_costs=False): + """ + Return a np.array which encodes the opimal order of mutiplications. + + The optimal order array is then used by `_multi_dot()` to do the + multiplication. + + Also return the cost matrix if `return_costs` is `True` + + The implementation CLOSELY follows Cormen, "Introduction to Algorithms", + Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices. + + cost[i, j] = min([ + cost[prefix] + cost[suffix] + cost_mult(prefix, suffix) + for k in range(i, j)]) + + """ + n = len(arrays) + # p stores the dimensions of the matrices + # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50] + p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]] + # m is a matrix of costs of the subproblems + # m[i,j]: min number of scalar multiplications needed to compute A_{i..j} + m = zeros((n, n), dtype=double) + # s is the actual ordering + # s[i, j] is the value of k at which we split the product A_i..A_j + s = empty((n, n), dtype=intp) + + for l in range(1, n): + for i in range(n - l): + j = i + l + m[i, j] = Inf + for k in range(i, j): + q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1] + if q < m[i, j]: + m[i, j] = q + s[i, j] = k # Note that Cormen uses 1-based index + + return (s, m) if return_costs else s + + +def _multi_dot(arrays, order, i, j): + """Actually do the multiplication with the given order.""" + if i == j: + return arrays[i] + else: + return dot(_multi_dot(arrays, order, i, order[i, j]), + _multi_dot(arrays, order, order[i, j] + 1, j)) diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py index 8edf36aa6..63baa590d 100644 --- a/numpy/linalg/tests/test_linalg.py +++ b/numpy/linalg/tests/test_linalg.py @@ -12,7 +12,8 @@ import numpy as np from numpy import array, single, double, csingle, cdouble, dot, identity from numpy import multiply, atleast_2d, inf, asarray, matrix from numpy import linalg -from numpy.linalg import matrix_power, norm, matrix_rank +from numpy.linalg import matrix_power, norm, matrix_rank, multi_dot +from numpy.linalg.linalg import _multi_dot_matrix_chain_order from numpy.testing import ( assert_, assert_equal, assert_raises, assert_array_equal, assert_almost_equal, assert_allclose, run_module_suite, @@ -207,7 +208,7 @@ for tgt, src in ((GENERALIZED_SQUARE_CASES, SQUARE_CASES), for case in src: if not isinstance(case.a, np.ndarray): continue - + a = np.array([case.a, 2*case.a, 3*case.a]) if case.b is None: b = None @@ -885,6 +886,48 @@ class _TestNorm(object): expected = [norm(B[:,:, k], ord=order) for k in range(B.shape[2])] assert_almost_equal(n, expected) + def test_keepdims(self): + A = np.arange(1,25, dtype=self.dt).reshape(2,3,4) + + allclose_err = 'order {0}, axis = {1}' + shape_err = 'Shape mismatch found {0}, expected {1}, order={2}, axis={3}' + + # check the order=None, axis=None case + expected = norm(A, ord=None, axis=None) + found = norm(A, ord=None, axis=None, keepdims=True) + assert_allclose(np.squeeze(found), expected, + err_msg=allclose_err.format(None,None)) + expected_shape = (1,1,1) + assert_(found.shape == expected_shape, + shape_err.format(found.shape, expected_shape, None, None)) + + # Vector norms. + for order in [None, -1, 0, 1, 2, 3, np.Inf, -np.Inf]: + for k in range(A.ndim): + expected = norm(A, ord=order, axis=k) + found = norm(A, ord=order, axis=k, keepdims=True) + assert_allclose(np.squeeze(found), expected, + err_msg=allclose_err.format(order,k)) + expected_shape = list(A.shape) + expected_shape[k] = 1 + expected_shape = tuple(expected_shape) + assert_(found.shape == expected_shape, + shape_err.format(found.shape, expected_shape, order, k)) + + # Matrix norms. + for order in [None, -2, 2, -1, 1, np.Inf, -np.Inf, 'fro']: + for k in itertools.permutations(range(A.ndim), 2): + expected = norm(A, ord=order, axis=k) + found = norm(A, ord=order, axis=k, keepdims=True) + assert_allclose(np.squeeze(found), expected, + err_msg=allclose_err.format(order,k)) + expected_shape = list(A.shape) + expected_shape[k[0]] = 1 + expected_shape[k[1]] = 1 + expected_shape = tuple(expected_shape) + assert_(found.shape == expected_shape, + shape_err.format(found.shape, expected_shape, order, k)) + def test_bad_args(self): # Check that bad arguments raise the appropriate exceptions. @@ -909,6 +952,7 @@ class _TestNorm(object): assert_raises(ValueError, norm, B, None, (2, 3)) assert_raises(ValueError, norm, B, None, (0, 1, 2)) +class TestNorm_NonSystematic(object): def test_longdouble_norm(self): # Non-regression test: p-norm of longdouble would previously raise # UnboundLocalError. @@ -1149,5 +1193,90 @@ def test_xerbla_override(): raise SkipTest('Numpy xerbla not linked in.') +class TestMultiDot(object): + def test_basic_function_with_three_arguments(self): + # multi_dot with three arguments uses a fast hand coded algorithm to + # determine the optimal order. Therefore test it separately. + A = np.random.random((6, 2)) + B = np.random.random((2, 6)) + C = np.random.random((6, 2)) + + assert_almost_equal(multi_dot([A, B, C]), A.dot(B).dot(C)) + assert_almost_equal(multi_dot([A, B, C]), np.dot(A, np.dot(B, C))) + + def test_basic_function_with_dynamic_programing_optimization(self): + # multi_dot with four or more arguments uses the dynamic programing + # optimization and therefore deserve a separate + A = np.random.random((6, 2)) + B = np.random.random((2, 6)) + C = np.random.random((6, 2)) + D = np.random.random((2, 1)) + assert_almost_equal(multi_dot([A, B, C, D]), A.dot(B).dot(C).dot(D)) + + def test_vector_as_first_argument(self): + # The first argument can be 1-D + A1d = np.random.random(2) # 1-D + B = np.random.random((2, 6)) + C = np.random.random((6, 2)) + D = np.random.random((2, 2)) + + # the result should be 1-D + assert_equal(multi_dot([A1d, B, C, D]).shape, (2,)) + + def test_vector_as_last_argument(self): + # The last argument can be 1-D + A = np.random.random((6, 2)) + B = np.random.random((2, 6)) + C = np.random.random((6, 2)) + D1d = np.random.random(2) # 1-D + + # the result should be 1-D + assert_equal(multi_dot([A, B, C, D1d]).shape, (6,)) + + def test_vector_as_first_and_last_argument(self): + # The first and last arguments can be 1-D + A1d = np.random.random(2) # 1-D + B = np.random.random((2, 6)) + C = np.random.random((6, 2)) + D1d = np.random.random(2) # 1-D + + # the result should be a scalar + assert_equal(multi_dot([A1d, B, C, D1d]).shape, ()) + + def test_dynamic_programming_logic(self): + # Test for the dynamic programming part + # This test is directly taken from Cormen page 376. + arrays = [np.random.random((30, 35)), + np.random.random((35, 15)), + np.random.random((15, 5)), + np.random.random((5, 10)), + np.random.random((10, 20)), + np.random.random((20, 25))] + m_expected = np.array([[0., 15750., 7875., 9375., 11875., 15125.], + [0., 0., 2625., 4375., 7125., 10500.], + [0., 0., 0., 750., 2500., 5375.], + [0., 0., 0., 0., 1000., 3500.], + [0., 0., 0., 0., 0., 5000.], + [0., 0., 0., 0., 0., 0.]]) + s_expected = np.array([[0, 1, 1, 3, 3, 3], + [0, 0, 2, 3, 3, 3], + [0, 0, 0, 3, 3, 3], + [0, 0, 0, 0, 4, 5], + [0, 0, 0, 0, 0, 5], + [0, 0, 0, 0, 0, 0]], dtype=np.int) + s_expected -= 1 # Cormen uses 1-based index, python does not. + + s, m = _multi_dot_matrix_chain_order(arrays, return_costs=True) + + # Only the upper triangular part (without the diagonal) is interesting. + assert_almost_equal(np.triu(s[:-1, 1:]), + np.triu(s_expected[:-1, 1:])) + assert_almost_equal(np.triu(m), np.triu(m_expected)) + + def test_too_few_input_arrays(self): + assert_raises(ValueError, multi_dot, []) + assert_raises(ValueError, multi_dot, [np.random.random((3, 3))]) + + if __name__ == "__main__": run_module_suite() diff --git a/numpy/ma/core.py b/numpy/ma/core.py index 00164b851..34e52d86e 100644 --- a/numpy/ma/core.py +++ b/numpy/ma/core.py @@ -780,6 +780,8 @@ class _DomainSafeDivide: # component of numpy's import time. if self.tolerance is None: self.tolerance = np.finfo(float).tiny + # don't call ma ufuncs from __array_wrap__ which would fail for scalars + a, b = np.asarray(a), np.asarray(b) return umath.absolute(a) * self.tolerance >= umath.absolute(b) @@ -2783,12 +2785,50 @@ class MaskedArray(ndarray): """ # Get main attributes ......... self._update_from(obj) + # We have to decide how to initialize self.mask, based on + # obj.mask. This is very difficult. There might be some + # correspondence between the elements in the array we are being + # created from (= obj) and us. Or... there might not. This method can + # be called in all kinds of places for all kinds of reasons -- could + # be empty_like, could be slicing, could be a ufunc, could be a view, + # ... The numpy subclassing interface simply doesn't give us any way + # to know, which means that at best this method will be based on + # guesswork and heuristics. To make things worse, there isn't even any + # clear consensus about what the desired behavior is. For instance, + # most users think that np.empty_like(marr) -- which goes via this + # method -- should return a masked array with an empty mask (see + # gh-3404 and linked discussions), but others disagree, and they have + # existing code which depends on empty_like returning an array that + # matches the input mask. + # + # Historically our algorithm was: if the template object mask had the + # same *number of elements* as us, then we used *it's mask object + # itself* as our mask, so that writes to us would also write to the + # original array. This is horribly broken in multiple ways. + # + # Now what we do instead is, if the template object mask has the same + # number of elements as us, and we do not have the same base pointer + # as the template object (b/c views like arr[...] should keep the same + # mask), then we make a copy of the template object mask and use + # that. This is also horribly broken but somewhat less so. Maybe. if isinstance(obj, ndarray): - odtype = obj.dtype - if odtype.names: - _mask = getattr(obj, '_mask', make_mask_none(obj.shape, odtype)) + # XX: This looks like a bug -- shouldn't it check self.dtype + # instead? + if obj.dtype.names: + _mask = getattr(obj, '_mask', + make_mask_none(obj.shape, obj.dtype)) else: _mask = getattr(obj, '_mask', nomask) + # If self and obj point to exactly the same data, then probably + # self is a simple view of obj (e.g., self = obj[...]), so they + # should share the same mask. (This isn't 100% reliable, e.g. self + # could be the first row of obj, or have strange strides, but as a + # heuristic it's not bad.) In all other cases, we make a copy of + # the mask, so that future modifications to 'self' do not end up + # side-effecting 'obj' as well. + if (obj.__array_interface__["data"][0] + != self.__array_interface__["data"][0]): + _mask = _mask.copy() else: _mask = nomask self._mask = _mask @@ -5073,7 +5113,11 @@ class MaskedArray(ndarray): return self if fill_value is None: if endwith: - filler = minimum_fill_value(self) + # nan > inf + if np.issubdtype(self.dtype, np.floating): + filler = np.nan + else: + filler = minimum_fill_value(self) else: filler = maximum_fill_value(self) else: @@ -6182,7 +6226,11 @@ def sort(a, axis= -1, kind='quicksort', order=None, endwith=True, fill_value=Non axis = 0 if fill_value is None: if endwith: - filler = minimum_fill_value(a) + # nan > inf + if np.issubdtype(a.dtype, np.floating): + filler = np.nan + else: + filler = minimum_fill_value(a) else: filler = maximum_fill_value(a) else: diff --git a/numpy/ma/extras.py b/numpy/ma/extras.py index a993fd05d..1849df72b 100644 --- a/numpy/ma/extras.py +++ b/numpy/ma/extras.py @@ -671,8 +671,8 @@ def median(a, axis=None, out=None, overwrite_input=False): """ if not hasattr(a, 'mask') or np.count_nonzero(a.mask) == 0: - return masked_array(np.median(a, axis=axis, out=out, - overwrite_input=overwrite_input), copy=False) + return masked_array(np.median(getattr(a, 'data', a), axis=axis, + out=out, overwrite_input=overwrite_input), copy=False) if overwrite_input: if axis is None: asorted = a.ravel() @@ -705,7 +705,14 @@ def median(a, axis=None, out=None, overwrite_input=False): low = high else: low[odd] = high[odd] - return np.ma.mean([low, high], axis=0, out=out) + + if np.issubdtype(asorted.dtype, np.inexact): + # avoid inf / x = masked + s = np.ma.sum([low, high], axis=0, out=out) + np.true_divide(s.data, 2., casting='unsafe', out=s.data) + else: + s = np.ma.mean([low, high], axis=0, out=out) + return s #.............................................................................. diff --git a/numpy/ma/tests/test_core.py b/numpy/ma/tests/test_core.py index 34951875d..4ac3465aa 100644 --- a/numpy/ma/tests/test_core.py +++ b/numpy/ma/tests/test_core.py @@ -2226,6 +2226,13 @@ class TestMaskedArrayMethods(TestCase): assert_equal(b.shape, a.shape) assert_equal(b.fill_value, a.fill_value) + # check empty_like mask handling + a = masked_array([1, 2, 3], mask=[False, True, False]) + b = empty_like(a) + assert_(not np.may_share_memory(a.mask, b.mask)) + b = a.view(masked_array) + assert_(np.may_share_memory(a.mask, b.mask)) + def test_put(self): # Tests put. d = arange(5) diff --git a/numpy/ma/tests/test_extras.py b/numpy/ma/tests/test_extras.py index 6ce1dc346..d7bc765a9 100644 --- a/numpy/ma/tests/test_extras.py +++ b/numpy/ma/tests/test_extras.py @@ -504,6 +504,9 @@ class TestApplyOverAxes(TestCase): class TestMedian(TestCase): + def test_pytype(self): + r = np.ma.median([[np.inf, np.inf], [np.inf, np.inf]], axis=-1) + assert_equal(r, np.inf) def test_2d(self): # Tests median w/ 2D diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py index 1fd49d774..1d3bef390 100644 --- a/numpy/polynomial/hermite.py +++ b/numpy/polynomial/hermite.py @@ -1577,13 +1577,13 @@ def hermcompanion(c): n = len(c) - 1 mat = np.zeros((n, n), dtype=c.dtype) - scl = np.hstack((1., np.sqrt(2.*np.arange(1, n)))) - scl = np.multiply.accumulate(scl) + scl = np.hstack((1., 1./np.sqrt(2.*np.arange(n - 1, 0, -1)))) + scl = np.multiply.accumulate(scl)[::-1] top = mat.reshape(-1)[1::n+1] bot = mat.reshape(-1)[n::n+1] top[...] = np.sqrt(.5*np.arange(1, n)) bot[...] = top - mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*.5 + mat[:, -1] -= scl*c[:-1]/(2.0*c[-1]) return mat @@ -1646,6 +1646,49 @@ def hermroots(c): return r +def _normed_hermite_n(x, n): + """ + Evaluate a normalized Hermite polynomial. + + Compute the value of the normalized Hermite polynomial of degree ``n`` + at the points ``x``. + + + Parameters + ---------- + x : ndarray of double. + Points at which to evaluate the function + n : int + Degree of the normalized Hermite function to be evaluated. + + Returns + ------- + values : ndarray + The shape of the return value is described above. + + Notes + ----- + .. versionadded:: 1.10.0 + + This function is needed for finding the Gauss points and integration + weights for high degrees. The values of the standard Hermite functions + overflow when n >= 207. + + """ + if n == 0: + return np.ones(x.shape)/np.sqrt(np.sqrt(np.pi)) + + c0 = 0. + c1 = 1./np.sqrt(np.sqrt(np.pi)) + nd = float(n) + for i in range(n - 1): + tmp = c0 + c0 = -c1*np.sqrt((nd - 1.)/nd) + c1 = tmp + c1*x*np.sqrt(2./nd) + nd = nd - 1.0 + return c0 + c1*x*np.sqrt(2) + + def hermgauss(deg): """ Gauss-Hermite quadrature. @@ -1688,22 +1731,21 @@ def hermgauss(deg): # first approximation of roots. We use the fact that the companion # matrix is symmetric in this case in order to obtain better zeros. - c = np.array([0]*deg + [1]) + c = np.array([0]*deg + [1], dtype=np.float64) m = hermcompanion(c) - x = la.eigvals(m) + x = la.eigvalsh(m) x.sort() # improve roots by one application of Newton - dy = hermval(x, c) - df = hermval(x, hermder(c)) + dy = _normed_hermite_n(x, ideg) + df = _normed_hermite_n(x, ideg - 1) * np.sqrt(2*ideg) x -= dy/df # compute the weights. We scale the factor to avoid possible numerical # overflow. - fm = hermval(x, c[1:]) + fm = _normed_hermite_n(x, ideg - 1) fm /= np.abs(fm).max() - df /= np.abs(df).max() - w = 1/(fm * df) + w = 1/(fm * fm) # for Hermite we can also symmetrize w = (w + w[::-1])/2 diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py index 6e33dc0bc..fce13a84e 100644 --- a/numpy/polynomial/hermite_e.py +++ b/numpy/polynomial/hermite_e.py @@ -1575,13 +1575,13 @@ def hermecompanion(c): n = len(c) - 1 mat = np.zeros((n, n), dtype=c.dtype) - scl = np.hstack((1., np.sqrt(np.arange(1, n)))) - scl = np.multiply.accumulate(scl) + scl = np.hstack((1., 1./np.sqrt(np.arange(n - 1, 0, -1)))) + scl = np.multiply.accumulate(scl)[::-1] top = mat.reshape(-1)[1::n+1] bot = mat.reshape(-1)[n::n+1] top[...] = np.sqrt(np.arange(1, n)) bot[...] = top - mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1]) + mat[:, -1] -= scl*c[:-1]/c[-1] return mat @@ -1644,6 +1644,49 @@ def hermeroots(c): return r +def _normed_hermite_e_n(x, n): + """ + Evaluate a normalized HermiteE polynomial. + + Compute the value of the normalized HermiteE polynomial of degree ``n`` + at the points ``x``. + + + Parameters + ---------- + x : ndarray of double. + Points at which to evaluate the function + n : int + Degree of the normalized HermiteE function to be evaluated. + + Returns + ------- + values : ndarray + The shape of the return value is described above. + + Notes + ----- + .. versionadded:: 1.10.0 + + This function is needed for finding the Gauss points and integration + weights for high degrees. The values of the standard HermiteE functions + overflow when n >= 207. + + """ + if n == 0: + return np.ones(x.shape)/np.sqrt(np.sqrt(2*np.pi)) + + c0 = 0. + c1 = 1./np.sqrt(np.sqrt(2*np.pi)) + nd = float(n) + for i in range(n - 1): + tmp = c0 + c0 = -c1*np.sqrt((nd - 1.)/nd) + c1 = tmp + c1*x*np.sqrt(1./nd) + nd = nd - 1.0 + return c0 + c1*x + + def hermegauss(deg): """ Gauss-HermiteE quadrature. @@ -1688,20 +1731,19 @@ def hermegauss(deg): # matrix is symmetric in this case in order to obtain better zeros. c = np.array([0]*deg + [1]) m = hermecompanion(c) - x = la.eigvals(m) + x = la.eigvalsh(m) x.sort() # improve roots by one application of Newton - dy = hermeval(x, c) - df = hermeval(x, hermeder(c)) + dy = _normed_hermite_e_n(x, ideg) + df = _normed_hermite_e_n(x, ideg - 1) * np.sqrt(ideg) x -= dy/df # compute the weights. We scale the factor to avoid possible numerical # overflow. - fm = hermeval(x, c[1:]) + fm = _normed_hermite_e_n(x, ideg - 1) fm /= np.abs(fm).max() - df /= np.abs(df).max() - w = 1/(fm * df) + w = 1/(fm * fm) # for Hermite_e we can also symmetrize w = (w + w[::-1])/2 diff --git a/numpy/random/mtrand/mtrand.pyx b/numpy/random/mtrand/mtrand.pyx index b089d7742..c03666527 100644 --- a/numpy/random/mtrand/mtrand.pyx +++ b/numpy/random/mtrand/mtrand.pyx @@ -2919,7 +2919,7 @@ cdef class RandomState: Plot Gaussian for comparison: - >>> g = (1/(scale * np.sqrt(2 * np.pi)) * + >>> g = (1/(scale * np.sqrt(2 * np.pi)) * ... np.exp(-(x - loc)**2 / (2 * scale**2))) >>> plt.plot(x,g) @@ -4220,8 +4220,8 @@ cdef class RandomState: mean : 1-D array_like, of length N Mean of the N-dimensional distribution. cov : 2-D array_like, of shape (N, N) - Covariance matrix of the distribution. Must be symmetric and - positive-semidefinite for "physically meaningful" results. + Covariance matrix of the distribution. It must be symmetric and + positive-semidefinite for proper sampling. size : int or tuple of ints, optional Given a shape of, for example, ``(m,n,k)``, ``m*n*k`` samples are generated, and packed in an `m`-by-`n`-by-`k` arrangement. Because @@ -4268,7 +4268,9 @@ cdef class RandomState: >>> x,y = np.random.multivariate_normal(mean,cov,5000).T >>> plt.plot(x,y,'x'); plt.axis('equal'); plt.show() - Note that the covariance matrix must be non-negative definite. + Note that the covariance matrix must be positive semidefinite (a.k.a. + nonnegative-definite). Otherwise, the behavior of this method is + undefined and backwards compatibility is not guaranteed. References ---------- diff --git a/numpy/testing/__init__.py b/numpy/testing/__init__.py index 258cbe928..dcc02ad57 100644 --- a/numpy/testing/__init__.py +++ b/numpy/testing/__init__.py @@ -10,7 +10,6 @@ from __future__ import division, absolute_import, print_function from unittest import TestCase from . import decorators as dec +from .nosetester import run_module_suite, NoseTester as Tester from .utils import * -from .nosetester import NoseTester as Tester -from .nosetester import run_module_suite test = Tester().test diff --git a/numpy/testing/utils.py b/numpy/testing/utils.py index dc5929753..3b20f9238 100644 --- a/numpy/testing/utils.py +++ b/numpy/testing/utils.py @@ -159,7 +159,7 @@ else: """ Return memory usage of running python. [Not implemented]""" raise NotImplementedError -if os.name=='nt' and sys.version[:3] > '2.3': +if os.name=='nt': # Code "stolen" from enthought/debug/memusage.py def GetPerformanceAttributes(object, counter, instance = None, inum=-1, format = None, machine=None): |