Initial version. Very dirty.

Author: sotte

.gitignore

@@ -0,0 +1,55 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+env/
+bin/
+build/
+develop-eggs/
+dist/
+eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+.tox/
+.coverage
+.cache
+nosetests.xml
+coverage.xml
+
+# Translations
+*.mo
+
+# Mr Developer
+.mr.developer.cfg
+.project
+.pydevproject
+
+# Rope
+.ropeproject
+
+# Django stuff:
+*.log
+*.pot
+
+# Sphinx documentation
+docs/_build/
+
+# ipython notebook
+.ipynb_checkpoints

README.md

@@ -0,0 +1,33 @@
+# Tired of numpy.dot?
+
+Note: this is pre alpha!
+
+TODO change that!
+
+TODO Ideally this thould be integrated in numpy.
+
+`mdot` chains multiplication calls and allows you to write
+```python
+mdot(A, B, C, D)
+```
+instead of
+```python
+np.dot(np.dot(np.dot(A, B), C), D)
+A.dot(B).dot(C).dot(D)
+```
+
+Did I mention that it automatically speeds up the multiplication by setting the
+parens in an optimal fashion:
+```python
+>>> %timeit np.dot(np.dot(np.dot(A, B), C), D)
+1 loops, best of 3: 694 ms per loop
+>>> %timeit mdot(A, B, C, D)
+100 loops, best of 3: 5.18 ms per loop
+```
+
+Still, not satisfied? Get red rid of the overhead for calculating the optimal
+parens once and then use the expression:
+```python
+>>> print_optimal(D, A, B, C, names=list("DABC"))
+"np.dot(np.dot(D, np.dot(A, B)), C)"
+```

mdot.py

@@ -0,0 +1,100 @@
+from __future__ import print_function
+
+import numpy as np
+import mytimer
+
+
+#@mytimer.timeit
+def chain_order_rec(args):
+    """
+    cost[i, k ] = min([cost[prefix] + cost[suffix] + cost_mult(prefix, suffix)
+                       for k in range(i, j)])
+    m[i, k ] = min([m[i, k] + m[k+1, j] + p[i] * p[k+1] * p[j+1]
+                    for k in range(i, j)])
+
+    """
+    # p is the list of the row length of all matrices plus the column of the
+    # last matrix
+    # example
+    # A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50]
+    # The cost for multipying AB is then: 10 * 100 * 5
+    p = [arg.shape[0] for arg in args]
+    p.append(args[-1].shape[1])
+
+    # determine the order of the multiplication using DP
+    n = len(p) - 1
+    # costs for subproblems
+    m = np.zeros((n, n))
+    # helper to actually multiply optimal solution
+    s = np.zeros((n, n))
+    for i in range(n):
+        for j in range(i+1, n):
+            cost, k = min((m[i, k] + m[k+1, j] + p[i] * p[k+1] * p[j+1], k)
+                          for k in range(i, j))
+            m[i, j] = cost
+            s[i, j] = k
+
+    return m, s
+
+
+#@mytimer.timeit
+def multiply_r(args, s, i, j):
+    if i == j:
+        return args[int(i)]
+    else:
+        return np.dot(multiply_r(args, s, i, s[i, j]),
+                      multiply_r(args, s, s[i, j] + 1, j))
+
+
+def _print_parens(args, s, i, j, names=None):
+    if i == j:
+        if names:
+            print(names[int(i)], end="")
+        else:
+            str_ = "M_{}".format(int(i))
+            print(str_, end="")
+    else:
+        print("np.dot(", end="")
+        _print_parens(args, s, i, s[i, j], names)
+        print(", ", end="")
+        _print_parens(args, s, s[i, j] + 1, j, names)
+        print(")", end="")
+
+
+def print_optimal(*args, **kwargs):
+    """Print the optimal chain of multiplications that minimizes the total
+    number of multiplications.
+
+    """
+    names = kwargs.get("names", None)
+    m, s = chain_order_rec(args)
+    _print_parens(args, s, 0, len(args) - 1, names=names)
+
+
+def mdot(*args, **kwargs):
+    """Multiply the given arrays.
+
+    `optimize` = True
+
+    TODO extend and document.
+
+    Minimize the number of required scalar multiplications for the given
+    matrices.
+
+    Example for the costs:
+    A_{10x100}, B_{100x5}, C_{5x50}
+
+    cost((AB)C) = 5000 + 2500 = 7500
+    cost(A(BC)) = 50000 + 25000 = 75000
+
+    """
+    if len(args) == 1:
+        return args[0]
+
+    optimize = kwargs.get("optimize", True)
+
+    if optimize:
+        m, s = chain_order_rec(args)
+        return multiply_r(args, s, 0, len(args) - 1)
+    else:
+        return reduce(np.dot, args)

mytimer.py

@@ -0,0 +1,33 @@
+import time
+
+
+def timeit(func=None, loops=1, verbose=False):
+    if func is not None:
+        def inner(*args, **kwargs):
+
+            sums = 0.0
+            mins = 1.7976931348623157e+308
+            maxs = 0.0
+            print '====%s Timing====' % func.__name__
+            for i in range(0, loops):
+                t0 = time.time()
+                result = func(*args, **kwargs)
+                dt = time.time() - t0
+                mins = dt if dt < mins else mins
+                maxs = dt if dt > maxs else maxs
+                sums += dt
+                if verbose is True:
+                    print "\t%r ran in %2.9f sec on run %s" % (
+                        func.__name__, dt, i)
+            print "%r min run time was %2.9f sec" % (func.__name__, mins)
+            print "%r max run time was %2.9f sec" % (func.__name__, maxs)
+            print "%r avg run time was %2.9f sec in %s runs" % (
+                func.__name__, sums/loops, loops)
+            print "==== end ===="
+            return result
+
+        return inner
+    else:
+        def partial_inner(func):
+            return timeit(func, loops, verbose)
+        return partial_inner

test_mdot.py

@@ -0,0 +1,35 @@
+from mdot import mdot
+import numpy as np
+
+
+##############################################################################
+def test_unoptimized_one_parameter():
+    I = np.eye(3, 3)
+    assert (mdot(I, optimize=False) == I).all()
+
+
+def test_unoptimized_multiple_parameters():
+    I = np.eye(3, 3)
+    assert (mdot(I, I, optimize=False) == I).all()
+    assert (mdot(I, I, I, optimize=False) == I).all()
+
+
+def test_unoptimized_fancy():
+    A = np.random.random((3, 3))
+    B = np.linalg.inv(A)
+    I = np.eye(3)
+    assert np.allclose(mdot(A, B, optimize=False), I)
+
+
+##############################################################################
+def test_optimized_general():
+    I = np.eye(3, 3)
+    assert np.allclose(mdot(I, I, I, I, I, I, I, I, optimize=True), I)
+
+
+def test_optimized_fancy():
+    A = np.random.random((3, 3))
+    B = np.linalg.inv(A)
+    I = np.eye(3)
+    assert np.allclose(mdot(A, B, I), I)
+    print mdot(A, B, I, optimize=True).shape