Skip to content

Commit

Permalink
Merge pull request Theano#381 from lamblin/intdiv_c_code
Browse files Browse the repository at this point in the history 
C code for integer division
  • Loading branch information
@nouiz
nouiz committed Jan 25, 2012
2 parents bc57896 + 96211a0 commit 0a74478
Show file tree
Hide file tree
Showing 2 changed files with 145 additions and 89 deletions.
185 changes: 121 additions & 64 deletions theano/scalar/basic.py
View file
Expand Up @@ -16,6 +16,7 @@
import warnings
from copy import copy
from itertools import imap
from textwrap import dedent

import numpy

Expand Down Expand Up @@ -1305,15 +1306,69 @@ def grad(self, (x, y), (gz, )):


class IntDiv(BinaryScalarOp):
complex_error = ComplexError(
"Theano does not support integer division (//) on "
"complex numbers, since numpy deprecated it.")

def impl(self, x, y):
return x // y

def c_support_code(self):
# We use a macro as python use % as a special string character,
# and the output of c_code may be run through another level
# of string formatting.
return "#define THEANO_MACRO_MOD(x,y) (x % y)"

def c_code(self, node, name, (x, y), (z,), sub):
raise NotImplementedError("For integer arguments the behavior of"
" division in C and in Python [can] differ"
" when the quotient is negative. C actually"
" does not even specify a correct behaviour"
" in this case, it is up to the chip.")
t = node.inputs[0].type.upcast(*[i.type for i in node.inputs[1:]])
if t in imap(str, discrete_types):
x_div_y_pp = '(%(x)s / %(y)s)' % locals()
x_div_y_mp = '((-%(x)s) / %(y)s)' % locals()
x_mod_y_mp = 'THEANO_MACRO_MOD((-%(x)s), %(y)s)' % locals()
x_div_y_pm = '(%(x)s / (-%(y)s))' % locals()
x_mod_y_pm = 'THEANO_MACRO_MOD(%(x)s, (-%(y)s))' % locals()
x_div_y_mm = '((-%(x)s) / (-%(y)s))' % locals()
elif t in imap(str, float_types):
# We need to call different functions of math.h
# depending on the type
if t == 'float32':
floor = 'floorf'
fmod = 'fmodf'
elif t == 'float64':
floor = 'floor'
fmod = 'fmod'
else:
raise NotImplementedError('type not supported', t)

x_div_y_pp = '%(floor)s(%(x)s / %(y)s)' % locals()
x_div_y_mp = '%(floor)s((-%(x)s) / %(y)s)' % locals()
x_mod_y_mp = '%(fmod)s((-%(x)s), %(y)s)' % locals()
x_div_y_pm = '%(floor)s(%(x)s / (-%(y)s))' % locals()
x_mod_y_pm = '%(fmod)s(%(x)s, (-%(y)s))' % locals()
x_div_y_mm = '%(floor)s((-%(x)s) / (-%(y)s))' % locals()
elif t in complex_types:
raise self.complex_error
else:
raise NotImplementedError('type not supported', t)

return dedent("""
if (%(x)s < 0) {
if (%(y)s < 0) {
%(z)s = %(x_div_y_mm)s;
} else {
%(z)s = - %(x_div_y_mp)s - ((%(x_mod_y_mp)s == 0) ? 0 : 1);
}
} else {
if (%(y)s < 0) {
%(z)s = - %(x_div_y_pm)s - ((%(x_mod_y_pm)s == 0) ? 0 : 1);
} else {
%(z)s = %(x_div_y_pp)s;
}
}
""") % locals()

def c_code_cache_version(self):
return (2,)

def grad(self, inputs, g_output):
return [None] * len(inputs)
Expand Down Expand Up @@ -1346,7 +1401,9 @@ def c_code_cache_version(self):
return (5,)

def c_support_code(self):
#We use a macro as python use % as a special string caractere.
# We use a macro as python use % as a special string character,
# and the output of c_code may be run through another level
# of string formatting.
return "#define THEANO_MACRO_MOD(x,y) (x % y)"

def c_code(self, node, name, (x, y), (z, ), sub):
Expand Down Expand Up @@ -1383,21 +1440,21 @@ def c_code(self, node, name, (x, y), (z, ), sub):
elif str(t) in imap(str, complex_types):
raise self.complex_error
else:
raise NotImplementedError('type not supported', type)

return """
if (%(x)s < 0){
if (%(y)s < 0){
%(z)s = -(%(x_mod_ymm)s);
}else{
%(z)s = - %(x_mod_ymp)s + (%(x_mod_ymp)s != 0 ? %(y)s : 0);
}
}else if (%(y)s < 0){
%(z)s = (%(x_mod_ypm)s) + (%(x_mod_ypm)s != 0 ? %(y)s : 0);
}else{
%(z)s = %(x_mod_y)s;
}
""" % locals()
raise NotImplementedError('type not supported', t)

return dedent("""
if (%(x)s < 0){
if (%(y)s < 0){
%(z)s = -(%(x_mod_ymm)s);
}else{
%(z)s = - %(x_mod_ymp)s + (%(x_mod_ymp)s != 0 ? %(y)s : 0);
}
}else if (%(y)s < 0){
%(z)s = (%(x_mod_ypm)s) + (%(x_mod_ypm)s != 0 ? %(y)s : 0);
}else{
%(z)s = %(x_mod_y)s;
}
""") % locals()

def grad(self, (x, y), (gz, )):
return None, None
Expand Down Expand Up @@ -1660,51 +1717,51 @@ def c_code___(self, node, name, (x, ), (z, ), sub):
if not node.outputs[0].type.dtype in ['float32', 'float64']:
Exception("The output should be float32 or float64")

return """
#ifndef ROUNDING_EPSILON
#define ROUNDING_EPSILON 0.0000001
#endif
if (%(x)s < 0.0){
// We implement the else part like that: -else( -%(x)s);
%(typ)s i;
std::modf( -%(x)s, &i );
// If %(x)s is exactly halfway between two integers
if ((-%(x)s -(i +0.5)) < epsilon){
// If 'i' is even then return 'i'
if (std::fmod( i, 2.0 ) < epsilon){
%(z)s = - i;
}else{
// Else return the nearest even integer
%(z)s = - ceil( i +0.5 );
}
}else{
// round to closest
%(z)s = - round(%(x)s+5);
}
}else{
%(typ)s i;
std::modf( %(x)s, &i );
// If %(x)s is exactly halfway between two integers
if ((%(x)s -(i +0.5)) < epsilon){
// If 'i' is even then return 'i'
if (std::fmod( i, 2.0 ) < epsilon){
%(z)s = i;
}else{
// Else return the nearest even integer
%(z)s = ceil( i +0.5 );
}
}else{
// round to closest
%(z)s = round(%(x)s+5);
}
}
return dedent("""
#ifndef ROUNDING_EPSILON
#define ROUNDING_EPSILON 0.0000001
#endif
if (%(x)s < 0.0){
// We implement the else part like that: -else( -%(x)s);
%(typ)s i;
std::modf( -%(x)s, &i );
// If %(x)s is exactly halfway between two integers
if ((-%(x)s -(i +0.5)) < epsilon){
// If 'i' is even then return 'i'
if (std::fmod( i, 2.0 ) < epsilon){
%(z)s = - i;
}else{
// Else return the nearest even integer
%(z)s = - ceil( i +0.5 );
}
}else{
// round to closest
%(z)s = - round(%(x)s+5);
}
}else{
%(typ)s i;
std::modf( %(x)s, &i );
// If %(x)s is exactly halfway between two integers
if ((%(x)s -(i +0.5)) < epsilon){
// If 'i' is even then return 'i'
if (std::fmod( i, 2.0 ) < epsilon){
%(z)s = i;
}else{
// Else return the nearest even integer
%(z)s = ceil( i +0.5 );
}
}else{
// round to closest
%(z)s = round(%(x)s+5);
}
}
#undef ROUNDING_EPSILON
#undef ROUNDING_EPSILON
"""
""")
round_half_to_even = RoundHalfToEven(same_out_float_only)


Expand Down
49 changes: 24 additions & 25 deletions theano/tensor/tests/test_basic.py
View file
Expand Up @@ -567,11 +567,9 @@ def copymod(dct, without=[], **kwargs):
column=(rand(2, 3), rand(2, 1)),
dtype_mixup_1=(rand(2, 3), randint_nonzero(2, 3)),
dtype_mixup_2=(randint_nonzero(2, 3), rand(2, 3)),
# Fix problem with integers and uintegers and add them.
# Then remove their specific addition to CeilIntDivTester tests.
# integer=(randint(2, 3), randint_nonzero(2, 3)),
# uinteger=(randint(2, 3).astype("uint8"),
# randint_nonzero(2, 3).astype("uint8")),
integer=(randint(2, 3), randint_nonzero(2, 3)),
uinteger=(randint(2, 3).astype("uint8"),
randint_nonzero(2, 3).astype("uint8")),
# This empty2 doesn't work for some tests. I don't remember why
#empty2=(numpy.asarray([0]), numpy.asarray([])),
)
Expand Down Expand Up @@ -610,31 +608,32 @@ def copymod(dct, without=[], **kwargs):
# float32.
# This is probably caused by our way of computing the gradient error.
div_grad_rtol=0.025
TrueDivTester = makeBroadcastTester(op = tensor.true_div,
expected = lambda x, y: check_floatX((x, y), x / y),
good = _good_broadcast_div_mod_normal_float,
# integers = (randint(2, 3), randint_nonzero(2, 3)),
# dtype_mixup_1 = (rand(2, 3), randint_nonzero(2, 3)),
# dtype_mixup_2 = (randint_nonzero(2, 3), rand(2, 3))),
grad = _grad_broadcast_div_mod_normal,
grad_rtol=div_grad_rtol,
)
TrueDivInplaceTester = makeBroadcastTester(op = inplace.true_div_inplace,
expected = lambda x, y: x / y,
good = _good_broadcast_div_mod_normal_float_inplace,
grad = _grad_broadcast_div_mod_normal,
grad_rtol=div_grad_rtol,
inplace = True)

TrueDivTester = makeBroadcastTester(
op=tensor.true_div,
expected = (lambda x, y:
check_floatX((x, y), numpy.true_divide(x, y))),
good=_good_broadcast_div_mod_normal_float,
grad=_grad_broadcast_div_mod_normal,
grad_rtol=div_grad_rtol,
)

TrueDivInplaceTester = makeBroadcastTester(
op=inplace.true_div_inplace,
expected=(lambda x, y: numpy.true_divide(x, y)),
good=copymod(
_good_broadcast_div_mod_normal_float_inplace,
# The output is now in float, we cannot work inplace on an int.
without=['integer', 'uinteger']),
grad=_grad_broadcast_div_mod_normal,
grad_rtol=div_grad_rtol,
inplace=True)


CeilIntDivTester = makeBroadcastTester(
op=tensor.ceil_intdiv,
expected=lambda x, y: check_floatX((x, y), (x // y) + ((x % y) != 0)),
good=copymod(_good_broadcast_div_mod_normal_float_no_complex,
integer=(randint(2, 3), randint_nonzero(2, 3)),
uinteger=(randint(2, 3).astype("uint8"),
randint_nonzero(2, 3).astype("uint8")),
),
good=_good_broadcast_div_mod_normal_float_no_complex,
name='CeilIntDiv',
# As we implement this function with neq, the gradient returned is always 0.
# grad=_grad_broadcast_div_mod_normal,
Expand Down

0 comments on commit 0a74478

Please sign in to comment.