From 0653f0487332beae5ecd0b05dccad37b62b8ebc8 Mon Sep 17 00:00:00 2001 From: Mark Florisson Date: Sun, 24 Jun 2012 13:08:57 +0100 Subject: [PATCH] Fix broadcasting & add omitted tests --- Cython/Compiler/Vector.py | 18 +- tests/array_expressions/broadcasting.pyx | 218 +++++++++++++++++++++++ tests/array_expressions/utils.pxi | 76 ++++++++ 3 files changed, 305 insertions(+), 7 deletions(-) create mode 100644 tests/array_expressions/broadcasting.pyx create mode 100644 tests/array_expressions/utils.pxi diff --git a/Cython/Compiler/Vector.py b/Cython/Compiler/Vector.py index f89826f1929..1bc053880ea 100644 --- a/Cython/Compiler/Vector.py +++ b/Cython/Compiler/Vector.py @@ -808,13 +808,16 @@ def advance_lhs_data_ptr(self, code): m3[:, :] = m1[:] * m2[:] m3[0, :] contains the data, which we need to broadcast over m3[1:, :] + + (This only works for one dimension at a time, so this would need a + wrapping loop. Currently disabled, it broadcasts itself with itself). """ lhs_offset, _ = offsets(self.lhs, self.rhs) lhs_r, rhs_r = self.lhs.result(), self.rhs.result() def advance(i): - code.putln("%s.data += %s.strides[%d];" % (lhs_r, lhs_r, i)) - code.putln("%s.shape[%d] -= 1;" % (lhs_r, i)) + #code.putln("%s.data += %s.strides[%d];" % (lhs_r, lhs_r, i)) + #code.putln("%s.shape[%d] -= 1;" % (lhs_r, i)) if not lhs_offset: self.final_broadcast.init_broadcast_flag(code, True) @@ -822,11 +825,12 @@ def advance(i): advance(i) self.final_broadcast.init_broadcast_flag(code, True) - for i in range(self.rhs.type.ndim): - code.putln("if (%s.shape[%d] > 1 && %s.shape[%d] == 1) {" % - (lhs_r, i + lhs_offset, rhs_r, i)) - advance(i + lhs_offset) - code.putln("}") + if not lhs_offset: + for i in range(self.rhs.type.ndim): + code.putln("if (%s.shape[%d] > 1 && %s.shape[%d] == 1) {" % + (lhs_r, i + lhs_offset, rhs_r, i)) + advance(i + lhs_offset) + code.putln("}") def verify_final_shape(self, code): call = "__pyx_verify_shapes(%s, %s, %d, %d)" % ( diff --git a/tests/array_expressions/broadcasting.pyx b/tests/array_expressions/broadcasting.pyx new file mode 100644 index 00000000000..98d7b86ebf6 --- /dev/null +++ b/tests/array_expressions/broadcasting.pyx @@ -0,0 +1,218 @@ +# tag: numpy +# mode: run + +include "utils.pxi" + +#@testcase +def test_broadcasting(fused_dtype_t[:] m1, fused_dtype_t[:, :] m2): + """ + >>> test_broadcasting(*operands('l')) + array([[ 0, 2, 4, 6, 8, 10, 12, 14, 16, 18], + [ 10, 12, 14, 16, 18, 20, 22, 24, 26, 28], + [ 20, 22, 24, 26, 28, 30, 32, 34, 36, 38], + [ 30, 32, 34, 36, 38, 40, 42, 44, 46, 48], + [ 40, 42, 44, 46, 48, 50, 52, 54, 56, 58], + [ 50, 52, 54, 56, 58, 60, 62, 64, 66, 68], + [ 60, 62, 64, 66, 68, 70, 72, 74, 76, 78], + [ 70, 72, 74, 76, 78, 80, 82, 84, 86, 88], + [ 80, 82, 84, 86, 88, 90, 92, 94, 96, 98], + [ 90, 92, 94, 96, 98, 100, 102, 104, 106, 108]]) + + >>> test_broadcasting(*operands(np.complex128)) + array([[ 0.+0.j, 2.+0.j, 4.+0.j, 6.+0.j, 8.+0.j, 10.+0.j, + 12.+0.j, 14.+0.j, 16.+0.j, 18.+0.j], + [ 10.+0.j, 12.+0.j, 14.+0.j, 16.+0.j, 18.+0.j, 20.+0.j, + 22.+0.j, 24.+0.j, 26.+0.j, 28.+0.j], + [ 20.+0.j, 22.+0.j, 24.+0.j, 26.+0.j, 28.+0.j, 30.+0.j, + 32.+0.j, 34.+0.j, 36.+0.j, 38.+0.j], + [ 30.+0.j, 32.+0.j, 34.+0.j, 36.+0.j, 38.+0.j, 40.+0.j, + 42.+0.j, 44.+0.j, 46.+0.j, 48.+0.j], + [ 40.+0.j, 42.+0.j, 44.+0.j, 46.+0.j, 48.+0.j, 50.+0.j, + 52.+0.j, 54.+0.j, 56.+0.j, 58.+0.j], + [ 50.+0.j, 52.+0.j, 54.+0.j, 56.+0.j, 58.+0.j, 60.+0.j, + 62.+0.j, 64.+0.j, 66.+0.j, 68.+0.j], + [ 60.+0.j, 62.+0.j, 64.+0.j, 66.+0.j, 68.+0.j, 70.+0.j, + 72.+0.j, 74.+0.j, 76.+0.j, 78.+0.j], + [ 70.+0.j, 72.+0.j, 74.+0.j, 76.+0.j, 78.+0.j, 80.+0.j, + 82.+0.j, 84.+0.j, 86.+0.j, 88.+0.j], + [ 80.+0.j, 82.+0.j, 84.+0.j, 86.+0.j, 88.+0.j, 90.+0.j, + 92.+0.j, 94.+0.j, 96.+0.j, 98.+0.j], + [ 90.+0.j, 92.+0.j, 94.+0.j, 96.+0.j, 98.+0.j, 100.+0.j, + 102.+0.j, 104.+0.j, 106.+0.j, 108.+0.j]]) + + >>> result1 = test_broadcasting(object_range(10), object_range(100, (10, 10))) + >>> result1 + array([[0, 2, 4, 6, 8, 10, 12, 14, 16, 18], + [10, 12, 14, 16, 18, 20, 22, 24, 26, 28], + [20, 22, 24, 26, 28, 30, 32, 34, 36, 38], + [30, 32, 34, 36, 38, 40, 42, 44, 46, 48], + [40, 42, 44, 46, 48, 50, 52, 54, 56, 58], + [50, 52, 54, 56, 58, 60, 62, 64, 66, 68], + [60, 62, 64, 66, 68, 70, 72, 74, 76, 78], + [70, 72, 74, 76, 78, 80, 82, 84, 86, 88], + [80, 82, 84, 86, 88, 90, 92, 94, 96, 98], + [90, 92, 94, 96, 98, 100, 102, 104, 106, 108]], dtype=object) + + >>> result2 = test_broadcasting(object_range(10, cls=UniqueObjectInplace), + ... object_range(100, (10, 10), cls=UniqueObjectInplace)) + >>> np.all(result1 == result2) + True + """ + m2[...] = m2 + m1 + return np.asarray(m2) + +#@testcase_like(test_broadcasting) +def test_broadcasting_c_contig(fused_dtype_t[::1] m1, fused_dtype_t[:, ::1] m2): + m2[...] = m2 + m1 + return np.asarray(m2) + +@testcase +def test_broadcasting_larger_rhs_ndim(double[:] m): + """ + >>> test_broadcasting_larger_rhs_ndim(np.arange(10, dtype=np.double)) + """ + m[:] = m[None, ::-1] + m[None, ::-1] + +@testcase +def test_broadcasting_larger_rhs_ndim_contig(double[::1] m): + """ + >>> test_broadcasting_larger_rhs_ndim_contig(np.arange(10, dtype=np.double)) + """ + m[:] = m[None, ::-1] + m[None, ::-1] + +@testcase +def test_broadcasting_runtime(double[:, :] m1, double[:, :] m2): + """ + >>> m1, m2 = operands(np.double) + >>> test_broadcasting_runtime(m1[None, :], m2[:2]) + array([[ 0., 2., 4., 6., 8., 10., 12., 14., 16., 18.], + [ 0., 2., 4., 6., 8., 10., 12., 14., 16., 18.]]) + >>> test_broadcasting_runtime(m1[:, None], m2[:, :2]) + array([[ 0., 0.], + [ 2., 2.], + [ 4., 4.], + [ 6., 6.], + [ 8., 8.], + [ 10., 10.], + [ 12., 12.], + [ 14., 14.], + [ 16., 16.], + [ 18., 18.]]) + >>> test_broadcasting_runtime(m1[:, None], m2[:, :4]) + array([[ 0., 0., 0., 0.], + [ 2., 2., 2., 2.], + [ 4., 4., 4., 4.], + [ 6., 6., 6., 6.], + [ 8., 8., 8., 8.], + [ 10., 10., 10., 10.], + [ 12., 12., 12., 12.], + [ 14., 14., 14., 14.], + [ 16., 16., 16., 16.], + [ 18., 18., 18., 18.]]) + >>> np.all(test_broadcasting_runtime(np.array([[1.0]], dtype=np.double), m2) == 2.0) + True + >>> np.all(test_broadcasting_runtime(np.array([1.0], dtype=np.double)[:, None], m2) == 2.0) + True + """ + m2[:, :] = m1 + m1 + return np.asarray(m2) + +@testcase +def test_broadcasting_multiple_dims(double[:, :, :] m1, double[:, :, :] m2): + """ + >>> m2 = np.empty((5, 5, 5), dtype=np.double) + >>> m1 = np.arange(5, dtype=np.double) + + >>> test_broadcasting_multiple_dims(m1[None, None, :], m2.copy()) + array([[[ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.]], + + [[ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.]], + + [[ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.]], + + [[ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.]], + + [[ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.], + [ 0., 2., 4., 6., 8.]]]) + >>> test_broadcasting_multiple_dims(m1[None, :, None], m2.copy()) + array([[[ 0., 0., 0., 0., 0.], + [ 2., 2., 2., 2., 2.], + [ 4., 4., 4., 4., 4.], + [ 6., 6., 6., 6., 6.], + [ 8., 8., 8., 8., 8.]], + + [[ 0., 0., 0., 0., 0.], + [ 2., 2., 2., 2., 2.], + [ 4., 4., 4., 4., 4.], + [ 6., 6., 6., 6., 6.], + [ 8., 8., 8., 8., 8.]], + + [[ 0., 0., 0., 0., 0.], + [ 2., 2., 2., 2., 2.], + [ 4., 4., 4., 4., 4.], + [ 6., 6., 6., 6., 6.], + [ 8., 8., 8., 8., 8.]], + + [[ 0., 0., 0., 0., 0.], + [ 2., 2., 2., 2., 2.], + [ 4., 4., 4., 4., 4.], + [ 6., 6., 6., 6., 6.], + [ 8., 8., 8., 8., 8.]], + + [[ 0., 0., 0., 0., 0.], + [ 2., 2., 2., 2., 2.], + [ 4., 4., 4., 4., 4.], + [ 6., 6., 6., 6., 6.], + [ 8., 8., 8., 8., 8.]]]) + >>> test_broadcasting_multiple_dims(m1[:, None, None], m2.copy()) + array([[[ 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0.]], + + [[ 2., 2., 2., 2., 2.], + [ 2., 2., 2., 2., 2.], + [ 2., 2., 2., 2., 2.], + [ 2., 2., 2., 2., 2.], + [ 2., 2., 2., 2., 2.]], + + [[ 4., 4., 4., 4., 4.], + [ 4., 4., 4., 4., 4.], + [ 4., 4., 4., 4., 4.], + [ 4., 4., 4., 4., 4.], + [ 4., 4., 4., 4., 4.]], + + [[ 6., 6., 6., 6., 6.], + [ 6., 6., 6., 6., 6.], + [ 6., 6., 6., 6., 6.], + [ 6., 6., 6., 6., 6.], + [ 6., 6., 6., 6., 6.]], + + [[ 8., 8., 8., 8., 8.], + [ 8., 8., 8., 8., 8.], + [ 8., 8., 8., 8., 8.], + [ 8., 8., 8., 8., 8.], + [ 8., 8., 8., 8., 8.]]]) + """ + m2[:] = m1 + m1 + return np.asarray(m2) diff --git a/tests/array_expressions/utils.pxi b/tests/array_expressions/utils.pxi new file mode 100644 index 00000000000..068a66a13ad --- /dev/null +++ b/tests/array_expressions/utils.pxi @@ -0,0 +1,76 @@ +cimport cython +from cython cimport view + +from cpython.object cimport Py_EQ +cimport numpy as np + +import sys +import numpy as np + + +__test__ = {} + +def testcase(f): + assert f.__doc__, f + __test__[f.__name__] = f.__doc__ + return f + +def testcase_like(similar_func, docstring_specializer=None): + def decorator(wrapper_func): + assert similar_func.__doc__ + doc = similar_func.__doc__.replace(similar_func.__name__, + wrapper_func.__name__) + if docstring_specializer: + doc = docstring_specializer(doc) + + wrapper_func.__doc__ = doc + return testcase(wrapper_func) + return decorator + +def testcase_like_sub(similar_func, sub_dict): + def docstring_specializer(s): + # not the most efficient implementation :) + for k, v in sub_dict.items(): + s = s.replace(k, v) + return s + return testcase_like(similar_func, docstring_specializer) + +cdef class UniqueObject(object): + cdef public object value + def __init__(self, value): + self.value = value + + def __add__(self, other): + return UniqueObject(self.value + other.value) + + def __richcmp__(self, other, int opid): + if opid == Py_EQ: + return self.value == other.value + return NotImplemented + + def __str__(self): + return str(self.value) + + def __dealloc__(self): + pass + # sys.stdout.write("dealloc %s " % self.value) + +class UniqueObjectInplace(UniqueObject): + def __add__(self, other): + self.value += other.value + return self + +def object_range(n, shape=None, cls=UniqueObject): + result = np.array([cls(i) for i in range(n)], dtype=np.object) + if shape: + result = result.reshape(shape) + return result + +def operands(dtype): + return np.arange(10, dtype=dtype), np.arange(100, dtype=dtype).reshape(10, 10) + +cdef fused fused_dtype_t: + long + long double + double complex + object