Do gradient via mutating and unmutating cell

Author: oxinabox

src/grad.jl

@@ -8,19 +8,25 @@ Approximate the gradient of `f` at `xs...` using `fdm`. Assumes that `f(xs...)`
 """
 function grad end
 
-function grad(fdm, f, x::AbstractArray{T}) where T <: Number
+function _grad(fdm, f, x::AbstractArray{T}) where T <: Number
+    # x must be mutable, we will mutate it and then mutate it back.
     dx = similar(x)
-    tmp = similar(x)
     for k in eachindex(x)
         dx[k] = fdm(zero(T)) do ϵ
-            tmp .= x
-            tmp[k] += ϵ
-            return f(tmp)
+            xk = x[k]
+            x[k] = xk +  ϵ
+            ret = f(x)
+            x[k] = xk  # Can't do `x[k] -= ϵ` as floating-point math is not associative
+            return ret
         end
     end
     return (dx, )
 end
 
+grad(fdm, f, x::Array{<:Number}) = _grad(fdm, f, x)
+# Fallback for when we don't know `x` will be mutable:
+grad(fdm, f, x::AbstractArray{<:Number}) = _grad(fdm, f, similar(x).=x)
+
 grad(fdm, f, x::Real) = (fdm(f, x), )
 grad(fdm, f, x::Tuple) = (grad(fdm, (xs...)->f(xs), x...), )