many changes motivated by theano_cuda_ndarray

theano/compile/debugmode.py

 """Provides `DebugMode`, an evaluation mode for debugging theano internals."""
 __docformat__ = "restructuredtext en"
 
-import time, copy, sys, copy_reg
+import time, copy, sys, copy_reg, gc
 from StringIO import StringIO
 
 import numpy
@@ -20,6 +20,22 @@ from ..compile.function_module import (FunctionMaker,
         Supervisor)
 from ..compile.mode import Mode, register_mode
 
+import logging
+_logger=logging.getLogger("theano.compile.debugmode")
+_logger.setLevel(logging.DEBUG)
+
+def error(*args):
+    #sys.stderr.write('ERROR:'+ ' '.join(str(a) for a in args)+'\n')
+    _logger.error("ERROR: "+' '.join(str(a) for a in args))
+def warning(*args):
+    #sys.stderr.write('WARNING:'+ ' '.join(str(a) for a in args)+'\n')
+    _logger.warning("WARNING: "+' '.join(str(a) for a in args))
+def info(*args):
+    #sys.stderr.write('INFO:'+ ' '.join(str(a) for a in args)+'\n')
+    _logger.info("INFO: "+' '.join(str(a) for a in args))
+def debug(*args):
+    #sys.stderr.write('DEBUG:'+ ' '.join(str(a) for a in args)+'\n')
+    _logger.debug("DEBUG: "+' '.join(str(a) for a in args))
 
 ########################
 #
@@ -116,8 +132,16 @@ class BadOptimization(DebugModeError):
         print >> sio, "  Variable: id", id(self.new_r), self.new_r 
         print >> sio, "  Op", self.new_r.owner
         print >> sio, "  Value Type:", type(self.new_r_val)
-        print >> sio, "  Old Value: ", self.old_r_val
-        print >> sio, "  New Value: ", self.new_r_val
+        str_old_r_val = str(self.old_r_val)
+        if len(str_old_r_val) > 80:
+            print >> sio, "  Old Value: ", str(self.old_r_val)[:80], '...'
+        else:
+            print >> sio, "  Old Value: ", str(self.old_r_val)
+        str_new_r_val = str(self.new_r_val)
+        if len(str_new_r_val) > 80:
+            print >> sio, "  New Value: ", str(self.new_r_val)[:80], '...'
+        else:
+            print >> sio, "  New Value: ", str(self.new_r_val)
         print >> sio, "  Reason: ", str(self.reason)
         print >> sio, "  Old Graph:"
         print >> sio, self.old_graph
@@ -135,6 +159,9 @@ class BadDestroyMap(DebugModeError):
         self.new_val = new_val
     
     def __str__(self):
+        npy_old_val = numpy.asarray(self.old_val)
+        npy_new_val = numpy.asarray(self.new_val)
+        try:
             sio = StringIO()
             print >> sio, "  node:", self.node
             print >> sio, "  node.inputs:", [(str(i), id(i)) for i in self.node.inputs]
@@ -143,15 +170,17 @@ class BadDestroyMap(DebugModeError):
             print >> sio, "  changed input type:", self.node.inputs[self.idx].type
             print >> sio, "  repr (old val):", repr(self.old_val)
             print >> sio, "  repr (new val):", repr(self.new_val)
-        print >> sio, "  value dtype (new <space> old):", self.new_val.dtype, self.old_val.dtype
-        print >> sio, "  value shape (new <space> old):", self.new_val.shape, self.old_val.shape
-        print >> sio, "  value min (new <space> old):", self.new_val.min(), self.old_val.min()
-        print >> sio, "  value max (new <space> old):", self.new_val.max(), self.old_val.max()
-        print >> sio, "  value min (new-old):", (self.new_val-self.old_val).min()
-        print >> sio, "  value max (new-old):", (self.new_val-self.old_val).max()
+            print >> sio, "  value dtype (new <space> old):", npy_new_val.dtype, npy_old_val.dtype
+            print >> sio, "  value shape (new <space> old):", npy_new_val.shape, npy_old_val.shape
+            print >> sio, "  value min (new <space> old):", npy_new_val.min(), npy_old_val.min()
+            print >> sio, "  value max (new <space> old):", npy_new_val.max(), npy_old_val.max()
+            print >> sio, "  value min (new-old):", (npy_new_val-npy_old_val).min()
+            print >> sio, "  value max (new-old):", (npy_new_val-npy_old_val).max()
             print >> sio, ""
             print >> sio, "  Hint: this can also be caused by a deficient values_eq_approx() or __eq__() implementation that compares node input values"
             return sio.getvalue()
+        except Exception, e:
+            return str(e)
 
 class BadViewMap(DebugModeError):
     """Exception: Some perform() or c_code() created a memory alias that wasn't in the view_map"""
@@ -191,15 +220,25 @@ class StochasticOrder(DebugModeError):
 
 class InvalidValueError(DebugModeError):
     """Exception: some Op an output value that is inconsistent with the Type of that output"""
-    def __init__(self, r, v):
+    def __init__(self, r, v, client_node=None):
         super(InvalidValueError, self).__init__()
         self.r = r
         self.v = v
+        self.client_node = client_node
 
     def __str__(self):
         r, v = self.r, self.v
-        return "InvalidValueError: Variable %s,  Type %s, type(Value) %s, Value %s"\
-                % (str(r), str(r.type), str(type(v)), str(v)[0:100])
+        type_r = type(r)
+        type_v = type(v)
+        v_val = str(v)[0:100]
+        client_node = self.client_node
+        return """InvalidValueError
+        type(variable) = %(type_r)s
+        variable       = %(r)s
+        type(value)    = %(type_v)s
+        value          = %(v_val)s
+        client_node    = %(client_node)s
+        """ % locals()
 
 ########################
 #
@@ -390,7 +429,7 @@ def _lessbroken_deepcopy(a):
 
     Returns a copy of `a` that shares no internal storage with the original.  A deep copy.
     This function handles numpy arrays specially to avoid some bug I had one time... (possibly
-    about copying 1-d arrays?)
+    about copying 0-d arrays?)
     """
     # this exists because numpy copies are broken
     if type(a) is numpy.ndarray:
@@ -806,6 +845,7 @@ class _Linker(gof.link.LocalLinker):
         # This is the function that runs when you evaluate the graph
         #####
         def f():
+            print >> sys.stderr, "starting f"
             for x in no_recycling:
                 x[0] = None
 
@@ -828,6 +868,7 @@ class _Linker(gof.link.LocalLinker):
                 assert len(thunks_py) == len(order)
 
                 # transfer the initial values from the storage_map to the r_vals
+                debug("DEBUGMODE: transfer initial values")
                 for r in storage_map:
                     if (r.owner is None):
                         if (storage_map[r][0] is None):
@@ -849,18 +890,27 @@ class _Linker(gof.link.LocalLinker):
                 for i, (thunk_py, thunk_c, node) in enumerate(zip(thunks_py, thunks_c, order)):
                     this_node_destroyed_variables = set()
 
+                    debug(i, "DEBUGMODE: starting node", i, node)
+                    debug(i, "DEBUGMODE: inputs broadcastable", [i.type.broadcastable for i in node.inputs])
+                    debug(i, "DEBUGMODE: outputs broadcastable", [i.type.broadcastable for i in node.outputs])
+
                     # put a copy of each input into the storage_map
                     # also, check that inputs have valid values
                     for r in node.inputs:
                         assert isinstance(r, gof.Variable)
                         assert r in r_vals
+                        # print >> sys.stderr,i,  "DEBUGMODE: deepcopy input ", r
                         storage_map[r][0] = _lessbroken_deepcopy(r_vals[r])
                         if not r.type.is_valid_value(storage_map[r][0]):
-                            raise InvalidValueError(r, storage_map[r][0])
+                            raise InvalidValueError(r, storage_map[r][0], client_node=node)
 
+                    debug(i, "DEBUGMODE running thunk_py")
                     if thunk_py:
+                        try:
                             thunk_py()
-
+                        except utils.MethodNotDefined:
+                            thunk_py = None #shouldn't have put it into the list in the first place
+                    if thunk_py:
                         # check output values for type-correctness
                         for r in node.outputs:
                             if not r.type.is_valid_value(storage_map[r][0]):
@@ -872,9 +922,15 @@ class _Linker(gof.link.LocalLinker):
 
                         _check_viewmap(node, storage_map)
 
+                        # print >> sys.stderr, i, "DEBUGMODE thunk_py %100s %50s %30s" % (node,
+                            #[(id(o), numpy.asarray(storage_map[o][0])[0,0]) for o in node.inputs],
+                            #[(id(o), numpy.asarray(storage_map[o][0])[0,0]) for o in node.outputs])
+                        sys.stdout.flush()
+
                         #retrieve each output from the storage_map
                         for r in node.outputs:
                             assert r not in r_vals
+                            # print >> sys.stderr, i, "DEBUGMODE storing reference output %x" % id(storage_map[r][0])
                             r_vals[r] = storage_map[r][0]
                             storage_map[r][0] = None #clear the storage_map of outputs for the thunk_c
 
@@ -882,8 +938,10 @@ class _Linker(gof.link.LocalLinker):
 
                         for r in node.inputs:
                             # TODO:  we only need to overwrite the non-destroyed inputs
+                            # print >> sys.stderr, i, "DEBUGMODE replacing input", r
                             storage_map[r][0] = _lessbroken_deepcopy(r_vals[r])
 
+                        debug(i, "DEBUGMODE running thunk_c")
                         thunk_c()
 
                         for r in node.outputs:
@@ -902,21 +960,35 @@ class _Linker(gof.link.LocalLinker):
 
                         _check_viewmap(node, storage_map)
 
+                        # print >> sys.stderr, i, "DEBUGMODE thunk_c  %100s %50s %30s" % (node,
+                            #[(id(o), numpy.asarray(storage_map[o][0])[0,0]) for o in node.inputs],
+                            #[(id(o), numpy.asarray(storage_map[o][0])[0,0]) for o in node.outputs])
+                        sys.stdout.flush()
+
                         for r in node.outputs:
                             if r in r_vals:
+                                #print >> sys.stderr, i, "DEBUGMODE clearing output", r
                                 # compares the version from thunk_py (in r_vals)
                                 # to the version produced by thunk_c (in storage_map)
                                 if not r.type.values_eq_approx(r_vals[r], storage_map[r][0]):
                                     raise BadCLinkerOutput(r, val_py=r_vals[r], val_c=storage_map[r][0])
                             else:
+                                #print >> sys.stderr, i, "DEBUGMODE storing reference output %x" % id(storage_map[r][0])
                                 #retrieve each output from the storage_map
                                 r_vals[r] = storage_map[r][0]
                             storage_map[r][0] = None #clear the storage_map for the thunk_c
 
+
+
                     # we're done with this thunk
                     # clear everything out of the storage_map
                     for r in node.inputs:
+                        #print >> sys.stderr, i, "DEBUGMODE clearing input", r
                         storage_map[r][0] = None
+                    debug("done with node")
+
+                    if True: 
+                        gc.collect()
 
                 #except:
                 #    raise_with_op(node)

theano/gof/cc.py

@@ -352,7 +352,6 @@ class CLinker(link.Linker):
         self.env = env
         self.fetch_variables()
         self.no_recycling = no_recycling
-        self.module_compile_str = cmodule.gcc_module_compile_str
         return self
 
     def fetch_variables(self):
@@ -392,6 +391,8 @@ class CLinker(link.Linker):
 
         self.consts = []
 
+        c_support_code_apply = []
+
         symbol = {}
 
         # (init_)tasks contains a list of pairs (Op/Variable, task_name)
@@ -472,7 +473,7 @@ class CLinker(link.Linker):
 
             id += 2
 
-        for node in self.node_order:
+        for node_num, node in enumerate(self.node_order):
 
             # We populate sub with a mapping from the variable names specified by the op's c_var_names
             # method to the actual variable names that we will use.
@@ -481,7 +482,7 @@ class CLinker(link.Linker):
 ##            for variable, vname in zip(op.inputs + op.outputs, ivnames + ovnames):
 ##                sub[vname] = symbol[variable]
 
-            name = "<invalid_c_thing>"
+            name = "node_%i" % node_num
             isyms, osyms = [symbol[r] for r in node.inputs], [symbol[r] for r in node.outputs]
 
             # c_validate_update is deprecated
@@ -493,6 +494,11 @@ class CLinker(link.Linker):
             sub['fail'] = failure_code(sub)
 
             op = node.op
+            # type-specific support code
+            try: c_support_code_apply.append(op.c_support_code_apply(node, name))
+            except utils.MethodNotDefined: pass
+
+            # emit c_code 
             try: behavior = op.c_code(node, name, isyms, osyms, sub)
             except utils.MethodNotDefined:
                 raise NotImplementedError("%s cannot produce C code" % op)
@@ -529,6 +535,7 @@ class CLinker(link.Linker):
         self.blocks = blocks
         self.tasks = tasks
         all = self.inputs + self.outputs + self.orphans
+        self.c_support_code_apply = c_support_code_apply
 
         if (self.init_tasks, self.tasks) != self.get_init_tasks():
             print >> sys.stderr, "init_tasks\n", self.init_tasks
@@ -551,6 +558,7 @@ class CLinker(link.Linker):
         This might contain duplicates.
         """
         ret = []
+        # generic support code
         for x in [y.type for y in self.variables] + [y.op for y in self.node_order]:
             try: ret.append(x.c_support_code())
             except utils.MethodNotDefined: pass
@@ -589,28 +597,70 @@ class CLinker(link.Linker):
     def headers(self):
         """WRITEME
         Returns a list of headers that are needed by one
-        or more Variables or Ops.
+        or more Types or Ops.
 
-        This might contain duplicates.
+        The return value will not contain duplicates.
         """
         ret = []
         for x in [y.type for y in self.variables] + [y.op for y in self.node_order]:
             try: ret += x.c_headers()
             except utils.MethodNotDefined: pass
-        return ret
+        return list(set(ret))
+
+    def c_compiler(self):
+        c_compiler = None
+        for x in [y.type for y in self.variables] + [y.op for y in self.node_order]:
+            if hasattr(x, 'c_compiler'):
+                x_compiler = x.c_compiler()
+            else:
+                continue
+
+            if c_compiler is None:
+                c_compiler = x_compiler
+            else:
+                if x_compiler and (x_compiler != c_compiler):
+                    raise Exception('Nodes have requested specific different compilers',
+                            (c_compiler, x_compiler))
+        return cmodule.gcc_module_compile_str if (c_compiler is None) else c_compiler
+
+    def header_dirs(self):
+        """WRITEME
+        Returns a list of lib directories that are needed by one
+        or more Types or Ops.
+
+        The return value will not contain duplicates.
+        """
+        ret = []
+        for x in [y.type for y in self.variables] + [y.op for y in self.node_order]:
+            try: ret += x.c_header_dirs()
+            except utils.MethodNotDefined: pass
+        return list(set(ret))
 
     def libraries(self):
         """WRITEME
         Returns a list of libraries that are needed by one
-        or more Variables or Ops.
+        or more Types or Ops.
 
-        This might contain duplicates.
+        The return value will not contain duplicates.
         """
         ret = []
         for x in [y.type for y in self.variables] + [y.op for y in self.node_order]:
             try: ret += x.c_libraries()
             except utils.MethodNotDefined: pass
-        return ret
+        return list(set(ret))
+
+    def lib_dirs(self):
+        """WRITEME
+        Returns a list of lib directories that are needed by one
+        or more Types or Ops.
+
+        The return value will not contain duplicates.
+        """
+        ret = []
+        for x in [y.type for y in self.variables] + [y.op for y in self.node_order]:
+            try: ret += x.c_lib_dirs()
+            except utils.MethodNotDefined: pass
+        return list(set(ret))
 
     def __compile__(self, input_storage = None, output_storage = None):
         """WRITEME
@@ -774,11 +824,13 @@ class CLinker(link.Linker):
         get_lock()
         try:
             debug("LOCATION", location)
-            module = self.module_compile_str(
+            c_compiler = self.c_compiler()
+            module = c_compiler(
                     module_name=mod.name,
                     src_code = mod.code(),
                     location=location,
-                    include_dirs=[],
+                    include_dirs=self.header_dirs(),
+                    lib_dirs=self.lib_dirs(),
                     libs=self.libraries(),
                     preargs=self.compile_args())
         finally:
@@ -821,7 +873,7 @@ class CLinker(link.Linker):
         """ % dict(struct_name = self.struct_name)
 
         # We add all the support code, compile args, headers and libs we need.
-        for support_code in self.support_code():
+        for support_code in self.support_code() + self.c_support_code_apply:
             mod.add_support_code(support_code)
         mod.add_support_code(self.struct_code)
         mod.add_support_code(static)

theano/gof/cmodule.py

@@ -93,7 +93,12 @@ class DynamicModule(object):
     def code(self):
         sio = StringIO.StringIO()
         for inc in self.includes:
+            if not inc:
+                continue
+            if inc[0] == '<' or inc[0] == '"':
                 print >> sio, "#include", inc
+            else:
+                print >> sio, '#include "%s"'%inc
 
         print  >> sio, "//////////////////////"
         print  >> sio, "////  Support Code"
@@ -464,16 +469,11 @@ def get_lib_extension():
     else:
         return 'so'
 
-def gcc_module_compile_str(module_name, src_code, location=None, include_dirs=[], lib_dirs=[], libs=[],
-        preargs=[], tmpdir=None):
-    #TODO: don't to the dlimport in this function
-    preargs= [] if preargs is None else list(preargs)
-    preargs.append('-fPIC')
-    no_opt = False
 
-    include_dirs = [distutils.sysconfig.get_python_inc()] + \
-                    numpy.distutils.misc_util.get_numpy_include_dirs()\
-                    + include_dirs
+def std_include_dirs():
+    return [distutils.sysconfig.get_python_inc()] + numpy.distutils.misc_util.get_numpy_include_dirs()
+
+def std_lib_dirs_and_libs():
     python_inc = distutils.sysconfig.get_python_inc()
     if sys.platform == 'win32':
         # Typical include directory: C:\Python26\include
@@ -481,15 +481,43 @@ def gcc_module_compile_str(module_name, src_code, location=None, include_dirs=[]
         # Also add directory containing the Python library to the library
         # directories.
         python_lib_dir = os.path.join(os.path.dirname(python_inc), 'libs')
-        lib_dirs = [python_lib_dir] + lib_dirs
+        lib_dirs = [python_lib_dir]
+        return [libname], [python_lib_dir]
     else:
         # Typical include directory: /usr/include/python2.6
         libname = os.path.basename(python_inc)
-    libs = [libname] + libs
+        return [libname], []
+
+def std_libs():
+    return std_lib_dirs_and_libs()[0]
+
+def std_lib_dirs():
+    return std_lib_dirs_and_libs()[1]
+
+
+def gcc_module_compile_str(module_name, src_code, location=None, include_dirs=[], lib_dirs=[], libs=[],
+        preargs=[]):
+    """
+    :param module_name: string (this has been embedded in the src_code
+    :param src_code: a complete c or c++ source listing for the module
+    :param location: a pre-existing filesystem directory where the cpp file and .so will be written
+    :param include_dirs: a list of include directory names (each gets prefixed with -I)
+    :param lib_dirs: a list of library search path directory names (each gets prefixed with -L)
+    :param libs: a list of libraries to link with (each gets prefixed with -l)
+    :param preargs: a list of extra compiler arguments
+
+    :returns: dynamically-imported python module of the compiled code.
+    """
+    #TODO: don't to the dlimport in this function
+    preargs= [] if preargs is None else list(preargs)
+    preargs.append('-fPIC')
+    no_opt = False
 
-    workdir = location
+    include_dirs = std_include_dirs() + include_dirs
+    libs = std_libs() + libs
+    lib_dirs = std_lib_dirs() + lib_dirs
 
-    cppfilename = os.path.join(workdir, 'mod.cpp')
+    cppfilename = os.path.join(location, 'mod.cpp')
     cppfile = file(cppfilename, 'w')
 
     debug('Writing module C++ code to', cppfilename)
@@ -499,7 +527,7 @@ def gcc_module_compile_str(module_name, src_code, location=None, include_dirs=[]
     cppfile.write(src_code)
     cppfile.close()
 
-    lib_filename = os.path.join(workdir, '%s.%s' %
+    lib_filename = os.path.join(location, '%s.%s' %
             (module_name, get_lib_extension()))
 
     debug('Generating shared lib', lib_filename)
@@ -519,83 +547,13 @@ def gcc_module_compile_str(module_name, src_code, location=None, include_dirs=[]
     status = p.wait()
 
     if status:
-        error('g++ return status', status)
-    else:
-        #touch the __init__ file
-        file(os.path.join(workdir, "__init__.py"),'w').close()      
-
-        rval = dlimport(lib_filename)
-    return rval
-
-
-def nvcc_module_compile_str(module_name, src_code, location=None, include_dirs=[], lib_dirs=[], libs=[],
-        preargs=[], tmpdir=None):
-    preargs= [] if preargs is None else list(preargs)
-    preargs.append('-fPIC')
-    no_opt = False
-
-
-    raise NotImplementedError()
-
-    #TODO: -O preargs should be passed globally, not to -Xcompiler
+        raise Exception('g++ return status', status)
 
-    #TODO: where to find these strings?  sys? distutils?
-    include_dirs = ['/usr/include/python2.6'] + include_dirs
-    libs = ['python2.6', 'cudart'] + libs
-    lib_dirs = ['/usr/local/cuda/lib']+lib_dirs
-
-    workdir = tempfile.mkdtemp(dir=location)
-
-    cppfilename = os.path.join(workdir, 'mod.cpp') #.cpp to use g++
-    cppfilename = os.path.join(workdir, 'mod.cu') #.cu to use nvopencc
-    cppfile = file(cppfilename, 'w')
-
-    debug('Writing module C++ code to', cppfilename)
-    ofiles = []
-    rval = None
-    try:
-        cppfile.write(src_code)
-        cppfile.close()
-        lib_filename = os.path.join(workdir, '%s.%s' %
-                (module_name, get_lib_extension()))
-
-        debug('Generating shared lib', lib_filename)
-        cmd = ['nvcc', '-shared', '-g']
-        cmd.extend(['-Xcompiler', ','.join(preargs)])
-        cmd.extend('-I%s'%idir for idir in include_dirs)
-        cmd.extend(['-o',lib_filename]) 
-        cmd.append(cppfilename)
-        cmd.extend(['-L%s'%ldir for ldir in lib_dirs])
-        cmd.extend(['-l%s'%l for l in libs])
-        debug('Running cmd', ' '.join(cmd))
-
-        p = subprocess.Popen(cmd)
-        status = p.wait()
-
-        if status:
-            warning('nvcc return status', status)
-        else:
     #touch the __init__ file
-            file(os.path.join(workdir, "__init__.py"),'w').close()      
-
-            #load the module
-            sys.path.insert(0, workdir)
-            try:
-                rval = __import__(module_name, {}, {}, [module_name])
-                if not rval:
-                    debug('__import__ failed')
-            finally:
-                del sys.path[0]
+    file(os.path.join(location, "__init__.py"),'w').close()      
 
-            assert pathcopy == sys.path
+    return dlimport(lib_filename)
 
-    finally:
-        warning("TODO: cleanup")
-        #os.remove(cppfilename)
-        for ofile in ofiles:
-            #os.remove(ofiles[0])
-            pass
-    return rval
 
 def icc_module_compile_str(*args):
     raise NotImplementedError()

theano/gof/op.py

@@ -11,8 +11,138 @@ __docformat__ = "restructuredtext en"
 import utils
 import traceback
 
+class CLinkerObject(object):
+    """Standard elements of an Op or Type used with the CLinker
+    """
+    def c_headers(self):
+        """Optional: Return a list of header files required by code returned by
+        this class.
+
+        For example: return ['<iostream>', '<math.h>', '/full/path/to/header.h']
+
+        These strings will be prefixed with "#include " and inserted at the beginning of the c
+        source code.
+
+        Strings in this list that start neither with '<' nor '"' will be enclosed in
+        double-quotes.
+
+        :Exceptions:
+         - `MethodNotDefined`: Subclass does not implement this method
+
+        """
+        raise utils.MethodNotDefined("c_headers", type(self), self.__class__.__name__)
+
+    def c_header_dirs(self):
+        """Optional: Return a list of header search paths required by code returned by
+        this class.
+
+        For example: return ['/usr/local/include', '/opt/weirdpath/src/include'].
+
+        Provide search paths for headers, in addition to those in any relevant environment
+        variables.
+        
+        Hint: for unix compilers, these are the things that get '-I' prefixed in the compiler
+        cmdline.
+        
+        :Exceptions:
+         - `MethodNotDefined`: Subclass does not implement this method
+
+        """
+        raise utils.MethodNotDefined("c_lib_dirs", type(self), self.__class__.__name__)
+
+    def c_libraries(self):
+        """Optional: Return a list of libraries required by code returned by
+        this class.
+
+        For example: return ['gsl', 'gslcblas', 'm', 'fftw3', 'g2c'].
+
+        The compiler will search the directories specified by the environment
+        variable LD_LIBRARY_PATH in addition to any returned by `c_lib_dirs`.
+
+        Hint: for unix compilers, these are the things that get '-l' prefixed in the compiler
+        cmdline.
+
+        :Exceptions:
+         - `MethodNotDefined`: Subclass does not implement this method
+
+        """
+        raise utils.MethodNotDefined("c_libraries", type(self), self.__class__.__name__)
+
+    def c_lib_dirs(self):
+        """Optional: Return a list of library search paths required by code returned by
+        this class.
+
+        For example: return ['/usr/local/lib', '/opt/weirdpath/build/libs'].
+
+        Provide search paths for libraries, in addition to those in any relevant environment
+        variables (e.g. LD_LIBRARY_PATH).
+        
+        Hint: for unix compilers, these are the things that get '-L' prefixed in the compiler
+        cmdline.
+        
+        :Exceptions:
+         - `MethodNotDefined`: Subclass does not implement this method
+
+        """
+        raise utils.MethodNotDefined("c_lib_dirs", type(self), self.__class__.__name__)
+
+    def c_support_code(self):
+        """Optional: Return utility code for use by a `Variable` or `Op` to be
+        included at global scope prior to the rest of the code for this class.
+
+        QUESTION: How many times will this support code be emitted for a graph
+        with many instances of the same type?
+
+        :Exceptions:
+         - `MethodNotDefined`: Subclass does not implement this method
+
+        """
+        raise utils.MethodNotDefined("c_support_code", type(self), self.__class__.__name__)
+
+    def c_code_cache_version(self):
+        """Return a tuple of integers indicating the version of this Op.
+
+        An empty tuple indicates an 'unversioned' Op that will not be cached between processes.
+
+        The cache mechanism may erase cached modules that have been superceded by newer
+        versions.  See `ModuleCache` for details.
+        """
+        return (1,)
+
+    def c_compile_args(self):
+        """Optional: Return a list of compile args recommended to compile the
+        code returned by other methods in this class.
+
+        Example: return ['-ffast-math']
+
+        Compiler arguments related to headers, libraries and search paths should be provided
+        via the functions `c_headers`, `c_libraries`, `c_header_dirs`, and `c_lib_dirs`.
+
+        :Exceptions:
+         - `MethodNotDefined`: Subclass does not implement this method
 
-class CLinkerOp(object):
+        """
+        raise utils.MethodNotDefined("c_compile_args", type(self), self.__class__.__name__)
+
+    def c_no_compile_args(self):
+        """Optional: Return a list of incompatible gcc compiler arguments.
+
+        We will remove those arguments from the command line of gcc. So if 
+        another Op adds a compile arg in the graph that is incompatible 
+        with this Op, the incompatible arg will not be used. 
+        Useful for instance to remove -ffast-math.
+
+        EXAMPLE
+
+        WRITEME
+
+        :Exceptions:
+         - `MethodNotDefined`: the subclass does not override this method
+
+        """
+        raise utils.MethodNotDefined("c_no_compile_args", type(self), self.__class__.__name__)
+
+class CLinkerOp(CLinkerObject):
     """
     Interface definition for `Op` subclasses compiled by `CLinker`.
 
@@ -83,94 +213,22 @@ class CLinkerOp(object):
         raise utils.MethodNotDefined('%s.c_code_cleanup' \
                 % self.__class__.__name__)
 
-    def c_compile_args(self):
-        """Optional: Return a list of recommended gcc compiler arguments.
-
-        QUESTION: is this function optional?
-        
-        This is only a hint.
-
-        EXAMPLE
-
-        WRITEME
-        """
-        raise utils.MethodNotDefined('%s.c_compile_args' \
-                % self.__class__.__name__)
-
-    def c_no_compile_args(self):
-        """Optional: Return a list of incompatible gcc compiler arguments.
-
-        We will remove those arguments from the command line of gcc. So if 
-        another Op adds a compile arg in the graph that is incompatible 
-        with this Op, the incompatible arg will not be used. 
-        Useful for instance to remove -ffast-math.
+    def c_support_code_apply(self, node, name):
+        """Optional: Return utility code for use by an `Op` that will be inserted at global
+        scope, that can be specialized for the support of a particular `Apply` node.
 
-        EXAMPLE
+        :param node: an Apply instance in the graph being compiled
 
-        WRITEME
+        :param node_id: a string or number that serves to uniquely identify this node.
+        Symbol names defined by this support code should include the node_id, so that they can
+        be called from the c_code, and so that they do not cause name collisions.
 
         :Exceptions:
-         - `MethodNotDefined`: the subclass does not override this method
+         - `MethodNotDefined`: Subclass does not implement this method
 
         """
-        raise utils.MethodNotDefined('%s.c_no_compile_args' \
-                % self.__class__.__name__)
-
-    def c_headers(self):
-        """Optional: Return a list of header files that must be included to compile the C code.
-
-        A subclass should override this method.
-
-        EXAMPLE
-
-        WRITEME
-
-        :Exceptions:
-         - `MethodNotDefined`: the subclass does not override this method
+        raise utils.MethodNotDefined("c_support_code_apply", type(self), self.__class__.__name__)
 
-        """
-        raise utils.MethodNotDefined('%s.c_headers' \
-                % self.__class__.__name__)
-
-    def c_libraries(self):
-        """Optional: Return a list of libraries to link against to manipulate this `Op`.
-
-        A subclass should override this method.
-
-        WRITEME
-
-        :Exceptions:
-         - `MethodNotDefined`: the subclass does not override this method
-
-        """
-        raise utils.MethodNotDefined('%s.c_libraries' \
-                % self.__class__.__name__)
-
-    def c_support_code(self):
-        """Optional: Return support code for use by the code that is returned by `c_code`.
-        
-        Support code is inserted into the C code at global scope.
-
-        A subclass should override this method.
-
-        WRITEME
-
-        :Exceptions:
-         - `MethodNotDefined`: the subclass does not override this method
-
-        """
-        raise utils.MethodNotDefined('%s.c_support_code' \
-                % self.__class__.__name__)
-
-    def c_code_cache_version(self):
-        """Return a tuple of integers indicating the version of this Op.
-
-        An empty tuple indicates an 'unversioned' Op that will not be cached between processes.
-
-        The cache mechanism may erase cached modules that have been superceded by newer
-        versions.  See `ModuleCache` for details.
-        """
-        return (1,)
 
 class PureOp(object):
     """

theano/gof/opt.py

@@ -745,6 +745,7 @@ class NavigatorOptimizer(Optimizer):
                 raise
         if replacements is False or replacements is None:
             return False
+        assert len(node.outputs) == len(replacements)
         repl_pairs = zip(node.outputs, replacements)
         try:
             env.replace_all_validate(repl_pairs, reason=lopt)

theano/gof/optdb.py

-
+import sys, StringIO
 from collections import defaultdict
 import opt
 
@@ -140,4 +140,14 @@ class SequenceDB(DB):
         opts.sort(key = lambda obj: self.__priority__[obj.name])
         return opt.SeqOptimizer(opts, failure_callback = self.failure_callback)
 
+    def print_summary(self, stream=sys.stdout):
+        print >> stream, "SequenceDB (id %i)"%id(self)
+        print >> stream, "  priority", self.__priority__
+        print >> stream, "  names", self._names
+        print >> stream, "  db", self.__db__
+    def __str__(self):
+        sio = StringIO.StringIO()
+        self.print_summary(sio)
+        return sio.getvalue()
+
 

theano/gof/type.py

@@ -12,8 +12,9 @@ import traceback
 ########
 # Type #
 ########
+from .op import CLinkerObject
 
-class CLinkerType(object):
+class CLinkerType(CLinkerObject):
     """Interface specification for Types that can be arguments to a `CLinkerOp`.
 
     A CLinkerType instance is mainly reponsible  for providing the C code that
@@ -176,89 +177,8 @@ class CLinkerType(object):
         """
         raise MethodNotDefined("c_sync", type(self), self.__class__.__name__)
 
-    def c_compile_args(self):
-        """Optional: Return a list of compile args recommended to compile the
-        code returned by other methods in this class.
 
-        WRITEME: example of formatting for -I, -L, -f args.
 
-        :Exceptions:
-         - `MethodNotDefined`: Subclass does not implement this method
-
-        """
-        raise MethodNotDefined("c_compile_args", type(self), self.__class__.__name__)
-
-    def c_no_compile_args(self):
-        """Optional: Return a list of incompatible gcc compiler arguments.
-
-        We will remove those arguments from the command line of gcc. So if 
-        another Op adds a compile arg in the graph that is incompatible 
-        with this Op, the incompatible arg will not be used. 
-        Useful for instance to remove -ffast-math.
-
-        EXAMPLE
-
-        WRITEME
-
-        :Exceptions:
-         - `MethodNotDefined`: the subclass does not override this method
-
-        """
-        raise MethodNotDefined("c_no_compile_args", type(self), self.__class__.__name__)
-
-    def c_headers(self):
-        """Optional: Return a list of header files required by code returned by
-        this class.
-
-        WRITEME: example of local file, standard file.
-
-        :Exceptions:
-         - `MethodNotDefined`: Subclass does not implement this method
-
-        """
-        raise MethodNotDefined("c_headers", type(self), self.__class__.__name__)
-
-    def c_libraries(self):
-        """Optional: Return a list of libraries required by code returned by
-        this class.
-
-        For example: return ['gsl', 'gslcblas', 'm', 'fftw3', 'g2c'].
-
-        The compiler will search the directories specified by the environment
-        variable LD_LIBRARY_PATH.  No option is provided for an Op to provide an
-        extra library directory because this would change the linking path for
-        other Ops in a potentially disasterous way.
-
-        QUESTION: What about via the c_compile_args? a -L option is allowed no?
-
-        :Exceptions:
-         - `MethodNotDefined`: Subclass does not implement this method
-
-        """
-        raise MethodNotDefined("c_libraries", type(self), self.__class__.__name__)
-
-    def c_support_code(self):
-        """Optional: Return utility code for use by a `Variable` or `Op` to be
-        included at global scope prior to the rest of the code for this class.
-
-        QUESTION: How many times will this support code be emitted for a graph
-        with many instances of the same type?
-
-        :Exceptions:
-         - `MethodNotDefined`: Subclass does not implement this method
-
-        """
-        raise MethodNotDefined("c_support_code", type(self), self.__class__.__name__)
-
-    def c_code_cache_version(self):
-        """Return a tuple of integers indicating the version of this Op.
-
-        An empty tuple indicates an 'unversioned' Op that will not be cached between processes.
-
-        The cache mechanism may erase cached modules that have been superceded by newer
-        versions.  See `ModuleCache` for details.
-        """
-        return (1,)
 
 class PureType(object):
     """Interface specification for variable type instances.

theano/gradient.py

@@ -9,9 +9,12 @@ import numpy #for numeric_grad
 from gof.python25 import all
 import gof.utils
 
-def warning(msg):
-    # replace this with logger.warning when adding logging support
-    print >> sys.stderr, 'WARNING', msg
+import logging
+_logger = logging.getLogger('theano.gradient')
+def warning(*msg):
+    _logger.warning('WARNING theano.gradient: '+' '.join(msg))
+def info(*msg):
+    _logger.info('INFO theano.gradient: '+' '.join(msg))
 
 _msg_retType = 'op.grad(...) returned a non-list'
 _msg_badlen = 'op.grad(...) returned wrong number of gradients'
@@ -103,7 +106,9 @@ def grad_sources_inputs(sources, graph_inputs, warn_type=True):
         for ii, (r, g_r) in enumerate(zip(node.inputs, g_inputs)):
             if warn_type:
                 if g_r and (getattr(r,'type',0) != getattr(g_r,'type', 1)):
-                    warning('%s.grad returned a different type for input %i: %s vs. %s'%(node.op, ii, r, g_r))
+                    r_type = getattr(r,'type', None)
+                    g_r_type = getattr(g_r,'type', None)
+                    info('%s.grad returned a different type for input %i: %s vs.  %s'%(node.op, ii, r_type, g_r_type))
             if g_r and len(sources) == 1 and sources[0][0].name and r.name:
                 g_r.name = "(d%s/d%s)" % (sources[0][0].name, r.name)
             if g_r is not None: 

theano/sandbox/conv.py

@@ -188,7 +188,6 @@ class ConvOp(Op):
             buf = N.zeros((batchsize,)+ self.imshp_logical, dtype=img2d.dtype)
             buf[:,:,::rstride, ::cstride] = img2d
             img2d = buf
-            print 'A'
             del buf, rstride, cstride
 
         if self.kshp != self.kshp_logical:
@@ -204,7 +203,6 @@ class ConvOp(Op):
                 assert coffset >= 0
             buf[:,:,roffset::rstride, coffset::cstride] = filtersflipped
             filtersflipped = buf
-            print 'B'
             del buf, rstride, cstride
 
         for b in range(batchsize):

theano/scalar/basic.py

@@ -31,10 +31,25 @@ def as_scalar(x, name = None):
 
 def constant(x):
     if isinstance(x, float):
-        return ScalarConstant(float64, x)
+        for dtype in ['float32', 'float64']:
+            x_ = numpy.asarray(x, dtype=dtype)
+            if numpy.all(x == x_):
+                break
+            x_ = None
+        assert x_ is not None
+        return ScalarConstant(Scalar(str(x_.dtype)), x)
     if isinstance(x, int):
-        return ScalarConstant(int64, x)
-    return ScalarConstant(float64, float(x))
+        for dtype in ['int8', 'int16', 'int32', 'int64']:
+            x_ = numpy.asarray(x, dtype=dtype)
+            if numpy.all(x == x_):
+                break
+            x_ = None
+        assert x_ is not None
+        return ScalarConstant(Scalar(str(x_.dtype)), x)
+    if isinstance(x, complex):
+        raise NotImplementedError()
+    raise TypeError(x)
+    #return ScalarConstant(float64, float(x))
 
 
 class Scalar(Type):
@@ -192,9 +207,9 @@ class _scalar_py_operators:
     def __neg__(self): return neg(self)
 
     #CASTS
-    def __int__(self): return AsInt(self).out
-    def __float__(self): return AsInt(self).out
-    def __complex__(self): return AsComplex(self).out
+    #def __int__(self): return AsInt(self).out
+    #def __float__(self): return AsDouble(self).out
+    #def __complex__(self): return AsComplex(self).out
 
     #BITWISE
     def __invert__(self): return invert(self) 

theano/tensor/basic.py

@@ -140,11 +140,26 @@ def constant_or_value(x, rtype, name=None, ndim=None):
      - `ValueError`: `x` could not be expanded to have ndim dimensions
 
     """
-    if isinstance(x, numpy.ndarray):
+    x_ = None
+    if rtype is TensorConstant and isinstance(x, int):
+        for dtype in ['int8', 'int16', 'int32', 'int64']:
+            x_ = numpy.asarray(x, dtype=dtype)
+            if numpy.all(x == x_):
+                break
+            x_ = None
+    elif rtype is TensorConstant and isinstance(x, float):
+        for dtype in ['float32', 'float64']:
+            x_ = numpy.asarray(x, dtype=dtype)
+            if numpy.all(x == x_):
+                break
+            x_ = None
+    elif isinstance(x, numpy.ndarray):
         x_ = x
     else:
         x_ = numpy.asarray(x)
 
+    assert type(x_) == numpy.ndarray
+
     bcastable = [d == 1 for d in x_.shape]
     if ndim is not None:
         if len(bcastable) < ndim:
@@ -261,7 +276,8 @@ class TensorType(Type):
         """Compare True iff other is the same kind of TensorType"""
         return type(self) == type(other) and other.dtype == self.dtype and other.broadcastable == self.broadcastable
 
-    def values_eq_approx(self, a, b):
+    @staticmethod
+    def values_eq_approx(a, b):
         if type(a) is numpy.ndarray and type(b) is numpy.ndarray:
             if a.shape != b.shape:
                 return False
@@ -653,18 +669,22 @@ class _tensor_py_operators:
         """
         return reshape(self, shape, ndim=ndim)
 
-    def dimshuffle(self, pattern):
+    def dimshuffle(self, *pattern):
         """Reorder the dimensions of this variable, optionally inserting broadcasted dimensions.
 
-        :param pattern: list of int mixed with 'x' for broadcastable dimensions
+        :param pattern: list/tuple of int mixed with 'x' for broadcastable dimensions
 
         For example, to create a 3D view of a [2D] matrix, call ``dimshuffle([0,'x',1])``.  This
         will create a 3D view such that the middle dimension is an implicit broadcasted
         dimension.  To do the same thing on the transpose of that matrix, call ``dimshuffle([1,
         'x', 0])``.
 
+        This function supports the pattern passed as a tuple, or as a variable-length argument (e.g. ``a.dimshuffle(pattern)`` is equivalent to ``a.dimshuffle(*pattern)`` where ``pattern`` is a list/tuple of ints mixed with 'x' characters).
+
         For more information, see `DimShuffle`.
         """
+        if (len(pattern) == 1) and (isinstance(pattern[0], (list, tuple))):
+            pattern = pattern[0]
         op = DimShuffle(list(self.type.broadcastable), pattern)
         return op(self)
 
@@ -683,7 +703,8 @@ class _tensor_py_operators:
         return Subtensor(args)(self, *Subtensor.collapse(args, lambda entry: isinstance(entry, Variable)))
     
     #COPYING
-    def copy(self): return tensor_copy(self)
+    def copy(self):
+        return tensor_copy(self)
 
     def __iter__(self): 
         try:
@@ -1542,6 +1563,7 @@ class Subtensor(Op):
         return type(self) == type(other) and self.idx_list == other.idx_list
 
     def __hash__(self):
+        #TODO: optimize by cache this hash value
         idx_list = tuple((entry.start, entry.stop, entry.step)
                          if isinstance(entry, slice)
                          else entry
@@ -2096,9 +2118,13 @@ class Reshape(Op):
         self.name = name
 
     def __eq__(self, other):
-        return (type(other) is Reshape) and (other.ndim == self.ndim) and self.name == other.name
+        # .name does not participate because it doesn't affect computations
+        return (type(other) == type(self)) and (other.ndim == self.ndim)
     def __hash__(self):
-        return hash(Reshape) ^ hash(self.ndim) ^ hash(self.name)
+        # .name does not participate because it doesn't affect computations
+        return hash(type(self)) ^ hash(self.ndim)
+    def __str__(self):
+        return '%s{%i}' % (self.__class__.__name__, self.ndim)
     def make_node(self, x, shp):
         x = as_tensor_variable(x)
         shp = as_tensor_variable(shp)

theano/tensor/blas.py

@@ -4,8 +4,8 @@ import os, sys, traceback
 import numpy
 
 from ..gof import (utils, Op, Apply, view_roots, PatternSub, DestroyHandler, 
-        SeqOptimizer, local_optimizer, LocalOptimizer, OpKeyOptimizer, 
-        InconsistencyError)
+        SeqOptimizer, local_optimizer, Optimizer, LocalOptimizer, OpKeyOptimizer, 
+        InconsistencyError, toolbox)
 from ..printing import pprint, FunctionPrinter
 from .opt import register_specialize, out2in, insert_inplace_optimizer
 # opt.py
@@ -369,30 +369,98 @@ def res_is_a(node, op, maxclients=None):
             and node.owner.op == op \
             and (len(node.clients) <= maxclients if maxclients is not None else True)
 
-class GemmLocalOptimizer(LocalOptimizer):
-    """This is a massive beast for recognizing all the ways that a subtraction or addition
-    could be replaced by a GEMM
+def _as_scalar(res):
+    """Return None or a TensorVariable whose type is in T.float_scalar_types"""
+    if res.owner and isinstance(res.owner.op, T.DimShuffle):
+        return _as_scalar(res.owner.inputs[0])
+    elif res.type in T.float_scalar_types:
+        return res
+    elif isinstance(res, T.Constant) and res.data.size == 1:
+        return res.data.flatten()[0]
+    else:
+        return None
 
-    It depends on `local_transposed_dot` to canonicalize the graph a bit by swapping
-    dot(a,b).T -> dot(b.T, a.T)
-    """
+def _is_real_matrix(res):
+    return res.type.dtype in ('float32', 'float64') \
+            and res.type.ndim == 2 \
+            and res.type.broadcastable[0] == False \
+            and res.type.broadcastable[1] == False #cope with tuple vs. list
 
-    def __init__(self):
-        super(LocalOptimizer, self).__init__()
+def _as_isolated_scalar_times_matrix(res):
+    if res_is_a(res, T.mul, 1):
+        if len(res.owner.inputs) == 2:
+            L, R = res.owner.inputs
+            sL = _as_scalar(L)
+            sR = _as_scalar(R)
+            if (sL is not None) and _is_real_matrix(R):
+                return (sL, R)
+            if (sR is not None) and _is_real_matrix(L):
+                return (sR, L)
+        else:
+            scalars = []
+            matrices = []
+            for input in res.owner.inputs:
+                scalar_input = _as_scalar(input)
+                if scalar_input is not None:
+                    scalars.append(scalar_input)
+                elif _is_real_matrix(input):
+                    matrices.append(input)
+                else:
+                    return None
+            if len(matrices) == 1:
+                rval = (T.mul(*scalars), matrices[0])
+                return rval
 
-    def op_key(self):
-        return [T.add, T.sub]
+def _beta_L_plus_alpha_M(beta, L, alpha, M, recurse_flip = True):
+    #print 'BETA L + ALPHA M', beta, L, alpha, M, recurse_flip
+    #EXPRESSION: (beta * L) + (alpha * M)
 
-    def add_requirements(self, env):
-        super(GemmLocalOptimizer,self).add_requirements(env)
-        env.extend(DestroyHandler())
+    if res_is_a(M, _dot22, 1):
+        Ml, Mr = M.owner.inputs
+        rval = [gemm(L, alpha, Ml, Mr, beta)]
+        #print 'GEMM 0', rval, beta, L, alpha, M
+        return rval
+
+    # this is False'd out because of inadequate testing.  
+    # TODO see ticket #237
+    if False and res_is_a(M, gemm, 1):
+        #EXPRESSION: (beta * L) + (alpha * (gemm(G, a, u, v, b)))
+        #EXPRESSION: (beta * L) + alpha * (b * G) + alpha * a * dot(u, v)
+        G, a, u, v, b = M.owner.inputs
+        #print 'GEMM', G, L
+
+        if res_is_a(G, _dot22, 1):
+            #EXPRESSION: (beta * L) + (alpha * (gemm(dot(x,y), a, u, v, b)))
+            x, y = G.owner.inputs
+
+            #EXPRESSION: (beta * L) + (alpha * ((b*dot(x,y) + (a * dot(u, v)))))
+            #EXPRESSION: (beta * L) + (alpha*b*dot(x,y)) + (alpha * a * dot(u, v))
+            rval = [gemm(gemm(L, alpha * b, x, y, beta), alpha * a, u, v, 1.0)]
+            print 'GEMM 1', rval
+            return rval
+        if (G is L):
+            #EXPRESSION: (beta * L) + (alpha*b*L) + (alpha * a * dot(u, v))
+            rval = [gemm(L, alpha*a, u, v, alpha * b + beta)]
+            print 'GEMM 2', rval
+            return rval
+        if (1.0 != alpha):
+            #at the very least, move the alpha inside the gemm
+            rval = [beta * L + gemm(G, alpha * a, u, v, alpha * b)]
+            print 'GEMM 3', rval
+            return rval
+
+    if recurse_flip:
+        return _beta_L_plus_alpha_M(alpha, M, beta, L, recurse_flip = False)
+    else:
+        return False
+
+def _gemm_from_node(node):
+    """
+    :todo: In many expressions, there are many ways to turn it into a gemm.  For example
+    dot(a,b) + c + d.  This function should return all of them, so that if one version of gemm
+    causes a cycle in the graph, then another application of gemm can be tried.
 
-    def transform(self, node):
-        _as_scalar, _is_real_matrix, _as_isolated_scalar_times_matrix, beta_L_plus_alpha_M\
-                = (GemmLocalOptimizer._as_scalar, 
-                        GemmLocalOptimizer._is_real_matrix, 
-                        GemmLocalOptimizer._as_isolated_scalar_times_matrix, 
-                        GemmLocalOptimizer.beta_L_plus_alpha_M)
+    """
     if node.op == T.sub:
         L, R = node.inputs
         if not _is_real_matrix(L):
@@ -411,7 +479,7 @@ class GemmLocalOptimizer(LocalOptimizer):
             sR, mR = tmp
         except:
             sR, mR = 1.0, R
-            rval = beta_L_plus_alpha_M(sL, mL, -sR, mR)
+        rval = _beta_L_plus_alpha_M(sL, mL, -sR, mR)
         return rval
     if node.op == T.add:
         # arguments of the form scalar * matrix
@@ -439,7 +507,7 @@ class GemmLocalOptimizer(LocalOptimizer):
 
         if len(sM_list) == 2:
             (sL, mL), (sR, mR) = sM_list
-                gemm_of_sM_list = beta_L_plus_alpha_M(sL, mL, sR, mR)
+            gemm_of_sM_list = _beta_L_plus_alpha_M(sL, mL, sR, mR)
             if gemm_of_sM_list: 
                 #we turned the two candidates into a gemm
                 # now we have to add the other_inputs and return the replacement graph
@@ -454,7 +522,7 @@ class GemmLocalOptimizer(LocalOptimizer):
                     assert i != j
                     sL, mL = sM_list[i]
                     sR, mR = sM_list[j]
-                        gemm_of_sM_list = beta_L_plus_alpha_M(sL, mL, sR, mR)
+                    gemm_of_sM_list = _beta_L_plus_alpha_M(sL, mL, sR, mR)
                     if gemm_of_sM_list:
                         assert len(gemm_of_sM_list) == 1
                         inputs_without_ij = [input for k, input in enumerate(sM_orig) if k not in (i,j)]
@@ -477,115 +545,43 @@ class GemmLocalOptimizer(LocalOptimizer):
                         return [T.add(*new_add_inputs)]
     return False
 
-    @staticmethod
-    def _as_scalar(res):
-        """Return None or a TensorVariable whose type is in T.float_scalar_types"""
-        if res.owner and isinstance(res.owner.op, T.DimShuffle):
-            return GemmLocalOptimizer._as_scalar(res.owner.inputs[0])
-        elif res.type in T.float_scalar_types:
-            return res
-        elif isinstance(res, T.Constant) and res.data.size == 1:
-            return res.data.flatten()[0]
-        else:
-            return None
-
-    @staticmethod
-    def _is_real_matrix(res):
-        return res.type in T.float_matrix_types \
-                and res.type.broadcastable[0] == False \
-                and res.type.broadcastable[1] == False #cope with tuple vs. list
-
-    @staticmethod
-    def _as_isolated_scalar_times_matrix(res):
-        _as_scalar, _is_real_matrix, _as_isolated_scalar_times_matrix, beta_L_plus_alpha_M\
-                = (GemmLocalOptimizer._as_scalar, 
-                        GemmLocalOptimizer._is_real_matrix, 
-                        GemmLocalOptimizer._as_isolated_scalar_times_matrix, 
-                        GemmLocalOptimizer.beta_L_plus_alpha_M)
-        if res_is_a(res, T.mul, 1):
-            if len(res.owner.inputs) == 2:
-                L, R = res.owner.inputs
-                sL = _as_scalar(L)
-                sR = _as_scalar(R)
-                if (sL is not None) and _is_real_matrix(R):
-                    return (sL, R)
-                if (sR is not None) and _is_real_matrix(L):
-                    return (sR, L)
-            else:
-                scalars = []
-                matrices = []
-                for input in res.owner.inputs:
-                    scalar_input = _as_scalar(input)
-                    if scalar_input is not None:
-                        scalars.append(scalar_input)
-                    elif _is_real_matrix(input):
-                        matrices.append(input)
-                    else:
-                        return None
-                if len(matrices) == 1:
-                    rval = (T.mul(*scalars), matrices[0])
-                    return rval
-
-    @staticmethod
-    def beta_L_plus_alpha_M(beta, L, alpha, M, recurse_flip = True):
-        #print 'BETA L + ALPHA M', beta, L, alpha, M, recurse_flip
-        #EXPRESSION: (beta * L) + (alpha * M)
-
-        if res_is_a(M, _dot22, 1):
-            Ml, Mr = M.owner.inputs
-            rval = [gemm(L, alpha, Ml, Mr, beta)]
-            #print 'GEMM 0', rval, beta, L, alpha, M
-            return rval
-
-        # this is False'd out because of inadequate testing.  
-        # TODO see ticket #237
-        if False and res_is_a(M, gemm, 1):
-            #EXPRESSION: (beta * L) + (alpha * (gemm(G, a, u, v, b)))
-            #EXPRESSION: (beta * L) + alpha * (b * G) + alpha * a * dot(u, v)
-            G, a, u, v, b = M.owner.inputs
-            #print 'GEMM', G, L
+class GemmOptimizer(Optimizer):
+    """Graph optimizer for inserting Gemm operations"""
+    def __init__(self):
+        Optimizer.__init__(self)
 
-            if res_is_a(G, _dot22, 1):
-                #EXPRESSION: (beta * L) + (alpha * (gemm(dot(x,y), a, u, v, b)))
-                x, y = G.owner.inputs
+    def add_requirements(self, env):
+        env.extend(toolbox.ReplaceValidate())
+        env.extend(DestroyHandler())
 
-                #EXPRESSION: (beta * L) + (alpha * ((b*dot(x,y) + (a * dot(u, v)))))
-                #EXPRESSION: (beta * L) + (alpha*b*dot(x,y)) + (alpha * a * dot(u, v))
-                rval = [gemm(gemm(L, alpha * b, x, y, beta), alpha * a, u, v, 1.0)]
-                print 'GEMM 1', rval
-                return rval
-            if (G is L):
-                #EXPRESSION: (beta * L) + (alpha*b*L) + (alpha * a * dot(u, v))
-                rval = [gemm(L, alpha*a, u, v, alpha * b + beta)]
-                print 'GEMM 2', rval
-                return rval
-            if (1.0 != alpha):
-                #at the very least, move the alpha inside the gemm
-                rval = [beta * L + gemm(G, alpha * a, u, v, alpha * b)]
-                print 'GEMM 3', rval
-                return rval
+    def apply(self, env):
+        did_something = True
+        while did_something:
+            nodelist = list(env.toposort())
+            did_something = False
+            for node in nodelist:
+                new_outputs = _gemm_from_node(node)
+                if new_outputs:
+                    assert len(new_outputs) == len(node.outputs)
+                    try:
+                        env.replace_all_validate(
+                                zip(node.outputs, new_outputs),
+                                reason = 'GemmOptimizer')
+                        did_something = True
+                        break
+                    except InconsistencyError, e:
+                        #TODO: retry other applications of gemm (see comment in _gemm_from_node
+                        pass
 
-        if recurse_flip:
-            return GemmLocalOptimizer.beta_L_plus_alpha_M(alpha, M, beta, L, recurse_flip = False)
-        else:
-            return False
+compile.optdb.register('inplace_gemm', GemmOptimizer(), 70.00, 'fast_run', 'inplace', 'gemm')
 
-#I think that three passes should suffice to catch all the GEMMs.
-# TODO: This could be an equilibriumOptmizer, but I don't know how to combine an OpKeyOptimizer and
-# an EquilibriumOptimizer.
-compile.optdb.register('inplace_gemm_0', OpKeyOptimizer(GemmLocalOptimizer(), 
-    failure_callback=OpKeyOptimizer.warn_inplace), 70.00, 'fast_run', 'inplace', 'gemm')
-compile.optdb.register('inplace_gemm_1', OpKeyOptimizer(GemmLocalOptimizer(), 
-    failure_callback=OpKeyOptimizer.warn_inplace), 70.01, 'fast_run', 'inplace', 'gemm')
-compile.optdb.register('inplace_gemm_2', OpKeyOptimizer(GemmLocalOptimizer(), 
-    failure_callback=OpKeyOptimizer.warn_inplace), 70.02, 'fast_run', 'inplace', 'gemm')
 
 class Dot22(GemmRelated):
     """Compute a matrix-matrix product.
     This is a specialization of the more general Dot()
     """
     def make_node(self, x, y):
-        assert GemmLocalOptimizer._is_real_matrix(x)
+        assert _is_real_matrix(x)
         assert y.type == x.type               #makes sure y is a matrix
         bz = [False, False]
         outputs = [T.tensor(x.type.dtype, bz)]
@@ -638,7 +634,7 @@ _dot22 = Dot22()
 def local_dot_to_dot22(node):
     if node.op == T.dot:
         x,y = node.inputs
-        if GemmLocalOptimizer._is_real_matrix(x) and y.type == x.type:
+        if _is_real_matrix(x) and y.type == x.type:
             return [_dot22(*node.inputs)]
     else:
         return False

theano/tensor/elemwise.py

@@ -624,7 +624,7 @@ class Elemwise(Op):
 
         task_code = self.scalar_op.c_code(Apply(self.scalar_op,
                                                 [Scalar(dtype = input.type.dtype)() for input in node.inputs],
-                                                [Scalar(dtype = output.type.dtype)() for input in node.outputs]),
+                                                [Scalar(dtype = output.type.dtype)() for output in node.outputs]),
                                           name + '_scalar_',
                                           ["%s_i" % s for s in _inames],
                                           ["%s_i" % s for s in onames],
@@ -847,11 +847,12 @@ class Sum(CAReduce):
         CAReduce.__init__(self, scalar.add, axis)
 
     def _output_dtype(self, idtype):
-        if idtype.startswith('int'):
-            return 'int64' #we want to protect against overflow
-        else:
-            return idtype
-
+        # we want to protect against overflow
+        return dict(
+                int8='int32',
+                int16='int32',
+                int32='int64',
+                ).get(idtype, idtype)
 
     def grad(self, (x, ), (gz, )):
         gz = as_tensor_variable(gz)

theano/tensor/nnet.py

@@ -133,6 +133,8 @@ class SoftmaxWithBias(gof.Op):
     def c_headers(self):
         return ['<iostream>','<cmath>']
 
+    def c_code_cache_version(self):
+        return ()
     @staticmethod
     def c_code_template():
         # this implementation was lifted from
@@ -157,14 +159,14 @@ class SoftmaxWithBias(gof.Op):
             PyErr_SetString(PyExc_ValueError, "b not 1d tensor");
             %(fail)s;
         }
-        if (%(x)s->descr->type_num != PyArray_DOUBLE)
+        if ((%(x)s->descr->type_num != PyArray_DOUBLE)&&(%(x)s->descr->type_num != PyArray_DOUBLE))
         {
-            PyErr_SetString(PyExc_TypeError, "a not float64");
+            PyErr_SetString(PyExc_TypeError, "a not float");
             %(fail)s;
         }
-        if (%(b)s->descr->type_num != PyArray_DOUBLE)
+        if ((%(b)s->descr->type_num != PyArray_DOUBLE) && (%(b)s->descr->type_num != PyArray_DOUBLE))
         {
-            PyErr_SetString(PyExc_TypeError, "b not float64");
+            PyErr_SetString(PyExc_TypeError, "b not float");
             %(fail)s;
         }
         if ((%(x)s->dimensions[1] != %(b)s->dimensions[0]))
@@ -193,22 +195,22 @@ class SoftmaxWithBias(gof.Op):
             double sum = 0.0;
             bool  discount_max = false;
 
-            const double* __restrict__ x_i = (double*)(%(x)s->data + %(x)s->strides[0] * i);
-            const double* __restrict__ b_i = (double*)(%(b)s->data);
-            double* __restrict__ sm_i = (double*)(%(sm)s->data + %(sm)s->strides[0] * i);
+            const REAL* __restrict__ x_i = (REAL*)(%(x)s->data + %(x)s->strides[0] * i);
+            const REAL* __restrict__ b_i = (REAL*)(%(b)s->data);
+            REAL* __restrict__ sm_i = (REAL*)(%(sm)s->data + %(sm)s->strides[0] * i);
         """
 
         inside_row_loop = """
-            npy_intp Sx = %(x)s->strides[1]/sizeof(double);
-            npy_intp Sb = %(b)s->strides[0]/sizeof(double);
-            npy_intp Ssm = %(sm)s->strides[1]/sizeof(double);
+            npy_intp Sx = %(x)s->strides[1]/sizeof(REAL);
+            npy_intp Sb = %(b)s->strides[0]/sizeof(REAL);
+            npy_intp Ssm = %(sm)s->strides[1]/sizeof(REAL);
 
             size_t row_max_j=0;
-            double row_max = x_i[0] + b_i[0];
+            REAL row_max = x_i[0] + b_i[0];
             // Get the maximum value of the row
             for (j = 0; j < Nx[1]; ++j)
             {
-                double row_ij = x_i[j * Sx] +  b_i[j * Sb];
+                REAL row_ij = x_i[j * Sx] +  b_i[j * Sb];
 //                std::cout << "1" << row_ij << "\\n";
                 row_max_j = (row_ij > row_max) ? j : row_max_j;
                 row_max   = (row_ij > row_max) ? row_ij : row_max;
@@ -216,9 +218,9 @@ class SoftmaxWithBias(gof.Op):
 
             for (j = 0; j < Nx[1]; ++j)
             {
-                double row_ij = x_i[j * Sx] +  b_i[j * Sb];
+                REAL row_ij = x_i[j * Sx] +  b_i[j * Sb];
 //                std::cout << "2" << row_ij << "\\n";
-                double sm_ij = exp(row_ij - row_max);
+                REAL sm_ij = exp(row_ij - row_max);
 //                std::cout << "3" << sm_ij << "\\n";
                 sum += sm_ij;
                 sm_i[j * Ssm] = sm_ij;
@@ -292,12 +294,24 @@ class SoftmaxGrad(gof.Op):
     def grad(self, *args):
         raise NotImplementedError()
 
+    def c_code_cache_version(self):
+        return ()
     def c_code(self, node, name, (dy, sm), (dx,), sub):
+        if node.inputs[1].type.dtype != node.inputs[0].type.dtype:
+            raise NotImplementedError()
+        if node.inputs[0].type.dtype == 'float32':
+            REAL = 'float'
+        else:
+            REAL = 'double'
         return '''
-        if ((%(dy)s->descr->type_num != PyArray_DOUBLE)
-            || (%(sm)s->descr->type_num != PyArray_DOUBLE))
+        if ((%(dy)s->descr->type_num != PyArray_DOUBLE) && ((%(dy)s->descr->type_num != PyArray_FLOAT)
+        {
+            PyErr_SetString(PyExc_TypeError, "types should be float or float64");
+            %(fail)s;
+        }
+        if ((%(sm)s->descr->type_num != PyArray_DOUBLE) && ((%(sm)s->descr->type_num != PyArray_FLOAT)
         {
-            PyErr_SetString(PyExc_TypeError, "types should be float64, float64");
+            PyErr_SetString(PyExc_TypeError, "types should be float or float64");
             %(fail)s;
         }
         if ((%(dy)s->nd != 2)
@@ -327,12 +341,12 @@ class SoftmaxGrad(gof.Op):
 
         for (size_t i = 0; i < %(dx)s->dimensions[0]; ++i)
         {
-            const double* __restrict__ dy_i = (double*) (%(dy)s->data + %(dy)s->strides[0] * i);
-            npy_intp Sdy = %(dy)s->strides[1]/sizeof(double);
-            const double* __restrict__ sm_i = (double*) (%(sm)s->data + %(sm)s->strides[0] * i);
-            npy_intp Ssm = %(sm)s->strides[1]/sizeof(double);
-            double* __restrict__ dx_i = (double*) (%(dx)s->data + %(dx)s->strides[0] * i);
-            npy_intp Sdx = %(dx)s->strides[1]/sizeof(double);
+            const REAL* __restrict__ dy_i = (REAL*) (%(dy)s->data + %(dy)s->strides[0] * i);
+            npy_intp Sdy = %(dy)s->strides[1]/sizeof(REAL);
+            const REAL* __restrict__ sm_i = (REAL*) (%(sm)s->data + %(sm)s->strides[0] * i);
+            npy_intp Ssm = %(sm)s->strides[1]/sizeof(REAL);
+            REAL* __restrict__ dx_i = (REAL*) (%(dx)s->data + %(dx)s->strides[0] * i);
+            npy_intp Sdx = %(dx)s->strides[1]/sizeof(REAL);
 
             double sum_dy_times_sm = 0.;
             for (size_t j = 0; j < %(dx)s->dimensions[1]; ++j)
@@ -587,7 +601,7 @@ class CrossentropySoftmaxArgmax1HotWithBias(gof.Op):
                 begin_row_loop,
                 """
             const %(y_idx_type)s y_i = ((%(y_idx_type)s*)(%(y_idx)s->data + %(y_idx)s->strides[0] * i))[0];
-            double* __restrict__ nll_i = (double*)(%(nll)s->data + %(nll)s->strides[0] * i);
+            REAL* __restrict__ nll_i = (REAL*)(%(nll)s->data + %(nll)s->strides[0] * i);
             %(am_type)s* __restrict__ am_i = (%(am_type)s*) (%(am)s->data + %(am)s->strides[0] * i);
                 """,
                 inside_row_loop,
@@ -610,6 +624,10 @@ class CrossentropySoftmaxArgmax1HotWithBias(gof.Op):
         y_idx_type = node.inputs[2].type.dtype_specs()[1]
         am_type = y_idx_type
         code_template = ''.join(self.c_code_template())
+        if node.inputs[0].type.dtype == 'float32':
+            REAL = 'float'
+        else:
+            REAL = 'double'
         return code_template % dict(locals(), **sub)
 
 class CrossentropySoftmax1HotWithBiasDx (gof.Op):
@@ -686,15 +704,15 @@ class CrossentropySoftmax1HotWithBiasDx (gof.Op):
 
         for (size_t i = 0; i < %(dx)s->dimensions[0]; ++i)
         {
-            const double dnll_i = ((double*)(%(dnll)s->data + %(dnll)s->strides[0] * i))[0];
+            const REAL dnll_i = ((REAL*)(%(dnll)s->data + %(dnll)s->strides[0] * i))[0];
 
             const %(y_idx_type)s y_i = ((%(y_idx_type)s*)(%(y_idx)s->data + %(y_idx)s->strides[0] * i))[0];
 
-            const double* __restrict__ sm_i = (double*)(%(sm)s->data + %(sm)s->strides[0] * i);
-            npy_intp Ssm = %(sm)s->strides[1]/sizeof(double);
+            const REAL* __restrict__ sm_i = (REAL*)(%(sm)s->data + %(sm)s->strides[0] * i);
+            npy_intp Ssm = %(sm)s->strides[1]/sizeof(REAL);
 
-            double* __restrict__ dx_i = (double*)(%(dx)s->data + %(dx)s->strides[0] * i);
-            npy_intp Sdx = %(dx)s->strides[1]/sizeof(double);
+            REAL* __restrict__ dx_i = (REAL*)(%(dx)s->data + %(dx)s->strides[0] * i);
+            npy_intp Sdx = %(dx)s->strides[1]/sizeof(REAL);
 
             for (size_t j = 0; j < %(dx)s->dimensions[1]; ++j)
             {
@@ -790,7 +808,8 @@ class CrossentropyCategorical1Hot(gof.Op):
                     '(got type: %s instead of: %s)' % (_true_one_of_n.type,
                         tensor.lvector))
 
-        return gof.Apply(self, [_coding_dist, _true_one_of_n], [tensor.dvector()])
+        return gof.Apply(self, [_coding_dist, _true_one_of_n],
+                [tensor.Tensor(dtype=_coding_dist.dtype, broadcastable=[False])()])
 
     def perform(self, node, (coding, one_of_n), (y_out,)):
         y = numpy.zeros_like(coding[:,0])

theano/tensor/opt.py

@@ -114,7 +114,7 @@ def local_dimshuffle_lift(node):
 
     input = node.inputs[0]
     inode = input.owner
-    if inode and isinstance(inode.op, Elemwise):
+    if inode and isinstance(inode.op, Elemwise) and (len(input.clients)==1):
         return inode.op.make_node(*[DimShuffle(input.type.broadcastable,
                                                op.new_order,
                                                op.inplace)(input) for input in inode.inputs]).outputs