Skip to content

Commit

Permalink
Fix broadcasting & add omitted tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@markflorisson
markflorisson committed Jul 23, 2012
1 parent ce10532 commit 0653f04
Show file tree
Hide file tree
Showing 3 changed files with 305 additions and 7 deletions.
18 changes: 11 additions & 7 deletions Cython/Compiler/Vector.py
View file
Expand Up @@ -808,25 +808,29 @@ 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)

for i in range(lhs_offset):
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)" % (
Expand Down
218 changes: 218 additions & 0 deletions tests/array_expressions/broadcasting.pyx
View file
@@ -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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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.]],
<BLANKLINE>
[[ 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)
76 changes: 76 additions & 0 deletions tests/array_expressions/utils.pxi
View file
@@ -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

0 comments on commit 0653f04

Please sign in to comment.