Merge pull request #151 from JuliaDiff/setinde

Piracy-free setindex

Author: ChrisRackauckas

src/FiniteDiff.jl

@@ -12,6 +12,26 @@ _mat(x::AbstractMatrix) = x
 _mat(x::StaticVector)   = reshape(x, (axes(x, 1),     SOneTo(1)))
 _mat(x::AbstractVector) = reshape(x, (axes(x, 1), Base.OneTo(1)))
 
+# Setindex overloads without piracy
+setindex(x...) = Base.setindex(x...)
+setindex(x::StaticArray, v, i::Int...) = StaticArrays.setindex(x, v, i...)
+
+function setindex(x::AbstractArray, v, i...)
+    _x = Base.copymutable(x)
+    _x[i...] = v
+    return _x
+end
+
+function setindex(x::AbstractVector, v, i::Int)
+    n = length(x)
+    x .* (i .!== 1:n) .+ v .* (i .== 1:n)
+end
+
+function setindex(x::AbstractMatrix, v, i::Int, j::Int)
+    n, m = Base.size(x)
+    x .* (i .!== 1:n) .* (j .!== i:m)' .+ v .* (i .== 1:n) .* (j .== i:m)'
+end
+
 include("iteration_utils.jl")
 include("epsilons.jl")
 include("derivatives.jl")

src/hessians.jl

@@ -75,8 +75,8 @@ function finite_difference_hessian!(H,f,x,
             ArrayInterfaceCore.allowed_setindex!(xpp,xi + epsilon,i)
             ArrayInterfaceCore.allowed_setindex!(xmm,xi - epsilon,i)
         else
-            _xpp = Base.setindex(xpp,xi + epsilon, i)
-            _xmm = Base.setindex(xmm,xi - epsilon, i)
+            _xpp = setindex(xpp,xi + epsilon, i)
+            _xmm = setindex(xmm,xi - epsilon, i)
         end
 
         ArrayInterfaceCore.allowed_setindex!(H,(f(_xpp) - 2*fx + f(_xmm)) / epsilon^2,i,i)
@@ -90,10 +90,10 @@ function finite_difference_hessian!(H,f,x,
             ArrayInterfaceCore.allowed_setindex!(xmp,xm,i)
             ArrayInterfaceCore.allowed_setindex!(xmm,xm,i)
         else
-            _xpp = Base.setindex(xpp,xp,i)
-            _xpm = Base.setindex(xpm,xp,i)
-            _xmp = Base.setindex(xmp,xm,i)
-            _xmm = Base.setindex(xmm,xm,i)
+            _xpp = setindex(xpp,xp,i)
+            _xpm = setindex(xpm,xp,i)
+            _xmp = setindex(xmp,xm,i)
+            _xmm = setindex(xmm,xm,i)
         end
 
         for j = i+1:n
@@ -108,10 +108,10 @@ function finite_difference_hessian!(H,f,x,
                 ArrayInterfaceCore.allowed_setindex!(xmp,xp,j)
                 ArrayInterfaceCore.allowed_setindex!(xmm,xm,j)
             else
-                _xpp = Base.setindex(_xpp,xp,j)
-                _xpm = Base.setindex(_xpm,xm,j)
-                _xmp = Base.setindex(_xmp,xp,j)
-                _xmm = Base.setindex(_xmm,xm,j)
+                _xpp = setindex(_xpp,xp,j)
+                _xpm = setindex(_xpm,xm,j)
+                _xmp = setindex(_xmp,xp,j)
+                _xmm = setindex(_xmm,xm,j)
             end
 
             ArrayInterfaceCore.allowed_setindex!(H,(f(_xpp) - f(_xpm) - f(_xmp) + f(_xmm))/(4*epsiloni*epsilonj),i,j)
@@ -122,10 +122,10 @@ function finite_difference_hessian!(H,f,x,
                 ArrayInterfaceCore.allowed_setindex!(xmp,xj,j)
                 ArrayInterfaceCore.allowed_setindex!(xmm,xj,j)
             else
-                _xpp = Base.setindex(_xpp,xj,j)
-                _xpm = Base.setindex(_xpm,xj,j)
-                _xmp = Base.setindex(_xmp,xj,j)
-                _xmm = Base.setindex(_xmm,xj,j)
+                _xpp = setindex(_xpp,xj,j)
+                _xpm = setindex(_xpm,xj,j)
+                _xmp = setindex(_xmp,xj,j)
+                _xmm = setindex(_xmm,xj,j)
             end
         end
 

src/jacobians.jl

@@ -193,7 +193,7 @@ function finite_difference_jacobian(
         function calculate_Ji_forward(i)
             x_save = ArrayInterfaceCore.allowed_getindex(vecx, i)
             epsilon = compute_epsilon(Val(:forward), x_save, relstep, absstep, dir)
-            _vecx1 = Base.setindex(vecx, x_save+epsilon, i)
+            _vecx1 = setindex(vecx, x_save+epsilon, i)
             _x1 = reshape(_vecx1, axes(x))
             vecfx1 = _vec(f(_x1))
             dx = (vecfx1-vecfx) / epsilon
@@ -226,8 +226,8 @@ function finite_difference_jacobian(
             x1_save = ArrayInterfaceCore.allowed_getindex(vecx1,i)
             x_save = ArrayInterfaceCore.allowed_getindex(vecx,i)
             epsilon = compute_epsilon(Val(:forward), x1_save, relstep, absstep, dir)
-            _vecx1 = Base.setindex(vecx1,x1_save+epsilon,i)
-            _vecx = Base.setindex(vecx,x_save-epsilon,i)
+            _vecx1 = setindex(vecx1,x1_save+epsilon,i)
+            _vecx = setindex(vecx,x_save-epsilon,i)
             _x1 = reshape(_vecx1, axes(x))
             _x = reshape(_vecx, axes(x))
             vecfx1 = _vec(f(_x1))
@@ -264,7 +264,7 @@ function finite_difference_jacobian(
 
         function calculate_Ji_complex(i)
             x_save = ArrayInterfaceCore.allowed_getindex(vecx,i)
-            _vecx = Base.setindex(complex.(vecx),x_save+im*epsilon,i)
+            _vecx = setindex(complex.(vecx),x_save+im*epsilon,i)
             _x = reshape(_vecx, axes(x))
             vecfx = _vec(f(_x))
             dx = imag(vecfx)/epsilon
@@ -325,18 +325,18 @@ function _findstructralnz(A::DenseMatrix)
     numnz = count(A .≠ 0)
     I = Vector{Int64}(undef, numnz)
     J = Vector{Int64}(undef, numnz)
-    idx = 1 
+    idx = 1
     for j in axes(A, 2)
         for i in axes(A, 1)
             if A[i, j] ≠ 0
-                I[idx] = i 
+                I[idx] = i
                 J[idx] = j
-                idx += 1 
+                idx += 1
             end
         end
     end
     I, J
-end 
+end
 
 function finite_difference_jacobian!(
     J,
@@ -361,7 +361,7 @@ function finite_difference_jacobian!(
     cols_index = nothing
     if _use_findstructralnz(sparsity)
         rows_index, cols_index = ArrayInterfaceCore.findstructralnz(sparsity)
-    elseif sparsity isa DenseMatrix 
+    elseif sparsity isa DenseMatrix
         rows_index, cols_index = FiniteDiff._findstructralnz(sparsity)
     end