diff options
| author | Travis Oliphant <oliphant@enthought.com> | 2009-08-28 15:36:42 +0000 |
|---|---|---|
| committer | Travis Oliphant <oliphant@enthought.com> | 2009-08-28 15:36:42 +0000 |
| commit | 2b01ee6b966b9ca298247e6717e3a5be16a92970 (patch) | |
| tree | 6bbb8ee8eebdfe2ef3eb26f13994193b313c6fe7 /numpy/core | |
| parent | c2191bc97da8a0879cec8d3e9a7a93fe9e66fcd8 (diff) | |
| parent | fddd4b9c3b8f18ba7cf386f766b70ec3328b1c69 (diff) | |
| download | numpy-2b01ee6b966b9ca298247e6717e3a5be16a92970.tar.gz | |
Re-base the date-time branch back to the trunk.
Diffstat (limited to 'numpy/core')
47 files changed, 4438 insertions, 1486 deletions
diff --git a/numpy/core/SConscript b/numpy/core/SConscript index f09b17618..fffbee5af 100644 --- a/numpy/core/SConscript +++ b/numpy/core/SConscript @@ -120,10 +120,9 @@ numpyconfig_sym.append(('NPY_NO_SMP', nosmp)) # Check whether we can use C99 printing formats #---------------------------------------------- if config.CheckDeclaration(('PRIdPTR'), includes = '#include <inttypes.h>'): - usec99 = 1 + numpyconfig_sym.append(('DEFINE_NPY_USE_C99_FORMATS', '#define NPY_USE_C99_FORMATS 1')) else: - usec99 = 0 -numpyconfig_sym.append(('USE_C99_FORMATS', usec99)) + numpyconfig_sym.append(('DEFINE_NPY_USE_C99_FORMATS', '')) #---------------------- # Checking the mathlib @@ -187,9 +186,10 @@ for f in ["isnan", "isinf", "signbit", "isfinite"]: """ st = config.CheckDeclaration(f, includes=includes) if st: - numpyconfig_sym.append(('DEFINE_NPY_HAVE_DECL_%s' % f.upper(), - '#define NPY_HAVE_DECL_%s' % f.upper())) - + numpyconfig_sym.append(('DEFINE_NPY_HAVE_DECL_%s' % f.upper(), + '#define NPY_HAVE_DECL_%s' % f.upper())) + else: + numpyconfig_sym.append(('DEFINE_NPY_HAVE_DECL_%s' % f.upper(), '')) inline = config.CheckInline() config.Define('inline', inline) @@ -198,7 +198,9 @@ numpyconfig_sym.append(('NPY_INLINE', inline)) if ENABLE_SEPARATE_COMPILATION: config.Define("ENABLE_SEPARATE_COMPILATION", 1) - numpyconfig_sym.append(('NPY_ENABLE_SEPARATE_COMPILATION', 1)) + numpyconfig_sym.append(('DEFINE_NPY_ENABLE_SEPARATE_COMPILATION', '#define NPY_ENABLE_SEPARATE_COMPILATION 1')) +else: + numpyconfig_sym.append(('DEFINE_NPY_ENABLE_SEPARATE_COMPILATION', '')) # Checking for visibility macro def visibility_define(): @@ -262,6 +264,10 @@ write_info(env) # Build #========== +# List of headers which need to be "installed " into the build directory for +# proper in-place build support +generated_headers = [] + #--------------------------------------- # Generate the public configuration file #--------------------------------------- @@ -272,7 +278,9 @@ for key, value in numpyconfig_sym: env['SUBST_DICT'] = config_dict include_dir = 'include/numpy' -env.SubstInFile(pjoin(include_dir, 'numpyconfig.h'), pjoin(include_dir, 'numpyconfig.h.in')) +target = env.SubstInFile(pjoin(include_dir, 'numpyconfig.h'), + pjoin(include_dir, 'numpyconfig.h.in')) +generated_headers.append(target[0]) env['CONFIG_H_GEN'] = numpyconfig_sym @@ -298,26 +306,45 @@ umathmodule_src = env.GenerateFromTemplate(pjoin('src', 'umath', 'umathmodule.c.src')) umath_tests_src = env.GenerateFromTemplate(pjoin('src', 'umath', 'umath_tests.c.src')) +multiarray_tests_src = env.GenerateFromTemplate(pjoin('src', 'multiarray', + 'multiarray_tests.c.src')) scalarmathmodule_src = env.GenerateFromTemplate( pjoin('src', 'scalarmathmodule.c.src')) umath = env.GenerateUmath('__umath_generated', pjoin('code_generators', 'generate_umath.py')) -multiarray_api = env.GenerateMultiarrayApi('multiarray_api', +multiarray_api = env.GenerateMultiarrayApi('include/numpy/multiarray_api', [ pjoin('code_generators', 'numpy_api_order.txt')]) +generated_headers.append(multiarray_api[0]) -ufunc_api = env.GenerateUfuncApi('ufunc_api', +ufunc_api = env.GenerateUfuncApi('include/numpy/ufunc_api', pjoin('code_generators', 'ufunc_api_order.txt')) +generated_headers.append(ufunc_api[0]) -env.Prepend(CPPPATH = ['include', '.']) +# include/numpy is added for compatibility reasons with distutils: this is +# needed for __multiarray_api.c and __ufunc_api.c included from multiarray and +# ufunc. +env.Prepend(CPPPATH = ['include', '.', 'include/numpy']) # npymath core lib -npymath_src = env.GenerateFromTemplate(pjoin('src', 'npy_math.c.src')) -env.DistutilsStaticExtLibrary("npymath", npymath_src) +npymath_src = env.GenerateFromTemplate(pjoin('src', 'npymath', 'npy_math.c.src')) +env.DistutilsInstalledStaticExtLibrary("npymath", npymath_src, install_dir='lib') env.Prepend(LIBS=["npymath"]) env.Prepend(LIBPATH=["."]) +subst_dict = {'@prefix@': '$distutils_install_prefix', + '@sep@': repr(os.path.sep)} +npymath_ini = env.SubstInFile(pjoin('lib', 'npy-pkg-config', 'npymath.ini'), + 'npymath.ini.in', SUBST_DICT=subst_dict) + +subst_dict = {'@posix_mathlib@': " ".join(['-l%s' % l for l in mlib]), + '@msvc_mathlib@': " ".join(['%s.mlib' % l for l in mlib])} +mlib_ini = env.SubstInFile(pjoin('lib', 'npy-pkg-config', 'mlib.ini'), + 'mlib.ini.in', SUBST_DICT=subst_dict) +env.Install('$distutils_installdir/lib/npy-pkg-config', mlib_ini) +env.Install('$distutils_installdir/lib/npy-pkg-config', npymath_ini) + #----------------- # Build multiarray #----------------- @@ -353,6 +380,7 @@ if ENABLE_SEPARATE_COMPILATION: else: multiarray_src = [pjoin('src', 'multiarray', 'multiarraymodule_onefile.c')] multiarray = env.DistutilsPythonExtension('multiarray', source = multiarray_src) +env.DistutilsPythonExtension('multiarray_tests', source=multiarray_tests_src) #------------------ # Build sort module @@ -395,3 +423,7 @@ if build_blasdot: dotblas_o = env.PythonObject('_dotblas', source = dotblas_src) env.Depends(dotblas_o, pjoin("blasdot", "cblas.h")) dotblas = env.DistutilsPythonExtension('_dotblas', dotblas_o) + +# "Install" the header in the build directory, so that in-place build works +for h in generated_headers: + env.Install(pjoin('$distutils_installdir', include_dir), h) diff --git a/numpy/core/arrayprint.py b/numpy/core/arrayprint.py index 6d3c52990..c8bc9438a 100644 --- a/numpy/core/arrayprint.py +++ b/numpy/core/arrayprint.py @@ -56,10 +56,14 @@ def set_printoptions(precision=None, threshold=None, edgeitems=None, suppress : bool, optional Whether or not suppress printing of small floating point values using scientific notation (default False). - nanstr : string, optional - String representation of floating point not-a-number (default nan). - infstr : string, optional - String representation of floating point infinity (default inf). + nanstr : str, optional + String representation of floating point not-a-number (default NaN). + infstr : str, optional + String representation of floating point infinity (default Inf). + + See Also + -------- + get_printoptions, set_string_function Examples -------- @@ -79,12 +83,9 @@ def set_printoptions(precision=None, threshold=None, edgeitems=None, >>> eps = np.finfo(float).eps >>> x = np.arange(4.) - >>> x**2 - (x + eps)**2 array([ -4.9304e-32, -4.4409e-16, 0.0000e+00, 0.0000e+00]) - >>> np.set_printoptions(suppress=True) - >>> x**2 - (x + eps)**2 array([-0., -0., 0., 0.]) diff --git a/numpy/core/code_generators/generate_numpy_api.py b/numpy/core/code_generators/generate_numpy_api.py index 509048471..69f8c2026 100644 --- a/numpy/core/code_generators/generate_numpy_api.py +++ b/numpy/core/code_generators/generate_numpy_api.py @@ -28,6 +28,7 @@ extern NPY_NO_EXPORT PyTypeObject PyArrayFlags_Type; extern NPY_NO_EXPORT PyTypeObject PyArrayIter_Type; extern NPY_NO_EXPORT PyTypeObject PyArrayMapIter_Type; extern NPY_NO_EXPORT PyTypeObject PyArrayMultiIter_Type; +extern NPY_NO_EXPORT PyTypeObject PyArrayNeighborhoodIter_Type; extern NPY_NO_EXPORT PyTypeObject PyBoolArrType_Type; extern NPY_NO_EXPORT PyBoolScalarObject _PyArrayScalar_BoolValues[2]; #else @@ -39,6 +40,7 @@ NPY_NO_EXPORT PyTypeObject PyArrayFlags_Type; NPY_NO_EXPORT PyTypeObject PyArrayIter_Type; NPY_NO_EXPORT PyTypeObject PyArrayMapIter_Type; NPY_NO_EXPORT PyTypeObject PyArrayMultiIter_Type; +NPY_NO_EXPORT PyTypeObject PyArrayNeighborhoodIter_Type; NPY_NO_EXPORT PyTypeObject PyBoolArrType_Type; NPY_NO_EXPORT PyBoolScalarObject _PyArrayScalar_BoolValues[2]; #endif @@ -113,13 +115,13 @@ _import_array(void) PyErr_Format(PyExc_RuntimeError, "FATAL: module compiled as unknown endian"); return -1; } -#ifdef NPY_BIG_ENDIAN +#if NPY_BYTE_ORDER ==NPY_BIG_ENDIAN if (st != NPY_CPU_BIG) { PyErr_Format(PyExc_RuntimeError, "FATAL: module compiled as "\ "big endian, but detected different endianness at runtime"); return -1; } -#elif defined(NPY_LITTLE_ENDIAN) +#elif NPY_BYTE_ORDER == NPY_LITTLE_ENDIAN if (st != NPY_CPU_LITTLE) { PyErr_Format(PyExc_RuntimeError, "FATAL: module compiled as "\ "little endian, but detected different endianness at runtime"); diff --git a/numpy/core/code_generators/generate_umath.py b/numpy/core/code_generators/generate_umath.py index 15e427f29..60a8b5a30 100644 --- a/numpy/core/code_generators/generate_umath.py +++ b/numpy/core/code_generators/generate_umath.py @@ -650,6 +650,11 @@ defdict = { docstrings.get('numpy.core.umath.signbit'), TD(flts, out='?'), ), +'copysign' : + Ufunc(2, 1, None, + docstrings.get('numpy.core.umath.copysign'), + TD(flts), + ), 'modf' : Ufunc(1, 2, None, docstrings.get('numpy.core.umath.modf'), diff --git a/numpy/core/code_generators/numpy_api_order.txt b/numpy/core/code_generators/numpy_api_order.txt index 6b7272459..72f8d5c82 100644 --- a/numpy/core/code_generators/numpy_api_order.txt +++ b/numpy/core/code_generators/numpy_api_order.txt @@ -173,9 +173,10 @@ PyArray_CompareString PyArray_MultiIterFromObjects PyArray_GetEndianness PyArray_GetNDArrayCFeatureVersion -PyArray_Acorrelate +PyArray_Correlate2 +PyArray_NeighborhoodIterNew PyArray_SetDatetimeParseFunction PyArray_DatetimeToDatetimeStruct PyArray_TimedeltaToTimedeltaStruct PyArray_DatetimeStructToDatetime -PyArray_TimedeltaStructToTimedelta
\ No newline at end of file +PyArray_TimedeltaStructToTimedelt diff --git a/numpy/core/code_generators/ufunc_docstrings.py b/numpy/core/code_generators/ufunc_docstrings.py index e880f21a2..7744bb4bf 100644 --- a/numpy/core/code_generators/ufunc_docstrings.py +++ b/numpy/core/code_generators/ufunc_docstrings.py @@ -93,6 +93,9 @@ add_newdoc('numpy.core.umath', 'arccos', `x`-coordinate on the unit circle. For real arguments, the domain is [-1, 1]. + out : ndarray, optional + Array to store results in. + Returns ------- angle : ndarray @@ -153,6 +156,8 @@ add_newdoc('numpy.core.umath', 'arccosh', ---------- x : array_like Input array. + out : ndarray, optional + Array of the same shape as `x`. Returns ------- @@ -715,16 +720,44 @@ add_newdoc('numpy.core.umath', 'cos', ---------- x : array_like Input array in radians. + out : ndarray, optional + Output array of same shape as `x`. Returns ------- - out : ndarray - Output array of same shape as `x`. + y : ndarray + The corresponding cosine values. + + Raises + ------ + ValueError: invalid return array shape + if `out` is provided and `out.shape` != `x.shape` (See Examples) + + Notes + ----- + If `out` is provided, the function writes the result into it, + and returns a reference to `out`. (See Examples) + + References + ---------- + M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions. + New York, NY: Dover, 1972. Examples -------- >>> np.cos(np.array([0, np.pi/2, np.pi])) array([ 1.00000000e+00, 6.12303177e-17, -1.00000000e+00]) + >>> + >>> # Example of providing the optional output parameter + >>> out2 = np.cos([0.1], out1) + >>> out2 is out1 + True + >>> + >>> # Example of ValueError due to provision of shape mis-matched `out` + >>> np.cos(np.zeros((3,3)),np.zeros((2,2))) + Traceback (most recent call last): + File "<stdin>", line 1, in <module> + ValueError: invalid return array shape """) @@ -903,7 +936,7 @@ add_newdoc('numpy.core.umath', 'equal', add_newdoc('numpy.core.umath', 'exp', """ - Calculate the exponential of the elements in the input array. + Calculate the exponential of all elements in the input array. Parameters ---------- @@ -913,7 +946,12 @@ add_newdoc('numpy.core.umath', 'exp', Returns ------- out : ndarray - Element-wise exponential of `x`. + Output array, element-wise exponential of `x`. + + See Also + -------- + expm1 : Calculate ``exp(x) - 1`` for all elements in the array. + exp2 : Calculate ``2**x`` for all elements in the array. Notes ----- @@ -968,20 +1006,34 @@ add_newdoc('numpy.core.umath', 'exp2', x : array_like Input values. + out : ndarray, optional + \tArray to insert results into. + Returns ------- out : ndarray Element-wise 2 to the power `x`. + See Also + -------- + exp : calculate x**p. + Notes ----- .. versionadded:: 1.3.0 + + + Examples + -------- + >>> np.exp2([2,3]) + array([4,9]) + """) add_newdoc('numpy.core.umath', 'expm1', """ - Compute ``exp(x) - 1`` for all elements in the array. + Calculate ``exp(x) - 1`` for all elements in the array. Parameters ---------- @@ -1621,7 +1673,7 @@ add_newdoc('numpy.core.umath', 'log', add_newdoc('numpy.core.umath', 'log10', """ - Compute the logarithm in base 10 element-wise. + Return the base 10 logarithm of the input array, element-wise. Parameters ---------- @@ -1631,7 +1683,8 @@ add_newdoc('numpy.core.umath', 'log10', Returns ------- y : ndarray - Base-10 logarithm of `x`. + The logarithm to the base 10 of `x`, element-wise. NaNs are + returned where x is negative. Notes ----- @@ -1656,7 +1709,7 @@ add_newdoc('numpy.core.umath', 'log10', Examples -------- - >>> np.log10([1.e-15,-3.]) + >>> np.log10([1e-15, -3.]) array([-15., NaN]) """) @@ -1687,71 +1740,91 @@ add_newdoc('numpy.core.umath', 'log2', add_newdoc('numpy.core.umath', 'logaddexp', """ - Logarithm of `exp(x) + exp(y)`. + Logarithm of the sum of exponentiations of the inputs. - This function is useful in statistics where the calculated probabilities of - events may be so small as to excede the range of normal floating point - numbers. In such cases the logarithm of the calculated probability is - stored. This function allows adding probabilities stored in such a fashion. + Calculates ``log(exp(x1) + exp(x2))``. This function is useful in + statistics where the calculated probabilities of events may be so small + as to exceed the range of normal floating point numbers. In such cases + the logarithm of the calculated probability is stored. This function + allows adding probabilities stored in such a fashion. Parameters ---------- - x : array_like - Input values. - y : array_like + x1, x2 : array_like Input values. - Returns ------- result : ndarray - Logarithm of `exp(x) + exp(y)`. + Logarithm of ``exp(x1) + exp(x2)``. See Also -------- - logaddexp2 + logaddexp2: Logarithm of the sum of exponentiations of inputs in base-2. Notes ----- .. versionadded:: 1.3.0 + Examples + -------- + >>> prob1 = np.log(1e-50) + >>> prob2 = np.log(2.5e-50) + >>> prob12 = np.logaddexp(prob1, prob2) + >>> prob12 + -113.87649168120691 + >>> np.exp(prob12) + 3.5000000000000057e-50 + """) add_newdoc('numpy.core.umath', 'logaddexp2', """ - Base-2 Logarithm of `2**x + 2**y`. + Logarithm of the sum of exponentiations of the inputs in base-2. - This function is useful in machine learning when the calculated probabilities of - events may be so small as to excede the range of normal floating point - numbers. In such cases the base-2 logarithm of the calculated probability - can be used instead. This function allows adding probabilities stored in such a fashion. + Calculates ``log2(2**x1 + 2**x2)``. This function is useful in machine + learning when the calculated probabilities of events may be so small + as to exceed the range of normal floating point numbers. In such cases + the base-2 logarithm of the calculated probability can be used instead. + This function allows adding probabilities stored in such a fashion. Parameters ---------- - x : array_like - Input values. - y : array_like + x1, x2 : array_like Input values. - + out : ndarray, optional + Array to store results in. Returns ------- result : ndarray - Base-2 logarithm of `2**x + 2**y`. + Base-2 logarithm of ``2**x1 + 2**x2``. See Also -------- - logaddexp + logaddexp: Logarithm of the sum of exponentiations of the inputs. Notes ----- .. versionadded:: 1.3.0 + Examples + -------- + >>> prob1 = np.log2(1e-50) + >>> prob2 = np.log2(2.5e-50) + >>> prob12 = np.logaddexp2(prob1, prob2) + >>> prob1, prob2, prob12 + (-166.09640474436813, -164.77447664948076, -164.28904982231052) + >>> 2**prob12 + 3.4999999999999914e-50 + """) add_newdoc('numpy.core.umath', 'log1p', """ - `log(1 + x)` in base `e`, elementwise. + Return the natural logarithm of one plus the input array, element-wise. + + Calculates ``log(1 + x)``. Parameters ---------- @@ -1761,7 +1834,11 @@ add_newdoc('numpy.core.umath', 'log1p', Returns ------- y : ndarray - Natural logarithm of `1 + x`, elementwise. + Natural logarithm of `1 + x`, element-wise. + + See Also + -------- + expm1 : ``exp(x) - 1``, the inverse of `log1p`. Notes ----- @@ -2022,8 +2099,6 @@ add_newdoc('numpy.core.umath', 'minimum', add_newdoc('numpy.core.umath', 'fmax', """ - fmax(x1, x2[, out]) - Element-wise maximum of array elements. Compare two arrays and returns a new array containing the element-wise @@ -2132,7 +2207,7 @@ add_newdoc('numpy.core.umath', 'fmin', add_newdoc('numpy.core.umath', 'modf', """ - Return the fractional and integral part of a number. + Return the fractional and integral parts of an array, element-wise. The fractional and integral parts are negative if the given number is negative. @@ -2140,7 +2215,7 @@ add_newdoc('numpy.core.umath', 'modf', Parameters ---------- x : array_like - Input number. + Input array. Returns ------- @@ -2149,33 +2224,37 @@ add_newdoc('numpy.core.umath', 'modf', y2 : ndarray Integral part of `x`. + Notes + ----- + For integer input the return values are floats. + Examples -------- - >>> np.modf(2.5) - (0.5, 2.0) - >>> np.modf(-.4) - (-0.40000000000000002, -0.0) + >>> np.modf([0, 3.5]) + (array([ 0. , 0.5]), array([ 0., 3.])) + >>> np.modf(-0.5) + (-0.5, -0) """) add_newdoc('numpy.core.umath', 'multiply', """ - Multiply arguments elementwise. + Multiply arguments element-wise. Parameters ---------- x1, x2 : array_like - The arrays to be multiplied. + Input arrays to be multiplied. Returns ------- y : ndarray - The product of `x1` and `x2`, elementwise. Returns a scalar if + The product of `x1` and `x2`, element-wise. Returns a scalar if both `x1` and `x2` are scalars. Notes ----- - Equivalent to `x1` * `x2` in terms of array-broadcasting. + Equivalent to `x1` * `x2` in terms of array broadcasting. Examples -------- @@ -2353,23 +2432,34 @@ add_newdoc('numpy.core.umath', 'deg2rad', add_newdoc('numpy.core.umath', 'reciprocal', """ - Return element-wise reciprocal. + Return the reciprocal of the argument, element-wise. + + Calculates ``1/x``. Parameters ---------- x : array_like - Input value. + Input array. Returns ------- y : ndarray - Return value. + Return array. + + Notes + ----- + .. note:: + This function is not designed to work with integers. + + For integer arguments with absolute value larger than 1 the result is + always zero because of the way Python handles integer division. + For integer zero the result is an overflow. Examples -------- - >>> reciprocal(2.) + >>> np.reciprocal(2.) 0.5 - >>> reciprocal([1, 2., 3.33]) + >>> np.reciprocal([1, 2., 3.33]) array([ 1. , 0.5 , 0.3003003]) """) @@ -2378,7 +2468,7 @@ add_newdoc('numpy.core.umath', 'remainder', """ Returns element-wise remainder of division. - Computes `x1 - floor(x1/x2)*x2`. + Computes ``x1 - floor(x1/x2)*x2``. Parameters ---------- @@ -2390,22 +2480,23 @@ add_newdoc('numpy.core.umath', 'remainder', Returns ------- y : ndarray - The remainder of the quotient `x1/x2`, element-wise. Returns a scalar + The remainder of the quotient ``x1/x2``, element-wise. Returns a scalar if both `x1` and `x2` are scalars. See Also -------- - divide - floor + divide, floor Notes ----- - Returns 0 when `x2` is 0. + Returns 0 when `x2` is 0 and both `x1` and `x2` are (arrays of) integers. Examples -------- - >>> np.remainder([4,7],[2,3]) + >>> np.remainder([4,7], [2,3]) array([0, 1]) + >>> np.remainder(np.arange(7), 5) + array([0, 1, 2, 3, 4, 0, 1]) """) @@ -2523,6 +2614,33 @@ add_newdoc('numpy.core.umath', 'signbit', """) +add_newdoc('numpy.core.umath', 'copysign', + """ + Change the sign of x to that of y element-wise. + + Parameters + ---------- + x: array_like + Values to change the sign of. + y: array_like + The sign of y is copied to x. + + Returns + ------- + out : array_like + values of x with the sign of y + + Examples + -------- + >>> np.copysign(1.3, -1) + -1.3 + >>> 1/np.copysign(0, 1) + inf + >>> 1/np.copysign(0, -1) + -inf + + """) + add_newdoc('numpy.core.umath', 'sin', """ Trigonometric sine, element-wise. @@ -2590,11 +2708,50 @@ add_newdoc('numpy.core.umath', 'sinh', ---------- x : array_like Input array. + out : ndarray, optional + Output array of same shape as `x`. Returns ------- - out : ndarray - Output array of same shape as `x`. + y : ndarray + The corresponding hyperbolic sine values. + + Raises + ------ + ValueError: invalid return array shape + if `out` is provided and `out.shape` != `x.shape` (See Examples) + + Notes + ----- + If `out` is provided, the function writes the result into it, + and returns a reference to `out`. (See Examples) + + References + ---------- + M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions. + New York, NY: Dover, 1972, pg. 83. + + Examples + -------- + >>> import numpy as np + >>> np.sinh(0) + 0.0 + >>> np.sinh(np.pi*1j/2) + 1j + >>> np.sinh(np.pi*1j) + 1.2246063538223773e-016j (exact value is 0) + >>> # Discrepancy due to vagaries of floating point arithmetic. + >>> + >>> # Example of providing the optional output parameter + >>> out2 = np.sinh([0.1], out1) + >>> out2 is out1 + True + >>> + >>> # Example of ValueError due to provision of shape mis-matched `out` + >>> np.sinh(np.zeros((3,3)),np.zeros((2,2))) + Traceback (most recent call last): + File "<stdin>", line 1, in <module> + ValueError: invalid return array shape """) @@ -2667,13 +2824,12 @@ add_newdoc('numpy.core.umath', 'square', add_newdoc('numpy.core.umath', 'subtract', """ - Subtract arguments element-wise. + Subtract arguments, element-wise. Parameters ---------- x1, x2 : array_like - The arrays to be subtracted from each other. If type is 'array_like' - the `x1` and `x2` shapes must be identical. + The arrays to be subtracted from each other. Returns ------- @@ -2683,7 +2839,7 @@ add_newdoc('numpy.core.umath', 'subtract', Notes ----- - Equivalent to `x1` - `x2` in terms of array-broadcasting. + Equivalent to ``x1 - x2`` in terms of array broadcasting. Examples -------- @@ -2703,39 +2859,107 @@ add_newdoc('numpy.core.umath', 'tan', """ Compute tangent element-wise. + Equivalent to ``np.sin(x)/np.cos(x)`` element-wise. + Parameters ---------- x : array_like - Angles in radians. + Input array. + out : ndarray, optional + Output array of same shape as `x`. Returns ------- y : ndarray The corresponding tangent values. + Raises + ------ + ValueError: invalid return array shape + if `out` is provided and `out.shape` != `x.shape` (See Examples) + + Notes + ----- + If `out` is provided, the function writes the result into it, + and returns a reference to `out`. (See Examples) + + References + ---------- + M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions. + New York, NY: Dover, 1972. Examples -------- >>> from math import pi >>> np.tan(np.array([-pi,pi/2,pi])) array([ 1.22460635e-16, 1.63317787e+16, -1.22460635e-16]) + >>> + >>> # Example of providing the optional output parameter illustrating + >>> # that what is returned is a reference to said parameter + >>> out2 = np.cos([0.1], out1) + >>> out2 is out1 + True + >>> + >>> # Example of ValueError due to provision of shape mis-matched `out` + >>> np.cos(np.zeros((3,3)),np.zeros((2,2))) + Traceback (most recent call last): + File "<stdin>", line 1, in <module> + ValueError: invalid return array shape """) add_newdoc('numpy.core.umath', 'tanh', """ - Hyperbolic tangent element-wise. + Compute hyperbolic tangent element-wise. + + Equivalent to ``np.sinh(x)/np.cosh(x)`` or + ``-1j * np.tan(1j*x)``. Parameters ---------- x : array_like Input array. + out : ndarray, optional + Output array of same shape as `x`. Returns ------- y : ndarray The corresponding hyperbolic tangent values. + Raises + ------ + ValueError: invalid return array shape + if `out` is provided and `out.shape` != `x.shape` (See Examples) + + Notes + ----- + If `out` is provided, the function writes the result into it, + and returns a reference to `out`. (See Examples) + + References + ---------- + M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions. + New York, NY: Dover, 1972, pg. 83. + + Examples + -------- + >>> import numpy as np + >>> np.tanh((0, np.pi*1j, np.pi*1j/2)) + array([ 0. +0.00000000e+00j, 0. -1.22460635e-16j, 0. +1.63317787e+16j]) + >>> + >>> # Example of providing the optional output parameter illustrating + >>> # that what is returned is a reference to said parameter + >>> out2 = np.tanh([0.1], out1) + >>> out2 is out1 + True + >>> + >>> # Example of ValueError due to provision of shape mis-matched `out` + >>> np.tanh(np.zeros((3,3)),np.zeros((2,2))) + Traceback (most recent call last): + File "<stdin>", line 1, in <module> + ValueError: invalid return array shape + """) add_newdoc('numpy.core.umath', 'true_divide', diff --git a/numpy/core/defmatrix.py b/numpy/core/defmatrix.py index d1636e8b5..354e40060 100644 --- a/numpy/core/defmatrix.py +++ b/numpy/core/defmatrix.py @@ -2,7 +2,7 @@ __all__ = ['matrix', 'bmat', 'mat', 'asmatrix'] import sys import numeric as N -from numeric import concatenate, isscalar, binary_repr, identity +from numeric import concatenate, isscalar, binary_repr, identity, asanyarray from numerictypes import issubdtype # make translation table @@ -115,6 +115,7 @@ def matrix_power(M,n): [ 0, -1]]) """ + M = asanyarray(M) if len(M.shape) != 2 or M.shape[0] != M.shape[1]: raise ValueError("input must be a square array") if not issubdtype(type(n),int): @@ -490,6 +491,26 @@ class matrix(N.ndarray): return N.ndarray.prod(self, axis, dtype, out)._align(axis) def any(self, axis=None, out=None): + """ + Test whether any array element along a given axis evaluates to True. + + Refer to `numpy.any` for full documentation. + + Parameters + ---------- + axis: int, optional + Axis along which logical OR is performed + out: ndarray, optional + Output to existing array instead of creating new one, must have + same shape as expected output + + Returns + ------- + any : bool, ndarray + Returns a single bool if `axis` is ``None``; otherwise, + returns `ndarray` + + """ return N.ndarray.any(self, axis, out)._align(axis) def all(self, axis=None, out=None): diff --git a/numpy/core/fromnumeric.py b/numpy/core/fromnumeric.py index 99b837ba2..f7f584d3d 100644 --- a/numpy/core/fromnumeric.py +++ b/numpy/core/fromnumeric.py @@ -256,15 +256,14 @@ def repeat(a, repeats, axis=None): def put(a, ind, v, mode='raise'): """ - Changes specific elements of one array by replacing from another array. + Replaces specified elements of an array with given values. - The indexing works on the flattened target array, `put` is roughly + The indexing works on the flattened target array. `put` is roughly equivalent to: :: - for i, val in zip(ind, v): - x.flat[i] = val + a.flat[ind] = v Parameters ---------- @@ -292,14 +291,14 @@ def put(a, ind, v, mode='raise'): Examples -------- - >>> x = np.arange(5) - >>> np.put(x, [0, 2], [-44, -55]) - >>> x + >>> a = np.arange(5) + >>> np.put(a, [0, 2], [-44, -55]) + >>> a array([-44, 1, -55, 3, 4]) - >>> x = np.arange(5) - >>> np.put(x, 22, -5, mode='clip') - >>> x + >>> a = np.arange(5) + >>> np.put(a, 22, -5, mode='clip') + >>> a array([ 0, 1, 2, 3, -5]) """ @@ -450,6 +449,22 @@ def sort(a, axis=-1, kind='quicksort', order=None): the last axis is faster and uses less space than sorting along any other axis. + The sort order for complex numbers is lexicographic. If both the real + and imaginary parts are non-nan then the order is determined by the + real parts except when they are equal, in which case the order is + determined by the imaginary parts. + + Previous to numpy 1.4.0 sorting real and complex arrays containing nan + values led to undefined behaviour. In numpy versions >= 1.4.0 nan + values are sorted to the end. The extended sort order is: + + Real: [R, nan] + Complex: [R + Rj, R + nanj, nan + Rj, nan + nanj] + + where R is a non-nan real value. Complex values with the same nan + placements are sorted according to the non-nan part if it exists. + Non-nan values are sorted as before. + Examples -------- >>> a = np.array([[1,4],[3,1]]) @@ -529,6 +544,9 @@ def argsort(a, axis=-1, kind='quicksort', order=None): ----- See `sort` for notes on the different sorting algorithms. + As of Numpy 1.4.0 argsort works with real/complex arrays containing + nan values. The enhanced sort order is documented in the numpy.sort. + Examples -------- One dimensional array: @@ -665,6 +683,9 @@ def searchsorted(a, v, side='left'): ----- Binary search is used to find the required insertion points. + As of Numpy 1.4.0 searchsorted works with real/complex arrays containing + nan values. The enhanced sort order is documented in the numpy.sort. + Examples -------- >>> np.searchsorted([1,2,3,4,5], 3) @@ -1086,10 +1107,14 @@ def compress(condition, a, axis=None, out=None): """ Return selected slices of an array along given axis. + When working along a given axis, a slice along that axis is returned in + `output` for each index where `condition` evaluates to True. When + working on a 1-D array, `compress` is equivalent to `extract`. + Parameters ---------- - condition : array_like - Boolean 1-D array selecting which entries to return. If len(condition) + condition : 1-D array of bools + Array that selects which entries to return. If len(condition) is less than the size of `a` along the given axis, then output is truncated to the length of the condition array. a : array_like @@ -1109,18 +1134,31 @@ def compress(condition, a, axis=None, out=None): See Also -------- - ndarray.compress: Equivalent method. + take, choose, diag, diagonal, select + ndarray.compress : Equivalent method. Examples -------- - >>> a = np.array([[1, 2], [3, 4]]) + >>> a = np.array([[1, 2], [3, 4], [5, 6]]) + >>> a + array([[1, 2], + [3, 4], + [5, 6]]) >>> np.compress([0, 1], a, axis=0) array([[3, 4]]) - >>> np.compress([1], a, axis=1) - array([[1], - [3]]) - >>> np.compress([0,1,1], a) - array([2, 3]) + >>> np.compress([False, True, True], a, axis=0) + array([[3, 4], + [5, 6]]) + >>> np.compress([False, True], a, axis=1) + array([[2], + [4], + [6]]) + + Working on the flattened array does not return slices along an axis but + selects elements. + + >>> np.compress([False, True], a) + array([2]) """ try: @@ -1306,6 +1344,8 @@ def any(a,axis=None, out=None): """ Test whether any array element along a given axis evaluates to True. + Returns single boolean unless `axis` is not ``None`` + Parameters ---------- a : array_like @@ -1322,8 +1362,8 @@ def any(a,axis=None, out=None): Returns ------- - any : ndarray, bool - A new boolean or array is returned unless `out` is + any : bool, ndarray + A new boolean or `ndarray` is returned unless `out` is specified, in which case a reference to `out` is returned. See Also @@ -1429,12 +1469,10 @@ def cumsum (a, axis=None, dtype=None, out=None): Parameters ---------- a : array_like - Input array or object that can be converted to an array. + Input array. axis : int, optional Axis along which the cumulative sum is computed. The default - (`axis` = `None`) is to compute the cumsum over the flattened - array. `axis` may be negative, in which case it counts from the - last to the first axis. + (None) is to compute the cumsum over the flattened array. dtype : dtype, optional Type of the returned array and of the accumulator in which the elements are summed. If `dtype` is not specified, it defaults @@ -1459,11 +1497,12 @@ def cumsum (a, axis=None, dtype=None, out=None): Examples -------- - >>> a = np.array([[1,2,3],[4,5,6]]) + >>> a = np.array([[1,2,3], [4,5,6]]) >>> np.cumsum(a) array([ 1, 3, 6, 10, 15, 21]) - >>> np.cumsum(a,dtype=float) # specifies type of output value(s) + >>> np.cumsum(a, dtype=float) # specifies type of output value(s) array([ 1., 3., 6., 10., 15., 21.]) + >>> np.cumsum(a,axis=0) # sum over rows for each of the 3 columns array([[1, 2, 3], [5, 7, 9]]) @@ -2122,14 +2161,13 @@ def std(a, axis=None, dtype=None, out=None, ddof=0): Returns ------- - standard_deviation : {ndarray, scalar}; see dtype parameter above. + standard_deviation : ndarray, see dtype parameter above. If `out` is None, return a new array containing the standard deviation, otherwise return a reference to the output array. See Also -------- - numpy.var : Variance - numpy.mean : Average + var, mean Notes ----- @@ -2145,7 +2183,7 @@ def std(a, axis=None, dtype=None, out=None, ddof=0): is the square root of the estimated variance, so even with ``ddof=1``, it will not be an unbiased estimate of the standard deviation per se. - Note that, for complex numbers, std takes the absolute + Note that, for complex numbers, `std` takes the absolute value before squaring, so that the result is always real and nonnegative. Examples @@ -2153,9 +2191,9 @@ def std(a, axis=None, dtype=None, out=None, ddof=0): >>> a = np.array([[1, 2], [3, 4]]) >>> np.std(a) 1.1180339887498949 - >>> np.std(a, 0) + >>> np.std(a, axis=0) array([ 1., 1.]) - >>> np.std(a, 1) + >>> np.std(a, axis=1) array([ 0.5, 0.5]) """ diff --git a/numpy/core/include/numpy/_neighborhood_iterator_imp.h b/numpy/core/include/numpy/_neighborhood_iterator_imp.h new file mode 100644 index 000000000..5a73784c1 --- /dev/null +++ b/numpy/core/include/numpy/_neighborhood_iterator_imp.h @@ -0,0 +1,90 @@ +#ifndef _NPY_INCLUDE_NEIGHBORHOOD_IMP +#error You should not include this header directly +#endif +/* + * Private API (here for inline) + */ +static NPY_INLINE int +_PyArrayNeighborhoodIter_IncrCoord(PyArrayNeighborhoodIterObject* iter); + +/* + * Update to next item of the iterator + * + * Note: this simply increment the coordinates vector, last dimension + * incremented first , i.e, for dimension 3 + * ... + * -1, -1, -1 + * -1, -1, 0 + * -1, -1, 1 + * .... + * -1, 0, -1 + * -1, 0, 0 + * .... + * 0, -1, -1 + * 0, -1, 0 + * .... + */ +#define _UPDATE_COORD_ITER(c) \ + wb = iter->coordinates[c] < iter->bounds[c][1]; \ + if (wb) { \ + iter->coordinates[c] += 1; \ + return 0; \ + } \ + else { \ + iter->coordinates[c] = iter->bounds[c][0]; \ + } + +static NPY_INLINE int +_PyArrayNeighborhoodIter_IncrCoord(PyArrayNeighborhoodIterObject* iter) +{ + int i, wb; + + for (i = iter->nd - 1; i >= 0; --i) { + _UPDATE_COORD_ITER(i) + } + + return 0; +} + +/* + * Version optimized for 2d arrays, manual loop unrolling + */ +static NPY_INLINE int +_PyArrayNeighborhoodIter_IncrCoord2D(PyArrayNeighborhoodIterObject* iter) +{ + int wb; + + _UPDATE_COORD_ITER(1) + _UPDATE_COORD_ITER(0) + + return 0; +} +#undef _UPDATE_COORD_ITER + +/* + * Advance to the next neighbour + */ +static NPY_INLINE int +PyArrayNeighborhoodIter_Next(PyArrayNeighborhoodIterObject* iter) +{ + _PyArrayNeighborhoodIter_IncrCoord (iter); + iter->dataptr = iter->translate((PyArrayIterObject*)iter, iter->coordinates); + + return 0; +} + +/* + * Reset functions + */ +static NPY_INLINE int +PyArrayNeighborhoodIter_Reset(PyArrayNeighborhoodIterObject* iter) +{ + int i; + + for (i = 0; i < iter->nd; ++i) { + iter->coordinates[i] = iter->bounds[i][0]; + } + iter->dataptr = iter->translate((PyArrayIterObject*)iter, iter->coordinates); + + return 0; +} diff --git a/numpy/core/include/numpy/ndarrayobject.h b/numpy/core/include/numpy/ndarrayobject.h index 8d0e444a2..34b080732 100644 --- a/numpy/core/include/numpy/ndarrayobject.h +++ b/numpy/core/include/numpy/ndarrayobject.h @@ -735,7 +735,18 @@ typedef int (PyArray_FinalizeFunc)(PyArrayObject *, PyObject *); #define NPY_DISABLE_C_API #endif -typedef struct { +/***************************** + * Basic iterator object + *****************************/ + +/* FWD declaration */ +typedef struct PyArrayIterObject_tag PyArrayIterObject; + +/* type of the function which translates a set of coordinates to a pointer to + * the data */ +typedef char* (*npy_iter_get_dataptr_t)(PyArrayIterObject* iter, npy_intp*); + +struct PyArrayIterObject_tag { PyObject_HEAD int nd_m1; /* number of dimensions - 1 */ npy_intp index, size; @@ -747,7 +758,12 @@ typedef struct { PyArrayObject *ao; char *dataptr; /* pointer to current item*/ npy_bool contiguous; -} PyArrayIterObject; + + npy_intp bounds[NPY_MAXDIMS][2]; + npy_intp limits[NPY_MAXDIMS][2]; + npy_intp limits_sizes[NPY_MAXDIMS]; + npy_iter_get_dataptr_t translate; +} ; /* Iterator API */ @@ -971,6 +987,72 @@ typedef struct { } PyArrayMapIterObject; +enum { + NPY_NEIGHBORHOOD_ITER_ZERO_PADDING, + NPY_NEIGHBORHOOD_ITER_ONE_PADDING, + NPY_NEIGHBORHOOD_ITER_CONSTANT_PADDING, + NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING, + NPY_NEIGHBORHOOD_ITER_MIRROR_PADDING, +}; + +typedef struct { + PyObject_HEAD + + /* + * PyArrayIterObject part: keep this in this exact order + */ + int nd_m1; /* number of dimensions - 1 */ + npy_intp index, size; + npy_intp coordinates[NPY_MAXDIMS];/* N-dimensional loop */ + npy_intp dims_m1[NPY_MAXDIMS]; /* ao->dimensions - 1 */ + npy_intp strides[NPY_MAXDIMS]; /* ao->strides or fake */ + npy_intp backstrides[NPY_MAXDIMS];/* how far to jump back */ + npy_intp factors[NPY_MAXDIMS]; /* shape factors */ + PyArrayObject *ao; + char *dataptr; /* pointer to current item*/ + npy_bool contiguous; + + npy_intp bounds[NPY_MAXDIMS][2]; + npy_intp limits[NPY_MAXDIMS][2]; + npy_intp limits_sizes[NPY_MAXDIMS]; + npy_iter_get_dataptr_t translate; + + /* + * New members + */ + npy_intp nd; + + /* Dimensions is the dimension of the array */ + npy_intp dimensions[NPY_MAXDIMS]; + + /* Neighborhood points coordinates are computed relatively to the point pointed + * by _internal_iter */ + PyArrayIterObject* _internal_iter; + /* To keep a reference to the representation of the constant value for + * constant padding */ + char* constant; + + int mode; +} PyArrayNeighborhoodIterObject; + +/* + * Neighborhood iterator API + */ + +/* General: those work for any mode */ +static NPY_INLINE int +PyArrayNeighborhoodIter_Reset(PyArrayNeighborhoodIterObject* iter); +static NPY_INLINE int +PyArrayNeighborhoodIter_Next(PyArrayNeighborhoodIterObject* iter); +// static NPY_INLINE int +// PyArrayNeighborhoodIter_Next2D(PyArrayNeighborhoodIterObject* iter); + +/* Include inline implementations - functions defined there are not considered + * public API */ +#define _NPY_INCLUDE_NEIGHBORHOOD_IMP +#include "_neighborhood_iterator_imp.h" +#undef _NPY_INCLUDE_NEIGHBORHOOD_IMP + /* The default array type */ #define NPY_DEFAULT_TYPE NPY_DOUBLE @@ -1112,7 +1194,7 @@ typedef struct { #define NPY_SWAP 's' #define NPY_IGNORE '|' -#ifdef NPY_BIG_ENDIAN +#if NPY_BYTE_ORDER == NPY_BIG_ENDIAN #define NPY_NATBYTE NPY_BIG #define NPY_OPPBYTE NPY_LITTLE #else diff --git a/numpy/core/include/numpy/npy_cpu.h b/numpy/core/include/numpy/npy_cpu.h index 106901e55..4d6579189 100644 --- a/numpy/core/include/numpy/npy_cpu.h +++ b/numpy/core/include/numpy/npy_cpu.h @@ -8,7 +8,12 @@ * NPY_CPU_SPARC * NPY_CPU_S390 * NPY_CPU_IA64 - * NPY_CPU_PARISC + * NPY_CPU_HPPA + * NPY_CPU_ALPHA + * NPY_CPU_ARMEL + * NPY_CPU_ARMEB + * NPY_CPU_SH_LE + * NPY_CPU_SH_BE */ #ifndef _NPY_CPUARCH_H_ #define _NPY_CPUARCH_H_ @@ -42,9 +47,18 @@ #define NPY_CPU_S390 #elif defined(__ia64) #define NPY_CPU_IA64 -#elif defined(__parisc__) - /* XXX: Not sure about this one... */ - #define NPY_CPU_PARISC +#elif defined(__hppa__) + #define NPY_CPU_HPPA +#elif defined(__alpha__) + #define NPY_CPU_ALPHA +#elif defined(__arm__) && defined(__ARMEL__) + #define NPY_CPU_ARMEL +#elif defined(__arm__) && defined(__ARMEB__) + #define NPY_CPU_ARMEB +#elif defined(__sh__) && defined(__LITTLE_ENDIAN__) + #define NPY_CPU_SH_LE +#elif defined(__sh__) && defined(__BIG_ENDIAN__) + #define NPY_CPU_SH_BE #else #error Unknown CPU, please report this to numpy maintainers with \ information about your platform (OS, CPU and compiler) diff --git a/numpy/core/include/numpy/npy_endian.h b/numpy/core/include/numpy/npy_endian.h index 0a5c05ef9..f3ae4dbff 100644 --- a/numpy/core/include/numpy/npy_endian.h +++ b/numpy/core/include/numpy/npy_endian.h @@ -9,26 +9,32 @@ #ifdef NPY_HAVE_ENDIAN_H /* Use endian.h if available */ #include <endian.h> + #define NPY_BYTE_ORDER __BYTE_ORDER - #if (__BYTE_ORDER == __LITTLE_ENDIAN) - #define NPY_LITTLE_ENDIAN - #elif (__BYTE_ORDER == __BIG_ENDIAN) - #define NPY_BIG_ENDIAN - #else - #error Unknown machine endianness detected. - #endif + #define NPY_LITTLE_ENDIAN __LITTLE_ENDIAN + #define NPY_BIG_ENDIAN __BIG_ENDIAN #else /* Set endianness info using target CPU */ #include "npy_cpu.h" - #if defined(NPY_CPU_X86) || defined(NPY_CPU_AMD64)\ - || defined(NPY_CPU_IA64) - #define NPY_LITTLE_ENDIAN - #define NPY_BYTE_ORDER 1234 - #elif defined(NPY_CPU_PPC) || defined(NPY_CPU_SPARC)\ - || defined(NPY_CPU_S390) || defined(NPY_CPU_PARISC) || defined(NPY_CPU_PPC64) - #define NPY_BIG_ENDIAN - #define NPY_BYTE_ORDER 4321 + #define NPY_LITTLE_ENDIAN 1234 + #define NPY_BIG_ENDIAN 4321 + + #if defined(NPY_CPU_X86) \ + || defined(NPY_CPU_AMD64) \ + || defined(NPY_CPU_IA64) \ + || defined(NPY_CPU_ALPHA) \ + || defined(NPY_CPU_ARMEL) \ + || defined(NPY_CPU_SH_LE) + #define NPY_BYTE_ORDER NPY_LITTLE_ENDIAN + #elif defined(NPY_CPU_PPC) \ + || defined(NPY_CPU_SPARC) \ + || defined(NPY_CPU_S390) \ + || defined(NPY_CPU_HPPA) \ + || defined(NPY_CPU_PPC64) \ + || defined(NPY_CPU_ARMEB) \ + || defined(NPY_CPU_SH_BE) + #define NPY_BYTE_ORDER NPY_BIG_ENDIAN #else #error Unknown CPU: can not set endianness #endif diff --git a/numpy/core/include/numpy/npy_math.h b/numpy/core/include/numpy/npy_math.h index 2a8ea182b..c219504e4 100644 --- a/numpy/core/include/numpy/npy_math.h +++ b/numpy/core/include/numpy/npy_math.h @@ -52,38 +52,47 @@ NPY_INLINE static float __npy_nzerof(void) /* * Useful constants */ -#define NPY_E 2.7182818284590452353602874713526625 /* e */ -#define NPY_LOG2E 1.4426950408889634073599246810018921 /* log_2 e */ -#define NPY_LOG10E 0.4342944819032518276511289189166051 /* log_10 e */ -#define NPY_LOGE2 0.6931471805599453094172321214581766 /* log_e 2 */ -#define NPY_LOGE10 2.3025850929940456840179914546843642 /* log_e 10 */ -#define NPY_PI 3.1415926535897932384626433832795029 /* pi */ -#define NPY_PI_2 1.5707963267948966192313216916397514 /* pi/2 */ -#define NPY_PI_4 0.7853981633974483096156608458198757 /* pi/4 */ -#define NPY_1_PI 0.3183098861837906715377675267450287 /* 1/pi */ -#define NPY_2_PI 0.6366197723675813430755350534900574 /* 2/pi */ - -#define NPY_Ef 2.7182818284590452353602874713526625F /* e */ -#define NPY_LOG2Ef 1.4426950408889634073599246810018921F /* log_2 e */ -#define NPY_LOG10Ef 0.4342944819032518276511289189166051F /* log_10 e */ -#define NPY_LOGE2f 0.6931471805599453094172321214581766F /* log_e 2 */ -#define NPY_LOGE10f 2.3025850929940456840179914546843642F /* log_e 10 */ -#define NPY_PIf 3.1415926535897932384626433832795029F /* pi */ -#define NPY_PI_2f 1.5707963267948966192313216916397514F /* pi/2 */ -#define NPY_PI_4f 0.7853981633974483096156608458198757F /* pi/4 */ -#define NPY_1_PIf 0.3183098861837906715377675267450287F /* 1/pi */ -#define NPY_2_PIf 0.6366197723675813430755350534900574F /* 2/pi */ - -#define NPY_El 2.7182818284590452353602874713526625L /* e */ -#define NPY_LOG2El 1.4426950408889634073599246810018921L /* log_2 e */ -#define NPY_LOG10El 0.4342944819032518276511289189166051L /* log_10 e */ -#define NPY_LOGE2l 0.6931471805599453094172321214581766L /* log_e 2 */ -#define NPY_LOGE10l 2.3025850929940456840179914546843642L /* log_e 10 */ -#define NPY_PIl 3.1415926535897932384626433832795029L /* pi */ -#define NPY_PI_2l 1.5707963267948966192313216916397514L /* pi/2 */ -#define NPY_PI_4l 0.7853981633974483096156608458198757L /* pi/4 */ -#define NPY_1_PIl 0.3183098861837906715377675267450287L /* 1/pi */ -#define NPY_2_PIl 0.6366197723675813430755350534900574L /* 2/pi */ +#define NPY_E 2.718281828459045235360287471352662498 /* e */ +#define NPY_LOG2E 1.442695040888963407359924681001892137 /* log_2 e */ +#define NPY_LOG10E 0.434294481903251827651128918916605082 /* log_10 e */ +#define NPY_LOGE2 0.693147180559945309417232121458176568 /* log_e 2 */ +#define NPY_LOGE10 2.302585092994045684017991454684364208 /* log_e 10 */ +#define NPY_PI 3.141592653589793238462643383279502884 /* pi */ +#define NPY_PI_2 1.570796326794896619231321691639751442 /* pi/2 */ +#define NPY_PI_4 0.785398163397448309615660845819875721 /* pi/4 */ +#define NPY_1_PI 0.318309886183790671537767526745028724 /* 1/pi */ +#define NPY_2_PI 0.636619772367581343075535053490057448 /* 2/pi */ +#define NPY_EULER 0.577215664901532860606512090082402431 /* Euler constant */ +#define NPY_SQRT2 1.414213562373095048801688724209698079 /* sqrt(2) */ +#define NPY_SQRT1_2 0.707106781186547524400844362104849039 /* 1/sqrt(2) */ + +#define NPY_Ef 2.718281828459045235360287471352662498F /* e */ +#define NPY_LOG2Ef 1.442695040888963407359924681001892137F /* log_2 e */ +#define NPY_LOG10Ef 0.434294481903251827651128918916605082F /* log_10 e */ +#define NPY_LOGE2f 0.693147180559945309417232121458176568F /* log_e 2 */ +#define NPY_LOGE10f 2.302585092994045684017991454684364208F /* log_e 10 */ +#define NPY_PIf 3.141592653589793238462643383279502884F /* pi */ +#define NPY_PI_2f 1.570796326794896619231321691639751442F /* pi/2 */ +#define NPY_PI_4f 0.785398163397448309615660845819875721F /* pi/4 */ +#define NPY_1_PIf 0.318309886183790671537767526745028724F /* 1/pi */ +#define NPY_2_PIf 0.636619772367581343075535053490057448F /* 2/pi */ +#define NPY_EULERf 0.577215664901532860606512090082402431F /* Euler constan*/ +#define NPY_SQRT2f 1.414213562373095048801688724209698079F /* sqrt(2) */ +#define NPY_SQRT1_2f 0.707106781186547524400844362104849039F /* 1/sqrt(2) */ + +#define NPY_El 2.718281828459045235360287471352662498L /* e */ +#define NPY_LOG2El 1.442695040888963407359924681001892137L /* log_2 e */ +#define NPY_LOG10El 0.434294481903251827651128918916605082L /* log_10 e */ +#define NPY_LOGE2l 0.693147180559945309417232121458176568L /* log_e 2 */ +#define NPY_LOGE10l 2.302585092994045684017991454684364208L /* log_e 10 */ +#define NPY_PIl 3.141592653589793238462643383279502884L /* pi */ +#define NPY_PI_2l 1.570796326794896619231321691639751442L /* pi/2 */ +#define NPY_PI_4l 0.785398163397448309615660845819875721L /* pi/4 */ +#define NPY_1_PIl 0.318309886183790671537767526745028724L /* 1/pi */ +#define NPY_2_PIl 0.636619772367581343075535053490057448L /* 2/pi */ +#define NPY_EULERl 0.577215664901532860606512090082402431L /* Euler constan*/ +#define NPY_SQRT2l 1.414213562373095048801688724209698079L /* sqrt(2) */ +#define NPY_SQRT1_2l 0.707106781186547524400844362104849039L /* 1/sqrt(2) */ /* * C99 double math funcs @@ -128,6 +137,8 @@ double npy_atan2(double x, double y); double npy_pow(double x, double y); double npy_modf(double x, double* y); +double npy_copysign(double x, double y); + /* * IEEE 754 fpu handling. Those are guaranteed to be macros */ @@ -198,6 +209,8 @@ float npy_fmodf(float x, float y); float npy_modff(float x, float* y); +float npy_copysignf(float x, float y); + /* * float C99 math functions */ @@ -235,6 +248,8 @@ npy_longdouble npy_fmodl(npy_longdouble x, npy_longdouble y); npy_longdouble npy_modfl(npy_longdouble x, npy_longdouble* y); +npy_longdouble npy_copysignl(npy_longdouble x, npy_longdouble y); + /* * Non standard functions */ diff --git a/numpy/core/include/numpy/numpyconfig.h.in b/numpy/core/include/numpy/numpyconfig.h.in index 9c3f40d17..b3c0d851d 100644 --- a/numpy/core/include/numpy/numpyconfig.h.in +++ b/numpy/core/include/numpy/numpyconfig.h.in @@ -21,10 +21,10 @@ @DEFINE_NPY_SIZEOF_PY_LONG_LONG@ #define NPY_INLINE @NPY_INLINE@ -#define NPY_ENABLE_SEPARATE_COMPILATION @NPY_ENABLE_SEPARATE_COMPILATION@ +@DEFINE_NPY_ENABLE_SEPARATE_COMPILATION@ #define NPY_VISIBILITY_HIDDEN @VISIBILITY_HIDDEN@ -#define NPY_USE_C99_FORMATS @USE_C99_FORMATS@ +@DEFINE_NPY_USE_C99_FORMATS@ #define NPY_ABI_VERSION @NPY_ABI_VERSION@ #define NPY_API_VERSION @NPY_API_VERSION@ diff --git a/numpy/core/memmap.py b/numpy/core/memmap.py index 65f4938fe..2392c3aa7 100644 --- a/numpy/core/memmap.py +++ b/numpy/core/memmap.py @@ -17,7 +17,7 @@ mode_equivalents = { class memmap(ndarray): """ - Create a memory-map to an array stored in a file on disk. + Create a memory-map to an array stored in a *binary* file on disk. Memory-mapped files are used for accessing small segments of large files on disk, without reading the entire file into memory. Numpy's @@ -58,7 +58,7 @@ class memmap(ndarray): order : {'C', 'F'}, optional Specify the order of the ndarray memory layout: C (row-major) or Fortran (column-major). This only has an effect if the shape is - greater than 1-D. The defaullt order is 'C'. + greater than 1-D. The default order is 'C'. Methods ------- diff --git a/numpy/core/mlib.ini.in b/numpy/core/mlib.ini.in new file mode 100644 index 000000000..badaa2ae9 --- /dev/null +++ b/numpy/core/mlib.ini.in @@ -0,0 +1,12 @@ +[meta] +Name = mlib +Description = Math library used with this version of numpy +Version = 1.0 + +[default] +Libs=@posix_mathlib@ +Cflags= + +[msvc] +Libs=@msvc_mathlib@ +Cflags= diff --git a/numpy/core/npymath.ini.in b/numpy/core/npymath.ini.in new file mode 100644 index 000000000..73379e47c --- /dev/null +++ b/numpy/core/npymath.ini.in @@ -0,0 +1,19 @@ +[meta] +Name=npymath +Description=Portable, core math library implementing C99 standard +Version=0.1 + +[variables] +prefix=@prefix@ +libdir=${prefix}@sep@lib +includedir=${prefix}@sep@include + +[default] +Libs=-L${libdir} -lnpymath +Cflags=-I${includedir} +Requires=mlib + +[msvc] +Libs=/LIBPATH:${libdir} npymath.lib +Cflags=/INCLUDE:${includedir} +Requires=mlib diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py index 50e3fd75b..c961bcc0f 100644 --- a/numpy/core/numeric.py +++ b/numpy/core/numeric.py @@ -6,7 +6,7 @@ __all__ = ['newaxis', 'ndarray', 'flatiter', 'ufunc', 'set_numeric_ops', 'can_cast', 'asarray', 'asanyarray', 'ascontiguousarray', 'asfortranarray', 'isfortran', 'empty_like', 'zeros_like', - 'correlate', 'acorrelate', 'convolve', 'inner', 'dot', 'outer', 'vdot', + 'correlate', 'convolve', 'inner', 'dot', 'outer', 'vdot', 'alterdot', 'restoredot', 'roll', 'rollaxis', 'cross', 'tensordot', 'array2string', 'get_printoptions', 'set_printoptions', 'array_repr', 'array_str', 'set_string_function', @@ -22,6 +22,7 @@ __all__ = ['newaxis', 'ndarray', 'flatiter', 'ufunc', 'CLIP', 'RAISE', 'WRAP', 'MAXDIMS', 'BUFSIZE', 'ALLOW_THREADS'] import sys +import warnings import multiarray import umath from umath import * @@ -362,24 +363,52 @@ def require(a, dtype=None, requirements=None): Parameters ---------- a : array_like - The object to be converted to a type-and-requirement satisfying array + The object to be converted to a type-and-requirement-satisfying array. dtype : data-type - The required data-type (None is the default data-type -- float64) - requirements : list of strings + The required data-type, the default data-type is float64). + requirements : str or list of str The requirements list can be any of the following - * 'ENSUREARRAY' ('E') - ensure that a base-class ndarray * 'F_CONTIGUOUS' ('F') - ensure a Fortran-contiguous array * 'C_CONTIGUOUS' ('C') - ensure a C-contiguous array * 'ALIGNED' ('A') - ensure a data-type aligned array - * 'WRITEABLE' ('W') - ensure a writeable array + * 'WRITEABLE' ('W') - ensure a writable array * 'OWNDATA' ('O') - ensure an array that owns its own data + See Also + -------- + asarray : Convert input to an ndarray. + asanyarray : Convert to an ndarray, but pass through ndarray subclasses. + ascontiguousarray : Convert input to a contiguous array. + asfortranarray : Convert input to an ndarray with column-major + memory order. + ndarray.flags : Information about the memory layout of the array. + Notes ----- The returned array will be guaranteed to have the listed requirements by making a copy if needed. + Examples + -------- + >>> x = np.arange(6).reshape(2,3) + >>> x.flags + C_CONTIGUOUS : True + F_CONTIGUOUS : False + OWNDATA : False + WRITEABLE : True + ALIGNED : True + UPDATEIFCOPY : False + + >>> y = np.require(x, dtype=np.float32, requirements=['A', 'O', 'W', 'F']) + >>> y.flags + C_CONTIGUOUS : False + F_CONTIGUOUS : True + OWNDATA : True + WRITEABLE : True + ALIGNED : True + UPDATEIFCOPY : False + """ if requirements is None: requirements = [] @@ -507,7 +536,7 @@ def argwhere(a): [1, 2]]) """ - return asarray(a.nonzero()).T + return transpose(asanyarray(a).nonzero()) def flatnonzero(a): """ @@ -557,7 +586,7 @@ def _mode_from_name(mode): return _mode_from_name_dict[mode.lower()[0]] return mode -def correlate(a,v,mode='valid'): +def correlate(a,v,mode='valid',old_behavior=True): """ Discrete, linear correlation of two 1-dimensional sequences. @@ -574,51 +603,48 @@ def correlate(a,v,mode='valid'): mode : {'valid', 'same', 'full'}, optional Refer to the `convolve` docstring. Note that the default is `valid`, unlike `convolve`, which uses `full`. + old_behavior : bool + If True, uses the old, numeric behavior (correlate(a,v) == correlate(v, + a), and the conjugate is not taken for complex arrays). If False, uses + the conventional signal processing definition (see note). See Also -------- convolve : Discrete, linear convolution of two one-dimensional sequences. - acorrelate: Discrete correlation following the usual signal processing - definition for complex arrays, and without assuming that - correlate(a, b) == correlate(b, a) - """ - mode = _mode_from_name(mode) - return multiarray.correlate(a,v,mode) - -def acorrelate(a, v, mode='valid'): - """ - Discrete, linear correlation of two 1-dimensional sequences. - This function computes the correlation as generally defined in signal - processing texts: + Note + ---- + If old_behavior is False, this function computes the correlation as + generally defined in signal processing texts:: - z[k] = sum_n a[n] * conj(v[n+k]) + z[k] = sum_n a[n] * conj(v[n+k]) with a and v sequences being zero-padded where necessary and conj being the conjugate. - Parameters - ---------- - a, v : array_like - Input sequences. - mode : {'valid', 'same', 'full'}, optional - Refer to the `convolve` docstring. Note that the default - is `valid`, unlike `convolve`, which uses `full`. - - Note - ---- - This is the function which corresponds to matlab xcorr. - - See Also + Examples -------- - convolve : Discrete, linear convolution of two - one-dimensional sequences. - correlate: Deprecated function to compute correlation + >>> np.correlate([1, 2, 3], [0, 1, 0.5]) + array([ 3.5]) + >>> np.correlate([1, 2, 3], [0, 1, 0.5], "same") + array([ 2. , 3.5, 3. ]) + >>> np.correlate([1, 2, 3], [0, 1, 0.5], "full") + array([ 0.5, 2. , 3.5, 3. , 0. ]) + """ mode = _mode_from_name(mode) - return multiarray.acorrelate(a, v, mode) - + if old_behavior: + warnings.warn(""" +The current behavior of correlate is deprecated for 1.4.0, and will be removed +for NumPy 1.5.0. + +The new behavior fits the conventional definition of correlation: inputs are +never swapped, and the second argument is conjugated for complex arrays.""", + DeprecationWarning) + return multiarray.correlate(a,v,mode) + else: + return multiarray.correlate2(a,v,mode) def convolve(a,v,mode='full'): """ @@ -1170,20 +1196,26 @@ def array_repr(arr, max_line_width=None, precision=None, suppress_small=None): Parameters ---------- arr : ndarray - Input array. - max_line_width : int - The maximum number of columns the string should span. Newline - characters splits the string appropriately after array elements. - precision : int - Floating point precision. - suppress_small : bool - Represent very small numbers as zero. + Input array. + max_line_width : int, optional + The maximum number of columns the string should span. Newline + characters split the string appropriately after array elements. + precision : int, optional + Floating point precision. Default is the current printing precision + (usually 8), which can be altered using `set_printoptions`. + suppress_small : bool, optional + Represent very small numbers as zero, default is False. Very small + is defined by `precision`, if the precision is 8 then + numbers smaller than 5e-9 are represented as zero. Returns ------- string : str The string representation of an array. + See Also + -------- + array_str, array2string, set_printoptions Examples -------- @@ -1194,6 +1226,10 @@ def array_repr(arr, max_line_width=None, precision=None, suppress_small=None): >>> np.array_repr(np.array([], np.int32)) 'array([], dtype=int32)' + >>> x = np.array([1e-6, 4e-7, 2, 3]) + >>> np.array_repr(x, precision=6, suppress_small=True) + 'array([ 0.000001, 0. , 2. , 3. ])' + """ if arr.size > 0 or arr.shape==(0,): lst = array2string(arr, max_line_width, precision, suppress_small, @@ -1221,7 +1257,11 @@ def array_repr(arr, max_line_width=None, precision=None, suppress_small=None): def array_str(a, max_line_width=None, precision=None, suppress_small=None): """ - Return a string representation of an array. + Return a string representation of the data in an array. + + The data in the array is returned as a single string. This function + is similar to `array_repr`, the difference is that `array_repr` also + returns information on the type of array and data type. Parameters ---------- @@ -1230,13 +1270,16 @@ def array_str(a, max_line_width=None, precision=None, suppress_small=None): max_line_width : int, optional Inserts newlines if text is longer than `max_line_width`. precision : int, optional - If `a` is float, `precision` sets floating point precision. - suppress_small : boolean, optional - Represent very small numbers as zero. + Floating point precision. Default is the current printing precision + (usually 8), which can be altered using set_printoptions. + suppress_small : bool, optional + Represent very small numbers as zero, default is False. Very small is + defined by precision, if the precision is 8 then numbers smaller than + 5e-9 are represented as zero. See Also -------- - array2string, array_repr + array2string, array_repr, set_printoptions Examples -------- @@ -1264,8 +1307,8 @@ def indices(dimensions, dtype=int): ---------- dimensions : sequence of ints The shape of the grid. - dtype : optional - Data_type of the result. + dtype : dtype, optional + Data type of the result. Returns ------- @@ -1291,7 +1334,7 @@ def indices(dimensions, dtype=int): Examples -------- - >>> grid = np.indices((2,3)) + >>> grid = np.indices((2, 3)) >>> grid.shape (2,2,3) >>> grid[0] # row indices @@ -1301,6 +1344,17 @@ def indices(dimensions, dtype=int): array([[0, 1, 2], [0, 1, 2]]) + The indices can be used as an index into an array. + + >>> x = np.arange(20).reshape(5, 4) + >>> row, col = np.indices((2, 3)) + >>> x[row, col] + array([[0, 1, 2], + [4, 5, 6]]) + + Note that it would be more straightforward in the above example to + extract the required elements directly with ``x[:2, :3]``. + """ dimensions = tuple(dimensions) N = len(dimensions) @@ -1638,15 +1692,9 @@ def identity(n, dtype=None): [ 0., 0., 1.]]) """ - a = array([1]+n*[0],dtype=dtype) - b = empty((n,n),dtype=dtype) - - # Note that this assignment depends on the convention that since the a - # array is shorter than the flattened b array, then the a array will - # be repeated until it is the appropriate size. Given a's construction, - # this nicely sets the diagonal to all ones. - b.flat = a - return b + a = zeros((n,n), dtype=dtype) + a.flat[::n+1] = 1 + return a def allclose(a, b, rtol=1.e-5, atol=1.e-8): """ @@ -1816,22 +1864,24 @@ def seterr(all=None, divide=None, over=None, under=None, invalid=None): Parameters ---------- - all : {'ignore', 'warn', 'raise', 'call'}, optional + all : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional Set treatment for all types of floating-point errors at once: - - ignore: Take no action when the exception occurs - - warn: Print a `RuntimeWarning` (via the Python `warnings` module) - - raise: Raise a `FloatingPointError` + - ignore: Take no action when the exception occurs. + - warn: Print a `RuntimeWarning` (via the Python `warnings` module). + - raise: Raise a `FloatingPointError`. - call: Call a function specified using the `seterrcall` function. + - print: Print a warning directly to ``stdout``. + - log: Record error in a Log object specified by `seterrcall`. The default is not to change the current behavior. - divide : {'ignore', 'warn', 'raise', 'call'}, optional + divide : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional Treatment for division by zero. - over : {'ignore', 'warn', 'raise', 'call'}, optional + over : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional Treatment for floating-point overflow. - under : {'ignore', 'warn', 'raise', 'call'}, optional + under : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional Treatment for floating-point underflow. - invalid : {'ignore', 'warn', 'raise', 'call'}, optional + invalid : {'ignore', 'warn', 'raise', 'call', 'print', 'log'}, optional Treatment for invalid floating-point operation. Returns @@ -1859,22 +1909,25 @@ def seterr(all=None, divide=None, over=None, under=None, invalid=None): Examples -------- - - Set mode: - - >>> seterr(over='raise') # doctest: +SKIP + >>> np.seterr(over='raise') {'over': 'ignore', 'divide': 'ignore', 'invalid': 'ignore', 'under': 'ignore'} + >>> np.seterr(all='ignore') # reset to default + {'over': 'raise', 'divide': 'warn', 'invalid': 'warn', 'under': 'warn'} - >>> old_settings = seterr(all='warn', over='raise') # doctest: +SKIP - - >>> int16(32000) * int16(3) # doctest: +SKIP + >>> np.int16(32000) * np.int16(3) + 30464 + >>> old_settings = np.seterr(all='warn', over='raise') + >>> np.int16(32000) * np.int16(3) Traceback (most recent call last): - File "<stdin>", line 1, in ? + File "<stdin>", line 1, in <module> FloatingPointError: overflow encountered in short_scalars - >>> seterr(all='ignore') # doctest: +SKIP - {'over': 'ignore', 'divide': 'ignore', 'invalid': 'ignore', - 'under': 'ignore'} + + >>> np.seterr(all='print') + {'over': 'print', 'divide': 'print', 'invalid': 'print', 'under': 'print'} + >>> np.int16(32000) * np.int16(3) + Warning: overflow encountered in short_scalars + 30464 """ @@ -1897,11 +1950,41 @@ def seterr(all=None, divide=None, over=None, under=None, invalid=None): def geterr(): - """Get the current way of handling floating-point errors. + """ + Get the current way of handling floating-point errors. + + Returns + ------- + res : dict + A dictionary with keys "divide", "over", "under", and "invalid", + whose values are from the strings "ignore", "print", "log", "warn", + "raise", and "call". The keys represent possible floating-point + exceptions, and the values define how these exceptions are handled. + + See Also + -------- + geterrcall, seterr, seterrcall + + Notes + ----- + For complete documentation of the types of floating-point exceptions and + treatment options, see `seterr`. + + Examples + -------- + >>> np.geterr() # default is all set to 'ignore' + {'over': 'ignore', 'divide': 'ignore', 'invalid': 'ignore', + 'under': 'ignore'} + >>> np.arange(3.) / np.arange(3.) + array([ NaN, 1., 1.]) + + >>> oldsettings = np.seterr(all='warn', over='raise') + >>> np.geterr() + {'over': 'raise', 'divide': 'warn', 'invalid': 'warn', 'under': 'warn'} + >>> np.arange(3.) / np.arange(3.) + __main__:1: RuntimeWarning: invalid value encountered in divide + array([ NaN, 1., 1.]) - Returns a dictionary with entries "divide", "over", "under", and - "invalid", whose values are from the strings - "ignore", "print", "log", "warn", "raise", and "call". """ maskvalue = umath.geterrobj()[1] mask = 7 @@ -1952,13 +2035,13 @@ def seterrcall(func): is to set the error-handler to 'call', using `seterr`. Then, set the function to call using this function. - The second is to set the error-handler to `log`, using `seterr`. + The second is to set the error-handler to 'log', using `seterr`. Floating-point errors then trigger a call to the 'write' method of the provided object. Parameters ---------- - log_func_or_obj : callable f(err, flag) or object with write method + func : callable f(err, flag) or object with write method Function to call upon floating-point errors ('call'-mode) or object whose 'write' method is used to log such message ('log'-mode). @@ -1971,7 +2054,7 @@ def seterrcall(func): In other words, ``flags = divide + 2*over + 4*under + 8*invalid``. - If an object is provided, it's write method should take one argument, + If an object is provided, its write method should take one argument, a string. Returns @@ -1979,6 +2062,10 @@ def seterrcall(func): h : callable or log instance The old error handler. + See Also + -------- + seterr, geterr, geterrcall + Examples -------- Callback upon error: @@ -2025,7 +2112,45 @@ def seterrcall(func): return old def geterrcall(): - """Return the current callback function used on floating-point errors. + """ + Return the current callback function used on floating-point errors. + + When the error handling for a floating-point error (one of "divide", + "over", "under", or "invalid") is set to 'call' or 'log', the function + that is called or the log instance that is written to is returned by + `geterrcall`. This function or log instance has been set with + `seterrcall`. + + Returns + ------- + errobj : callable, log instance or None + The current error handler. If no handler was set through `seterrcall`, + ``None`` is returned. + + See Also + -------- + seterrcall, seterr, geterr + + Notes + ----- + For complete documentation of the types of floating-point exceptions and + treatment options, see `seterr`. + + Examples + -------- + >>> np.geterrcall() # we did not yet set a handler, returns None + + >>> oldsettings = np.seterr(all='call') + >>> def err_handler(type, flag): + ... print "Floating point error (%s), with flag %s" % (type, flag) + >>> oldhandler = np.seterrcall(err_handler) + >>> np.array([1,2,3])/0.0 + Floating point error (divide by zero), with flag 1 + array([ Inf, Inf, Inf]) + >>> cur_handler = np.geterrcall() + >>> cur_handler is err_handler + True + """ return umath.geterrobj()[2] diff --git a/numpy/core/numerictypes.py b/numpy/core/numerictypes.py index c72cc122a..70315f1e0 100644 --- a/numpy/core/numerictypes.py +++ b/numpy/core/numerictypes.py @@ -521,6 +521,9 @@ def issubdtype(arg1, arg2): arg2 : dtype_like dtype or string representing a typecode. + Returns + ------- + out : bool See Also -------- @@ -660,14 +663,24 @@ def _find_common_coerce(a, b): thisind = __test_types.index(a.char) except ValueError: return None + return _can_coerce_all([a,b], start=thisind) + +# Find a data-type that all data-types in a list can be coerced to +def _can_coerce_all(dtypelist, start=0): + N = len(dtypelist) + if N == 0: + return None + if N == 1: + return dtypelist[0] + thisind = start while thisind < __len_test_types: newdtype = dtype(__test_types[thisind]) - if newdtype >= b and newdtype >= a: + numcoerce = len([x for x in dtypelist if newdtype >= x]) + if numcoerce == N: return newdtype thisind += 1 return None - def find_common_type(array_types, scalar_types): """ Determine common type following standard coercion rules @@ -696,16 +709,14 @@ def find_common_type(array_types, scalar_types): array_types = [dtype(x) for x in array_types] scalar_types = [dtype(x) for x in scalar_types] - if len(scalar_types) == 0: - if len(array_types) == 0: - return None - else: - return max(array_types) - if len(array_types) == 0: - return max(scalar_types) + maxa = _can_coerce_all(array_types) + maxsc = _can_coerce_all(scalar_types) - maxa = max(array_types) - maxsc = max(scalar_types) + if maxa is None: + return maxsc + + if maxsc is None: + return maxa try: index_a = _kind_list.index(maxa.kind) diff --git a/numpy/core/setup.py b/numpy/core/setup.py index dca1787a9..6c66512b0 100644 --- a/numpy/core/setup.py +++ b/numpy/core/setup.py @@ -1,11 +1,13 @@ import imp import os import sys +import shutil from os.path import join from numpy.distutils import log from distutils.dep_util import newer from distutils.sysconfig import get_config_var import warnings +import re from setup_common import * @@ -142,7 +144,7 @@ def check_math_capabilities(config, moredefs, mathlibs): # config.h in the public namespace, so we have a clash for the common # functions we test. We remove every function tested by python's # autoconf, hoping their own test are correct - if sys.version_info[:2] >= (2, 6): + if sys.version_info[:2] >= (2, 5): for f in OPTIONAL_STDFUNCS_MAYBE: if config.check_decl(fname2def(f), headers=["Python.h", "math.h"]): @@ -321,6 +323,7 @@ def configuration(parent_package='',top_path=None): d = os.path.dirname(target) if not os.path.exists(d): os.makedirs(d) + if newer(__file__,target): config_cmd = config.get_config_cmd() log.info('Generating %s',target) @@ -434,8 +437,6 @@ def configuration(parent_package='',top_path=None): # Check wether we can use inttypes (C99) formats if config_cmd.check_decl('PRIdPTR', headers = ['inttypes.h']): moredefs.append(('NPY_USE_C99_FORMATS', 1)) - else: - moredefs.append(('NPY_USE_C99_FORMATS', 0)) # Inline check inline = config_cmd.check_inline() @@ -548,12 +549,13 @@ def configuration(parent_package='',top_path=None): return [] config.add_data_files('include/numpy/*.h') + config.add_include_dirs(join('src', 'npymath')) config.add_include_dirs(join('src', 'multiarray')) config.add_include_dirs(join('src', 'umath')) config.numpy_include_dirs.extend(config.paths('include')) - deps = [join('src','_signbit.c'), + deps = [join('src','npymath','_signbit.c'), join('include','numpy','*object.h'), 'include/numpy/fenv/fenv.c', 'include/numpy/fenv/fenv.h', @@ -578,10 +580,28 @@ def configuration(parent_package='',top_path=None): # (don't ask). Because clib are generated before extensions, we have to # explicitly add an extension which has generate_config_h and # generate_numpyconfig_h as sources *before* adding npymath. - config.add_library('npymath', - sources=[join('src', 'npy_math.c.src')], - depends=[]) - + + subst_dict = dict([("sep", os.path.sep)]) + def get_mathlib_info(*args): + # Another ugly hack: the mathlib info is known once build_src is run, + # but we cannot use add_installed_pkg_config here either, so we only + # updated the substition dictionary during npymath build + config_cmd = config.get_config_cmd() + mlibs = check_mathlib(config_cmd) + + posix_mlib = ' '.join(['-l%s' % l for l in mlibs]) + msvc_mlib = ' '.join(['%s.lib' % l for l in mlibs]) + subst_dict["posix_mathlib"] = posix_mlib + subst_dict["msvc_mathlib"] = msvc_mlib + + config.add_installed_library('npymath', + sources=[join('src', 'npymath', 'npy_math.c.src'), get_mathlib_info], + install_dir='lib') + config.add_npy_pkg_config("npymath.ini.in", "lib/npy-pkg-config", + subst_dict) + config.add_npy_pkg_config("mlib.ini.in", "lib/npy-pkg-config", + subst_dict) + multiarray_deps = [ join('src', 'multiarray', 'arrayobject.h'), join('src', 'multiarray', 'arraytypes.h'), @@ -706,6 +726,9 @@ def configuration(parent_package='',top_path=None): config.add_extension('umath_tests', sources = [join('src','umath', 'umath_tests.c.src')]) + config.add_extension('multiarray_tests', + sources = [join('src', 'multiarray', 'multiarray_tests.c.src')]) + config.add_data_dir('tests') config.add_data_dir('tests/data') diff --git a/numpy/core/setup_common.py b/numpy/core/setup_common.py index 37023e8db..ab801fc6d 100644 --- a/numpy/core/setup_common.py +++ b/numpy/core/setup_common.py @@ -76,8 +76,8 @@ def check_api_version(apiversion, codegen_dir): "with checksum %s, but recorded checksum for C API version %d in " \ "codegen_dir/cversions.txt is %s. If functions were added in the " \ "C API, you have to update C_API_VERSION in %s." - warnings.warn(msg % (apiversion, curapi_hash, apiversion, api_hash, - __file__), + warnings.warn(msg % (apiversion, curapi_hash, apiversion, api_hash, + __file__), MismatchCAPIWarning) # Mandatory functions: if not found, fail the build MANDATORY_FUNCS = ["sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs", @@ -87,17 +87,19 @@ MANDATORY_FUNCS = ["sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs", # Standard functions which may not be available and for which we have a # replacement implementation. Note that some of these are C99 functions. OPTIONAL_STDFUNCS = ["expm1", "log1p", "acosh", "asinh", "atanh", - "rint", "trunc", "exp2", "log2"] + "rint", "trunc", "exp2", "log2", "hypot", "atan2", "pow", + "copysign"] # Subset of OPTIONAL_STDFUNCS which may alreay have HAVE_* defined by Python.h -OPTIONAL_STDFUNCS_MAYBE = ["expm1", "log1p", "acosh", "atanh", "asinh"] +OPTIONAL_STDFUNCS_MAYBE = ["expm1", "log1p", "acosh", "atanh", "asinh", "hypot", + "copysign"] # C99 functions: float and long double versions C99_FUNCS = ["sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs", "floor", "ceil", "rint", "trunc", "sqrt", "log10", "log", "log1p", "exp", "expm1", "asin", "acos", "atan", "asinh", "acosh", "atanh", "hypot", "atan2", "pow", "fmod", "modf", 'frexp', 'ldexp', - "exp2", "log2"] + "exp2", "log2", "copysign"] C99_FUNCS_SINGLE = [f + 'f' for f in C99_FUNCS] C99_FUNCS_EXTENDED = [f + 'l' for f in C99_FUNCS] diff --git a/numpy/core/setupscons.py b/numpy/core/setupscons.py index be42246ad..3dfaa48d9 100644 --- a/numpy/core/setupscons.py +++ b/numpy/core/setupscons.py @@ -63,8 +63,10 @@ def configuration(parent_package='',top_path=None): # XXX: I really have to think about how to communicate path info # between scons and distutils, and set the options at one single # location. - h_file = join(get_scons_pkg_build_dir(config.name), '__multiarray_api.h') - t_file = join(get_scons_pkg_build_dir(config.name), 'multiarray_api.txt') + h_file = join(get_scons_pkg_build_dir(config.name), + 'include/numpy/__multiarray_api.h') + t_file = join(get_scons_pkg_build_dir(config.name), + 'include/numpy/multiarray_api.txt') config.add_data_files((header_dir, h_file), (header_dir, t_file)) @@ -73,8 +75,10 @@ def configuration(parent_package='',top_path=None): # XXX: I really have to think about how to communicate path info # between scons and distutils, and set the options at one single # location. - h_file = join(get_scons_pkg_build_dir(config.name), '__ufunc_api.h') - t_file = join(get_scons_pkg_build_dir(config.name), 'ufunc_api.txt') + h_file = join(get_scons_pkg_build_dir(config.name), + 'include/numpy/__ufunc_api.h') + t_file = join(get_scons_pkg_build_dir(config.name), + 'include/numpy/ufunc_api.txt') config.add_data_files((header_dir, h_file), (header_dir, t_file)) @@ -87,6 +91,7 @@ def configuration(parent_package='',top_path=None): config.add_sconscript('SConstruct', post_hook = add_generated_files, source_files = source_files) + config.add_scons_installed_library('npymath', 'lib') config.add_data_files('include/numpy/*.h') config.add_include_dirs('src') diff --git a/numpy/core/src/_sortmodule.c.src b/numpy/core/src/_sortmodule.c.src index 51c5feb41..28299c1a7 100644 --- a/numpy/core/src/_sortmodule.c.src +++ b/numpy/core/src/_sortmodule.c.src @@ -1,85 +1,257 @@ /* -*- c -*- */ -/* The purpose of this module is to add faster sort functions - that are type-specific. This is done by altering the - function table for the builtin descriptors. - - These sorting functions are copied almost directly from numarray - with a few modifications (complex comparisons compare the imaginary - part if the real parts are equal, for example), and the names - are changed. - - The original sorting code is due to Charles R. Harris who wrote - it for numarray. -*/ - -/* Quick sort is usually the fastest, but the worst case scenario can - be slower than the merge and heap sorts. The merge sort requires - extra memory and so for large arrays may not be useful. - - The merge sort is *stable*, meaning that equal components - are unmoved from their entry versions, so it can be used to - implement lexigraphic sorting on multiple keys. +/* + * The purpose of this module is to add faster sort functions + * that are type-specific. This is done by altering the + * function table for the builtin descriptors. + * + * These sorting functions are copied almost directly from numarray + * with a few modifications (complex comparisons compare the imaginary + * part if the real parts are equal, for example), and the names + * are changed. + * + * The original sorting code is due to Charles R. Harris who wrote + * it for numarray. + */ - The heap sort is included for completeness. -*/ +/* + * Quick sort is usually the fastest, but the worst case scenario can + * be slower than the merge and heap sorts. The merge sort requires + * extra memory and so for large arrays may not be useful. + * + * The merge sort is *stable*, meaning that equal components + * are unmoved from their entry versions, so it can be used to + * implement lexigraphic sorting on multiple keys. + * + * The heap sort is included for completeness. + */ #include "Python.h" #include "numpy/noprefix.h" +#include "numpy/npy_math.h" +#define NOT_USED NPY_UNUSED(unused) #define PYA_QS_STACK 100 #define SMALL_QUICKSORT 15 #define SMALL_MERGESORT 20 #define SMALL_STRING 16 -#define SWAP(a,b) {SWAP_temp = (b); (b)=(a); (a) = SWAP_temp;} -#define STDC_LT(a,b) ((a) < (b)) -#define STDC_LE(a,b) ((a) <= (b)) -#define STDC_EQ(a,b) ((a) == (b)) -#define NUMC_LT(p,q) ((((p).real==(q).real) ? ((p).imag < (q).imag): ((p).real < (q).real))) -#define NUMC_LE(p,q) ((((p).real==(q).real) ? ((p).imag <= (q).imag): ((p).real <= (q).real))) -#define NUMC_EQ(p,q) (((p).real==(q).real) && ((p).imag == (q).imag)) -#define STRING_LT(pa, pb, len) (compare_string(pa, pb, len) < 0) -#define STRING_LE(pa, pb, len) (compare_string(pa, pb, len) <= 0) -#define STRING_EQ(pa, pb, len) (compare_string(pa, pb, len) == 0) -#define UNICODE_LT(pa, pb, len) (compare_ucs4(pa, pb, len) < 0) -#define UNICODE_LE(pa, pb, len) (compare_ucs4(pa, pb, len) <= 0) -#define UNICODE_EQ(pa, pb, len) (compare_ucs4(pa, pb, len) == 0) + +/* + ***************************************************************************** + ** SWAP MACROS ** + ***************************************************************************** + */ + +/**begin repeat + * + * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG, + * LONGLONG, ULONGLONG, FLOAT, DOUBLE, LONGDOUBLE, CFLOAT, + * CDOUBLE,CLONGDOUBLE, INTP# + * #type = npy_bool, npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, + * npy_uint, npy_long, npy_ulong, npy_longlong, npy_ulonglong, + * npy_float, npy_double, npy_longdouble, npy_cfloat, npy_cdouble, + * npy_clongdouble, npy_intp# + */ +#define @TYPE@_SWAP(a,b) {@type@ tmp = (b); (b)=(a); (a) = tmp;} + +/**end repeat**/ + +/* + ***************************************************************************** + ** COMPARISON FUNCTIONS ** + ***************************************************************************** + */ + +/**begin repeat + * + * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG, + * LONGLONG, ULONGLONG# + * #type = Bool, byte, ubyte, short, ushort, int, uint, long, ulong, + * longlong, ulonglong# + */ +NPY_INLINE static int +@TYPE@_LT(@type@ a, @type@ b) +{ + return a < b; +} +/**end repeat**/ + + +/**begin repeat + * + * #TYPE = FLOAT, DOUBLE, LONGDOUBLE# + * #type = float, double, longdouble# + */ +NPY_INLINE static int +@TYPE@_LT(@type@ a, @type@ b) +{ + return a < b || (b != b && a == a); +} +/**end repeat**/ + + +/* + * For inline functions SUN recommends not using a return in the then part + * of an if statement. It's a SUN compiler thing, so assign the return value + * to a variable instead. + */ + +/**begin repeat + * + * #TYPE = CFLOAT, CDOUBLE, CLONGDOUBLE# + * #type = cfloat, cdouble, clongdouble# + */ +NPY_INLINE static int +@TYPE@_LT(@type@ a, @type@ b) +{ + int ret; + + if (a.real < b.real) { + ret = a.imag == a.imag || b.imag != b.imag; + } + else if (a.real > b.real) { + ret = b.imag != b.imag && a.imag == a.imag; + } + else if (a.real == b.real || (a.real != a.real && b.real != b.real)) { + ret = a.imag < b.imag || (b.imag != b.imag && a.imag == a.imag); + } + else { + ret = b.real != b.real; + } + + return ret; +} +/**end repeat**/ + + +/* The PyObject functions are stubs for later use */ +NPY_INLINE static int +PyObject_LT(PyObject *pa, PyObject *pb) +{ + return 0; +} + + +NPY_INLINE static void +STRING_COPY(char *s1, char *s2, size_t len) +{ + memcpy(s1, s2, len); +} + + +NPY_INLINE static void +STRING_SWAP(char *s1, char *s2, size_t len) +{ + while(len--) { + const char t = *s1; + *s1++ = *s2; + *s2++ = t; + } +} + + +NPY_INLINE static int +STRING_LT(char *s1, char *s2, size_t len) +{ + const unsigned char *c1 = (unsigned char *)s1; + const unsigned char *c2 = (unsigned char *)s2; + size_t i; + int ret = 0; + + for (i = 0; i < len; ++i) { + if (c1[i] != c2[i]) { + ret = c1[i] < c2[i]; + break; + } + } + return ret; +} + + +NPY_INLINE static void +UNICODE_COPY(npy_ucs4 *s1, npy_ucs4 *s2, size_t len) +{ + while(len--) { + *s1++ = *s2++; + } +} + + +NPY_INLINE static void +UNICODE_SWAP(npy_ucs4 *s1, npy_ucs4 *s2, size_t len) +{ + while(len--) { + const npy_ucs4 t = *s1; + *s1++ = *s2; + *s2++ = t; + } +} + + +NPY_INLINE static int +UNICODE_LT(npy_ucs4 *s1, npy_ucs4 *s2, size_t len) +{ + size_t i; + int ret = 0; + + for (i = 0; i < len; ++i) { + if (s1[i] != s2[i]) { + ret = s1[i] < s2[i]; + break; + } + } + return ret; +} + + +/* + ***************************************************************************** + ** NUMERIC SORTS ** + ***************************************************************************** + */ /**begin repeat - #TYPE=BOOL,BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE,CFLOAT,CDOUBLE,CLONGDOUBLE,DATETIME,TIMEDELTA# - #type=Bool,byte,ubyte,short,ushort,int,uint,long,ulong,longlong,ulonglong,float,double,longdouble,cfloat,cdouble,clongdouble,datetime,timedelta# - #lessthan=STDC_LT*14,NUMC_LT*3,STDC_LT*2# - #lessequal=STDC_LE*14,NUMC_LE*3,STDC_LE*2# -**/ + * + * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG, + * LONGLONG, ULONGLONG, FLOAT, DOUBLE, LONGDOUBLE, + * CFLOAT, CDOUBLE, CLONGDOUBLE# + * #type = Bool, byte, ubyte, short, ushort, int, uint, long, ulong, + * longlong, ulonglong, float, double, longdouble, + * cfloat, cdouble, clongdouble# + */ + + static int -@TYPE@_quicksort(@type@ *start, intp num, void * NPY_UNUSED(unused)) +@TYPE@_quicksort(@type@ *start, npy_intp num, void *NOT_USED) { @type@ *pl = start; @type@ *pr = start + num - 1; - @type@ vp, SWAP_temp; + @type@ vp; @type@ *stack[PYA_QS_STACK], **sptr = stack, *pm, *pi, *pj, *pk; - for(;;) { + for (;;) { while ((pr - pl) > SMALL_QUICKSORT) { /* quicksort partition */ pm = pl + ((pr - pl) >> 1); - if (@lessthan@(*pm, *pl)) SWAP(*pm, *pl); - if (@lessthan@(*pr, *pm)) SWAP(*pr, *pm); - if (@lessthan@(*pm, *pl)) SWAP(*pm, *pl); + if (@TYPE@_LT(*pm, *pl)) @TYPE@_SWAP(*pm, *pl); + if (@TYPE@_LT(*pr, *pm)) @TYPE@_SWAP(*pr, *pm); + if (@TYPE@_LT(*pm, *pl)) @TYPE@_SWAP(*pm, *pl); vp = *pm; pi = pl; pj = pr - 1; - SWAP(*pm, *pj); - for(;;) { - do ++pi; while (@lessthan@(*pi, vp)); - do --pj; while (@lessthan@(vp, *pj)); - if (pi >= pj) break; - SWAP(*pi,*pj); + @TYPE@_SWAP(*pm, *pj); + for (;;) { + do ++pi; while (@TYPE@_LT(*pi, vp)); + do --pj; while (@TYPE@_LT(vp, *pj)); + if (pi >= pj) { + break; + } + @TYPE@_SWAP(*pi,*pj); } pk = pr - 1; - SWAP(*pi, *pk); + @TYPE@_SWAP(*pi, *pk); /* push largest partition on stack */ if (pi - pl < pr - pi) { *sptr++ = pi + 1; @@ -94,16 +266,18 @@ static int } /* insertion sort */ - for(pi = pl + 1; pi <= pr; ++pi) { + for (pi = pl + 1; pi <= pr; ++pi) { vp = *pi; pj = pi; pk = pi - 1; - while (pj > pl && @lessthan@(vp, *pk)) { + while (pj > pl && @TYPE@_LT(vp, *pk)) { *pj-- = *pk--; } *pj = vp; } - if (sptr == stack) break; + if (sptr == stack) { + break; + } pr = *(--sptr); pl = *(--sptr); } @@ -112,34 +286,36 @@ static int } static int -@TYPE@_aquicksort(@type@ *v, intp* tosort, intp num, void *NPY_UNUSED(unused)) +@TYPE@_aquicksort(@type@ *v, npy_intp* tosort, npy_intp num, void *NOT_USED) { @type@ vp; - intp *pl, *pr, SWAP_temp; - intp *stack[PYA_QS_STACK], **sptr=stack, *pm, *pi, *pj, *pk, vi; + npy_intp *pl, *pr; + npy_intp *stack[PYA_QS_STACK], **sptr=stack, *pm, *pi, *pj, *pk, vi; pl = tosort; pr = tosort + num - 1; - for(;;) { + for (;;) { while ((pr - pl) > SMALL_QUICKSORT) { /* quicksort partition */ pm = pl + ((pr - pl) >> 1); - if (@lessthan@(v[*pm],v[*pl])) SWAP(*pm,*pl); - if (@lessthan@(v[*pr],v[*pm])) SWAP(*pr,*pm); - if (@lessthan@(v[*pm],v[*pl])) SWAP(*pm,*pl); + if (@TYPE@_LT(v[*pm],v[*pl])) INTP_SWAP(*pm,*pl); + if (@TYPE@_LT(v[*pr],v[*pm])) INTP_SWAP(*pr,*pm); + if (@TYPE@_LT(v[*pm],v[*pl])) INTP_SWAP(*pm,*pl); vp = v[*pm]; pi = pl; pj = pr - 1; - SWAP(*pm,*pj); - for(;;) { - do ++pi; while (@lessthan@(v[*pi],vp)); - do --pj; while (@lessthan@(vp,v[*pj])); - if (pi >= pj) break; - SWAP(*pi,*pj); - } - pk = pr - 1; - SWAP(*pi,*pk); + INTP_SWAP(*pm,*pj); + for (;;) { + do ++pi; while (@TYPE@_LT(v[*pi],vp)); + do --pj; while (@TYPE@_LT(vp,v[*pj])); + if (pi >= pj) { + break; + } + INTP_SWAP(*pi,*pj); + } + pk = pr - 1; + INTP_SWAP(*pi,*pk); /* push largest partition on stack */ if (pi - pl < pr - pi) { *sptr++ = pi + 1; @@ -154,17 +330,19 @@ static int } /* insertion sort */ - for(pi = pl + 1; pi <= pr; ++pi) { + for (pi = pl + 1; pi <= pr; ++pi) { vi = *pi; vp = v[vi]; pj = pi; pk = pi - 1; - while (pj > pl && @lessthan@(vp, v[*pk])) { + while (pj > pl && @TYPE@_LT(vp, v[*pk])) { *pj-- = *pk--; } *pj = vi; } - if (sptr == stack) break; + if (sptr == stack) { + break; + } pr = *(--sptr); pl = *(--sptr); } @@ -174,10 +352,10 @@ static int static int -@TYPE@_heapsort(@type@ *start, intp n, void *NPY_UNUSED(unused)) +@TYPE@_heapsort(@type@ *start, npy_intp n, void *NOT_USED) { @type@ tmp, *a; - intp i,j,l; + npy_intp i,j,l; /* The array needs to be offset by one for heapsort indexing */ a = start - 1; @@ -185,15 +363,17 @@ static int for (l = n>>1; l > 0; --l) { tmp = a[l]; for (i = l, j = l<<1; j <= n;) { - if (j < n && @lessthan@(a[j], a[j+1])) + if (j < n && @TYPE@_LT(a[j], a[j+1])) { j += 1; - if (@lessthan@(tmp, a[j])) { + } + if (@TYPE@_LT(tmp, a[j])) { a[i] = a[j]; i = j; j += j; } - else + else { break; + } } a[i] = tmp; } @@ -203,15 +383,17 @@ static int a[n] = a[1]; n -= 1; for (i = 1, j = 2; j <= n;) { - if (j < n && @lessthan@(a[j], a[j+1])) + if (j < n && @TYPE@_LT(a[j], a[j+1])) { j++; - if (@lessthan@(tmp, a[j])) { + } + if (@TYPE@_LT(tmp, a[j])) { a[i] = a[j]; i = j; j += j; } - else + else { break; + } } a[i] = tmp; } @@ -220,24 +402,26 @@ static int } static int -@TYPE@_aheapsort(@type@ *v, intp *tosort, intp n, void *NPY_UNUSED(unused)) +@TYPE@_aheapsort(@type@ *v, npy_intp *tosort, npy_intp n, void *NOT_USED) { - intp *a, i,j,l, tmp; + npy_intp *a, i,j,l, tmp; /* The arrays need to be offset by one for heapsort indexing */ a = tosort - 1; for (l = n>>1; l > 0; --l) { tmp = a[l]; for (i = l, j = l<<1; j <= n;) { - if (j < n && @lessthan@(v[a[j]], v[a[j+1]])) + if (j < n && @TYPE@_LT(v[a[j]], v[a[j+1]])) { j += 1; - if (@lessthan@(v[tmp], v[a[j]])) { + } + if (@TYPE@_LT(v[tmp], v[a[j]])) { a[i] = a[j]; i = j; j += j; } - else + else { break; + } } a[i] = tmp; } @@ -247,15 +431,17 @@ static int a[n] = a[1]; n -= 1; for (i = 1, j = 2; j <= n;) { - if (j < n && @lessthan@(v[a[j]], v[a[j+1]])) + if (j < n && @TYPE@_LT(v[a[j]], v[a[j+1]])) { j++; - if (@lessthan@(v[tmp], v[a[j]])) { + } + if (@TYPE@_LT(v[tmp], v[a[j]])) { a[i] = a[j]; i = j; j += j; } - else + else { break; + } } a[i] = tmp; } @@ -273,17 +459,17 @@ static void pm = pl + ((pr - pl) >> 1); @TYPE@_mergesort0(pl, pm, pw); @TYPE@_mergesort0(pm, pr, pw); - for(pi = pw, pj = pl; pj < pm;) { + for (pi = pw, pj = pl; pj < pm;) { *pi++ = *pj++; } pj = pw; pk = pl; while (pj < pi && pm < pr) { - if (@lessequal@(*pj,*pm)) { - *pk = *pj++; + if (@TYPE@_LT(*pm,*pj)) { + *pk = *pm++; } else { - *pk = *pm++; + *pk = *pj++; } pk++; } @@ -293,11 +479,11 @@ static void } else { /* insertion sort */ - for(pi = pl + 1; pi < pr; ++pi) { + for (pi = pl + 1; pi < pr; ++pi) { vp = *pi; pj = pi; pk = pi -1; - while (pj > pl && @lessthan@(vp, *pk)) { + while (pj > pl && @TYPE@_LT(vp, *pk)) { *pj-- = *pk--; } *pj = vp; @@ -306,7 +492,7 @@ static void } static int -@TYPE@_mergesort(@type@ *start, intp num, void *NPY_UNUSED(unused)) +@TYPE@_mergesort(@type@ *start, npy_intp num, void *NOT_USED) { @type@ *pl, *pr, *pw; @@ -324,39 +510,39 @@ static int } static void -@TYPE@_amergesort0(intp *pl, intp *pr, @type@ *v, intp *pw) +@TYPE@_amergesort0(npy_intp *pl, npy_intp *pr, @type@ *v, npy_intp *pw) { @type@ vp; - intp vi, *pi, *pj, *pk, *pm; + npy_intp vi, *pi, *pj, *pk, *pm; if (pr - pl > SMALL_MERGESORT) { /* merge sort */ pm = pl + ((pr - pl + 1)>>1); @TYPE@_amergesort0(pl,pm-1,v,pw); @TYPE@_amergesort0(pm,pr,v,pw); - for(pi = pw, pj = pl; pj < pm; ++pi, ++pj) { + for (pi = pw, pj = pl; pj < pm; ++pi, ++pj) { *pi = *pj; } - for(pk = pw, pm = pl; pk < pi && pj <= pr; ++pm) { - if (@lessequal@(v[*pk],v[*pj])) { - *pm = *pk; - ++pk; - } - else { + for (pk = pw, pm = pl; pk < pi && pj <= pr; ++pm) { + if (@TYPE@_LT(v[*pj],v[*pk])) { *pm = *pj; ++pj; } + else { + *pm = *pk; + ++pk; + } } - for(; pk < pi; ++pm, ++pk) { + for (; pk < pi; ++pm, ++pk) { *pm = *pk; } } else { /* insertion sort */ - for(pi = pl + 1; pi <= pr; ++pi) { + for (pi = pl + 1; pi <= pr; ++pi) { vi = *pi; vp = v[vi]; - for(pj = pi, pk = pi - 1; pj > pl && @lessthan@(vp, v[*pk]); --pj, --pk) { + for (pj = pi, pk = pi - 1; pj > pl && @TYPE@_LT(vp, v[*pk]); --pj, --pk) { *pj = *pk; } *pj = vi; @@ -365,9 +551,9 @@ static void } static int -@TYPE@_amergesort(@type@ *v, intp *tosort, intp num, void *NPY_UNUSED(unused)) +@TYPE@_amergesort(@type@ *v, npy_intp *tosort, npy_intp num, void *NOT_USED) { - intp *pl, *pr, *pw; + npy_intp *pl, *pr, *pw; pl = tosort; pr = pl + num - 1; pw = PyDimMem_NEW((1+num/2)); @@ -382,85 +568,22 @@ static int return 0; } + + /**end repeat**/ /* - * Subroutines that will hopefully be inlined when the code - * for strings and unicode is compiled with proper flags. + ***************************************************************************** + ** STRING SORTS ** + ***************************************************************************** */ -#define copy_string memcpy - - -static void -swap_string(char *s1, char *s2, size_t len) -{ - while(len--) { - const char t = *s1; - *s1++ = *s2; - *s2++ = t; - } -} - - -static int -compare_string(char *s1, char *s2, size_t len) -{ - const unsigned char *c1 = (unsigned char *)s1; - const unsigned char *c2 = (unsigned char *)s2; - size_t i; - - for(i = 0; i < len; ++i) { - if (c1[i] != c2[i]) { - return (c1[i] > c2[i]) ? 1 : -1; - } - } - return 0; -} - - -static void -copy_ucs4(npy_ucs4 *s1, npy_ucs4 *s2, size_t len) -{ - while(len--) { - *s1++ = *s2++; - } -} - - -static void -swap_ucs4(npy_ucs4 *s1, npy_ucs4 *s2, size_t len) -{ - while(len--) { - const npy_ucs4 t = *s1; - *s1++ = *s2; - *s2++ = t; - } -} - - -static int -compare_ucs4(npy_ucs4 *s1, npy_ucs4 *s2, size_t len) -{ - size_t i; - - for(i = 0; i < len; ++i) { - if (s1[i] != s2[i]) { - return (s1[i] > s2[i]) ? 1 : -1; - } - } - return 0; -} - /**begin repeat - #TYPE=STRING, UNICODE# - #type=char, PyArray_UCS4# - #lessthan=STRING_LT, UNICODE_LT# - #lessequal=STRING_LE, UNICODE_LE# - #swap=swap_string, swap_ucs4# - #copy=copy_string, copy_ucs4# -**/ + * + * #TYPE = STRING, UNICODE# + * #type = char, PyArray_UCS4# + */ static void @TYPE@_mergesort0(@type@ *pl, @type@ *pr, @type@ *pw, @type@ *vp, size_t len) @@ -472,41 +595,41 @@ static void pm = pl + (((pr - pl)/len) >> 1)*len; @TYPE@_mergesort0(pl, pm, pw, vp, len); @TYPE@_mergesort0(pm, pr, pw, vp, len); - @copy@(pw, pl, pm - pl); + @TYPE@_COPY(pw, pl, pm - pl); pi = pw + (pm - pl); pj = pw; pk = pl; while (pj < pi && pm < pr) { - if (@lessequal@(pj, pm, len)) { - @copy@(pk, pj, len); - pj += len; + if (@TYPE@_LT(pm, pj, len)) { + @TYPE@_COPY(pk, pm, len); + pm += len; } else { - @copy@(pk, pm, len); - pm += len; + @TYPE@_COPY(pk, pj, len); + pj += len; } pk += len; } - @copy@(pk, pj, pi - pj); + @TYPE@_COPY(pk, pj, pi - pj); } else { /* insertion sort */ - for(pi = pl + len; pi < pr; pi += len) { - @copy@(vp, pi, len); + for (pi = pl + len; pi < pr; pi += len) { + @TYPE@_COPY(vp, pi, len); pj = pi; pk = pi - len; - while (pj > pl && @lessthan@(vp, pk, len)) { - @copy@(pj, pk, len); + while (pj > pl && @TYPE@_LT(vp, pk, len)) { + @TYPE@_COPY(pj, pk, len); pj -= len; pk -= len; } - @copy@(pj, vp, len); + @TYPE@_COPY(pj, vp, len); } } } static int -@TYPE@_mergesort(@type@ *start, intp num, PyArrayObject *arr) +@TYPE@_mergesort(@type@ *start, npy_intp num, PyArrayObject *arr) { const size_t elsize = arr->descr->elsize; const size_t len = elsize / sizeof(@type@); @@ -537,7 +660,7 @@ fail_0: } static int -@TYPE@_quicksort(@type@ *start, intp num, PyArrayObject *arr) +@TYPE@_quicksort(@type@ *start, npy_intp num, PyArrayObject *arr) { const size_t len = arr->descr->elsize/sizeof(@type@); @type@ *vp = malloc(arr->descr->elsize); @@ -545,25 +668,27 @@ static int @type@ *pr = start + (num - 1)*len; @type@ *stack[PYA_QS_STACK], **sptr = stack, *pm, *pi, *pj, *pk; - for(;;) { + for (;;) { while ((size_t)(pr - pl) > SMALL_QUICKSORT*len) { /* quicksort partition */ pm = pl + (((pr - pl)/len) >> 1)*len; - if (@lessthan@(pm, pl, len)) @swap@(pm, pl, len); - if (@lessthan@(pr, pm, len)) @swap@(pr, pm, len); - if (@lessthan@(pm, pl, len)) @swap@(pm, pl, len); - @copy@(vp, pm, len); + if (@TYPE@_LT(pm, pl, len)) @TYPE@_SWAP(pm, pl, len); + if (@TYPE@_LT(pr, pm, len)) @TYPE@_SWAP(pr, pm, len); + if (@TYPE@_LT(pm, pl, len)) @TYPE@_SWAP(pm, pl, len); + @TYPE@_COPY(vp, pm, len); pi = pl; pj = pr - len; - @swap@(pm, pj, len); - for(;;) { - do pi += len; while (@lessthan@(pi, vp, len)); - do pj -= len; while (@lessthan@(vp, pj, len)); - if (pi >= pj) break; - @swap@(pi, pj, len); + @TYPE@_SWAP(pm, pj, len); + for (;;) { + do pi += len; while (@TYPE@_LT(pi, vp, len)); + do pj -= len; while (@TYPE@_LT(vp, pj, len)); + if (pi >= pj) { + break; + } + @TYPE@_SWAP(pi, pj, len); } pk = pr - len; - @swap@(pi, pk, len); + @TYPE@_SWAP(pi, pk, len); /* push largest partition on stack */ if (pi - pl < pr - pi) { *sptr++ = pi + len; @@ -578,18 +703,20 @@ static int } /* insertion sort */ - for(pi = pl + len; pi <= pr; pi += len) { - @copy@(vp, pi, len); + for (pi = pl + len; pi <= pr; pi += len) { + @TYPE@_COPY(vp, pi, len); pj = pi; pk = pi - len; - while (pj > pl && @lessthan@(vp, pk, len)) { - @copy@(pj, pk, len); + while (pj > pl && @TYPE@_LT(vp, pk, len)) { + @TYPE@_COPY(pj, pk, len); pj -= len; pk -= len; } - @copy@(pj, vp, len); + @TYPE@_COPY(pj, vp, len); + } + if (sptr == stack) { + break; } - if (sptr == stack) break; pr = *(--sptr); pl = *(--sptr); } @@ -600,45 +727,47 @@ static int static int -@TYPE@_heapsort(@type@ *start, intp n, PyArrayObject *arr) +@TYPE@_heapsort(@type@ *start, npy_intp n, PyArrayObject *arr) { size_t len = arr->descr->elsize/sizeof(@type@); @type@ *tmp = malloc(arr->descr->elsize); @type@ *a = start - len; - intp i,j,l; + npy_intp i,j,l; for (l = n>>1; l > 0; --l) { - @copy@(tmp, a + l*len, len); + @TYPE@_COPY(tmp, a + l*len, len); for (i = l, j = l<<1; j <= n;) { - if (j < n && @lessthan@(a + j*len, a + (j+1)*len, len)) + if (j < n && @TYPE@_LT(a + j*len, a + (j+1)*len, len)) j += 1; - if (@lessthan@(tmp, a + j*len, len)) { - @copy@(a + i*len, a + j*len, len); + if (@TYPE@_LT(tmp, a + j*len, len)) { + @TYPE@_COPY(a + i*len, a + j*len, len); i = j; j += j; } - else + else { break; + } } - @copy@(a + i*len, tmp, len); + @TYPE@_COPY(a + i*len, tmp, len); } for (; n > 1;) { - @copy@(tmp, a + n*len, len); - @copy@(a + n*len, a + len, len); + @TYPE@_COPY(tmp, a + n*len, len); + @TYPE@_COPY(a + n*len, a + len, len); n -= 1; for (i = 1, j = 2; j <= n;) { - if (j < n && @lessthan@(a + j*len, a + (j+1)*len, len)) + if (j < n && @TYPE@_LT(a + j*len, a + (j+1)*len, len)) j++; - if (@lessthan@(tmp, a + j*len, len)) { - @copy@(a + i*len, a + j*len, len); + if (@TYPE@_LT(tmp, a + j*len, len)) { + @TYPE@_COPY(a + i*len, a + j*len, len); i = j; j += j; } - else + else { break; + } } - @copy@(a + i*len, tmp, len); + @TYPE@_COPY(a + i*len, tmp, len); } free(tmp); @@ -647,10 +776,10 @@ static int static int -@TYPE@_aheapsort(@type@ *v, intp *tosort, intp n, PyArrayObject *arr) +@TYPE@_aheapsort(@type@ *v, npy_intp *tosort, npy_intp n, PyArrayObject *arr) { size_t len = arr->descr->elsize/sizeof(@type@); - intp *a, i,j,l, tmp; + npy_intp *a, i,j,l, tmp; /* The array needs to be offset by one for heapsort indexing */ a = tosort - 1; @@ -658,15 +787,16 @@ static int for (l = n>>1; l > 0; --l) { tmp = a[l]; for (i = l, j = l<<1; j <= n;) { - if (j < n && @lessthan@(v + a[j]*len, v + a[j+1]*len, len)) + if (j < n && @TYPE@_LT(v + a[j]*len, v + a[j+1]*len, len)) j += 1; - if (@lessthan@(v + tmp*len, v + a[j]*len, len)) { + if (@TYPE@_LT(v + tmp*len, v + a[j]*len, len)) { a[i] = a[j]; i = j; j += j; } - else + else { break; + } } a[i] = tmp; } @@ -676,15 +806,16 @@ static int a[n] = a[1]; n -= 1; for (i = 1, j = 2; j <= n;) { - if (j < n && @lessthan@(v + a[j]*len, v + a[j+1]*len, len)) + if (j < n && @TYPE@_LT(v + a[j]*len, v + a[j+1]*len, len)) j++; - if (@lessthan@(v + tmp*len, v + a[j]*len, len)) { + if (@TYPE@_LT(v + tmp*len, v + a[j]*len, len)) { a[i] = a[j]; i = j; j += j; } - else + else { break; + } } a[i] = tmp; } @@ -694,35 +825,37 @@ static int static int -@TYPE@_aquicksort(@type@ *v, intp* tosort, intp num, PyArrayObject *arr) +@TYPE@_aquicksort(@type@ *v, npy_intp* tosort, npy_intp num, PyArrayObject *arr) { size_t len = arr->descr->elsize/sizeof(@type@); @type@ *vp; - intp *pl = tosort; - intp *pr = tosort + num - 1; - intp *stack[PYA_QS_STACK]; - intp **sptr=stack; - intp *pm, *pi, *pj, *pk, vi, SWAP_temp; + npy_intp *pl = tosort; + npy_intp *pr = tosort + num - 1; + npy_intp *stack[PYA_QS_STACK]; + npy_intp **sptr=stack; + npy_intp *pm, *pi, *pj, *pk, vi; - for(;;) { + for (;;) { while ((pr - pl) > SMALL_QUICKSORT) { /* quicksort partition */ pm = pl + ((pr - pl) >> 1); - if (@lessthan@(v + (*pm)*len, v + (*pl)*len, len)) SWAP(*pm, *pl); - if (@lessthan@(v + (*pr)*len, v + (*pm)*len, len)) SWAP(*pr, *pm); - if (@lessthan@(v + (*pm)*len, v + (*pl)*len, len)) SWAP(*pm, *pl); + if (@TYPE@_LT(v + (*pm)*len, v + (*pl)*len, len)) INTP_SWAP(*pm, *pl); + if (@TYPE@_LT(v + (*pr)*len, v + (*pm)*len, len)) INTP_SWAP(*pr, *pm); + if (@TYPE@_LT(v + (*pm)*len, v + (*pl)*len, len)) INTP_SWAP(*pm, *pl); vp = v + (*pm)*len; pi = pl; pj = pr - 1; - SWAP(*pm,*pj); - for(;;) { - do ++pi; while (@lessthan@(v + (*pi)*len, vp, len)); - do --pj; while (@lessthan@(vp, v + (*pj)*len, len)); - if (pi >= pj) break; - SWAP(*pi,*pj); + INTP_SWAP(*pm,*pj); + for (;;) { + do ++pi; while (@TYPE@_LT(v + (*pi)*len, vp, len)); + do --pj; while (@TYPE@_LT(vp, v + (*pj)*len, len)); + if (pi >= pj) { + break; + } + INTP_SWAP(*pi,*pj); } pk = pr - 1; - SWAP(*pi,*pk); + INTP_SWAP(*pi,*pk); /* push largest partition on stack */ if (pi - pl < pr - pi) { *sptr++ = pi + 1; @@ -737,17 +870,19 @@ static int } /* insertion sort */ - for(pi = pl + 1; pi <= pr; ++pi) { + for (pi = pl + 1; pi <= pr; ++pi) { vi = *pi; vp = v + vi*len; pj = pi; pk = pi - 1; - while (pj > pl && @lessthan@(vp, v + (*pk)*len, len)) { + while (pj > pl && @TYPE@_LT(vp, v + (*pk)*len, len)) { *pj-- = *pk--; } *pj = vi; } - if (sptr == stack) break; + if (sptr == stack) { + break; + } pr = *(--sptr); pl = *(--sptr); } @@ -757,26 +892,26 @@ static int static void -@TYPE@_amergesort0(intp *pl, intp *pr, @type@ *v, intp *pw, int len) +@TYPE@_amergesort0(npy_intp *pl, npy_intp *pr, @type@ *v, npy_intp *pw, int len) { @type@ *vp; - intp vi, *pi, *pj, *pk, *pm; + npy_intp vi, *pi, *pj, *pk, *pm; if (pr - pl > SMALL_MERGESORT) { /* merge sort */ pm = pl + ((pr - pl) >> 1); @TYPE@_amergesort0(pl,pm,v,pw,len); @TYPE@_amergesort0(pm,pr,v,pw,len); - for(pi = pw, pj = pl; pj < pm;) { + for (pi = pw, pj = pl; pj < pm;) { *pi++ = *pj++; } pj = pw; pk = pl; while (pj < pi && pm < pr) { - if (@lessequal@(v + (*pj)*len, v + (*pm)*len, len)) { - *pk = *pj++; - } else { + if (@TYPE@_LT(v + (*pm)*len, v + (*pj)*len, len)) { *pk = *pm++; + } else { + *pk = *pj++; } pk++; } @@ -785,12 +920,12 @@ static void } } else { /* insertion sort */ - for(pi = pl + 1; pi < pr; ++pi) { + for (pi = pl + 1; pi < pr; ++pi) { vi = *pi; vp = v + vi*len; pj = pi; pk = pi -1; - while (pj > pl && @lessthan@(vp, v + (*pk)*len, len)) { + while (pj > pl && @TYPE@_LT(vp, v + (*pk)*len, len)) { *pj-- = *pk--; } *pj = vi; @@ -800,11 +935,11 @@ static void static int -@TYPE@_amergesort(@type@ *v, intp *tosort, intp num, PyArrayObject *arr) +@TYPE@_amergesort(@type@ *v, npy_intp *tosort, npy_intp num, PyArrayObject *arr) { const size_t elsize = arr->descr->elsize; const size_t len = elsize / sizeof(@type@); - intp *pl, *pr, *pw; + npy_intp *pl, *pr, *pw; pl = tosort; pr = pl + num; @@ -826,20 +961,23 @@ add_sortfuncs(void) PyArray_Descr *descr; /**begin repeat - #TYPE=BOOL,BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE,CFLOAT,CDOUBLE,CLONGDOUBLE,STRING,UNICODE,DATETIME,TIMEDELTA# - **/ + * + * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG, + * LONGLONG, ULONGLONG, FLOAT, DOUBLE, LONGDOUBLE, + * CFLOAT, CDOUBLE, CLONGDOUBLE, STRING, UNICODE# + */ descr = PyArray_DescrFromType(PyArray_@TYPE@); - descr->f->sort[PyArray_QUICKSORT] = \ + descr->f->sort[PyArray_QUICKSORT] = (PyArray_SortFunc *)@TYPE@_quicksort; - descr->f->argsort[PyArray_QUICKSORT] = \ + descr->f->argsort[PyArray_QUICKSORT] = (PyArray_ArgSortFunc *)@TYPE@_aquicksort; - descr->f->sort[PyArray_HEAPSORT] = \ + descr->f->sort[PyArray_HEAPSORT] = (PyArray_SortFunc *)@TYPE@_heapsort; - descr->f->argsort[PyArray_HEAPSORT] = \ + descr->f->argsort[PyArray_HEAPSORT] = (PyArray_ArgSortFunc *)@TYPE@_aheapsort; - descr->f->sort[PyArray_MERGESORT] = \ + descr->f->sort[PyArray_MERGESORT] = (PyArray_SortFunc *)@TYPE@_mergesort; - descr->f->argsort[PyArray_MERGESORT] = \ + descr->f->argsort[PyArray_MERGESORT] = (PyArray_ArgSortFunc *)@TYPE@_amergesort; /**end repeat**/ diff --git a/numpy/core/src/multiarray/arrayobject.c b/numpy/core/src/multiarray/arrayobject.c index dccd4984f..2f49e03a5 100644 --- a/numpy/core/src/multiarray/arrayobject.c +++ b/numpy/core/src/multiarray/arrayobject.c @@ -1099,16 +1099,15 @@ PyArray_CheckStrides(int elsize, int nd, intp numbytes, intp offset, static PyObject * array_new(PyTypeObject *subtype, PyObject *args, PyObject *kwds) { - static char *kwlist[] = {"shape", "dtype", "buffer", - "offset", "strides", + static char *kwlist[] = {"shape", "dtype", "buffer", "offset", "strides", "order", NULL}; - PyArray_Descr *descr=NULL; + PyArray_Descr *descr = NULL; int itemsize; PyArray_Dims dims = {NULL, 0}; PyArray_Dims strides = {NULL, 0}; PyArray_Chunk buffer; - longlong offset=0; - NPY_ORDER order=PyArray_CORDER; + longlong offset = 0; + NPY_ORDER order = PyArray_CORDER; int fortran = 0; PyArrayObject *ret; @@ -1268,73 +1267,65 @@ array_alloc(PyTypeObject *type, Py_ssize_t NPY_UNUSED(nitems)) NPY_NO_EXPORT PyTypeObject PyArray_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size */ - "numpy.ndarray", /* tp_name */ - sizeof(PyArrayObject), /* tp_basicsize */ - 0, /* tp_itemsize */ + 0, /* ob_size */ +#endif + "numpy.ndarray", /* tp_name */ + sizeof(PyArrayObject), /* tp_basicsize */ + 0, /* tp_itemsize */ /* methods */ - (destructor)array_dealloc, /* tp_dealloc */ - (printfunc)NULL, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - (cmpfunc)0, /* tp_compare */ - (reprfunc)array_repr, /* tp_repr */ - &array_as_number, /* tp_as_number */ - &array_as_sequence, /* tp_as_sequence */ - &array_as_mapping, /* tp_as_mapping */ - (hashfunc)0, /* tp_hash */ - (ternaryfunc)0, /* tp_call */ - (reprfunc)array_str, /* tp_str */ - (getattrofunc)0, /* tp_getattro */ - (setattrofunc)0, /* tp_setattro */ - &array_as_buffer, /* tp_as_buffer */ + (destructor)array_dealloc, /* tp_dealloc */ + (printfunc)NULL, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ +#endif + (reprfunc)array_repr, /* tp_repr */ + &array_as_number, /* tp_as_number */ + &array_as_sequence, /* tp_as_sequence */ + &array_as_mapping, /* tp_as_mapping */ + (hashfunc)0, /* tp_hash */ + (ternaryfunc)0, /* tp_call */ + (reprfunc)array_str, /* tp_str */ + (getattrofunc)0, /* tp_getattro */ + (setattrofunc)0, /* tp_setattro */ + &array_as_buffer, /* tp_as_buffer */ (Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE - | Py_TPFLAGS_CHECKTYPES), /* tp_flags */ - /*Documentation string */ - 0, /* tp_doc */ - - (traverseproc)0, /* tp_traverse */ - (inquiry)0, /* tp_clear */ - (richcmpfunc)array_richcompare, /* tp_richcompare */ - offsetof(PyArrayObject, weakreflist), /* tp_weaklistoffset */ - - /* Iterator support (use standard) */ - - (getiterfunc)array_iter, /* tp_iter */ - (iternextfunc)0, /* tp_iternext */ - - /* Sub-classing (new-style object) support */ - - array_methods, /* tp_methods */ - 0, /* tp_members */ - array_getsetlist, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - (initproc)0, /* tp_init */ - array_alloc, /* tp_alloc */ - (newfunc)array_new, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ - -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ -#endif + | Py_TPFLAGS_CHECKTYPES), /* tp_flags */ + 0, /* tp_doc */ + + (traverseproc)0, /* tp_traverse */ + (inquiry)0, /* tp_clear */ + (richcmpfunc)array_richcompare, /* tp_richcompare */ + offsetof(PyArrayObject, weakreflist), /* tp_weaklistoffset */ + (getiterfunc)array_iter, /* tp_iter */ + (iternextfunc)0, /* tp_iternext */ + array_methods, /* tp_methods */ + 0, /* tp_members */ + array_getsetlist, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + (initproc)0, /* tp_init */ + array_alloc, /* tp_alloc */ + (newfunc)array_new, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; diff --git a/numpy/core/src/multiarray/arraytypes.c.src b/numpy/core/src/multiarray/arraytypes.c.src index 7080fbe7a..420fcea7d 100644 --- a/numpy/core/src/multiarray/arraytypes.c.src +++ b/numpy/core/src/multiarray/arraytypes.c.src @@ -2147,7 +2147,13 @@ VOID_nonzero (char *ip, PyArrayObject *ap) #undef __ALIGNED -/****************** compare **********************************/ +/* + ***************************************************************************** + ** COMPARISON FUNCTIONS ** + ***************************************************************************** + */ + +/* boolean type */ static int BOOL_compare(Bool *ip1, Bool *ip2, PyArrayObject *NPY_UNUSED(ap)) @@ -2155,48 +2161,144 @@ BOOL_compare(Bool *ip1, Bool *ip2, PyArrayObject *NPY_UNUSED(ap)) return (*ip1 ? (*ip2 ? 0 : 1) : (*ip2 ? -1 : 0)); } + +/* integer types */ + /**begin repeat -#fname=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE,DATETIME,TIMEDELTA# -#type=byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble, datetime, timedelta# -*/ + * #TYPE = BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG, + * LONGLONG, ULONGLONG, DATETIME, TIMEDELTA# + * #type = byte, ubyte, short, ushort, int, uint, long, ulong, + * longlong, ulonglong, datetime, timedelta# + */ static int -@fname@_compare (@type@ *ip1, @type@ *ip2, PyArrayObject *NPY_UNUSED(ap)) +@TYPE@_compare (@type@ *pa, @type@ *pb, PyArrayObject *NPY_UNUSED(ap)) { - return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; + const @type@ a = *pa; + const @type@ b = *pb; + + return a < b ? -1 : a == b ? 0 : 1; } /**end repeat**/ -/* compare imaginary part first, then complex if equal imaginary */ + +/* float types */ + +/* + * The real/complex comparison functions are compatible with the new sort + * order for nans introduced in numpy 1.4.0. All nan values now compare + * larger than non-nan values and are sorted to the end. The comparison + * order is: + * + * Real: [R, nan] + * Complex: [R + Rj, R + nanj, nan + Rj, nan + nanj] + * + * where complex values with the same nan placements are sorted according + * to the non-nan part if it exists. If both the real and imaginary parts + * of complex types are non-nan the order is the same as the real parts + * unless they happen to be equal, in which case the order is that of the + * imaginary parts. + */ + /**begin repeat -#fname=CFLOAT, CDOUBLE, CLONGDOUBLE# -#type= float, double, longdouble# -*/ + * #TYPE = FLOAT, DOUBLE, LONGDOUBLE# + * #type = float, double, longdouble# + */ + +#define LT(a,b) ((a) < (b) || ((b) != (b) && (a) ==(a))) static int -@fname@_compare (@type@ *ip1, @type@ *ip2, PyArrayObject *NPY_UNUSED(ap)) +@TYPE@_compare(@type@ *pa, @type@ *pb) { - if (*ip1 == *ip2) { - return ip1[1]<ip2[1] ? -1 : (ip1[1] == ip2[1] ? 0 : 1); + const @type@ a = *pa; + const @type@ b = *pb; + int ret; + + if (LT(a,b)) { + ret = -1; + } + else if (LT(b,a)) { + ret = 1; + } + else { + ret = 0; + } + return ret; +} + + +static int +C@TYPE@_compare(@type@ *pa, @type@ *pb) +{ + const @type@ ar = pa[0]; + const @type@ ai = pa[1]; + const @type@ br = pb[0]; + const @type@ bi = pb[1]; + int ret; + + if (ar < br) { + if (ai == ai || bi != bi) { + ret = -1; + } + else { + ret = 1; + } + } + else if (br < ar) { + if (bi == bi || ai != ai) { + ret = 1; + } + else { + ret = -1; + } + } + else if (ar == br || (ar != ar && br != br)) { + if (LT(ai,bi)) { + ret = -1; + } + else if (LT(bi,ai)) { + ret = 1; + } + else { + ret = 0; + } + } + else if (ar == ar) { + ret = -1; } else { - return *ip1 < *ip2 ? -1 : 1; + ret = 1; } + + return ret; } - /**end repeat**/ + +#undef LT + +/**end repeat**/ + + +/* object type */ static int OBJECT_compare(PyObject **ip1, PyObject **ip2, PyArrayObject *NPY_UNUSED(ap)) { if ((*ip1 == NULL) || (*ip2 == NULL)) { - if (ip1 == ip2) return 1; - if (ip1 == NULL) return -1; + if (ip1 == ip2) { + return 1; + } + if (ip1 == NULL) { + return -1; + } return 1; } return PyObject_Compare(*ip1, *ip2); } + +/* string type */ + static int STRING_compare(char *ip1, char *ip2, PyArrayObject *ap) { @@ -2213,33 +2315,38 @@ STRING_compare(char *ip1, char *ip2, PyArrayObject *ap) return 0; } -/* taken from Python */ + +/* unicode type */ + static int UNICODE_compare(PyArray_UCS4 *ip1, PyArray_UCS4 *ip2, PyArrayObject *ap) { - int itemsize=ap->descr->elsize; - PyArray_UCS4 c1, c2; - - if (itemsize < 0) return 0; + int itemsize = ap->descr->elsize; + if (itemsize < 0) { + return 0; + } while(itemsize-- > 0) { - c1 = *ip1++; - c2 = *ip2++; - - if (c1 != c2) + PyArray_UCS4 c1 = *ip1++; + PyArray_UCS4 c2 = *ip2++; + if (c1 != c2) { return (c1 < c2) ? -1 : 1; + } } return 0; } -/* If fields are defined, then compare on first field and if equal - compare on second field. Continue until done or comparison results - in not_equal. - Must align data passed on to sub-comparisons. -*/ +/* void type */ +/* + * If fields are defined, then compare on first field and if equal + * compare on second field. Continue until done or comparison results + * in not_equal. + * + * Must align data passed on to sub-comparisons. + */ static int VOID_compare(char *ip1, char *ip2, PyArrayObject *ap) { @@ -2247,17 +2354,18 @@ VOID_compare(char *ip1, char *ip2, PyArrayObject *ap) PyObject *names, *key; PyObject *tup, *title; char *nip1, *nip2; - int i, offset, res=0; + int i, offset, res = 0; - if (!PyArray_HASFIELDS(ap)) + if (!PyArray_HASFIELDS(ap)) { return STRING_compare(ip1, ip2, ap); - + } descr = ap->descr; - /* Compare on the first-field. If equal, then - compare on the second-field, etc. + /* + * Compare on the first-field. If equal, then + * compare on the second-field, etc. */ names = descr->names; - for (i=0; i<PyTuple_GET_SIZE(names); i++) { + for (i = 0; i < PyTuple_GET_SIZE(names); i++) { key = PyTuple_GET_ITEM(names, i); tup = PyDict_GetItem(descr->fields, key); if (!PyArg_ParseTuple(tup, "Oi|O", &new, &offset, @@ -2271,15 +2379,18 @@ VOID_compare(char *ip1, char *ip2, PyArrayObject *ap) if (((intp)(nip1) % new->alignment) != 0) { /* create buffer and copy */ nip1 = _pya_malloc(new->elsize); - if (nip1 == NULL) goto finish; + if (nip1 == NULL) { + goto finish; + } memcpy(nip1, ip1+offset, new->elsize); } if (((intp)(nip2) % new->alignment) != 0) { /* copy data to a buffer */ nip2 = _pya_malloc(new->elsize); if (nip2 == NULL) { - if (nip1 != ip1+offset) + if (nip1 != ip1+offset) { _pya_free(nip1); + } goto finish; } memcpy(nip2, ip2+offset, new->elsize); @@ -2294,7 +2405,9 @@ VOID_compare(char *ip1, char *ip2, PyArrayObject *ap) _pya_free(nip2); } } - if (res != 0) break; + if (res != 0) { + break; + } } finish: diff --git a/numpy/core/src/multiarray/ctors.c b/numpy/core/src/multiarray/ctors.c index 92723dbf1..e4873913a 100644 --- a/numpy/core/src/multiarray/ctors.c +++ b/numpy/core/src/multiarray/ctors.c @@ -2978,16 +2978,16 @@ PyArray_FromFile(FILE *fp, PyArray_Descr *dtype, intp num, char *sep) { PyArrayObject *ret; size_t nread = 0; - char *tmp; if (PyDataType_REFCHK(dtype)) { PyErr_SetString(PyExc_ValueError, - "cannot read into object array"); + "Cannot read into object array"); Py_DECREF(dtype); return NULL; } if (dtype->elsize == 0) { - PyErr_SetString(PyExc_ValueError, "0-sized elements."); + PyErr_SetString(PyExc_ValueError, + "The elements are 0-sized."); Py_DECREF(dtype); return NULL; } @@ -2997,28 +2997,20 @@ PyArray_FromFile(FILE *fp, PyArray_Descr *dtype, intp num, char *sep) else { if (dtype->f->scanfunc == NULL) { PyErr_SetString(PyExc_ValueError, - "don't know how to read " \ - "character files with that " \ - "array type"); + "Unable to read character files of that array type"); Py_DECREF(dtype); return NULL; } - ret = array_from_text(dtype, num, sep, &nread, - fp, - (next_element) fromfile_next_element, - (skip_separator) fromfile_skip_separator, - NULL); + ret = array_from_text(dtype, num, sep, &nread, fp, + (next_element) fromfile_next_element, + (skip_separator) fromfile_skip_separator, NULL); } if (((intp) nread) < num) { - fprintf(stderr, "%ld items requested but only %ld read\n", - (long) num, (long) nread); - /* Make sure realloc is > 0 */ - tmp = PyDataMem_RENEW(ret->data, - NPY_MAX(nread,1) * ret->descr->elsize); - /* FIXME: This should not raise a memory error when nread == 0 - We should return an empty array or at least raise an EOF Error. - */ - if ((tmp == NULL) || (nread == 0)) { + /* Realloc memory for smaller number of elements */ + const size_t nsize = NPY_MAX(nread,1)*ret->descr->elsize; + char *tmp; + + if((tmp = PyDataMem_RENEW(ret->data, nsize)) == NULL) { Py_DECREF(ret); return PyErr_NoMemory(); } diff --git a/numpy/core/src/multiarray/descriptor.c b/numpy/core/src/multiarray/descriptor.c index 909af2243..edcab8b09 100644 --- a/numpy/core/src/multiarray/descriptor.c +++ b/numpy/core/src/multiarray/descriptor.c @@ -2626,62 +2626,62 @@ static PyMappingMethods descr_as_mapping = { }; /****************** End of Mapping Protocol ******************************/ + NPY_NO_EXPORT PyTypeObject PyArrayDescr_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size */ - "numpy.dtype", /* tp_name */ - sizeof(PyArray_Descr), /* tp_basicsize */ - 0, /* tp_itemsize */ + 0, /* ob_size */ +#endif + "numpy.dtype", /* tp_name */ + sizeof(PyArray_Descr), /* tp_basicsize */ + 0, /* tp_itemsize */ /* methods */ - (destructor)arraydescr_dealloc, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_compare */ - (reprfunc)arraydescr_repr, /* tp_repr */ - 0, /* tp_as_number */ - &descr_as_sequence, /* tp_as_sequence */ - &descr_as_mapping, /* tp_as_mapping */ - 0, /* tp_hash */ - 0, /* tp_call */ - (reprfunc)arraydescr_str, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - Py_TPFLAGS_DEFAULT, /* tp_flags */ - 0, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - (richcmpfunc)arraydescr_richcompare, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - arraydescr_methods, /* tp_methods */ - arraydescr_members, /* tp_members */ - arraydescr_getsets, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - arraydescr_new, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ - -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + (destructor)arraydescr_dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + (void *)0, /* tp_reserved */ +#else + 0, /* tp_compare */ #endif + (reprfunc)arraydescr_repr, /* tp_repr */ + 0, /* tp_as_number */ + &descr_as_sequence, /* tp_as_sequence */ + &descr_as_mapping, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + (reprfunc)arraydescr_str, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + Py_TPFLAGS_DEFAULT, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + (richcmpfunc)arraydescr_richcompare, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + arraydescr_methods, /* tp_methods */ + arraydescr_members, /* tp_members */ + arraydescr_getsets, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + arraydescr_new, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; diff --git a/numpy/core/src/multiarray/flagsobject.c b/numpy/core/src/multiarray/flagsobject.c index 869ed613e..38121f910 100644 --- a/numpy/core/src/multiarray/flagsobject.c +++ b/numpy/core/src/multiarray/flagsobject.c @@ -553,61 +553,60 @@ arrayflags_new(PyTypeObject *NPY_UNUSED(self), PyObject *args, PyObject *NPY_UNU } NPY_NO_EXPORT PyTypeObject PyArrayFlags_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, + 0, /* ob_size */ +#endif "numpy.flagsobj", sizeof(PyArrayFlagsObject), - 0, /* tp_itemsize */ + 0, /* tp_itemsize */ /* methods */ - (destructor)arrayflags_dealloc, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - (cmpfunc)arrayflags_compare, /* tp_compare */ - (reprfunc)arrayflags_print, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - &arrayflags_as_mapping, /* tp_as_mapping */ - 0, /* tp_hash */ - 0, /* tp_call */ - (reprfunc)arrayflags_print, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - Py_TPFLAGS_DEFAULT, /* tp_flags */ - 0, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - 0, /* tp_methods */ - 0, /* tp_members */ - arrayflags_getsets, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - arrayflags_new, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ - -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + (destructor)arrayflags_dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + (cmpfunc)arrayflags_compare, /* tp_compare */ #endif + (reprfunc)arrayflags_print, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + &arrayflags_as_mapping, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + (reprfunc)arrayflags_print, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + Py_TPFLAGS_DEFAULT, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + arrayflags_getsets, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + arrayflags_new, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; diff --git a/numpy/core/src/multiarray/global.c b/numpy/core/src/multiarray/global.c deleted file mode 100644 index 22306da23..000000000 --- a/numpy/core/src/multiarray/global.c +++ /dev/null @@ -1,3 +0,0 @@ -#include "config.h" - -NPY_NO_EXPORT int NPY_NUMUSERTYPES = 0; diff --git a/numpy/core/src/multiarray/iterators.c b/numpy/core/src/multiarray/iterators.c index 6f25d9432..d6eac84ef 100644 --- a/numpy/core/src/multiarray/iterators.c +++ b/numpy/core/src/multiarray/iterators.c @@ -266,27 +266,33 @@ slice_GetIndices(PySliceObject *r, intp length, /* Aided by Peter J. Verveer's nd_image package and numpy's arraymap ****/ /* and Python's array iterator ***/ -/*NUMPY_API - * Get Iterator. - */ -NPY_NO_EXPORT PyObject * -PyArray_IterNew(PyObject *obj) +/* get the dataptr from its current coordinates for simple iterator */ +static char* +get_ptr_simple(PyArrayIterObject* iter, npy_intp *coordinates) { - PyArrayIterObject *it; - int i, nd; - PyArrayObject *ao = (PyArrayObject *)obj; + npy_intp i; + char *ret; - if (!PyArray_Check(ao)) { - PyErr_BadInternalCall(); - return NULL; - } + ret = iter->ao->data; - it = (PyArrayIterObject *)_pya_malloc(sizeof(PyArrayIterObject)); - PyObject_Init((PyObject *)it, &PyArrayIter_Type); - /* it = PyObject_New(PyArrayIterObject, &PyArrayIter_Type);*/ - if (it == NULL) { - return NULL; + for(i = 0; i < iter->ao->nd; ++i) { + ret += coordinates[i] * iter->strides[i]; } + + return ret; +} + +/* + * This is common initialization code between PyArrayIterObject and + * PyArrayNeighborhoodIterObject + * + * Increase ao refcount + */ +static PyObject * +array_iter_base_init(PyArrayIterObject *it, PyArrayObject *ao) +{ + int nd, i; + nd = ao->nd; PyArray_UpdateFlags(ao, CONTIGUOUS); if (PyArray_ISCONTIGUOUS(ao)) { @@ -307,12 +313,50 @@ PyArray_IterNew(PyObject *obj) if (i > 0) { it->factors[nd-i-1] = it->factors[nd-i] * ao->dimensions[nd-i]; } + it->bounds[i][0] = 0; + it->bounds[i][1] = ao->dimensions[i] - 1; + it->limits[i][0] = 0; + it->limits[i][1] = ao->dimensions[i] - 1; + it->limits_sizes[i] = it->limits[i][1] - it->limits[i][0] + 1; } + + it->translate = &get_ptr_simple; PyArray_ITER_RESET(it); return (PyObject *)it; } +static void +array_iter_base_dealloc(PyArrayIterObject *it) +{ + Py_XDECREF(it->ao); +} + +/*NUMPY_API + * Get Iterator. + */ +NPY_NO_EXPORT PyObject * +PyArray_IterNew(PyObject *obj) +{ + PyArrayIterObject *it; + PyArrayObject *ao = (PyArrayObject *)obj; + + if (!PyArray_Check(ao)) { + PyErr_BadInternalCall(); + return NULL; + } + + it = (PyArrayIterObject *)_pya_malloc(sizeof(PyArrayIterObject)); + PyObject_Init((PyObject *)it, &PyArrayIter_Type); + /* it = PyObject_New(PyArrayIterObject, &PyArrayIter_Type);*/ + if (it == NULL) { + return NULL; + } + + array_iter_base_init(it, ao); + return (PyObject *)it; +} + /*NUMPY_API * Get Iterator broadcast to a particular shape */ @@ -502,7 +546,7 @@ arrayiter_next(PyArrayIterObject *it) static void arrayiter_dealloc(PyArrayIterObject *it) { - Py_XDECREF(it->ao); + array_iter_base_dealloc(it); _pya_free(it); } @@ -1174,8 +1218,12 @@ iter_coords_get(PyArrayIterObject *self) int i; val = self->index; for (i = 0; i < nd; i++) { - self->coordinates[i] = val / self->factors[i]; - val = val % self->factors[i]; + if (self->factors[i] != 0) { + self->coordinates[i] = val / self->factors[i]; + val = val % self->factors[i]; + } else { + self->coordinates[i] = 0; + } } } return PyArray_IntTupleFromIntp(nd, self->coordinates); @@ -1189,63 +1237,62 @@ static PyGetSetDef iter_getsets[] = { }; NPY_NO_EXPORT PyTypeObject PyArrayIter_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size */ - "numpy.flatiter", /* tp_name */ - sizeof(PyArrayIterObject), /* tp_basicsize */ - 0, /* tp_itemsize */ + 0, /* ob_size */ +#endif + "numpy.flatiter", /* tp_name */ + sizeof(PyArrayIterObject), /* tp_basicsize */ + 0, /* tp_itemsize */ /* methods */ - (destructor)arrayiter_dealloc, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_compare */ - 0, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - &iter_as_mapping, /* tp_as_mapping */ - 0, /* tp_hash */ - 0, /* tp_call */ - 0, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - Py_TPFLAGS_DEFAULT, /* tp_flags */ - 0, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - (richcmpfunc)iter_richcompare, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - (iternextfunc)arrayiter_next, /* tp_iternext */ - iter_methods, /* tp_methods */ - iter_members, /* tp_members */ - iter_getsets, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - 0, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + (destructor)arrayiter_dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ #endif - + 0, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + &iter_as_mapping, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + Py_TPFLAGS_DEFAULT, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + (richcmpfunc)iter_richcompare, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + (iternextfunc)arrayiter_next, /* tp_iternext */ + iter_methods, /* tp_methods */ + iter_members, /* tp_members */ + iter_getsets, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; /** END of Array Iterator **/ @@ -1652,61 +1699,387 @@ static PyMethodDef arraymultiter_methods[] = { }; NPY_NO_EXPORT PyTypeObject PyArrayMultiIter_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size */ - "numpy.broadcast", /* tp_name */ - sizeof(PyArrayMultiIterObject), /* tp_basicsize */ - 0, /* tp_itemsize */ + 0, /* ob_size */ +#endif + "numpy.broadcast", /* tp_name */ + sizeof(PyArrayMultiIterObject), /* tp_basicsize */ + 0, /* tp_itemsize */ /* methods */ - (destructor)arraymultiter_dealloc, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_compare */ - 0, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - 0, /* tp_as_mapping */ - 0, /* tp_hash */ - 0, /* tp_call */ - 0, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - Py_TPFLAGS_DEFAULT, /* tp_flags */ - 0, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - (iternextfunc)arraymultiter_next, /* tp_iternext */ - arraymultiter_methods, /* tp_methods */ - arraymultiter_members, /* tp_members */ - arraymultiter_getsetlist, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - (initproc)0, /* tp_init */ - 0, /* tp_alloc */ - arraymultiter_new, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ - -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + (destructor)arraymultiter_dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ +#endif + 0, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + Py_TPFLAGS_DEFAULT, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + (iternextfunc)arraymultiter_next, /* tp_iternext */ + arraymultiter_methods, /* tp_methods */ + arraymultiter_members, /* tp_members */ + arraymultiter_getsetlist, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + (initproc)0, /* tp_init */ + 0, /* tp_alloc */ + arraymultiter_new, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ +}; + +/*========================= Neighborhood iterator ======================*/ + +static void neighiter_dealloc(PyArrayNeighborhoodIterObject* iter); + +static char* _set_constant(PyArrayNeighborhoodIterObject* iter, + PyArrayObject *fill) +{ + char *ret; + PyArrayIterObject *ar = iter->_internal_iter; + int storeflags, st; + + ret = PyDataMem_NEW(ar->ao->descr->elsize); + if (ret == NULL) { + PyErr_SetNone(PyExc_MemoryError); + return NULL; + } + + if (PyArray_ISOBJECT(ar->ao)) { + memcpy(ret, fill->data, sizeof(PyObject*)); + Py_INCREF(*(PyObject**)ret); + } else { + /* Non-object types */ + + storeflags = ar->ao->flags; + ar->ao->flags |= BEHAVED; + st = ar->ao->descr->f->setitem((PyObject*)fill, ret, ar->ao); + ar->ao->flags = storeflags; + + if (st < 0) { + PyDataMem_FREE(ret); + return NULL; + } + } + + return ret; +} + +#define _INF_SET_PTR(c) \ + bd = coordinates[c] + p->coordinates[c]; \ + if (bd < p->limits[c][0] || bd > p->limits[c][1]) { \ + return niter->constant; \ + } \ + _coordinates[c] = bd; + +/* set the dataptr from its current coordinates */ +static char* +get_ptr_constant(PyArrayIterObject* _iter, npy_intp *coordinates) +{ + int i; + npy_intp bd, _coordinates[NPY_MAXDIMS]; + PyArrayNeighborhoodIterObject *niter = (PyArrayNeighborhoodIterObject*)_iter; + PyArrayIterObject *p = niter->_internal_iter; + + for(i = 0; i < niter->nd; ++i) { + _INF_SET_PTR(i) + } + + return p->translate(p, _coordinates); +} +#undef _INF_SET_PTR + +#define _NPY_IS_EVEN(x) ((x) % 2 == 0) + +/* For an array x of dimension n, and given index i, returns j, 0 <= j < n + * such as x[i] = x[j], with x assumed to be mirrored. For example, for x = + * {1, 2, 3} (n = 3) + * + * index -5 -4 -3 -2 -1 0 1 2 3 4 5 6 + * value 2 3 3 2 1 1 2 3 3 2 1 1 + * + * _npy_pos_index_mirror(4, 3) will return 1, because x[4] = x[1]*/ +static inline npy_intp +__npy_pos_remainder(npy_intp i, npy_intp n) +{ + npy_intp k, l, j; + + /* Mirror i such as it is guaranteed to be positive */ + if (i < 0) { + i = - i - 1; + } + + /* compute k and l such as i = k * n + l, 0 <= l < k */ + k = i / n; + l = i - k * n; + + if (_NPY_IS_EVEN(k)) { + j = l; + } else { + j = n - 1 - l; + } + return j; +} +#undef _NPY_IS_EVEN + +#define _INF_SET_PTR_MIRROR(c) \ + lb = p->limits[c][0]; \ + bd = coordinates[c] + p->coordinates[c] - lb; \ + _coordinates[c] = lb + __npy_pos_remainder(bd, p->limits_sizes[c]); + +/* set the dataptr from its current coordinates */ +static char* +get_ptr_mirror(PyArrayIterObject* _iter, npy_intp *coordinates) +{ + int i; + npy_intp bd, _coordinates[NPY_MAXDIMS], lb; + PyArrayNeighborhoodIterObject *niter = (PyArrayNeighborhoodIterObject*)_iter; + PyArrayIterObject *p = niter->_internal_iter; + + for(i = 0; i < niter->nd; ++i) { + _INF_SET_PTR_MIRROR(i) + } + + return p->translate(p, _coordinates); +} +#undef _INF_SET_PTR_MIRROR + +/* compute l such as i = k * n + l, 0 <= l < |k| */ +static inline npy_intp +__npy_euclidean_division(npy_intp i, npy_intp n) +{ + npy_intp l; + + l = i % n; + if (l < 0) { + l += n; + } + return l; +} + +#define _INF_SET_PTR_CIRCULAR(c) \ + lb = p->limits[c][0]; \ + bd = coordinates[c] + p->coordinates[c] - lb; \ + _coordinates[c] = lb + __npy_euclidean_division(bd, p->limits_sizes[c]); + +static char* +get_ptr_circular(PyArrayIterObject* _iter, npy_intp *coordinates) +{ + int i; + npy_intp bd, _coordinates[NPY_MAXDIMS], lb; + PyArrayNeighborhoodIterObject *niter = (PyArrayNeighborhoodIterObject*)_iter; + PyArrayIterObject *p = niter->_internal_iter; + + for(i = 0; i < niter->nd; ++i) { + _INF_SET_PTR_CIRCULAR(i) + } + return p->translate(p, _coordinates); +} + +#undef _INF_SET_PTR_CIRCULAR + +/* + * fill and x->ao should have equivalent types + */ +/*NUMPY_API*/ +NPY_NO_EXPORT PyObject* +PyArray_NeighborhoodIterNew(PyArrayIterObject *x, intp *bounds, + int mode, PyArrayObject* fill) +{ + int i; + PyArrayNeighborhoodIterObject *ret; + + ret = _pya_malloc(sizeof(*ret)); + if (ret == NULL) { + return NULL; + } + PyObject_Init((PyObject *)ret, &PyArrayNeighborhoodIter_Type); + + array_iter_base_init((PyArrayIterObject*)ret, x->ao); + Py_INCREF(x); + ret->_internal_iter = x; + + ret->nd = x->ao->nd; + + for (i = 0; i < ret->nd; ++i) { + ret->dimensions[i] = x->ao->dimensions[i]; + } + + /* Compute the neighborhood size and copy the shape */ + ret->size = 1; + for (i = 0; i < ret->nd; ++i) { + ret->bounds[i][0] = bounds[2 * i]; + ret->bounds[i][1] = bounds[2 * i + 1]; + ret->size *= (ret->bounds[i][1] - ret->bounds[i][0]) + 1; + + /* limits keep track of valid ranges for the neighborhood: if a bound + * of the neighborhood is outside the array, then limits is the same as + * boundaries. On the contrary, if a bound is strictly inside the + * array, then limits correspond to the array range. For example, for + * an array [1, 2, 3], if bounds are [-1, 3], limits will be [-1, 3], + * but if bounds are [1, 2], then limits will be [0, 2]. + * + * This is used by neighborhood iterators stacked on top of this one */ + ret->limits[i][0] = ret->bounds[i][0] < 0 ? ret->bounds[i][0] : 0; + ret->limits[i][1] = ret->bounds[i][1] >= ret->dimensions[i] - 1 ? + ret->bounds[i][1] : + ret->dimensions[i] - 1; + ret->limits_sizes[i] = (ret->limits[i][1] - ret->limits[i][0]) + 1; + } + + switch (mode) { + case NPY_NEIGHBORHOOD_ITER_ZERO_PADDING: + ret->constant = PyArray_Zero(x->ao); + ret->mode = mode; + ret->translate = &get_ptr_constant; + break; + case NPY_NEIGHBORHOOD_ITER_ONE_PADDING: + ret->constant = PyArray_One(x->ao); + ret->mode = mode; + ret->translate = &get_ptr_constant; + break; + case NPY_NEIGHBORHOOD_ITER_CONSTANT_PADDING: + /* New reference in returned value of _set_constant if array + * object */ + assert(PyArray_EquivArrTypes(x->ao, fill) == NPY_TRUE); + ret->constant = _set_constant(ret, fill); + if (ret->constant == NULL) { + goto clean_x; + } + ret->mode = mode; + ret->translate = &get_ptr_constant; + break; + case NPY_NEIGHBORHOOD_ITER_MIRROR_PADDING: + ret->mode = mode; + ret->constant = NULL; + ret->translate = &get_ptr_mirror; + break; + case NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING: + ret->mode = mode; + ret->constant = NULL; + ret->translate = &get_ptr_circular; + break; + default: + PyErr_SetString(PyExc_ValueError, "Unsupported padding mode"); + goto clean_x; + } + + /* + * XXX: we force x iterator to be non contiguous because we need + * coordinates... Modifying the iterator here is not great + */ + x->contiguous = 0; + + PyArrayNeighborhoodIter_Reset(ret); + + return (PyObject*)ret; + +clean_x: + Py_DECREF(ret->_internal_iter); + array_iter_base_dealloc((PyArrayIterObject*)ret); + _pya_free((PyArrayObject*)ret); + return NULL; +} + +static void neighiter_dealloc(PyArrayNeighborhoodIterObject* iter) +{ + if (iter->mode == NPY_NEIGHBORHOOD_ITER_CONSTANT_PADDING) { + if (PyArray_ISOBJECT(iter->_internal_iter->ao)) { + Py_DECREF(*(PyObject**)iter->constant); + } + } + if (iter->constant != NULL) { + PyDataMem_FREE(iter->constant); + } + Py_DECREF(iter->_internal_iter); + + array_iter_base_dealloc((PyArrayIterObject*)iter); + _pya_free((PyArrayObject*)iter); +} + +NPY_NO_EXPORT PyTypeObject PyArrayNeighborhoodIter_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else + PyObject_HEAD_INIT(NULL) + 0, /* ob_size */ +#endif + "numpy.neigh_internal_iter", /* tp_name*/ + sizeof(PyArrayNeighborhoodIterObject), /* tp_basicsize*/ + 0, /* tp_itemsize*/ + (destructor)neighiter_dealloc, /* tp_dealloc*/ + 0, /* tp_print*/ + 0, /* tp_getattr*/ + 0, /* tp_setattr*/ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ #endif + 0, /* tp_repr*/ + 0, /* tp_as_number*/ + 0, /* tp_as_sequence*/ + 0, /* tp_as_mapping*/ + 0, /* tp_hash */ + 0, /* tp_call*/ + 0, /* tp_str*/ + 0, /* tp_getattro*/ + 0, /* tp_setattro*/ + 0, /* tp_as_buffer*/ + Py_TPFLAGS_DEFAULT, /* tp_flags*/ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + (iternextfunc)0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + (initproc)0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; diff --git a/numpy/core/src/multiarray/mapping.c b/numpy/core/src/multiarray/mapping.c index a1c787a97..389adf02f 100644 --- a/numpy/core/src/multiarray/mapping.c +++ b/numpy/core/src/multiarray/mapping.c @@ -22,9 +22,9 @@ static PyObject * array_subscript_simple(PyArrayObject *self, PyObject *op); -/************************************************************************* - **************** Implement Mapping Protocol *************************** - *************************************************************************/ +/****************************************************************************** + *** IMPLEMENT MAPPING PROTOCOL *** + *****************************************************************************/ NPY_NO_EXPORT Py_ssize_t array_length(PyArrayObject *self) @@ -1612,63 +1612,62 @@ arraymapiter_dealloc(PyArrayMapIterObject *mit) * slice syntax. */ NPY_NO_EXPORT PyTypeObject PyArrayMapIter_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size */ - "numpy.mapiter", /* tp_name */ - sizeof(PyArrayIterObject), /* tp_basicsize */ - 0, /* tp_itemsize */ + 0, /* ob_size */ +#endif + "numpy.mapiter", /* tp_name */ + sizeof(PyArrayIterObject), /* tp_basicsize */ + 0, /* tp_itemsize */ /* methods */ - (destructor)arraymapiter_dealloc, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_compare */ - 0, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - 0, /* tp_as_mapping */ - 0, /* tp_hash */ - 0, /* tp_call */ - 0, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - Py_TPFLAGS_DEFAULT, /* tp_flags */ - 0, /* tp_doc */ - (traverseproc)0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - (iternextfunc)0, /* tp_iternext */ - 0, /* tp_methods */ - 0, /* tp_members */ - 0, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - (initproc)0, /* tp_init */ - 0, /* tp_alloc */ - 0, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ - -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + (destructor)arraymapiter_dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ #endif + 0, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + Py_TPFLAGS_DEFAULT, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; /** END of Subscript Iterator **/ diff --git a/numpy/core/src/multiarray/methods.c b/numpy/core/src/multiarray/methods.c index 8e9bf24e4..de99ca137 100644 --- a/numpy/core/src/multiarray/methods.c +++ b/numpy/core/src/multiarray/methods.c @@ -767,6 +767,49 @@ array_wraparray(PyArrayObject *self, PyObject *args) return NULL; } arr = PyTuple_GET_ITEM(args, 0); + if (arr == NULL) { + return NULL; + } + if (!PyArray_Check(arr)) { + PyErr_SetString(PyExc_TypeError, + "can only be called with ndarray object"); + return NULL; + } + + if (self->ob_type != arr->ob_type){ + Py_INCREF(PyArray_DESCR(arr)); + ret = PyArray_NewFromDescr(self->ob_type, + PyArray_DESCR(arr), + PyArray_NDIM(arr), + PyArray_DIMS(arr), + PyArray_STRIDES(arr), PyArray_DATA(arr), + PyArray_FLAGS(arr), (PyObject *)self); + if (ret == NULL) { + return NULL; + } + Py_INCREF(arr); + PyArray_BASE(ret) = arr; + return ret; + } else { + /*The type was set in __array_prepare__*/ + Py_INCREF(arr); + return arr; + } +} + + +static PyObject * +array_preparearray(PyArrayObject *self, PyObject *args) +{ + PyObject *arr; + PyObject *ret; + + if (PyTuple_Size(args) < 1) { + PyErr_SetString(PyExc_TypeError, + "only accepts 1 argument"); + return NULL; + } + arr = PyTuple_GET_ITEM(args, 0); if (!PyArray_Check(arr)) { PyErr_SetString(PyExc_TypeError, "can only be called with ndarray object"); @@ -2031,6 +2074,8 @@ NPY_NO_EXPORT PyMethodDef array_methods[] = { /* for subtypes */ {"__array__", (PyCFunction)array_getarray, METH_VARARGS, NULL}, + {"__array_prepare__", (PyCFunction)array_preparearray, + METH_VARARGS, NULL}, {"__array_wrap__", (PyCFunction)array_wraparray, METH_VARARGS, NULL}, diff --git a/numpy/core/src/multiarray/multiarray_tests.c.src b/numpy/core/src/multiarray/multiarray_tests.c.src new file mode 100644 index 000000000..c09ccbc9d --- /dev/null +++ b/numpy/core/src/multiarray/multiarray_tests.c.src @@ -0,0 +1,390 @@ +#include <Python.h> +#include "numpy/ndarrayobject.h" + +/* + * TODO: + * - Handle mode + */ + +/**begin repeat + * #type = double, int# + * #typenum = NPY_DOUBLE, NPY_INT# + */ +static int copy_@type@(PyArrayIterObject *itx, PyArrayNeighborhoodIterObject *niterx, + npy_intp *bounds, + PyObject **out) +{ + npy_intp i, j; + @type@ *ptr; + npy_intp odims[NPY_MAXDIMS]; + PyArrayObject *aout; + + /* + * For each point in itx, copy the current neighborhood into an array which + * is appended at the output list + */ + for (i = 0; i < itx->size; ++i) { + PyArrayNeighborhoodIter_Reset(niterx); + + for (j = 0; j < itx->ao->nd; ++j) { + odims[j] = bounds[2 * j + 1] - bounds[2 * j] + 1; + } + aout = (PyArrayObject*)PyArray_SimpleNew(itx->ao->nd, odims, @typenum@); + if (aout == NULL) { + return -1; + } + + ptr = (@type@*)aout->data; + + for (j = 0; j < niterx->size; ++j) { + *ptr = *((@type@*)niterx->dataptr); + PyArrayNeighborhoodIter_Next(niterx); + ptr += 1; + } + + Py_INCREF(aout); + PyList_Append(*out, (PyObject*)aout); + Py_DECREF(aout); + PyArray_ITER_NEXT(itx); + } + + return 0; +} +/**end repeat**/ + +static int copy_object(PyArrayIterObject *itx, PyArrayNeighborhoodIterObject *niterx, + npy_intp *bounds, + PyObject **out) +{ + npy_intp i, j; + npy_intp odims[NPY_MAXDIMS]; + PyArrayObject *aout; + PyArray_CopySwapFunc *copyswap = itx->ao->descr->f->copyswap; + npy_int itemsize = PyArray_ITEMSIZE(itx->ao); + + /* + * For each point in itx, copy the current neighborhood into an array which + * is appended at the output list + */ + for (i = 0; i < itx->size; ++i) { + PyArrayNeighborhoodIter_Reset(niterx); + + for (j = 0; j < itx->ao->nd; ++j) { + odims[j] = bounds[2 * j + 1] - bounds[2 * j] + 1; + } + aout = (PyArrayObject*)PyArray_SimpleNew(itx->ao->nd, odims, NPY_OBJECT); + if (aout == NULL) { + return -1; + } + + for (j = 0; j < niterx->size; ++j) { + copyswap(aout->data + j * itemsize, niterx->dataptr, 0, NULL); + PyArrayNeighborhoodIter_Next(niterx); + } + + Py_INCREF(aout); + PyList_Append(*out, (PyObject*)aout); + Py_DECREF(aout); + PyArray_ITER_NEXT(itx); + } + + return 0; +} + +static PyObject* +test_neighborhood_iterator(PyObject* NPY_UNUSED(self), PyObject* args) +{ + PyObject *x, *fill, *out, *b; + PyArrayObject *ax, *afill; + PyArrayIterObject *itx; + int i, typenum, mode, st; + npy_intp bounds[NPY_MAXDIMS*2]; + PyArrayNeighborhoodIterObject *niterx; + + if (!PyArg_ParseTuple(args, "OOOi", &x, &b, &fill, &mode)) { + return NULL; + } + + if (!PySequence_Check(b)) { + return NULL; + } + + typenum = PyArray_ObjectType(x, 0); + typenum = PyArray_ObjectType(fill, typenum); + + ax = (PyArrayObject*)PyArray_FromObject(x, typenum, 1, 10); + if (ax == NULL) { + return NULL; + } + if (PySequence_Size(b) != 2 * ax->nd) { + PyErr_SetString(PyExc_ValueError, + "bounds sequence size not compatible with x input"); + goto clean_ax; + } + + out = PyList_New(0); + if (out == NULL) { + goto clean_ax; + } + + itx = (PyArrayIterObject*)PyArray_IterNew(x); + if (itx == NULL) { + goto clean_out; + } + + /* Compute boundaries for the neighborhood iterator */ + for (i = 0; i < 2 * ax->nd; ++i) { + PyObject* bound; + bound = PySequence_GetItem(b, i); + if (bounds == NULL) { + goto clean_itx; + } + if (!PyInt_Check(bound)) { + PyErr_SetString(PyExc_ValueError, "bound not long"); + Py_DECREF(bound); + goto clean_itx; + } + bounds[i] = PyInt_AsLong(bound); + Py_DECREF(bound); + } + + /* Create the neighborhood iterator */ + afill = NULL; + if (mode == NPY_NEIGHBORHOOD_ITER_CONSTANT_PADDING) { + afill = (PyArrayObject *)PyArray_FromObject(fill, typenum, 0, 0); + if (afill == NULL) { + goto clean_itx; + } + } + + niterx = (PyArrayNeighborhoodIterObject*)PyArray_NeighborhoodIterNew( + (PyArrayIterObject*)itx, bounds, mode, afill); + if (niterx == NULL) { + goto clean_afill; + } + + switch (typenum) { + case NPY_OBJECT: + st = copy_object(itx, niterx, bounds, &out); + break; + case NPY_INT: + st = copy_int(itx, niterx, bounds, &out); + break; + case NPY_DOUBLE: + st = copy_double(itx, niterx, bounds, &out); + break; + default: + PyErr_SetString(PyExc_ValueError, "Type not supported"); + goto clean_niterx; + } + + if (st) { + goto clean_niterx; + } + + Py_DECREF(niterx); + Py_XDECREF(afill); + Py_DECREF(itx); + + Py_DECREF(ax); + + return out; + +clean_niterx: + Py_DECREF(niterx); +clean_afill: + Py_XDECREF(afill); +clean_itx: + Py_DECREF(itx); +clean_out: + Py_DECREF(out); +clean_ax: + Py_DECREF(ax); + return NULL; +} + +static int +copy_double_double(PyArrayNeighborhoodIterObject *itx, + PyArrayNeighborhoodIterObject *niterx, + npy_intp *bounds, + PyObject **out) +{ + npy_intp i, j; + double *ptr; + npy_intp odims[NPY_MAXDIMS]; + PyArrayObject *aout; + + /* + * For each point in itx, copy the current neighborhood into an array which + * is appended at the output list + */ + PyArrayNeighborhoodIter_Reset(itx); + for (i = 0; i < itx->size; ++i) { + for (j = 0; j < itx->ao->nd; ++j) { + odims[j] = bounds[2 * j + 1] - bounds[2 * j] + 1; + } + aout = (PyArrayObject*)PyArray_SimpleNew(itx->ao->nd, odims, NPY_DOUBLE); + if (aout == NULL) { + return -1; + } + + ptr = (double*)aout->data; + + PyArrayNeighborhoodIter_Reset(niterx); + for (j = 0; j < niterx->size; ++j) { + *ptr = *((double*)niterx->dataptr); + ptr += 1; + PyArrayNeighborhoodIter_Next(niterx); + } + Py_INCREF(aout); + PyList_Append(*out, (PyObject*)aout); + Py_DECREF(aout); + PyArrayNeighborhoodIter_Next(itx); + } + return 0; +} + +static PyObject* +test_neighborhood_iterator_oob(PyObject* NPY_UNUSED(self), PyObject* args) +{ + PyObject *x, *out, *b1, *b2; + PyArrayObject *ax; + PyArrayIterObject *itx; + int i, typenum, mode1, mode2, st; + npy_intp bounds[NPY_MAXDIMS*2]; + PyArrayNeighborhoodIterObject *niterx1, *niterx2; + + if (!PyArg_ParseTuple(args, "OOiOi", &x, &b1, &mode1, &b2, &mode2)) { + return NULL; + } + + if (!PySequence_Check(b1) || !PySequence_Check(b2)) { + return NULL; + } + + typenum = PyArray_ObjectType(x, 0); + + ax = (PyArrayObject*)PyArray_FromObject(x, typenum, 1, 10); + if (ax == NULL) { + return NULL; + } + if (PySequence_Size(b1) != 2 * ax->nd) { + PyErr_SetString(PyExc_ValueError, + "bounds sequence 1 size not compatible with x input"); + goto clean_ax; + } + if (PySequence_Size(b2) != 2 * ax->nd) { + PyErr_SetString(PyExc_ValueError, + "bounds sequence 2 size not compatible with x input"); + goto clean_ax; + } + + out = PyList_New(0); + if (out == NULL) { + goto clean_ax; + } + + itx = (PyArrayIterObject*)PyArray_IterNew(x); + if (itx == NULL) { + goto clean_out; + } + + /* Compute boundaries for the neighborhood iterator */ + for (i = 0; i < 2 * ax->nd; ++i) { + PyObject* bound; + bound = PySequence_GetItem(b1, i); + if (bounds == NULL) { + goto clean_itx; + } + if (!PyInt_Check(bound)) { + PyErr_SetString(PyExc_ValueError, "bound not long"); + Py_DECREF(bound); + goto clean_itx; + } + bounds[i] = PyInt_AsLong(bound); + Py_DECREF(bound); + } + + /* Create the neighborhood iterator */ + niterx1 = (PyArrayNeighborhoodIterObject*)PyArray_NeighborhoodIterNew( + (PyArrayIterObject*)itx, bounds, + mode1, NULL); + if (niterx1 == NULL) { + goto clean_out; + } + + for (i = 0; i < 2 * ax->nd; ++i) { + PyObject* bound; + bound = PySequence_GetItem(b2, i); + if (bounds == NULL) { + goto clean_itx; + } + if (!PyInt_Check(bound)) { + PyErr_SetString(PyExc_ValueError, "bound not long"); + Py_DECREF(bound); + goto clean_itx; + } + bounds[i] = PyInt_AsLong(bound); + Py_DECREF(bound); + } + + niterx2 = (PyArrayNeighborhoodIterObject*)PyArray_NeighborhoodIterNew( + (PyArrayIterObject*)niterx1, bounds, + mode2, NULL); + if (niterx1 == NULL) { + goto clean_niterx1; + } + + switch (typenum) { + case NPY_DOUBLE: + st = copy_double_double(niterx1, niterx2, bounds, &out); + break; + default: + PyErr_SetString(PyExc_ValueError, "Type not supported"); + goto clean_niterx2; + } + + if (st) { + goto clean_niterx2; + } + + Py_DECREF(niterx2); + Py_DECREF(niterx1); + Py_DECREF(itx); + Py_DECREF(ax); + return out; + +clean_niterx2: + Py_DECREF(niterx2); +clean_niterx1: + Py_DECREF(niterx1); +clean_itx: + Py_DECREF(itx); +clean_out: + Py_DECREF(out); +clean_ax: + Py_DECREF(ax); + return NULL; +} + +static PyMethodDef Multiarray_TestsMethods[] = { + {"test_neighborhood_iterator", test_neighborhood_iterator, METH_VARARGS, NULL}, + {"test_neighborhood_iterator_oob", test_neighborhood_iterator_oob, METH_VARARGS, NULL}, + {NULL, NULL, 0, NULL} /* Sentinel */ +}; + +PyMODINIT_FUNC +initmultiarray_tests(void) +{ + PyObject *m; + + m = Py_InitModule("multiarray_tests", Multiarray_TestsMethods); + if (m == NULL) { + return; + } + import_array(); + if (PyErr_Occurred()) { + PyErr_SetString(PyExc_RuntimeError, + "cannot load umath_tests module."); + } +} diff --git a/numpy/core/src/multiarray/multiarraymodule.c b/numpy/core/src/multiarray/multiarraymodule.c index 6f9c34313..a793e0364 100644 --- a/numpy/core/src/multiarray/multiarraymodule.c +++ b/numpy/core/src/multiarray/multiarraymodule.c @@ -27,7 +27,7 @@ #include "config.h" -#include "global.c" +NPY_NO_EXPORT int NPY_NUMUSERTYPES = 0; #define PyAO PyArrayObject @@ -906,7 +906,7 @@ PyArray_CopyAndTranspose(PyObject *op) } /* - * Implementation which is common between PyArray_Correlate and PyArray_Acorrelate + * Implementation which is common between PyArray_Correlate and PyArray_Correlate2 * * inverted is set to 1 if computed correlate(ap2, ap1), 0 otherwise */ @@ -1065,13 +1065,13 @@ _pyarray_revert(PyArrayObject *ret) } /*NUMPY_API - * acorrelate(a1,a2,mode) + * correlate(a1,a2,mode) * - * This function computes the usual correlation (acorrelate(a1, a2) != - * accorrelate(a2, a1), and conjugate the second argument for complex inputs + * This function computes the usual correlation (correlate(a1, a2) != + * correlate(a2, a1), and conjugate the second argument for complex inputs */ NPY_NO_EXPORT PyObject * -PyArray_Acorrelate(PyObject *op1, PyObject *op2, int mode) +PyArray_Correlate2(PyObject *op1, PyObject *op2, int mode) { PyArrayObject *ap1, *ap2, *ret = NULL; int typenum; @@ -1796,7 +1796,7 @@ static PyObject *array_correlate(PyObject *NPY_UNUSED(dummy), PyObject *args, Py } static PyObject* -array_acorrelate(PyObject *NPY_UNUSED(dummy), PyObject *args, PyObject *kwds) +array_correlate2(PyObject *NPY_UNUSED(dummy), PyObject *args, PyObject *kwds) { PyObject *shape, *a0; int mode = 0; @@ -1806,7 +1806,7 @@ array_acorrelate(PyObject *NPY_UNUSED(dummy), PyObject *args, PyObject *kwds) &a0, &shape, &mode)) { return NULL; } - return PyArray_Acorrelate(a0, shape, mode); + return PyArray_Correlate2(a0, shape, mode); } static PyObject * @@ -2448,8 +2448,8 @@ static struct PyMethodDef array_module_methods[] = { {"correlate", (PyCFunction)array_correlate, METH_VARARGS | METH_KEYWORDS, NULL}, - {"acorrelate", - (PyCFunction)array_acorrelate, + {"correlate2", + (PyCFunction)array_correlate2, METH_VARARGS | METH_KEYWORDS, NULL}, {"frombuffer", (PyCFunction)array_frombuffer, @@ -2710,6 +2710,11 @@ PyMODINIT_FUNC initmultiarray(void) { if (PyType_Ready(&PyArrayMultiIter_Type) < 0) { return; } + PyArrayNeighborhoodIter_Type.tp_new = PyType_GenericNew; + if (PyType_Ready(&PyArrayNeighborhoodIter_Type) < 0) { + return; + } + PyArrayDescr_Type.tp_hash = PyArray_DescrHash; if (PyType_Ready(&PyArrayDescr_Type) < 0) { return; diff --git a/numpy/core/src/multiarray/scalartypes.c.src b/numpy/core/src/multiarray/scalartypes.c.src index 28ba7f47a..e50106866 100644 --- a/numpy/core/src/multiarray/scalartypes.c.src +++ b/numpy/core/src/multiarray/scalartypes.c.src @@ -37,63 +37,62 @@ NPY_NO_EXPORT PyBoolScalarObject _PyArrayScalar_BoolValues[] = { * Floating, ComplexFloating, Flexible, Character, TimeInteger# */ NPY_NO_EXPORT PyTypeObject Py@NAME@ArrType_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size*/ - "numpy.@name@", /* tp_name*/ - sizeof(PyObject), /* tp_basicsize*/ - 0, /* tp_itemsize */ + 0, /* ob_size */ +#endif + "numpy.@name@", /* tp_name*/ + sizeof(PyObject), /* tp_basicsize*/ + 0, /* tp_itemsize */ /* methods */ - 0, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_compare */ - 0, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - 0, /* tp_as_mapping */ - 0, /* tp_hash */ - 0, /* tp_call */ - 0, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - 0, /* tp_flags */ - 0, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - 0, /* tp_methods */ - 0, /* tp_members */ - 0, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - 0, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ - -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + 0, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ #endif + 0, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + 0, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; /**end repeat**/ @@ -1805,63 +1804,62 @@ static PyBufferProcs gentype_as_buffer = { #define LEAFFLAGS Py_TPFLAGS_DEFAULT | Py_TPFLAGS_CHECKTYPES NPY_NO_EXPORT PyTypeObject PyGenericArrType_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size*/ - "numpy.generic", /* tp_name*/ - sizeof(PyObject), /* tp_basicsize*/ - 0, /* tp_itemsize */ + 0, /* ob_size */ +#endif + "numpy.generic", /* tp_name*/ + sizeof(PyObject), /* tp_basicsize*/ + 0, /* tp_itemsize */ /* methods */ - 0, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_compare */ - 0, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - 0, /* tp_as_mapping */ - 0, /* tp_hash */ - 0, /* tp_call */ - 0, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - 0, /* tp_flags */ - 0, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - 0, /* tp_methods */ - 0, /* tp_members */ - 0, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - 0, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ - -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + 0, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ #endif + 0, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + 0, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; static void @@ -2610,61 +2608,61 @@ object_arrtype_call(PyObjectScalarObject *obj, PyObject *args, PyObject *kwds) } NPY_NO_EXPORT PyTypeObject PyObjectArrType_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size*/ - "numpy.object_", /* tp_name*/ - sizeof(PyObjectScalarObject), /* tp_basicsize*/ - 0, /* tp_itemsize */ - (destructor)object_arrtype_dealloc, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_compare */ - 0, /* tp_repr */ - 0, /* tp_as_number */ - &object_arrtype_as_sequence, /* tp_as_sequence */ - &object_arrtype_as_mapping, /* tp_as_mapping */ - 0, /* tp_hash */ - (ternaryfunc)object_arrtype_call, /* tp_call */ - 0, /* tp_str */ - (getattrofunc)object_arrtype_getattro, /* tp_getattro */ - (setattrofunc)object_arrtype_setattro, /* tp_setattro */ - &object_arrtype_as_buffer, /* tp_as_buffer */ - 0, /* tp_flags */ - 0, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - 0, /* tp_methods */ - 0, /* tp_members */ - 0, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - 0, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + 0, /* ob_size */ #endif + "numpy.object_", /* tp_name*/ + sizeof(PyObjectScalarObject), /* tp_basicsize*/ + 0, /* tp_itemsize */ + (destructor)object_arrtype_dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ +#endif + 0, /* tp_repr */ + 0, /* tp_as_number */ + &object_arrtype_as_sequence, /* tp_as_sequence */ + &object_arrtype_as_mapping, /* tp_as_mapping */ + 0, /* tp_hash */ + (ternaryfunc)object_arrtype_call, /* tp_call */ + 0, /* tp_str */ + (getattrofunc)object_arrtype_getattro, /* tp_getattro */ + (setattrofunc)object_arrtype_setattro, /* tp_setattro */ + &object_arrtype_as_buffer, /* tp_as_buffer */ + 0, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; static PyObject * @@ -2714,61 +2712,61 @@ gen_arrtype_subscript(PyObject *self, PyObject *key) * #ex = _,_,_,# */ NPY_NO_EXPORT PyTypeObject Py@NAME@ArrType_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size*/ - "numpy.@name@@ex@", /* tp_name*/ - sizeof(Py@NAME@ScalarObject), /* tp_basicsize*/ - 0, /* tp_itemsize */ - 0, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_compare */ - 0, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - 0, /* tp_as_mapping */ - 0, /* tp_hash */ - 0, /* tp_call */ - 0, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - 0, /* tp_flags */ - 0, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - 0, /* tp_methods */ - 0, /* tp_members */ - 0, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - 0, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + 0, /* ob_size */ #endif + "numpy.@name@@ex@", /* tp_name*/ + sizeof(Py@NAME@ScalarObject), /* tp_basicsize*/ + 0, /* tp_itemsize */ + 0, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ +#endif + 0, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + 0, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; /**end repeat**/ @@ -2797,61 +2795,61 @@ NPY_NO_EXPORT PyTypeObject Py@NAME@ArrType_Type = { #define _THIS_SIZE "256" #endif NPY_NO_EXPORT PyTypeObject Py@NAME@ArrType_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else PyObject_HEAD_INIT(NULL) - 0, /* ob_size*/ - "numpy.@name@" _THIS_SIZE, /* tp_name*/ - sizeof(Py@NAME@ScalarObject), /* tp_basicsize*/ - 0, /* tp_itemsize */ - 0, /* tp_dealloc */ - 0, /* tp_print */ - 0, /* tp_getattr */ - 0, /* tp_setattr */ - 0, /* tp_compare */ - 0, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - 0, /* tp_as_mapping */ - 0, /* tp_hash */ - 0, /* tp_call */ - 0, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - 0, /* tp_flags */ - 0, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - 0, /* tp_methods */ - 0, /* tp_members */ - 0, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - 0, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + 0, /* ob_size */ +#endif + "numpy.@name@" _THIS_SIZE, /* tp_name*/ + sizeof(Py@NAME@ScalarObject), /* tp_basicsize*/ + 0, /* tp_itemsize */ + 0, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ #endif + 0, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + 0, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; #undef _THIS_SIZE @@ -2895,62 +2893,62 @@ static PyMappingMethods gentype_as_mapping = { #define _THIS_DOC "Composed of two " _THIS_SIZE2 " bit floats" - NPY_NO_EXPORT PyTypeObject Py@NAME@ArrType_Type = { - PyObject_HEAD_INIT(NULL) - 0, /* ob_size*/ - "numpy.@name@" _THIS_SIZE1, /* tp_name*/ - sizeof(Py@NAME@ScalarObject), /* tp_basicsize*/ - 0, /* tp_itemsize*/ - 0, /* tp_dealloc*/ - 0, /* tp_print*/ - 0, /* tp_getattr*/ - 0, /* tp_setattr*/ - 0, /* tp_compare*/ - 0, /* tp_repr*/ - 0, /* tp_as_number*/ - 0, /* tp_as_sequence*/ - 0, /* tp_as_mapping*/ - 0, /* tp_hash */ - 0, /* tp_call*/ - 0, /* tp_str*/ - 0, /* tp_getattro*/ - 0, /* tp_setattro*/ - 0, /* tp_as_buffer*/ - Py_TPFLAGS_DEFAULT, /* tp_flags*/ - _THIS_DOC, /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - 0, /* tp_methods */ - 0, /* tp_members */ - 0, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - 0, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ +NPY_NO_EXPORT PyTypeObject Py@NAME@ArrType_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(0, 0) +#else + PyObject_HEAD_INIT(0) + 0, /* ob_size */ +#endif + "numpy.@name@" _THIS_SIZE1, /* tp_name*/ + sizeof(Py@NAME@ScalarObject), /* tp_basicsize*/ + 0, /* tp_itemsize*/ + 0, /* tp_dealloc*/ + 0, /* tp_print*/ + 0, /* tp_getattr*/ + 0, /* tp_setattr*/ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ #endif + 0, /* tp_repr*/ + 0, /* tp_as_number*/ + 0, /* tp_as_sequence*/ + 0, /* tp_as_mapping*/ + 0, /* tp_hash */ + 0, /* tp_call*/ + 0, /* tp_str*/ + 0, /* tp_getattro*/ + 0, /* tp_setattro*/ + 0, /* tp_as_buffer*/ + Py_TPFLAGS_DEFAULT, /* tp_flags*/ + _THIS_DOC, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; #undef _THIS_SIZE1 #undef _THIS_SIZE2 diff --git a/numpy/core/src/_signbit.c b/numpy/core/src/npymath/_signbit.c index a2ad38162..a2ad38162 100644 --- a/numpy/core/src/_signbit.c +++ b/numpy/core/src/npymath/_signbit.c diff --git a/numpy/core/src/npy_math.c.src b/numpy/core/src/npymath/npy_math.c.src index 21fc7d427..3fde802a2 100644 --- a/numpy/core/src/npy_math.c.src +++ b/numpy/core/src/npymath/npy_math.c.src @@ -40,6 +40,18 @@ * #ifdef SYMBOL_DEFINED_WEIRD_PLATFORM * double exp(double); * #endif + * + * Some of the code is taken from msun library in FreeBSD, with the following + * notice: + * + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== */ #include <Python.h> @@ -48,6 +60,8 @@ #include "config.h" #include "numpy/npy_math.h" +#include "npy_math_private.h" + /* ***************************************************************************** ** BASIC MATH FUNCTIONS ** @@ -56,61 +70,158 @@ /* Original code by Konrad Hinsen. */ #ifndef HAVE_EXPM1 -static double expm1(double x) +double npy_expm1(double x) { - double u = exp(x); + double u = npy_exp(x); if (u == 1.0) { return x; } else if (u-1.0 == -1.0) { return -1; } else { - return (u-1.0) * x/log(u); + return (u-1.0) * x/npy_log(u); } } #endif #ifndef HAVE_LOG1P -static double log1p(double x) +double npy_log1p(double x) { double u = 1. + x; if (u == 1.0) { return x; } else { - return log(u) * x / (u - 1); + return npy_log(u) * x / (u - 1); } } #endif +/* Taken from FreeBSD mlib, adapted for numpy + * + * XXX: we could be a bit faster by reusing high/low words for inf/nan + * classification instead of calling npy_isinf/npy_isnan: we should have some + * macros for this, though, instead of doing it manually + */ +#ifndef HAVE_ATAN2 +/* XXX: we should have this in npy_math.h */ +#define NPY_DBL_EPSILON 1.2246467991473531772E-16 +double npy_atan2(double y, double x) +{ + npy_int32 k, m, iy, ix, hx, hy; + npy_uint32 lx,ly; + double z; + + EXTRACT_WORDS(hx, lx, x); + ix = hx & 0x7fffffff; + EXTRACT_WORDS(hy, ly, y); + iy = hy & 0x7fffffff; + + /* if x or y is nan, return nan */ + if (npy_isnan(x * y)) { + return x + y; + } + + if (x == 1.0) { + return npy_atan(y); + } + + m = 2 * npy_signbit(x) + npy_signbit(y); + if (y == 0.0) { + switch(m) { + case 0: + case 1: return y; /* atan(+-0,+anything)=+-0 */ + case 2: return NPY_PI;/* atan(+0,-anything) = pi */ + case 3: return -NPY_PI;/* atan(-0,-anything) =-pi */ + } + } + + if (x == 0.0) { + return y > 0 ? NPY_PI_2 : -NPY_PI_2; + } + + if (npy_isinf(x)) { + if (npy_isinf(y)) { + switch(m) { + case 0: return NPY_PI_4;/* atan(+INF,+INF) */ + case 1: return -NPY_PI_4;/* atan(-INF,+INF) */ + case 2: return 3.0*NPY_PI_4;/*atan(+INF,-INF)*/ + case 3: return -3.0*NPY_PI_4;/*atan(-INF,-INF)*/ + } + } else { + switch(m) { + case 0: return NPY_PZERO; /* atan(+...,+INF) */ + case 1: return NPY_NZERO; /* atan(-...,+INF) */ + case 2: return NPY_PI; /* atan(+...,-INF) */ + case 3: return -NPY_PI; /* atan(-...,-INF) */ + } + } + } + + if (npy_isinf(y)) { + return y > 0 ? NPY_PI_2 : -NPY_PI_2; + } + + /* compute y/x */ + k = (iy - ix)>>20; + if(k > 60) { /* |y/x| > 2**60 */ + z = NPY_PI_2 + 0.5 * NPY_DBL_EPSILON; + m &= 1; + } else if(hx < 0 && k < -60) { + z = 0.0; /* 0 > |y|/x > -2**-60 */ + } else { + z = npy_atan(npy_fabs(y/x)); /* safe to do y/x */ + } + + switch (m) { + case 0: return z ; /* atan(+,+) */ + case 1: return -z ; /* atan(-,+) */ + case 2: return NPY_PI - (z - NPY_DBL_EPSILON);/* atan(+,-) */ + default: /* case 3 */ + return (z - NPY_DBL_EPSILON) - NPY_PI;/* atan(-,-) */ + } +} + +#endif + #ifndef HAVE_HYPOT -static double hypot(double x, double y) +double npy_hypot(double x, double y) { double yx; - x = fabs(x); - y = fabs(y); + /* Handle the case where x or y is a NaN */ + if (npy_isnan(x * y)) { + if (npy_isinf(x) || npy_isinf(y)) { + return NPY_INFINITY; + } else { + return NPY_NAN; + } + } + + x = npy_fabs(x); + y = npy_fabs(y); if (x < y) { double temp = x; x = y; y = temp; } - if (x == 0.) + if (x == 0.) { return 0.; + } else { yx = y/x; - return x*sqrt(1.+yx*yx); + return x*npy_sqrt(1.+yx*yx); } } #endif #ifndef HAVE_ACOSH -static double acosh(double x) +double npy_acosh(double x) { - return 2*log(sqrt((x+1.0)/2)+sqrt((x-1.0)/2)); + return 2*npy_log(npy_sqrt((x+1.0)/2)+npy_sqrt((x-1.0)/2)); } #endif #ifndef HAVE_ASINH -static double asinh(double xx) +double npy_asinh(double xx) { double x, d; int sign; @@ -125,37 +236,37 @@ static double asinh(double xx) if (x > 1e8) { d = x; } else { - d = sqrt(x*x + 1); + d = npy_sqrt(x*x + 1); } - return sign*log1p(x*(1.0 + x/(d+1))); + return sign*npy_log1p(x*(1.0 + x/(d+1))); } #endif #ifndef HAVE_ATANH -static double atanh(double x) +double npy_atanh(double x) { if (x > 0) { - return -0.5*log1p(-2.0*x/(1.0 + x)); + return -0.5*npy_log1p(-2.0*x/(1.0 + x)); } else { - return 0.5*log1p(2.0*x/(1.0 - x)); + return 0.5*npy_log1p(2.0*x/(1.0 - x)); } } #endif #ifndef HAVE_RINT -static double rint(double x) +double npy_rint(double x) { double y, r; - y = floor(x); + y = npy_floor(x); r = x - y; if (r > 0.5) goto rndup; /* Round to nearest even */ if (r==0.5) { - r = y - 2.0*floor(0.5*y); + r = y - 2.0*npy_floor(0.5*y); if (r==1.0) { rndup: y+=1.0; @@ -166,30 +277,41 @@ static double rint(double x) #endif #ifndef HAVE_TRUNC -static double trunc(double x) +double npy_trunc(double x) { - return x < 0 ? ceil(x) : floor(x); + return x < 0 ? npy_ceil(x) : npy_floor(x); } #endif #ifndef HAVE_EXP2 #define LOG2 0.69314718055994530943 -static double exp2(double x) +double npy_exp2(double x) { - return exp(LOG2*x); + return npy_exp(LOG2*x); } #undef LOG2 #endif #ifndef HAVE_LOG2 #define INVLOG2 1.4426950408889634074 -static double log2(double x) +double npy_log2(double x) { - return INVLOG2*log(x); + return INVLOG2*npy_log(x); } #undef INVLOG2 #endif +#ifndef HAVE_COPYSIGN +double npy_copysign(double x, double y) +{ + npy_uint32 hx,hy; + GET_HIGH_WORD(hx,x); + GET_HIGH_WORD(hy,y); + SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000)); + return x; +} +#endif + /* ***************************************************************************** ** IEEE 754 FPU HANDLING ** @@ -247,25 +369,25 @@ int _npy_signbit_ld (long double x) #undef @kind@@c@ #endif #ifndef HAVE_@KIND@@C@ -static @type@ @kind@@c@(@type@ x) +@type@ npy_@kind@@c@(@type@ x) { - return (@type@) @kind@((double)x); + return (@type@) npy_@kind@((double)x); } #endif /**end repeat1**/ /**begin repeat1 - * #kind = atan2,hypot,pow,fmod# - * #KIND = ATAN2,HYPOT,POW,FMOD# + * #kind = atan2,hypot,pow,fmod,copysign# + * #KIND = ATAN2,HYPOT,POW,FMOD,COPYSIGN# */ #ifdef @kind@@c@ #undef @kind@@c@ #endif #ifndef HAVE_@KIND@@C@ -static @type@ @kind@@c@(@type@ x, @type@ y) +@type@ npy_@kind@@c@(@type@ x, @type@ y) { - return (@type@) @kind@((double)x, (double) y); + return (@type@) npy_@kind@((double)x, (double) y); } #endif /**end repeat1**/ @@ -274,10 +396,10 @@ static @type@ @kind@@c@(@type@ x, @type@ y) #undef modf@c@ #endif #ifndef HAVE_MODF@C@ -static @type@ modf@c@(@type@ x, @type@ *iptr) +@type@ npy_modf@c@(@type@ x, @type@ *iptr) { double niptr; - double y = modf((double)x, &niptr); + double y = npy_modf((double)x, &niptr); *iptr = (@type@) niptr; return (@type@) y; } @@ -285,27 +407,6 @@ static @type@ modf@c@(@type@ x, @type@ *iptr) /**end repeat**/ -/* - * Useful constants in three precisions: - * XXX: those should really be in the header - */ - -/**begin repeat - * #c = f, ,l# - * #C = F, ,L# - */ -#define NPY_E@c@ 2.7182818284590452353602874713526625@C@ /* e */ -#define NPY_LOG2E@c@ 1.4426950408889634073599246810018921@C@ /* log_2 e */ -#define NPY_LOG10E@c@ 0.4342944819032518276511289189166051@C@ /* log_10 e */ -#define NPY_LOGE2@c@ 0.6931471805599453094172321214581766@C@ /* log_e 2 */ -#define NPY_LOGE10@c@ 2.3025850929940456840179914546843642@C@ /* log_e 10 */ -#define NPY_PI@c@ 3.1415926535897932384626433832795029@C@ /* pi */ -#define NPY_PI_2@c@ 1.5707963267948966192313216916397514@C@ /* pi/2 */ -#define NPY_PI_4@c@ 0.7853981633974483096156608458198757@C@ /* pi/4 */ -#define NPY_1_PI@c@ 0.3183098861837906715377675267450287@C@ /* 1/pi */ -#define NPY_2_PI@c@ 0.6366197723675813430755350534900574@C@ /* 2/pi */ -/**end repeat**/ - /* * Non standard functions */ @@ -321,17 +422,17 @@ static @type@ modf@c@(@type@ x, @type@ *iptr) #define RAD2DEG (180.0@c@/NPY_PI@c@) #define DEG2RAD (NPY_PI@c@/180.0@c@) -static @type@ rad2deg@c@(@type@ x) +@type@ npy_rad2deg@c@(@type@ x) { return x*RAD2DEG; } -static @type@ deg2rad@c@(@type@ x) +@type@ npy_deg2rad@c@(@type@ x) { return x*DEG2RAD; } -static @type@ log2_1p@c@(@type@ x) +@type@ npy_log2_1p@c@(@type@ x) { @type@ u = 1 + x; if (u == 1) { @@ -341,9 +442,9 @@ static @type@ log2_1p@c@(@type@ x) } } -static @type@ exp2_1m@c@(@type@ x) +@type@ npy_exp2_1m@c@(@type@ x) { - @type@ u = exp@c@(x); + @type@ u = npy_exp@c@(x); if (u == 1.0) { return LOGE2*x; } else if (u - 1 == -1) { @@ -353,31 +454,36 @@ static @type@ exp2_1m@c@(@type@ x) } } -static @type@ logaddexp@c@(@type@ x, @type@ y) +@type@ npy_logaddexp@c@(@type@ x, @type@ y) { const @type@ tmp = x - y; if (tmp > 0) { return x + npy_log1p@c@(npy_exp@c@(-tmp)); } - else { + else if (tmp <= 0) { return y + npy_log1p@c@(npy_exp@c@(tmp)); } + else { + /* NaNs, or infinities of the same sign involved */ + return x + y; + } } -static @type@ logaddexp2@c@(@type@ x, @type@ y) +@type@ npy_logaddexp2@c@(@type@ x, @type@ y) { const @type@ tmp = x - y; if (tmp > 0) { - return x + log2_1p@c@(npy_exp2@c@(-tmp)); + return x + npy_log2_1p@c@(npy_exp2@c@(-tmp)); + } + else if (tmp <= 0) { + return y + npy_log2_1p@c@(npy_exp2@c@(tmp)); } else { - return y + log2_1p@c@(npy_exp2@c@(tmp)); + /* NaNs, or infinities of the same sign involved */ + return x + y; } } -#define degrees@c@ rad2deg@c@ -#define radians@c@ deg2rad@c@ - #undef LOGE2 #undef LOG2E #undef RAD2DEG @@ -386,38 +492,46 @@ static @type@ logaddexp2@c@(@type@ x, @type@ y) /**end repeat**/ /* - * Decorate all the functions: those are the public ones + * Decorate all the math functions which are available on the current platform */ /**begin repeat * #type = npy_longdouble,double,float# * #c = l,,f# + * #C = L,,F# */ /**begin repeat1 * #kind = sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,rint,trunc,sqrt,log10, - * log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2, - * rad2deg,deg2rad,exp2_1m# + * log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2# + * #KIND = SIN,COS,TAN,SINH,COSH,TANH,FABS,FLOOR,CEIL,RINT,TRUNC,SQRT,LOG10, + * LOG,EXP,EXPM1,ASIN,ACOS,ATAN,ASINH,ACOSH,ATANH,LOG1P,EXP2,LOG2# */ - +#ifdef HAVE_@KIND@@C@ @type@ npy_@kind@@c@(@type@ x) { return @kind@@c@(x); } +#endif /**end repeat1**/ /**begin repeat1 - * #kind = atan2,hypot,pow,fmod,logaddexp,logaddexp2# + * #kind = atan2,hypot,pow,fmod,copysign# + * #KIND = ATAN2,HYPOT,POW,FMOD,COPYSIGN# */ +#ifdef HAVE_@KIND@@C@ @type@ npy_@kind@@c@(@type@ x, @type@ y) { return @kind@@c@(x, y); } +#endif /**end repeat1**/ +#ifdef HAVE_MODF@C@ @type@ npy_modf@c@(@type@ x, @type@ *iptr) { return modf@c@(x, iptr); } +#endif /**end repeat**/ diff --git a/numpy/core/src/npymath/npy_math_private.h b/numpy/core/src/npymath/npy_math_private.h new file mode 100644 index 000000000..ea7c47fe8 --- /dev/null +++ b/numpy/core/src/npymath/npy_math_private.h @@ -0,0 +1,121 @@ +/* + * + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * from: @(#)fdlibm.h 5.1 93/09/24 + * $FreeBSD$ + */ + +#ifndef _NPY_MATH_PRIVATE_H_ +#define _NPY_MATH_PRIVATE_H_ + +#include <numpy/npy_endian.h> + +/* + * The original fdlibm code used statements like: + * n0 = ((*(int*)&one)>>29)^1; * index of high word * + * ix0 = *(n0+(int*)&x); * high word of x * + * ix1 = *((1-n0)+(int*)&x); * low word of x * + * to dig two 32 bit words out of the 64 bit IEEE floating point + * value. That is non-ANSI, and, moreover, the gcc instruction + * scheduler gets it wrong. We instead use the following macros. + * Unlike the original code, we determine the endianness at compile + * time, not at run time; I don't see much benefit to selecting + * endianness at run time. + */ + +/* + * A union which permits us to convert between a double and two 32 bit + * ints. + */ + +/* XXX: not really, but we already make this assumption elsewhere. Will have to + * fix this at some point */ +#define IEEE_WORD_ORDER NPY_BYTE_ORDER + +#if IEEE_WORD_ORDER == NPY_BIG_ENDIAN + +typedef union +{ + double value; + struct + { + npy_uint32 msw; + npy_uint32 lsw; + } parts; +} ieee_double_shape_type; + +#endif + +#if IEEE_WORD_ORDER == NPY_LITTLE_ENDIAN + +typedef union +{ + double value; + struct + { + npy_uint32 lsw; + npy_uint32 msw; + } parts; +} ieee_double_shape_type; + +#endif + +/* Get two 32 bit ints from a double. */ + +#define EXTRACT_WORDS(ix0,ix1,d) \ +do { \ + ieee_double_shape_type ew_u; \ + ew_u.value = (d); \ + (ix0) = ew_u.parts.msw; \ + (ix1) = ew_u.parts.lsw; \ +} while (0) + +/* Get the more significant 32 bit int from a double. */ + +#define GET_HIGH_WORD(i,d) \ +do { \ + ieee_double_shape_type gh_u; \ + gh_u.value = (d); \ + (i) = gh_u.parts.msw; \ +} while (0) + +/* Get the less significant 32 bit int from a double. */ + +#define GET_LOW_WORD(i,d) \ +do { \ + ieee_double_shape_type gl_u; \ + gl_u.value = (d); \ + (i) = gl_u.parts.lsw; \ +} while (0) + +/* Set the more significant 32 bits of a double from an int. */ + +#define SET_HIGH_WORD(d,v) \ +do { \ + ieee_double_shape_type sh_u; \ + sh_u.value = (d); \ + sh_u.parts.msw = (v); \ + (d) = sh_u.value; \ +} while (0) + +/* Set the less significant 32 bits of a double from an int. */ + +#define SET_LOW_WORD(d,v) \ +do { \ + ieee_double_shape_type sl_u; \ + sl_u.value = (d); \ + sl_u.parts.lsw = (v); \ + (d) = sl_u.value; \ +} while (0) + +#endif /* !_NPY_MATH_PRIVATE_H_ */ diff --git a/numpy/core/src/py3k_notes.txt b/numpy/core/src/py3k_notes.txt new file mode 100644 index 000000000..e31755012 --- /dev/null +++ b/numpy/core/src/py3k_notes.txt @@ -0,0 +1,197 @@ +Notes on making the transition to python 3.x +============================================ + +PyTypeObject +------------ + +The PyTypeObject of py3k is binary compatible with the py2k version and the +old initializers should work. However, there are several considerations to +keep in mind. + +1) Because the first three slots are now part of a struct some compilers issue +warnings if they are initialized in the old way. + +2) The compare slot has been made reserved in order to preserve binary +compatibily while the tp_compare function went away. The tp_richcompare +function has replaced it and we need to use that slot instead. This will +likely require modifications in the searchsorted functions and generic sorts +that currently use the compare function. + +3) The previous numpy practice of initializing the COUNT_ALLOCS slots was +bogus. They are not supposed to be explicitly initialized and were out of +place in any case because an extra base slot was added in python 2.6. + +Because of these facts it was thought better to use #ifdefs to bring the old +initializers up to py3k snuff rather than just fill the tp_richcompare slot. +They also serve to mark the places where changes have been made. The new form +is shown below. Note that explicit initialization can stop once none of the +remaining entries are non-zero, because zero is the default value that +variables with non-local linkage receive. + + +NPY_NO_EXPORT PyTypeObject Foo_Type = { +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(0,0) +#else + PyObject_HEAD_INIT(0) + 0, /* ob_size */ +#endif + "numpy.foo" /* tp_name */ + 0, /* tp_basicsize */ + 0, /* tp_itemsize */ + /* methods */ + 0, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + (void *)0, /* tp_reserved */ +#else + 0, /* tp_compare */ +#endif + 0, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + 0, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + 0, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ + 0 /* tp_version_tag (2.6) */ +}; + +checklist of types having tp_compare but no tp_richcompare + +1) multiarray/flagsobject.c + +PyNumberMethods +--------------- + +Types with tp_as_number defined + +1) multiarray/arrayobject.c + +The PyNumberMethods struct has changed enough that it looks easiest to just +have an alternate version. Note that np_divide, np_long, np_oct, np_hex, and +np_inplace_divide have gone away. The slot np_int is what np_long used to be, +tp_divide is now tp_floor_divide, and np_inplace_divide is now +np_inplace_floor_divide. We will also have to make sure the *_true_divide +variants are defined. This should also be done for python < 3.x, but that +introduces a requirement for the Py_TPFLAGS_HAVE_CLASS in the type flag. + +/* + * Number implementations must check *both* arguments for proper type and + * implement the necessary conversions in the slot functions themselves. +*/ +PyNumberMethods foo_number_methods = { + (binaryfunc)0, /* nb_add */ + (binaryfunc)0, /* nb_subtract */ + (binaryfunc)0, /* nb_multiply */ + (binaryfunc)0, /* nb_remainder */ + (binaryfunc)0, /* nb_divmod */ + (ternaryfunc)0, /* nb_power */ + (unaryfunc)0, /* nb_negative */ + (unaryfunc)0, /* nb_positive */ + (unaryfunc)0, /* nb_absolute */ + (inquiry)0, /* nb_bool, nee nb_nonzero */ + (unaryfunc)0, /* nb_invert */ + (binaryfunc)0, /* nb_lshift */ + (binaryfunc)0, /* nb_rshift */ + (binaryfunc)0, /* nb_and */ + (binaryfunc)0, /* nb_xor */ + (binaryfunc)0, /* nb_or */ + (unaryfunc)0, /* nb_int */ + (void *)0, /* nb_reserved, nee nb_long */ + (unaryfunc)0, /* nb_float */ + (binaryfunc)0, /* nb_inplace_add */ + (binaryfunc)0, /* nb_inplace_subtract */ + (binaryfunc)0, /* nb_inplace_multiply */ + (binaryfunc)0, /* nb_inplace_remainder */ + (ternaryfunc)0, /* nb_inplace_power */ + (binaryfunc)0, /* nb_inplace_lshift */ + (binaryfunc)0, /* nb_inplace_rshift */ + (binaryfunc)0, /* nb_inplace_and */ + (binaryfunc)0, /* nb_inplace_xor */ + (binaryfunc)0, /* nb_inplace_or */ + (binaryfunc)0, /* nb_floor_divide */ + (binaryfunc)0, /* nb_true_divide */ + (binaryfunc)0, /* nb_inplace_floor_divide */ + (binaryfunc)0, /* nb_inplace_true_divide */ + (unaryfunc)0 /* nb_index */ +}; + +PySequenceMethods +----------------- + +Types with tp_as_sequence defined + +1) multiarray/descriptor.c +2) multiarray/scalartypes.c.src +3) multiarray/arrayobject.c + +PySequenceMethods in py3k are binary compatible with py2k, but some of the +slots have gone away. I suspect this means some functions need redefining so +the semantics of the slots needs to be checked. + +PySequenceMethods foo_sequence_methods = { + (lenfunc)0, /* sq_length */ + (binaryfunc)0, /* sq_concat */ + (ssizeargfunc)0, /* sq_repeat */ + (ssizeargfunc)0, /* sq_item */ + (void *)0, /* nee sq_slice */ + (ssizeobjargproc)0, /* sq_ass_item */ + (void *)0, /* nee sq_ass_slice */ + (objobjproc)0, /* sq_contains */ + (binaryfunc)0, /* sq_inplace_concat */ + (ssizeargfunc)0 /* sq_inplace_repeat */ +}; + +PyMappingMethods +---------------- + +Types with tp_as_mapping defined + +1) multiarray/descriptor.c +2) multiarray/iterators.c +3) multiarray/scalartypes.c.src +4) multiarray/flagsobject.c +5) multiarray/arrayobject.c + +PyMappingMethods in py3k look to be the same as in py2k. The semantics +of the slots needs to be checked. + +PyMappingMethods foo_mapping_methods = { + (lenfunc)0, /* mp_length */ + (binaryfunc)0, /* mp_subscript */ + (objobjargproc)0 /* mp_ass_subscript */ +}; + diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src index 25e15bca1..10cfa8716 100644 --- a/numpy/core/src/umath/loops.c.src +++ b/numpy/core/src/umath/loops.c.src @@ -1044,6 +1044,16 @@ NPY_NO_EXPORT void } /**end repeat1**/ +NPY_NO_EXPORT void +@TYPE@_copysign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) +{ + BINARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + const @type@ in2 = *(@type@ *)ip2; + *((@type@ *)op1)= npy_copysign@c@(in1, in2); + } +} + /**begin repeat1 * #kind = maximum, minimum# * #OP = >=, <=# @@ -1265,9 +1275,18 @@ C@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; - @type@ d = in2r*in2r + in2i*in2i; - ((@type@ *)op1)[0] = (in1r*in2r + in1i*in2i)/d; - ((@type@ *)op1)[1] = (in1i*in2r - in1r*in2i)/d; + if (npy_fabs@c@(in2r) >= npy_fabs@c@(in2i)) { + const @type@ rat = in2i/in2r; + const @type@ scl = 1.0@c@/(in2r + in2i*rat); + ((@type@ *)op1)[0] = (in1r + in1i*rat)*scl; + ((@type@ *)op1)[1] = (in1i - in1r*rat)*scl; + } + else { + const @type@ rat = in2r/in2i; + const @type@ scl = 1.0@c@/(in2i + in2r*rat); + ((@type@ *)op1)[0] = (in1r*rat + in1i)*scl; + ((@type@ *)op1)[1] = (in1i*rat - in1r)*scl; + } } } @@ -1279,9 +1298,16 @@ C@TYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSE const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; - @type@ d = in2r*in2r + in2i*in2i; - ((@type@ *)op1)[0] = npy_floor@c@((in1r*in2r + in1i*in2i)/d); - ((@type@ *)op1)[1] = 0; + if (npy_fabs@c@(in2r) >= npy_fabs@c@(in2i)) { + const @type@ rat = in2i/in2r; + ((@type@ *)op1)[0] = npy_floor@c@((in1r + in1i*rat)/(in2r + in2i*rat)); + ((@type@ *)op1)[1] = 0; + } + else { + const @type@ rat = in2r/in2i; + ((@type@ *)op1)[0] = npy_floor@c@((in1r*rat + in1i)/(in2i + in2r*rat)); + ((@type@ *)op1)[1] = 0; + } } } diff --git a/numpy/core/src/umath/loops.h b/numpy/core/src/umath/loops.h index 9de4c5893..bf33ea88c 100644 --- a/numpy/core/src/umath/loops.h +++ b/numpy/core/src/umath/loops.h @@ -1527,6 +1527,9 @@ NPY_NO_EXPORT void FLOAT_signbit(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); +NPY_NO_EXPORT void +FLOAT_copysign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); + NPY_NO_EXPORT void FLOAT_maximum(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); @@ -1668,6 +1671,9 @@ NPY_NO_EXPORT void DOUBLE_signbit(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); +NPY_NO_EXPORT void +DOUBLE_copysign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); + NPY_NO_EXPORT void DOUBLE_maximum(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); @@ -1809,6 +1815,9 @@ NPY_NO_EXPORT void LONGDOUBLE_signbit(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); +NPY_NO_EXPORT void +LONGDOUBLE_copysign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); + NPY_NO_EXPORT void LONGDOUBLE_maximum(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)); diff --git a/numpy/core/src/umath/ufunc_object.c b/numpy/core/src/umath/ufunc_object.c index 97fcc8124..94152ccb7 100644 --- a/numpy/core/src/umath/ufunc_object.c +++ b/numpy/core/src/umath/ufunc_object.c @@ -239,6 +239,122 @@ static char *_types_msg = "function not supported for these types, " \ "and can't coerce safely to supported types"; /* + * This function analyzes the input arguments + * and determines an appropriate __array_prepare__ function to call + * for the outputs. + * + * If an output argument is provided, then it is wrapped + * with its own __array_prepare__ not with the one determined by + * the input arguments. + * + * if the provided output argument is already an ndarray, + * the wrapping function is None (which means no wrapping will + * be done --- not even PyArray_Return). + * + * A NULL is placed in output_wrap for outputs that + * should just have PyArray_Return called. + */ +static void +_find_array_prepare(PyObject *args, PyObject **output_wrap, int nin, int nout) +{ + Py_ssize_t nargs; + int i; + int np = 0; + PyObject *with_wrap[NPY_MAXARGS], *wraps[NPY_MAXARGS]; + PyObject *obj, *wrap = NULL; + + nargs = PyTuple_GET_SIZE(args); + for (i = 0; i < nin; i++) { + obj = PyTuple_GET_ITEM(args, i); + if (PyArray_CheckExact(obj) || PyArray_IsAnyScalar(obj)) { + continue; + } + wrap = PyObject_GetAttrString(obj, "__array_prepare__"); + if (wrap) { + if (PyCallable_Check(wrap)) { + with_wrap[np] = obj; + wraps[np] = wrap; + ++np; + } + else { + Py_DECREF(wrap); + wrap = NULL; + } + } + else { + PyErr_Clear(); + } + } + if (np > 0) { + /* If we have some wraps defined, find the one of highest priority */ + wrap = wraps[0]; + if (np > 1) { + double maxpriority = PyArray_GetPriority(with_wrap[0], + PyArray_SUBTYPE_PRIORITY); + for (i = 1; i < np; ++i) { + double priority = PyArray_GetPriority(with_wrap[i], + PyArray_SUBTYPE_PRIORITY); + if (priority > maxpriority) { + maxpriority = priority; + Py_DECREF(wrap); + wrap = wraps[i]; + } + else { + Py_DECREF(wraps[i]); + } + } + } + } + + /* + * Here wrap is the wrapping function determined from the + * input arrays (could be NULL). + * + * For all the output arrays decide what to do. + * + * 1) Use the wrap function determined from the input arrays + * This is the default if the output array is not + * passed in. + * + * 2) Use the __array_prepare__ method of the output object. + * This is special cased for + * exact ndarray so that no PyArray_Return is + * done in that case. + */ + for (i = 0; i < nout; i++) { + int j = nin + i; + int incref = 1; + output_wrap[i] = wrap; + if (j < nargs) { + obj = PyTuple_GET_ITEM(args, j); + if (obj == Py_None) { + continue; + } + if (PyArray_CheckExact(obj)) { + output_wrap[i] = Py_None; + } + else { + PyObject *owrap = PyObject_GetAttrString(obj, + "__array_prepare__"); + incref = 0; + if (!(owrap) || !(PyCallable_Check(owrap))) { + Py_XDECREF(owrap); + owrap = wrap; + incref = 1; + PyErr_Clear(); + } + output_wrap[i] = owrap; + } + } + if (incref) { + Py_XINCREF(output_wrap[i]); + } + } + Py_XDECREF(wrap); + return; +} + +/* * Called for non-NULL user-defined functions. * The object should be a CObject pointing to a linked-list of functions * storing the function, data, and signature of all user-defined functions. @@ -1059,6 +1175,7 @@ construct_arrays(PyUFuncLoopObject *loop, PyObject *args, PyArrayObject **mps, npy_intp temp_dims[NPY_MAXDIMS]; npy_intp *out_dims; int out_nd; + PyObject *wraparr[NPY_MAXARGS]; /* Check number of arguments */ nargs = PyTuple_Size(args); @@ -1337,16 +1454,60 @@ construct_arrays(PyUFuncLoopObject *loop, PyObject *args, PyArrayObject **mps, return -1; } - /* Recover mps[i]. */ - if (self->core_enabled) { - PyArrayObject *ao = mps[i]; - mps[i] = (PyArrayObject *)mps[i]->base; - Py_DECREF(ao); - } + /* Recover mps[i]. */ + if (self->core_enabled) { + PyArrayObject *ao = mps[i]; + mps[i] = (PyArrayObject *)mps[i]->base; + Py_DECREF(ao); + } } /* + * Use __array_prepare__ on all outputs + * if present on one of the input arguments. + * If present for multiple inputs: + * use __array_prepare__ of input object with largest + * __array_priority__ (default = 0.0) + * + * Exception: we should not wrap outputs for items already + * passed in as output-arguments. These items should either + * be left unwrapped or wrapped by calling their own __array_prepare__ + * routine. + * + * For each output argument, wrap will be either + * NULL --- call PyArray_Return() -- default if no output arguments given + * None --- array-object passed in don't call PyArray_Return + * method --- the __array_prepare__ method to call. + */ + _find_array_prepare(args, wraparr, loop->ufunc->nin, loop->ufunc->nout); + + /* wrap outputs */ + for (i = 0; i < loop->ufunc->nout; i++) { + int j = loop->ufunc->nin+i; + PyObject *wrap; + wrap = wraparr[i]; + if (wrap != NULL) { + if (wrap == Py_None) { + Py_DECREF(wrap); + continue; + } + PyObject *res = PyObject_CallFunction(wrap, "O(OOi)", + mps[j], loop->ufunc, args, i); + Py_DECREF(wrap); + if ((res == NULL) || (res == Py_None)) { + if (!PyErr_Occurred()){ + PyErr_SetString(PyExc_TypeError, + "__array_prepare__ must return an ndarray or subclass thereof"); + } + return -1; + } + Py_DECREF(mps[j]); + mps[j] = (PyArrayObject *)res; + } + } + + /* * If any of different type, or misaligned or swapped * then must use buffers */ @@ -3827,7 +3988,10 @@ ufunc_repr(PyUFuncObject *self) } -/* -------------------------------------------------------- */ +/****************************************************************************** + *** UFUNC METHODS *** + *****************************************************************************/ + /* * op.outer(a,b) is equivalent to op(a[:,NewAxis,NewAxis,etc.],b) @@ -3962,6 +4126,10 @@ static struct PyMethodDef ufunc_methods[] = { }; +/****************************************************************************** + *** UFUNC GETSET *** + *****************************************************************************/ + /* construct the string y1,y2,...,yn */ static PyObject * @@ -4000,7 +4168,8 @@ _typecharfromnum(int num) { static PyObject * ufunc_get_doc(PyUFuncObject *self) { - /* Put docstring first or FindMethod finds it... could so some + /* + * Put docstring first or FindMethod finds it... could so some * introspection on name and nin + nout to automate the first part * of it the doc string shouldn't need the calling convention * construct name(x1, x2, ...,[ out1, out2, ...]) __doc__ @@ -4148,65 +4317,68 @@ static PyGetSetDef ufunc_getset[] = { {NULL, NULL, NULL, NULL, NULL}, /* Sentinel */ }; + +/****************************************************************************** + *** UFUNC TYPE OBJECT *** + *****************************************************************************/ + NPY_NO_EXPORT PyTypeObject PyUFunc_Type = { - PyObject_HEAD_INIT(0) - 0, /* ob_size */ - "numpy.ufunc", /* tp_name */ - sizeof(PyUFuncObject), /* tp_basicsize */ - 0, /* tp_itemsize */ +#if defined(NPY_PY3K) + PyVarObject_HEAD_INIT(NULL, 0) +#else + PyObject_HEAD_INIT(NULL) + 0, /* ob_size */ +#endif + "numpy.ufunc", /* tp_name */ + sizeof(PyUFuncObject), /* tp_basicsize */ + 0, /* tp_itemsize */ /* methods */ - (destructor)ufunc_dealloc, /* tp_dealloc */ - (printfunc)0, /* tp_print */ - (getattrfunc)0, /* tp_getattr */ - (setattrfunc)0, /* tp_setattr */ - (cmpfunc)0, /* tp_compare */ - (reprfunc)ufunc_repr, /* tp_repr */ - 0, /* tp_as_number */ - 0, /* tp_as_sequence */ - 0, /* tp_as_mapping */ - (hashfunc)0, /* tp_hash */ - (ternaryfunc)ufunc_generic_call, /* tp_call */ - (reprfunc)ufunc_repr, /* tp_str */ - 0, /* tp_getattro */ - 0, /* tp_setattro */ - 0, /* tp_as_buffer */ - Py_TPFLAGS_DEFAULT, /* tp_flags */ - NULL, /* tp_doc */ /* was Ufunctype__doc__ */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - ufunc_methods, /* tp_methods */ - 0, /* tp_members */ - ufunc_getset, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - 0, /* tp_init */ - 0, /* tp_alloc */ - 0, /* tp_new */ - 0, /* tp_free */ - 0, /* tp_is_gc */ - 0, /* tp_bases */ - 0, /* tp_mro */ - 0, /* tp_cache */ - 0, /* tp_subclasses */ - 0, /* tp_weaklist */ - 0, /* tp_del */ - -#ifdef COUNT_ALLOCS - /* these must be last and never explicitly initialized */ - 0, /* tp_allocs */ - 0, /* tp_frees */ - 0, /* tp_maxalloc */ - 0, /* tp_prev */ - 0, /* *tp_next */ + (destructor)ufunc_dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ +#if defined(NPY_PY3K) + 0, /* tp_reserved */ +#else + 0, /* tp_compare */ #endif + (reprfunc)ufunc_repr, /* tp_repr */ + 0, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + (ternaryfunc)ufunc_generic_call, /* tp_call */ + (reprfunc)ufunc_repr, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ + Py_TPFLAGS_DEFAULT, /* tp_flags */ + 0, /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + ufunc_methods, /* tp_methods */ + 0, /* tp_members */ + ufunc_getset, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + 0, /* tp_init */ + 0, /* tp_alloc */ + 0, /* tp_new */ + 0, /* tp_free */ + 0, /* tp_is_gc */ + 0, /* tp_bases */ + 0, /* tp_mro */ + 0, /* tp_cache */ + 0, /* tp_subclasses */ + 0, /* tp_weaklist */ + 0, /* tp_del */ }; /* End of code for ufunc objects */ -/* -------------------------------------------------------- */ diff --git a/numpy/core/tests/test_defmatrix.py b/numpy/core/tests/test_defmatrix.py index e9f0d9a7f..40728bd29 100644 --- a/numpy/core/tests/test_defmatrix.py +++ b/numpy/core/tests/test_defmatrix.py @@ -1,5 +1,6 @@ from numpy.testing import * from numpy.core import * +from numpy.core.defmatrix import matrix_power import numpy as np class TestCtor(TestCase): @@ -358,6 +359,15 @@ class TestNewScalarIndexing(TestCase): assert_array_equal(x[:,[1,0]],x[:,::-1]) assert_array_equal(x[[2,1,0],:],x[::-1,:]) +class TestPower(TestCase): + def test_returntype(self): + a = array([[0,1],[0,0]]) + assert type(matrix_power(a, 2)) is ndarray + a = mat(a) + assert type(matrix_power(a, 2)) is matrix + + def test_list(self): + assert_array_equal(matrix_power([[0, 1], [0, 0]], 2), [[0, 0], [0, 0]]) if __name__ == "__main__": run_module_suite() diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py index 7022ef14d..57f1bd4c6 100644 --- a/numpy/core/tests/test_multiarray.py +++ b/numpy/core/tests/test_multiarray.py @@ -4,6 +4,7 @@ import os import numpy as np from numpy.testing import * from numpy.core import * +from numpy.core.multiarray_tests import test_neighborhood_iterator, test_neighborhood_iterator_oob from test_print import in_foreign_locale @@ -279,6 +280,24 @@ class TestMethods(TestCase): self.failUnlessRaises(ValueError, lambda: a.transpose(0,1,2)) def test_sort(self): + # test ordering for floats and complex containing nans. It is only + # necessary to check the lessthan comparison, so sorts that + # only follow the insertion sort path are sufficient. We only + # test doubles and complex doubles as the logic is the same. + + # check doubles + msg = "Test real sort order with nans" + a = np.array([np.nan, 1, 0]) + b = sort(a) + assert_equal(b, a[::-1], msg) + # check complex + msg = "Test complex sort order with nans" + a = np.zeros(9, dtype=np.complex128) + a.real += [np.nan, np.nan, np.nan, 1, 0, 1, 1, 0, 0] + a.imag += [np.nan, 1, 0, np.nan, np.nan, 1, 0, 1, 0] + b = sort(a) + assert_equal(b, a[::-1], msg) + # all c scalar sorts use the same code with different types # so it suffices to run a quick check with one type. The number # of sorted items must be greater than ~50 to check the actual @@ -466,6 +485,33 @@ class TestMethods(TestCase): a = np.array(['aaaaaaaaa' for i in range(100)], dtype=np.unicode) assert_equal(a.argsort(kind='m'), r) + def test_searchsorted(self): + # test for floats and complex containing nans. The logic is the + # same for all float types so only test double types for now. + # The search sorted routines use the compare functions for the + # array type, so this checks if that is consistent with the sort + # order. + + # check double + a = np.array([np.nan, 1, 0]) + a = np.array([0, 1, np.nan]) + msg = "Test real searchsorted with nans, side='l'" + b = a.searchsorted(a, side='l') + assert_equal(b, np.arange(3), msg) + msg = "Test real searchsorted with nans, side='r'" + b = a.searchsorted(a, side='r') + assert_equal(b, np.arange(1,4), msg) + # check double complex + a = np.zeros(9, dtype=np.complex128) + a.real += [0, 0, 1, 1, 0, 1, np.nan, np.nan, np.nan] + a.imag += [0, 1, 0, 1, np.nan, np.nan, 0, 1, np.nan] + msg = "Test complex searchsorted with nans, side='l'" + b = a.searchsorted(a, side='l') + assert_equal(b, np.arange(9), msg) + msg = "Test complex searchsorted with nans, side='r'" + b = a.searchsorted(a, side='r') + assert_equal(b, np.arange(1,10), msg) + def test_flatten(self): x0 = np.array([[1,2,3],[4,5,6]], np.int32) x1 = np.array([[[1,2],[3,4]],[[5,6],[7,8]]], np.int32) @@ -1058,6 +1104,282 @@ class TestChoose(TestCase): A = np.choose(self.ind, (self.x, self.y2)) assert_equal(A, [[2,2,3],[2,2,3]]) - +def can_use_decimal(): + try: + from decimal import Decimal + return True + except ImportError: + return False + +# TODO: test for multidimensional +NEIGH_MODE = {'zero': 0, 'one': 1, 'constant': 2, 'circular': 3, 'mirror': 4} +class TestNeighborhoodIter(TestCase): + # Simple, 2d tests + def _test_simple2d(self, dt): + # Test zero and one padding for simple data type + x = np.array([[0, 1], [2, 3]], dtype=dt) + r = [np.array([[0, 0, 0], [0, 0, 1]], dtype=dt), + np.array([[0, 0, 0], [0, 1, 0]], dtype=dt), + np.array([[0, 0, 1], [0, 2, 3]], dtype=dt), + np.array([[0, 1, 0], [2, 3, 0]], dtype=dt)] + l = test_neighborhood_iterator(x, [-1, 0, -1, 1], x[0], NEIGH_MODE['zero']) + assert_array_equal(l, r) + + r = [np.array([[1, 1, 1], [1, 0, 1]], dtype=dt), + np.array([[1, 1, 1], [0, 1, 1]], dtype=dt), + np.array([[1, 0, 1], [1, 2, 3]], dtype=dt), + np.array([[0, 1, 1], [2, 3, 1]], dtype=dt)] + l = test_neighborhood_iterator(x, [-1, 0, -1, 1], x[0], NEIGH_MODE['one']) + assert_array_equal(l, r) + + r = [np.array([[4, 4, 4], [4, 0, 1]], dtype=dt), + np.array([[4, 4, 4], [0, 1, 4]], dtype=dt), + np.array([[4, 0, 1], [4, 2, 3]], dtype=dt), + np.array([[0, 1, 4], [2, 3, 4]], dtype=dt)] + l = test_neighborhood_iterator(x, [-1, 0, -1, 1], 4, NEIGH_MODE['constant']) + assert_array_equal(l, r) + + def test_simple2d(self): + self._test_simple2d(np.float) + + @dec.skipif(not can_use_decimal(), + "Skip neighborhood iterator tests for decimal objects " \ + "(decimal module not available") + def test_simple2d_object(self): + from decimal import Decimal + self._test_simple2d(Decimal) + + def _test_mirror2d(self, dt): + x = np.array([[0, 1], [2, 3]], dtype=dt) + r = [np.array([[0, 0, 1], [0, 0, 1]], dtype=dt), + np.array([[0, 1, 1], [0, 1, 1]], dtype=dt), + np.array([[0, 0, 1], [2, 2, 3]], dtype=dt), + np.array([[0, 1, 1], [2, 3, 3]], dtype=dt)] + l = test_neighborhood_iterator(x, [-1, 0, -1, 1], x[0], NEIGH_MODE['mirror']) + assert_array_equal(l, r) + + def test_mirror2d(self): + self._test_mirror2d(np.float) + + @dec.skipif(not can_use_decimal(), + "Skip neighborhood iterator tests for decimal objects " \ + "(decimal module not available") + def test_mirror2d_object(self): + from decimal import Decimal + self._test_mirror2d(Decimal) + + # Simple, 1d tests + def _test_simple(self, dt): + # Test padding with constant values + x = np.linspace(1, 5, 5).astype(dt) + r = [[0, 1, 2], [1, 2, 3], [2, 3, 4], [3, 4, 5], [4, 5, 0]] + l = test_neighborhood_iterator(x, [-1, 1], x[0], NEIGH_MODE['zero']) + assert_array_equal(l, r) + + r = [[1, 1, 2], [1, 2, 3], [2, 3, 4], [3, 4, 5], [4, 5, 1]] + l = test_neighborhood_iterator(x, [-1, 1], x[0], NEIGH_MODE['one']) + assert_array_equal(l, r) + + r = [[x[4], 1, 2], [1, 2, 3], [2, 3, 4], [3, 4, 5], [4, 5, x[4]]] + l = test_neighborhood_iterator(x, [-1, 1], x[4], NEIGH_MODE['constant']) + assert_array_equal(l, r) + + def test_simple_float(self): + self._test_simple(np.float) + + @dec.skipif(not can_use_decimal(), + "Skip neighborhood iterator tests for decimal objects " \ + "(decimal module not available") + def test_simple_object(self): + from decimal import Decimal + self._test_simple(Decimal) + + # Test mirror modes + def _test_mirror(self, dt): + x = np.linspace(1, 5, 5).astype(dt) + r = np.array([[2, 1, 1, 2, 3], [1, 1, 2, 3, 4], [1, 2, 3, 4, 5], + [2, 3, 4, 5, 5], [3, 4, 5, 5, 4]], dtype=dt) + l = test_neighborhood_iterator(x, [-2, 2], x[1], NEIGH_MODE['mirror']) + self.failUnless([i.dtype == dt for i in l]) + assert_array_equal(l, r) + + def test_mirror(self): + self._test_mirror(np.float) + + @dec.skipif(not can_use_decimal(), + "Skip neighborhood iterator tests for decimal objects " \ + "(decimal module not available") + def test_mirror_object(self): + from decimal import Decimal + self._test_mirror(Decimal) + + # Circular mode + def _test_circular(self, dt): + x = np.linspace(1, 5, 5).astype(dt) + r = np.array([[4, 5, 1, 2, 3], [5, 1, 2, 3, 4], [1, 2, 3, 4, 5], + [2, 3, 4, 5, 1], [3, 4, 5, 1, 2]], dtype=dt) + l = test_neighborhood_iterator(x, [-2, 2], x[0], NEIGH_MODE['circular']) + assert_array_equal(l, r) + + def test_circular(self): + self._test_circular(np.float) + + @dec.skipif(not can_use_decimal(), + "Skip neighborhood iterator tests for decimal objects " \ + "(decimal module not available") + def test_circular_object(self): + from decimal import Decimal + self._test_circular(Decimal) + +# Test stacking neighborhood iterators +class TestStackedNeighborhoodIter(TestCase): + # Simple, 1d test: stacking 2 constant-padded neigh iterators + def test_simple_const(self): + dt = np.float64 + # Test zero and one padding for simple data type + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([0], dtype=dt), + np.array([0], dtype=dt), + np.array([1], dtype=dt), + np.array([2], dtype=dt), + np.array([3], dtype=dt), + np.array([0], dtype=dt), + np.array([0], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-2, 4], NEIGH_MODE['zero'], + [0, 0], NEIGH_MODE['zero']) + assert_array_equal(l, r) + + r = [np.array([1, 0, 1], dtype=dt), + np.array([0, 1, 2], dtype=dt), + np.array([1, 2, 3], dtype=dt), + np.array([2, 3, 0], dtype=dt), + np.array([3, 0, 1], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-1, 3], NEIGH_MODE['zero'], + [-1, 1], NEIGH_MODE['one']) + assert_array_equal(l, r) + + # 2nd simple, 1d test: stacking 2 neigh iterators, mixing const padding and + # mirror padding + def test_simple_mirror(self): + dt = np.float64 + # Stacking zero on top of mirror + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([0, 1, 1], dtype=dt), + np.array([1, 1, 2], dtype=dt), + np.array([1, 2, 3], dtype=dt), + np.array([2, 3, 3], dtype=dt), + np.array([3, 3, 0], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-1, 3], NEIGH_MODE['mirror'], + [-1, 1], NEIGH_MODE['zero']) + assert_array_equal(l, r) + + # Stacking mirror on top of zero + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([1, 0, 0], dtype=dt), + np.array([0, 0, 1], dtype=dt), + np.array([0, 1, 2], dtype=dt), + np.array([1, 2, 3], dtype=dt), + np.array([2, 3, 0], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-1, 3], NEIGH_MODE['zero'], + [-2, 0], NEIGH_MODE['mirror']) + assert_array_equal(l, r) + + # Stacking mirror on top of zero: 2nd + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([0, 1, 2], dtype=dt), + np.array([1, 2, 3], dtype=dt), + np.array([2, 3, 0], dtype=dt), + np.array([3, 0, 0], dtype=dt), + np.array([0, 0, 3], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-1, 3], NEIGH_MODE['zero'], + [0, 2], NEIGH_MODE['mirror']) + assert_array_equal(l, r) + + # Stacking mirror on top of zero: 3rd + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([1, 0, 0, 1, 2], dtype=dt), + np.array([0, 0, 1, 2, 3], dtype=dt), + np.array([0, 1, 2, 3, 0], dtype=dt), + np.array([1, 2, 3, 0, 0], dtype=dt), + np.array([2, 3, 0, 0, 3], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-1, 3], NEIGH_MODE['zero'], + [-2, 2], NEIGH_MODE['mirror']) + assert_array_equal(l, r) + + # 3rd simple, 1d test: stacking 2 neigh iterators, mixing const padding and + # circular padding + def test_simple_circular(self): + dt = np.float64 + # Stacking zero on top of mirror + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([0, 3, 1], dtype=dt), + np.array([3, 1, 2], dtype=dt), + np.array([1, 2, 3], dtype=dt), + np.array([2, 3, 1], dtype=dt), + np.array([3, 1, 0], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-1, 3], NEIGH_MODE['circular'], + [-1, 1], NEIGH_MODE['zero']) + assert_array_equal(l, r) + + # Stacking mirror on top of zero + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([3, 0, 0], dtype=dt), + np.array([0, 0, 1], dtype=dt), + np.array([0, 1, 2], dtype=dt), + np.array([1, 2, 3], dtype=dt), + np.array([2, 3, 0], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-1, 3], NEIGH_MODE['zero'], + [-2, 0], NEIGH_MODE['circular']) + assert_array_equal(l, r) + + # Stacking mirror on top of zero: 2nd + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([0, 1, 2], dtype=dt), + np.array([1, 2, 3], dtype=dt), + np.array([2, 3, 0], dtype=dt), + np.array([3, 0, 0], dtype=dt), + np.array([0, 0, 1], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-1, 3], NEIGH_MODE['zero'], + [0, 2], NEIGH_MODE['circular']) + assert_array_equal(l, r) + + # Stacking mirror on top of zero: 3rd + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([3, 0, 0, 1, 2], dtype=dt), + np.array([0, 0, 1, 2, 3], dtype=dt), + np.array([0, 1, 2, 3, 0], dtype=dt), + np.array([1, 2, 3, 0, 0], dtype=dt), + np.array([2, 3, 0, 0, 1], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [-1, 3], NEIGH_MODE['zero'], + [-2, 2], NEIGH_MODE['circular']) + assert_array_equal(l, r) + + # 4th simple, 1d test: stacking 2 neigh iterators, but with lower iterator + # being strictly within the array + def test_simple_strict_within(self): + dt = np.float64 + # Stacking zero on top of zero, first neighborhood strictly inside the + # array + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([1, 2, 3, 0], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [1, 1], NEIGH_MODE['zero'], + [-1, 2], NEIGH_MODE['zero']) + assert_array_equal(l, r) + + # Stacking mirror on top of zero, first neighborhood strictly inside the + # array + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([1, 2, 3, 3], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [1, 1], NEIGH_MODE['zero'], + [-1, 2], NEIGH_MODE['mirror']) + assert_array_equal(l, r) + + # Stacking mirror on top of zero, first neighborhood strictly inside the + # array + x = np.array([1, 2, 3], dtype=dt) + r = [np.array([1, 2, 3, 1], dtype=dt)] + l = test_neighborhood_iterator_oob(x, [1, 1], NEIGH_MODE['zero'], + [-1, 2], NEIGH_MODE['circular']) + assert_array_equal(l, r) if __name__ == "__main__": run_module_suite() diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py index 832b2893f..206c06e66 100644 --- a/numpy/core/tests/test_numeric.py +++ b/numpy/core/tests/test_numeric.py @@ -896,19 +896,20 @@ class _TestCorrelate(TestCase): def test_float(self): self._setup(np.float) - z = np.correlate(self.x, self.y, 'full') - assert_array_almost_equal(z, self.z1) - z = np.correlate(self.y, self.x, 'full') + z = np.correlate(self.x, self.y, 'full', old_behavior=self.old_behavior) assert_array_almost_equal(z, self.z1) + z = np.correlate(self.y, self.x, 'full', old_behavior=self.old_behavior) + assert_array_almost_equal(z, self.z2) def test_object(self): self._setup(Decimal) - z = np.correlate(self.x, self.y, 'full') - assert_array_almost_equal(z, self.z1) - z = np.correlate(self.y, self.x, 'full') + z = np.correlate(self.x, self.y, 'full', old_behavior=self.old_behavior) assert_array_almost_equal(z, self.z1) + z = np.correlate(self.y, self.x, 'full', old_behavior=self.old_behavior) + assert_array_almost_equal(z, self.z2) class TestCorrelate(_TestCorrelate): + old_behavior = True def _setup(self, dt): # correlate uses an unconventional definition so that correlate(a, b) # == correlate(b, a), so force the corresponding outputs to be the same @@ -916,6 +917,7 @@ class TestCorrelate(_TestCorrelate): _TestCorrelate._setup(self, dt) self.z2 = self.z1 + @dec.deprecated() def test_complex(self): x = np.array([1, 2, 3, 4+1j], dtype=np.complex) y = np.array([-1, -2j, 3+1j], dtype=np.complex) @@ -923,7 +925,16 @@ class TestCorrelate(_TestCorrelate): z = np.correlate(x, y, 'full') assert_array_almost_equal(z, r_z) -class TestAcorrelate(_TestCorrelate): + @dec.deprecated() + def test_float(self): + _TestCorrelate.test_float(self) + + @dec.deprecated() + def test_object(self): + _TestCorrelate.test_object(self) + +class TestCorrelateNew(_TestCorrelate): + old_behavior = False def test_complex(self): x = np.array([1, 2, 3, 4+1j], dtype=np.complex) y = np.array([-1, -2j, 3+1j], dtype=np.complex) @@ -932,8 +943,24 @@ class TestAcorrelate(_TestCorrelate): #assert_array_almost_equal(z, r_z) r_z = r_z[::-1].conjugate() - z = np.acorrelate(y, x, 'full') + z = np.correlate(y, x, 'full', old_behavior=self.old_behavior) assert_array_almost_equal(z, r_z) +class TestArgwhere: + def test_2D(self): + x = np.arange(6).reshape((2, 3)) + assert_array_equal(np.argwhere(x > 1), + [[0, 2], + [1, 0], + [1, 1], + [1, 2]]) + + def test_list(self): + assert_equal(np.argwhere([4, 0, 2, 1, 3]), [[0], [2], [3], [4]]) + + def test_masked_array(self): + a = np.ma.array([0, 1, 2, 3], mask=[0, 0, 1, 0]) + assert_equal(np.argwhere(a), [[1], [3]]) + if __name__ == "__main__": run_module_suite() diff --git a/numpy/core/tests/test_numerictypes.py b/numpy/core/tests/test_numerictypes.py index 4e0bb462b..56ed4dbb1 100644 --- a/numpy/core/tests/test_numerictypes.py +++ b/numpy/core/tests/test_numerictypes.py @@ -338,13 +338,13 @@ class TestEmptyField(TestCase): class TestCommonType(TestCase): def test_scalar_loses1(self): - res = np.find_common_type(['f4','f4','i4'],['f8']) + res = np.find_common_type(['f4','f4','i2'],['f8']) assert(res == 'f4') def test_scalar_loses2(self): res = np.find_common_type(['f4','f4'],['i8']) assert(res == 'f4') def test_scalar_wins(self): - res = np.find_common_type(['f4','f4','i4'],['c8']) + res = np.find_common_type(['f4','f4','i2'],['c8']) assert(res == 'c8') def test_scalar_wins2(self): res = np.find_common_type(['u4','i4','i4'],['f4']) diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py index 94b10edb1..abea0a222 100644 --- a/numpy/core/tests/test_umath.py +++ b/numpy/core/tests/test_umath.py @@ -10,6 +10,29 @@ class TestDivision(TestCase): assert_equal(x // 100, [0, 0, 0, 1, -1, -1, -1, -1, -2]) assert_equal(x % 100, [5, 10, 90, 0, 95, 90, 10, 0, 80]) + def test_division_complex(self): + # check that implementation is correct + msg = "Complex division implementation check" + x = np.array([1. + 1.*1j, 1. + .5*1j, 1. + 2.*1j], dtype=np.complex128) + assert_almost_equal(x**2/x, x, err_msg=msg) + # check overflow, underflow + msg = "Complex division overflow/underflow check" + x = np.array([1.e+110, 1.e-110], dtype=np.complex128) + y = x**2/x + assert_almost_equal(y/x, [1, 1], err_msg=msg) + + def test_floor_division_complex(self): + # check that implementation is correct + msg = "Complex floor division implementation check" + x = np.array([.9 + 1j, -.1 + 1j, .9 + .5*1j, .9 + 2.*1j], dtype=np.complex128) + y = np.array([0., -1., 0., 0.], dtype=np.complex128) + assert_equal(np.floor_divide(x**2,x), y, err_msg=msg) + # check overflow, underflow + msg = "Complex floor division overflow/underflow check" + x = np.array([1.e+110, 1.e-110], dtype=np.complex128) + y = np.floor_divide(x**2, x) + assert_equal(y, [1.e+110, 0], err_msg=msg) + class TestPower(TestCase): def test_power_float(self): x = np.array([1., 2., 3.]) @@ -42,7 +65,7 @@ class TestPower(TestCase): def assert_complex_equal(x, y): assert_array_equal(x.real, y.real) assert_array_equal(x.imag, y.imag) - + for z in [complex(0, np.inf), complex(1, np.inf)]: z = np.array([z], dtype=np.complex_) assert_complex_equal(z**1, z) @@ -87,7 +110,25 @@ class TestLogAddExp2(object): logxf = np.array(x, dtype=dt) logyf = np.array(y, dtype=dt) logzf = np.array(z, dtype=dt) - assert_almost_equal(np.logaddexp(logxf, logyf), logzf) + assert_almost_equal(np.logaddexp2(logxf, logyf), logzf) + + def test_inf(self) : + inf = np.inf + x = [inf, -inf, inf, -inf, inf, 1, -inf, 1] + y = [inf, inf, -inf, -inf, 1, inf, 1, -inf] + z = [inf, inf, inf, -inf, inf, inf, 1, 1] + for dt in ['f','d','g'] : + logxf = np.array(x, dtype=dt) + logyf = np.array(y, dtype=dt) + logzf = np.array(z, dtype=dt) + assert_equal(np.logaddexp2(logxf, logyf), logzf) + + def test_nan(self): + assert np.isnan(np.logaddexp2(np.nan, np.inf)) + assert np.isnan(np.logaddexp2(np.inf, np.nan)) + assert np.isnan(np.logaddexp2(np.nan, 0)) + assert np.isnan(np.logaddexp2(0, np.nan)) + assert np.isnan(np.logaddexp2(np.nan, np.nan)) class TestLog(TestCase): def test_log_values(self) : @@ -130,6 +171,24 @@ class TestLogAddExp(object): logzf = np.array(z, dtype=dt) assert_almost_equal(np.logaddexp(logxf, logyf), logzf) + def test_inf(self) : + inf = np.inf + x = [inf, -inf, inf, -inf, inf, 1, -inf, 1] + y = [inf, inf, -inf, -inf, 1, inf, 1, -inf] + z = [inf, inf, inf, -inf, inf, inf, 1, 1] + for dt in ['f','d','g'] : + logxf = np.array(x, dtype=dt) + logyf = np.array(y, dtype=dt) + logzf = np.array(z, dtype=dt) + assert_equal(np.logaddexp(logxf, logyf), logzf) + + def test_nan(self): + assert np.isnan(np.logaddexp(np.nan, np.inf)) + assert np.isnan(np.logaddexp(np.inf, np.nan)) + assert np.isnan(np.logaddexp(np.nan, 0)) + assert np.isnan(np.logaddexp(0, np.nan)) + assert np.isnan(np.logaddexp(np.nan, np.nan)) + class TestLog1p(TestCase): def test_log1p(self): assert_almost_equal(ncu.log1p(0.2), ncu.log(1.2)) @@ -140,6 +199,94 @@ class TestExpm1(TestCase): assert_almost_equal(ncu.expm1(0.2), ncu.exp(0.2)-1) assert_almost_equal(ncu.expm1(1e-6), ncu.exp(1e-6)-1) +class TestHypot(TestCase, object): + def test_simple(self): + assert_almost_equal(ncu.hypot(1, 1), ncu.sqrt(2)) + assert_almost_equal(ncu.hypot(0, 0), 0) + +def assert_hypot_isnan(x, y): + assert np.isnan(ncu.hypot(x, y)) + +def assert_hypot_isinf(x, y): + assert np.isinf(ncu.hypot(x, y)) + +def test_hypot_special_values(): + yield assert_hypot_isnan, np.nan, np.nan + yield assert_hypot_isnan, np.nan, 1 + yield assert_hypot_isinf, np.nan, np.inf + yield assert_hypot_isinf, np.inf, np.nan + yield assert_hypot_isinf, np.inf, 0 + yield assert_hypot_isinf, 0, np.inf + +def test_arctan2_special_values(): + def assert_arctan2_isnan(x, y): + assert np.isnan(ncu.arctan2(x, y)) + + def assert_arctan2_ispinf(x, y): + assert np.isinf(ncu.arctan2(x, y)) and ncu.arctan2(x, y) > 0 + + def assert_arctan2_isninf(x, y): + assert np.isinf(ncu.arctan2(x, y)) and ncu.arctan2(x, y) < 0 + + def assert_arctan2_ispzero(x, y): + assert ncu.arctan2(x, y) == 0 and not np.signbit(ncu.arctan2(x, y)) + + def assert_arctan2_isnzero(x, y): + assert ncu.arctan2(x, y) == 0 and np.signbit(ncu.arctan2(x, y)) + + # atan2(1, 1) returns pi/4. + yield assert_almost_equal, ncu.arctan2(1, 1), 0.25 * np.pi + yield assert_almost_equal, ncu.arctan2(-1, 1), -0.25 * np.pi + yield assert_almost_equal, ncu.arctan2(1, -1), 0.75 * np.pi + + # atan2(+-0, -0) returns +-pi. + yield assert_almost_equal, ncu.arctan2(np.PZERO, np.NZERO), np.pi + yield assert_almost_equal, ncu.arctan2(np.NZERO, np.NZERO), -np.pi + # atan2(+-0, +0) returns +-0. + yield assert_arctan2_ispzero, np.PZERO, np.PZERO + yield assert_arctan2_isnzero, np.NZERO, np.PZERO + + # atan2(+-0, x) returns +-pi for x < 0. + yield assert_almost_equal, ncu.arctan2(np.PZERO, -1), np.pi + yield assert_almost_equal, ncu.arctan2(np.NZERO, -1), -np.pi + + # atan2(+-0, x) returns +-0 for x > 0. + yield assert_arctan2_ispzero, np.PZERO, 1 + yield assert_arctan2_isnzero, np.NZERO, 1 + + # atan2(y, +-0) returns +pi/2 for y > 0. + yield assert_almost_equal, ncu.arctan2(1, np.PZERO), 0.5 * np.pi + yield assert_almost_equal, ncu.arctan2(1, np.NZERO), 0.5 * np.pi + + # atan2(y, +-0) returns -pi/2 for y < 0. + yield assert_almost_equal, ncu.arctan2(-1, np.PZERO), -0.5 * np.pi + yield assert_almost_equal, ncu.arctan2(-1, np.NZERO), -0.5 * np.pi + + # atan2(+-y, -infinity) returns +-pi for finite y > 0. + yield assert_almost_equal, ncu.arctan2(1, np.NINF), np.pi + yield assert_almost_equal, ncu.arctan2(-1, np.NINF), -np.pi + + # atan2(+-y, +infinity) returns +-0 for finite y > 0. + yield assert_arctan2_ispzero, 1, np.inf + yield assert_arctan2_isnzero, -1, np.inf + + # atan2(+-infinity, x) returns +-pi/2 for finite x. + yield assert_almost_equal, ncu.arctan2( np.inf, 1), 0.5 * np.pi + yield assert_almost_equal, ncu.arctan2(-np.inf, 1), -0.5 * np.pi + + # atan2(+-infinity, -infinity) returns +-3*pi/4. + yield assert_almost_equal, ncu.arctan2( np.inf, -np.inf), 0.75 * np.pi + yield assert_almost_equal, ncu.arctan2(-np.inf, -np.inf), -0.75 * np.pi + + # atan2(+-infinity, +infinity) returns +-pi/4. + yield assert_almost_equal, ncu.arctan2( np.inf, np.inf), 0.25 * np.pi + yield assert_almost_equal, ncu.arctan2(-np.inf, np.inf), -0.25 * np.pi + + # atan2(nan, x) returns nan for any x, including inf + yield assert_arctan2_isnan, np.nan, np.inf + yield assert_arctan2_isnan, np.inf, np.nan + yield assert_arctan2_isnan, np.nan, np.nan + class TestMaximum(TestCase): def test_reduce_complex(self): assert_equal(np.maximum.reduce([1,2j]),1) @@ -337,6 +484,38 @@ class TestSpecialMethods(TestCase): a = A() self.failUnlessRaises(RuntimeError, ncu.maximum, a, a) + def test_default_prepare(self): + class with_wrap(object): + __array_priority__ = 10 + def __array__(self): + return np.zeros(1) + def __array_wrap__(self, arr, context): + return arr + a = with_wrap() + x = ncu.minimum(a, a) + assert_equal(x, np.zeros(1)) + assert_equal(type(x), np.ndarray) + + def test_prepare(self): + class with_prepare(np.ndarray): + __array_priority__ = 10 + def __array_prepare__(self, arr, context): + # make sure we can return a new + return np.array(arr).view(type=with_prepare) + a = np.array(1).view(type=with_prepare) + x = np.add(a, a) + assert_equal(x, np.array(2)) + assert_equal(type(x), with_prepare) + + def test_failing_prepare(self): + class A(object): + def __array__(self): + return np.zeros(1) + def __array_prepare__(self, arr, context=None): + raise RuntimeError + a = A() + self.failUnlessRaises(RuntimeError, ncu.maximum, a, a) + def test_array_with_context(self): class A(object): def __array__(self, dtype=None, context=None): @@ -637,6 +816,13 @@ def _check_branch_cut(f, x0, dx, re_sign=1, im_sign=-1, sig_zero_ok=False, assert np.all(np.absolute(y0.real - ym.real*re_sign) < atol), (y0, ym) assert np.all(np.absolute(y0.imag - ym.imag*im_sign) < atol), (y0, ym) +def test_copysign(): + assert np.copysign(1, -1) == -1 + assert 1 / np.copysign(0, -1) < 0 + assert 1 / np.copysign(0, 1) > 0 + assert np.signbit(np.copysign(np.nan, -1)) + assert not np.signbit(np.copysign(np.nan, 1)) + def test_pos_nan(): """Check np.nan is a positive nan.""" assert np.signbit(np.nan) == 0 |
