Merge pull request #351 from opesci/new-tti-again

New tti again

devito/__init__.py

@@ -11,6 +11,7 @@
 from devito.logger import error, warning, info  # noqa
 from devito.parameters import (configuration, init_configuration,  # noqa
                                env_vars_mapper)
+from devito.tools import *  # noqa
 
 
 def clear_cache():

devito/dle/blocking_utils.py

@@ -40,9 +40,13 @@ def fold_blockable_tree(node, exclude_innermost=False):
             pairwise_folds = list(zip(*reversed(trees)))
             if any(not is_foldable(j) for j in pairwise_folds):
                 continue
-            # Perform folding
+            # Maybe heuristically exclude innermost Iteration
             if exclude_innermost is True:
                 pairwise_folds = pairwise_folds[:-1]
+            # Perhaps there's nothing to fold
+            if len(pairwise_folds) == 1:
+                continue
+            # Perform folding
             for j in pairwise_folds:
                 root, remainder = j[0], j[1:]
                 folds = [(tuple(y-x for x, y in zip(i.offsets, root.offsets)), i.nodes)

devito/dse/backends/advanced.py

@@ -9,7 +9,7 @@
 from devito.dse.clusterizer import clusterize
 from devito.dse.inspection import estimate_cost
 from devito.dse.manipulation import (common_subexprs_elimination, collect_nested,
-                                     xreplace_constrained)
+                                     xreplace_constrained, compact_temporaries)
 from devito.dse.queries import iq_timeinvariant
 from devito.interfaces import Indexed, ScalarFunction, TensorFunction
 
@@ -34,14 +34,18 @@ def _extract_time_invariants(self, cluster, template, with_cse=True,
         rule = iq_timeinvariant(cluster.trace)
         costmodel = costmodel or (lambda e: estimate_cost(e) > 0)
         processed, found = xreplace_constrained(cluster.exprs, make, rule, costmodel)
-        leaves = [i for i in processed if i not in found]
 
-        # Search for common sub-expressions amongst them (and only them)
         if with_cse:
+            leaves = [i for i in processed if i not in found]
+
+            # Search for common sub-expressions amongst them (and only them)
             make = lambda i: ScalarFunction(name=template(i + len(found))).indexify()
             found = common_subexprs_elimination(found, make)
 
-        return cluster.reschedule(found + leaves)
+            # Some temporaries may be droppable at this point
+            processed = compact_temporaries(found + leaves)
+
+        return cluster.reschedule(processed)
 
     @dse_pass
     def _factorize(self, cluster, *args, **kwargs):

devito/dse/clusterizer.py

@@ -105,7 +105,8 @@ def optimize(clusters):
         mapper = {}
         temporaries = []
         for k, v in c1.trace.items():
-            if v.is_tensor and not any(v.function in c2.unknown for c2 in clusters):
+            if v.function.is_TensorFunction and\
+                    not any(v.function in c2.unknown for c2 in clusters):
                 for i in c1.tensors[v.function]:
                     # LHS scalarization
                     scalarized = ScalarFunction(name='s%d' % len(mapper)).indexify()

devito/dse/graph.py

@@ -92,10 +92,6 @@ def is_tensor(self):
     def is_scalar(self):
         return not self.is_tensor
 
-    @property
-    def is_dead(self):
-        return self.is_scalar and self.is_terminal and len(self.reads) == 1
-
     def construct(self, rule):
         """
         Create a new temporary starting from ``self`` replacing symbols in

devito/dse/manipulation.py

@@ -10,7 +10,7 @@
 from devito.dse.extended_sympy import Add, Eq, Mul
 from devito.dse.inspection import count, estimate_cost, retrieve_indexed
 from devito.dse.graph import temporaries_graph
-from devito.dse.queries import q_indexed, q_op
+from devito.dse.queries import q_indexed, q_op, q_leaf
 from devito.interfaces import Indexed, TensorFunction
 from devito.tools import as_tuple
 
@@ -105,7 +105,7 @@ def run(expr):
     return run(expr)[0]
 
 
-def xreplace_constrained(exprs, make, rule, costmodel=lambda e: True, repeat=False):
+def xreplace_constrained(exprs, make, rule=None, costmodel=lambda e: True, repeat=False):
     """
     Unlike ``xreplace``, which replaces all objects specified in a mapper,
     this function replaces all objects satisfying two criteria: ::
@@ -122,28 +122,38 @@ def xreplace_constrained(exprs, make, rule, costmodel=lambda e: True, repeat=Fal
     ``a + b`` satisfies the cost model.
 
     :param exprs: The target SymPy expression, or a collection of SymPy expressions.
-    :param make: A function to construct symbols used for replacement.
-                 The function takes as input an integer ID; ID is computed internally
-                 and used as a unique identifier for the constructed symbols.
-    :param rule: The matching rule (a lambda function).
+    :param make: Either a mapper M: K -> V, indicating how to replace an expression
+                 in K with a symbol in V, or a function, used to construct new, unique
+                 symbols. Such a function should take as input a parameter, used to
+                 enumerate the new symbols.
+    :param rule: The matching rule (a lambda function). May be left unspecified if
+                 ``make`` is a mapper.
     :param costmodel: The cost model (a lambda function, optional).
     :param repeat: Repeatedly apply ``xreplace`` until no more replacements are
                    possible (optional, defaults to False).
     """
-
     found = OrderedDict()
     rebuilt = []
 
-    def replace(expr):
-        temporary = found.get(expr)
-        if temporary:
-            return temporary
-        else:
-            temporary = make(replace.c)
-            found[expr] = temporary
-            replace.c += 1
-            return temporary
-    replace.c = 0  # Unique identifier for new temporaries
+    # Define /replace()/ based on the user-provided /make/
+    if isinstance(make, dict):
+        rule = rule if rule is not None else (lambda i: i in make)
+        replace = lambda i: make[i]
+    else:
+        assert callable(make) and callable(rule)
+
+        def replace(expr):
+            if isinstance(make, dict):
+                return make[expr]
+            temporary = found.get(expr)
+            if temporary:
+                return temporary
+            else:
+                temporary = make(replace.c)
+                found[expr] = temporary
+                replace.c += 1
+                return temporary
+        replace.c = 0  # Unique identifier for new temporaries
 
     def run(expr):
         if expr.is_Atom or q_indexed(expr):
@@ -254,9 +264,6 @@ def common_subexprs_elimination(exprs, make, mode='default'):
     mapper = {i.lhs: j.lhs for i, j in zip(mapped, reversed(mapped))}
     processed = [e.xreplace(mapper) for e in processed]
 
-    # Some temporaries may be droppable at this point
-    processed = compact_temporaries(processed)
-
     return processed
 
 
@@ -266,14 +273,16 @@ def compact_temporaries(exprs):
     """
     g = temporaries_graph(exprs)
 
-    mapper = {list(v.reads)[0]: k for k, v in g.items() if v.is_dead}
+    mapper = {k: v.rhs for k, v in g.items()
+              if v.is_scalar and (q_leaf(v.rhs) or v.rhs.is_Function)}
 
     processed = []
     for k, v in g.items():
-        if k in mapper:
-            processed.append(Eq(mapper[k], v.rhs))
-        elif not v.is_dead:
-            processed.append(v.xreplace(mapper))
+        if k not in mapper:
+            # The temporary /v/ is retained, and substitutions may be applied
+            handle, _ = xreplace_constrained(v, mapper, repeat=True)
+            assert len(handle) == 1
+            processed.extend(handle)
 
     return processed
 

devito/tools.py

@@ -1,6 +1,7 @@
 import os
 import ctypes
-from collections import Callable, Iterable, OrderedDict
+from collections import Callable, Iterable, OrderedDict, Hashable
+from functools import partial
 try:
     from itertools import izip_longest as zip_longest
 except ImportError:
@@ -9,6 +10,8 @@
 
 import numpy as np
 
+__all__ = ['memoized']
+
 
 def as_tuple(item, type=None, length=None):
     """
@@ -210,3 +213,39 @@ def __enter__(self):
 
     def __exit__(self, etype, value, traceback):
         os.chdir(self.savedPath)
+
+
+class memoized(object):
+    """
+    Decorator. Caches a function's return value each time it is called.
+    If called later with the same arguments, the cached value is returned
+    (not reevaluated).
+
+    Adapted from: ::
+
+        https://wiki.python.org/moin/PythonDecoratorLibrary#Memoize
+    """
+
+    def __init__(self, func):
+        self.func = func
+        self.cache = {}
+
+    def __call__(self, *args):
+        if not isinstance(args, Hashable):
+            # Uncacheable, a list, for instance.
+            # Better to not cache than blow up.
+            return self.func(*args)
+        if args in self.cache:
+            return self.cache[args]
+        else:
+            value = self.func(*args)
+            self.cache[args] = value
+            return value
+
+    def __repr__(self):
+        """Return the function's docstring."""
+        return self.func.__doc__
+
+    def __get__(self, obj, objtype):
+        """Support instance methods."""
+        return partial(self.__call__, obj)

examples/seismic/acoustic/acoustic_example.py

@@ -21,10 +21,10 @@ def smooth10(vel, shape):
     return out
 
 
-def acoustic_setup(dimensions=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
+def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
                    tn=500., time_order=2, space_order=4, nbpml=10, **kwargs):
-    nrec = dimensions[0]
-    model = demo_model('layers', ratio=3, shape=dimensions,
+    nrec = shape[0]
+    model = demo_model('layers-isotropic', ratio=3, shape=shape,
                        spacing=spacing, nbpml=nbpml)
 
     # Derive timestepping from model spacing
@@ -50,11 +50,11 @@ def acoustic_setup(dimensions=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
     return solver
 
 
-def run(dimensions=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
+def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
         time_order=2, space_order=4, nbpml=40, full_run=False,
         autotune=False, **kwargs):
 
-    solver = acoustic_setup(dimensions=dimensions, spacing=spacing,
+    solver = acoustic_setup(dimensions=shape, spacing=spacing,
                             nbpml=nbpml, tn=tn, space_order=space_order,
                             time_order=time_order, **kwargs)
 
@@ -89,22 +89,25 @@ def run(dimensions=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
                         type=int, help="Space order of the simulation")
     parser.add_argument("--nbpml", default=40,
                         type=int, help="Number of PML layers around the domain")
-    parser.add_argument("-dse", default='advanced',
-                        type=str, help="DSE backend choice")
-    parser.add_argument("-dle", default='advanced',
-                        type=str, help="DLE backend choice")
+    parser.add_argument("-dse", "-dse", default="advanced",
+                        choices=["noop", "basic", "advanced",
+                                 "speculative", "aggressive"],
+                        help="Devito symbolic engine (DSE) mode")
+    parser.add_argument("-dle", default="advanced",
+                        choices=["noop", "advanced", "speculative"],
+                        help="Devito loop engine (DSE) mode")
     args = parser.parse_args()
 
     # 3D preset parameters
     if args.dim2:
-        dimensions = (150, 150)
+        shape = (150, 150)
         spacing = (15.0, 15.0)
         tn = 750.0
     else:
-        dimensions = (50, 50, 50)
+        shape = (50, 50, 50)
         spacing = (20.0, 20.0, 20.0)
         tn = 250.0
 
-    run(dimensions=dimensions, spacing=spacing, nbpml=args.nbpml, tn=tn,
+    run(shape=shape, spacing=spacing, nbpml=args.nbpml, tn=tn,
         space_order=args.space_order, time_order=args.time_order,
         autotune=args.autotune, dse=args.dse, dle=args.dle, full_run=args.full)

examples/seismic/benchmark.py

@@ -30,7 +30,7 @@
                         choices=["acoustic", "tti"], help="Problem")
 
     simulation = parser.add_argument_group("Simulation")
-    simulation.add_argument("-d", "--dimensions", nargs=3, default=[50, 50, 50],
+    simulation.add_argument("-d", "--shape", nargs=3, default=[50, 50, 50],
                             type=int, help="Number of grid points along each axis",
                             metavar=("dim1", "dim2", "dim3"))
     simulation.add_argument("-s", "--spacing", nargs=3, default=[20.0, 20.0, 20.0],
@@ -91,7 +91,7 @@
     del parameters["max_flops"]
     del parameters["point_runtime"]
 
-    parameters["dimensions"] = tuple(parameters["dimensions"])
+    parameters["shape"] = tuple(parameters["shape"])
     parameters["spacing"] = tuple(parameters["spacing"])
 
     if args.execmode == "run":