Merge pull request odlgroup#1225 from kohr-h/issue-1181__op_aliasing_bug

BUG: fix in-out aliasing issue in OperatorSum, closes odlgroup#1181

odl/operator/operator.py

@@ -1154,8 +1154,10 @@ def _call(self, x, out=None):
         else:
             tmp = (self.__tmp_ran if self.__tmp_ran is not None
                    else self.range.element())
-            self.left(x, out=out)
-            self.right(x, out=tmp)
+            # Write to `tmp` first, otherwise aliased `x` and `out` lead
+            # to wrong result
+            self.left(x, out=tmp)
+            self.right(x, out=out)
             out += tmp
 
     def derivative(self, x):
@@ -1492,8 +1494,10 @@ def _call(self, x, out=None):
             return self.left(x) * self.right(x)
         else:
             tmp = self.right.range.element()
-            self.left(x, out=out)
-            self.right(x, out=tmp)
+            # Write to `tmp` first, otherwise aliased `x` and `out` lead
+            # to wrong result
+            self.left(x, out=tmp)
+            self.right(x, out=out)
             out *= tmp
 
     def derivative(self, x):

odl/operator/tensor_ops.py

@@ -858,6 +858,8 @@ def inverse(self):
 
     def _call(self, x, out=None):
         """Return ``self(x[, out])``."""
+        from pkg_resources import parse_version
+
         if out is None:
             return self.range.element(self.matrix.dot(x))
         else:
@@ -866,8 +868,14 @@ def _call(self, x, out=None):
                 # sparse matrices
                 out[:] = self.matrix.dot(x)
             else:
-                with writable_array(out) as out_arr:
-                    self.matrix.dot(x, out=out_arr)
+                if (parse_version(np.__version__) < parse_version('1.13.0') and
+                        x is out):
+                    # Workaround for bug in Numpy < 1.13 with aliased in and
+                    # out in np.dot
+                    out[:] = self.matrix.dot(x)
+                else:
+                    with writable_array(out) as out_arr:
+                        self.matrix.dot(x, out=out_arr)
 
     def __repr__(self):
         """Return ``repr(self)``."""