Added new behaviors for data type casts (enabled through config.cast_policy)

doc/library/config.txt

@@ -144,7 +144,7 @@ import theano and print the config variable, as in:
 
 .. attribute:: floatX
 
-    String value: either 'float64' or 'float32'.
+    String value: either 'float64' or 'float32'
 
     Default: 'float64'
 
@@ -152,6 +152,23 @@ import theano and print the config variable, as in:
     and similar functions.  It also sets the default theano bit width for
     arguments passed as Python floating-point numbers.
 
+.. attribute:: cast_policy
+
+    String value: either 'numpy+floatX', 'numpy' or 'custom'
+
+    Default: 'custom'
+
+    This specifies how data types are implicitly figured out in Theano, e.g.  for
+    constants or in the result of arithmetic operations. The recommended value is
+    'numpy+floatX', that mimics numpy's behavior except for floats when
+    ``config.floatX`` is set to 'float32', for which we use float32 instead of
+    float64 unless the user is explicitly using data typed as float64. When
+    'numpy' is used, this specific floatX behavior is discarded. The current
+    default value is 'custom' for backward compatibility reason, and corresponds
+    to a set of custom rules originally used in Theano (which can be partially
+    customized, see e.g. the in-code help of ``tensor.NumpyAutocaster``). The
+    'custom' option will be deprecated in a future release of Theano.
+
 .. attribute:: mode
 
     String value: 'Mode', 'ProfileMode', 'DebugMode', 'FAST_RUN', 'FAST_COMPILE'

theano/configdefaults.py

@@ -15,11 +15,10 @@ AddConfigVar('floatX',
         EnumStr('float64', 'float32'),
         )
 
-# TODO Work-in-progress
-#AddConfigVar('casting_policy',
-#        "Rules for implicit casts of constants in arithmetic operations",
-#        EnumStr('theano_0.3', 'numpy'),
-#        )
+AddConfigVar('cast_policy',
+        "Rules for implicit type casting.",
+        EnumStr('custom', 'numpy+floatX', 'numpy'),
+        )
 
 #gpu mean let the driver select the gpu. Needed in case of gpu in exclusive mode.
 #gpuX mean use the gpu number X.

theano/scalar/basic.py

@@ -26,11 +26,28 @@ builtin_complex = complex
 builtin_int = int
 builtin_float = float
 
+
 def upcast(dtype, *dtypes):
-    z = numpy.zeros((), dtype = dtype)
-    for dtype in dtypes:
-        z = z + numpy.zeros((), dtype = dtype)
-    return str(z.dtype)
+    # Should we try to keep float32 instead of float64? This is used so that
+    # for instance mixing int64 with float32 yields float32 instead of float64.
+    # Note that we store this boolean as a one-element list so that it can be
+    # modified within `make_array`.
+    keep_float32 = [(config.cast_policy == 'numpy+floatX' and
+                     config.floatX == 'float32')]
+    def make_array(dt):
+        if dt == 'float64':
+            # There is an explicit float64 dtype: we cannot keep float32.
+            keep_float32[0] = False
+        return numpy.zeros((), dtype=dt)
+    z = make_array(dtype)
+    for dt in dtypes:
+        z = z + make_array(dt=dt)
+    rval = str(z.dtype)
+    if rval == 'float64' and keep_float32[0]:
+        return 'float32'
+    else:
+        return rval
+
 
 def as_scalar(x, name = None):
     if isinstance(x, gof.Apply):
@@ -47,6 +64,7 @@ def as_scalar(x, name = None):
     except TypeError:
         raise TypeError("Cannot convert %s to Scalar" % x, type(x))
 
+
 def constant(x):
     # pass through numpy scalars, since they are already typed on purpose typically.
     if hasattr(x,'dtype'):
@@ -383,6 +401,7 @@ uint_types = uint8, uint16, uint32, uint64
 float_types = float32, float64
 complex_types = complex64, complex128
 
+discrete_types = int_types + uint_types
 continuous_types = float_types + complex_types
 
 class _scalar_py_operators:
@@ -416,6 +435,7 @@ class _scalar_py_operators:
     def __sub__(self,other): return sub(self,other)
     def __mul__(self,other): return mul(self,other)
     def __div__(self,other): return div_proxy(self,other)
+    def __floordiv__(self,other): return int_div(self,other)
     def __mod__(self,other): return mod(self,other)
     def __pow__(self,other): return pow(self,other)
 
@@ -995,32 +1015,48 @@ class Sub(BinaryScalarOp):
         return first_part, second_part
 sub = Sub(upcast_out, name = 'sub')
 
+
 def div_proxy(x, y):
-    """Proxy for either true_div or int_div, depending on types of x, y.
     """
-    if as_scalar(x).type.dtype.startswith('int') and as_scalar(y).type.dtype.startswith('int'):
-        return int_div(x, y)
+    Currently used as a check to ensure we are not trying to divide integers.
+
+    In 0.4 we will get rid of this function to always use true_div:
+        http://trac-hg.assembla.com/theano/ticket/669
+    """
+    if (as_scalar(x).type in discrete_types and
+        as_scalar(y).type in discrete_types):
+        # Following discussion on theano-dev ("Inconsistent behavior in integer
+        # division"), we will change the semantics of "/" on integer types in
+        # Theano 0.4. Until then, it is forbidden to use "/" on integers.
+        raise NotImplementedError(
+                "Dividing two integers with '/' is forbidden until Theano v0.4"
+                " is released (where the result will be a floating point "
+                "number). In the meantime, please either use '//' for integer "
+                "division, or cast one of the arguments to a floating point "
+                "type for float division.")
     else:
         return true_div(x, y)
 
+
 class TrueDiv(BinaryScalarOp):
     def output_types(self, types):
-        if all(t.dtype.startswith('int') for t in types):
-            return [float64]
+        if all(t in discrete_types for t in types):
+            return [Scalar(config.floatX)]
         else:
             return super(TrueDiv, self).output_types(types)
     def impl(self, x, y):
         x = numpy.asarray(x)
         y = numpy.asarray(y)
-        if str(x.dtype).startswith('int') and str(y.dtype).startswith('int'):
-            return float(x) / y
+        if all(a.dtype in discrete_types for a in (x, y)):
+            return numpy.array(float(x) / y, dtype=config.floatX)
         else:
             return x / y
     def c_code(self, node, name, (x, y), (z, ), sub):
         #we generate good c code only when both are complex!
         if sum([node.inputs[0].type in complex_types, node.inputs[1].type in complex_types])==1:
             raise NotImplementedError('type not supported', type)
-        if node.inputs[0].type in int_types and node.inputs[1].type in int_types:
+        if (node.inputs[0].type in discrete_types and
+            node.inputs[1].type in discrete_types):
             return "%(z)s = ((double)%(x)s) / %(y)s;" % locals()
         return "%(z)s = %(x)s / %(y)s;" % locals()
     def grad(self, (x, y), (gz, )):
@@ -1029,11 +1065,15 @@ class TrueDiv(BinaryScalarOp):
         if x.type in float_types:
             first_part = cast(gz / y, x.type.dtype)
         else:
+            assert x.type in discrete_types
             first_part = None
 
+        if y.type in complex_types:
+            raise NotImplementedError()
         if y.type in float_types:
             second_part = cast(-(gz * x) / (y * y), y.type.dtype)
         else:
+            assert y.type in discrete_types
             second_part = None
         return first_part, second_part
 true_div = TrueDiv(upcast_out, name = 'true_div')

theano/tensor/elemwise.py

@@ -454,7 +454,49 @@ class Elemwise(Op):
         """
 
         inputs = map(as_tensor_variable, inputs)
-        shadow = self.scalar_op.make_node(*[Scalar(dtype = t.type.dtype)() for t in inputs])
+        input_dtypes = [i.dtype for i in inputs]
+
+        scalar_inputs = []
+        array_inputs = []
+        for input_idx, input in enumerate(inputs):
+            if input.ndim == 0:
+                scalar_inputs.append((input_idx, input))
+            else:
+                array_inputs.append((input_idx, input))
+
+        if (scalar_inputs and
+            array_inputs and
+            theano.config.cast_policy in ('numpy', 'numpy+floatX')):
+            # We need to make sure that scalars do not upcast arrays unless
+            # they are fundamentally different. This is specified in
+            #   http://docs.scipy.org/doc/numpy/reference/ufuncs.html
+            # in the 'casting rules' section.
+            array_dtype = scalar.upcast(*[a[1].dtype for a in array_inputs])
+            for input_idx, input in scalar_inputs:
+                # Replace this scalar input's type with the one that numpy
+                # would use when adding this scalar to the array.
+                # Note that currently numpy's behavior is not consistent, which
+                # is a bug (will be fixed in 1.6). See for details
+                #   http://projects.scipy.org/numpy/ticket/1827
+                # As a result, we pick the highest precision data type that
+                # numpy may decide to use (although we may prefer float32 over
+                # float64).
+                n_inputs = [
+                        numpy.array(0, dtype=input_dtypes[input_idx]),
+                        numpy.array([0], dtype=array_dtype)]
+                n_types = [(n_inputs[0] + n_inputs[1]).dtype,
+                           (n_inputs[1] + n_inputs[0]).dtype]
+                n_highest_type = scalar.upcast(*map(str, n_types))
+                if (n_highest_type == 'float64' and
+                    theano.config.cast_policy == 'numpy+floatX' and
+                    theano.config.floatX == 'float32' and
+                    array_dtype != 'float64' and
+                    input_dtypes[input_idx] != 'float64'):
+                    # Prefer float 32 instead.
+                    n_highest_type = 'float32'
+                input_dtypes[input_idx] = n_highest_type
+
+        shadow = self.scalar_op.make_node(*[Scalar(dtype=dtype)() for dtype in input_dtypes])
 
         target_length = max([input.type.ndim for input in inputs])