Skip to content

Commit ae07e58

Browse files
ChrisRackauckasChrisRackauckas
committed
Use ArrayInterfaceCore
1 parent 0b90bf4 commit ae07e58

8 files changed

Lines changed: 93 additions & 82 deletions

File tree

‎Project.toml‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -3,24 +3,25 @@ uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
33
version = "2.11.1"
44

55
[deps]
6-
ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
6+
ArrayInterfaceCore = "30b0a656-2188-435a-8636-2ec0e6a096e2"
77
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
88
Requires = "ae029012-a4dd-5104-9daa-d747884805df"
99
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
1010
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
1111

1212
[compat]
13-
ArrayInterface = "1.1, 2.0, 3.0, 4, 5"
13+
ArrayInterfaceCore = "0.1.1"
1414
Requires = "0.5, 1.0"
1515
StaticArrays = "0.10, 0.11, 0.12, 1.0"
16-
julia = "1.2"
16+
julia = "1.6"
1717

1818
[extras]
19+
ArrayInterfaceBlockBandedMatrices = "5331f1e9-51c7-46b0-a9b0-df4434785e0a"
1920
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
2021
BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
2122
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
2223
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
2324
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
2425

2526
[targets]
26-
test = ["Test", "BlockBandedMatrices", "BandedMatrices", "Pkg", "SafeTestsets"]
27+
test = ["Test", "ArrayInterfaceBlockBandedMatrices", "BlockBandedMatrices", "BandedMatrices", "Pkg", "SafeTestsets"]

‎README.md‎

Lines changed: 14 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -180,7 +180,7 @@ fcalls #4
180180
Handling dense matrices? Easy. Handling sparse matrices? Cool stuff. Automatically
181181
specializing on the exact structure of a matrix? Even better. FiniteDiff can
182182
specialize on types which implement the
183-
[ArrayInterface.jl](https://github.com/JuliaDiffEq/ArrayInterface.jl) interface.
183+
[ArrayInterfaceCore.jl](https://github.com/JuliaDiffEq/ArrayInterfaceCore.jl) interface.
184184
This includes:
185185

186186
- Diagonal
@@ -194,7 +194,7 @@ Our previous example had a Tridiagonal Jacobian, so let's use this. If we just
194194
do
195195

196196
```julia
197-
using ArrayInterface, LinearAlgebra
197+
using ArrayInterfaceCore, LinearAlgebra
198198
tridiagjac = Tridiagonal(output)
199199
colors = matrix_colors(jac)
200200
```
@@ -235,7 +235,7 @@ let's specialize the Jacobian calculation on this fact:
235235
```julia
236236
using FillArrays, BlockBandedMatrices
237237
Jbbb = BandedBlockBandedMatrix(Ones(10000, 10000), fill(100, 100), fill(100, 100), (1, 1), (1, 1))
238-
colorsbbb = ArrayInterface.matrix_colors(Jbbb)
238+
colorsbbb = ArrayInterfaceCore.matrix_colors(Jbbb)
239239
bbbcache = FiniteDiff.JacobianCache(x,colorvec=colorsbbb,sparsity=Jbbb)
240240
FiniteDiff.finite_difference_jacobian!(Jbbb, pde, x, bbbcache)
241241
```
@@ -403,7 +403,7 @@ jacobian will return a similar type as `jac_prototype` if it is not a `nothing`.
403403
Jacobians, a cache which specifies the vector `fx` is required.
404404

405405
For sparse differentiation, pass a `colorvec` of matrix colors. `sparsity` should be a sparse
406-
or structured matrix (`Tridiagonal`, `Banded`, etc. according to the ArrayInterface.jl specs)
406+
or structured matrix (`Tridiagonal`, `Banded`, etc. according to the ArrayInterfaceCore.jl specs)
407407
to allow for decompression, otherwise the result will be the colorvec compressed Jacobian.
408408

409409
### Differencing Calls
@@ -430,7 +430,7 @@ finite_difference_jacobian!(J::AbstractMatrix,
430430
relstep=default_relstep(fdtype, eltype(x)),
431431
absstep=relstep,
432432
colorvec = 1:length(x),
433-
sparsity = ArrayInterface.has_sparsestruct(J) ? J : nothing)
433+
sparsity = ArrayInterfaceCore.has_sparsestruct(J) ? J : nothing)
434434

435435
# Cached
436436
FiniteDiff.finite_difference_jacobian(
@@ -528,3 +528,12 @@ HessianCache(xpp,xpm,xmp,xmm,
528528
fdtype::Type{T1}=Val{:hcentral},
529529
inplace::Type{Val{T2}} = x isa StaticArray ? Val{true} : Val{false})
530530
```
531+
532+
# Note about sparse differentiation of BandedMatrices and BlockBandedMatrices
533+
534+
These two matrix types need the dependencies ArrayInterfaceBandedMatrices.jl and
535+
ArrayInterfaceBlockBandedMatrices.jl to basically work with any functionality
536+
(anywhere). For now, the right thing to do is to add these libraries and do
537+
`import` on them if you are using BandedMatrices.jl or BlockBandedMatrices.jl
538+
for sparsity patterns. In the future, those two packages should just depend on
539+
ArrayInterface.jl and remove this issue entirely from the user space.

‎src/FiniteDiff.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
module FiniteDiff
22

33

4-
using LinearAlgebra, SparseArrays, StaticArrays, ArrayInterface, Requires
4+
using LinearAlgebra, SparseArrays, StaticArrays, ArrayInterfaceCore, Requires
55

66
import Base: resize!
77

‎src/gradients.jl‎

Lines changed: 26 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -136,20 +136,20 @@ function finite_difference_gradient!(
136136
# NOTE: in this case epsilon is a vector, we need two arrays for epsilon and x1
137137
# c1 denotes x1, c2 is epsilon
138138
fx, c1, c2, c3 = cache.fx, cache.c1, cache.c2, cache.c3
139-
if fdtype != Val(:complex) && ArrayInterface.fast_scalar_indexing(c2)
139+
if fdtype != Val(:complex) && ArrayInterfaceCore.fast_scalar_indexing(c2)
140140
@. c2 = compute_epsilon(fdtype, x, relstep, absstep, dir)
141141
copyto!(c1,x)
142142
end
143143
copyto!(c3,x)
144144
if fdtype == Val(:forward)
145145
@inbounds for i ∈ eachindex(x)
146-
if ArrayInterface.fast_scalar_indexing(c2)
147-
epsilon = ArrayInterface.allowed_getindex(c2,i)*dir
146+
if ArrayInterfaceCore.fast_scalar_indexing(c2)
147+
epsilon = ArrayInterfaceCore.allowed_getindex(c2,i)*dir
148148
else
149149
epsilon = compute_epsilon(fdtype, x, relstep, absstep, dir)*dir
150150
end
151-
c1_old = ArrayInterface.allowed_getindex(c1,i)
152-
ArrayInterface.allowed_setindex!(c1,c1_old + epsilon,i)
151+
c1_old = ArrayInterfaceCore.allowed_getindex(c1,i)
152+
ArrayInterfaceCore.allowed_setindex!(c1,c1_old + epsilon,i)
153153
if typeof(fx) != Nothing
154154
dfi = (f(c1) - fx) / epsilon
155155
else
@@ -158,51 +158,51 @@ function finite_difference_gradient!(
158158
end
159159
df_tmp = real(dfi)
160160
if eltype(df)<:Complex
161-
ArrayInterface.allowed_setindex!(c1,c1_old + im * epsilon,i)
161+
ArrayInterfaceCore.allowed_setindex!(c1,c1_old + im * epsilon,i)
162162
if typeof(fx) != Nothing
163163
dfi = (f(c1) - fx) / (im*epsilon)
164164
else
165165
dfi = (f(c1) - fx0) / (im*epsilon)
166166
end
167-
ArrayInterface.allowed_setindex!(c1,c1_old,i)
168-
ArrayInterface.allowed_setindex!(df, df_tmp - im * imag(dfi), i)
167+
ArrayInterfaceCore.allowed_setindex!(c1,c1_old,i)
168+
ArrayInterfaceCore.allowed_setindex!(df, df_tmp - im * imag(dfi), i)
169169
else
170-
ArrayInterface.allowed_setindex!(df, df_tmp, i)
171-
ArrayInterface.allowed_setindex!(c1,c1_old,i)
170+
ArrayInterfaceCore.allowed_setindex!(df, df_tmp, i)
171+
ArrayInterfaceCore.allowed_setindex!(c1,c1_old,i)
172172
end
173173
end
174174
elseif fdtype == Val(:central)
175175
@inbounds for i ∈ eachindex(x)
176-
if ArrayInterface.fast_scalar_indexing(c2)
177-
epsilon = ArrayInterface.allowed_getindex(c2,i)*dir
176+
if ArrayInterfaceCore.fast_scalar_indexing(c2)
177+
epsilon = ArrayInterfaceCore.allowed_getindex(c2,i)*dir
178178
else
179179
epsilon = compute_epsilon(fdtype, x, relstep, absstep, dir)*dir
180180
end
181-
c1_old = ArrayInterface.allowed_getindex(c1,i)
182-
ArrayInterface.allowed_setindex!(c1,c1_old + epsilon, i)
183-
x_old = ArrayInterface.allowed_getindex(x,i)
184-
ArrayInterface.allowed_setindex!(c3,x_old - epsilon,i)
181+
c1_old = ArrayInterfaceCore.allowed_getindex(c1,i)
182+
ArrayInterfaceCore.allowed_setindex!(c1,c1_old + epsilon, i)
183+
x_old = ArrayInterfaceCore.allowed_getindex(x,i)
184+
ArrayInterfaceCore.allowed_setindex!(c3,x_old - epsilon,i)
185185
df_tmp = real((f(c1) - f(c3)) / (2*epsilon))
186186
if eltype(df)<:Complex
187-
ArrayInterface.allowed_setindex!(c1,c1_old + im*epsilon,i)
188-
ArrayInterface.allowed_setindex!(c3,x_old - im*epsilon,i)
187+
ArrayInterfaceCore.allowed_setindex!(c1,c1_old + im*epsilon,i)
188+
ArrayInterfaceCore.allowed_setindex!(c3,x_old - im*epsilon,i)
189189
df_tmp2 = im*imag( (f(c1) - f(c3)) / (2*im*epsilon) )
190-
ArrayInterface.allowed_setindex!(df,df_tmp-df_tmp2,i)
190+
ArrayInterfaceCore.allowed_setindex!(df,df_tmp-df_tmp2,i)
191191
else
192-
ArrayInterface.allowed_setindex!(df,df_tmp,i)
192+
ArrayInterfaceCore.allowed_setindex!(df,df_tmp,i)
193193
end
194-
ArrayInterface.allowed_setindex!(c1,c1_old, i)
195-
ArrayInterface.allowed_setindex!(c3,x_old,i)
194+
ArrayInterfaceCore.allowed_setindex!(c1,c1_old, i)
195+
ArrayInterfaceCore.allowed_setindex!(c3,x_old,i)
196196
end
197197
elseif fdtype == Val(:complex) && returntype <: Real
198198
copyto!(c1,x)
199199
epsilon_complex = eps(real(eltype(x)))
200200
# we use c1 here to avoid typing issues with x
201201
@inbounds for i ∈ eachindex(x)
202-
c1_old = ArrayInterface.allowed_getindex(c1,i)
203-
ArrayInterface.allowed_setindex!(c1,c1_old+im*epsilon_complex,i)
204-
ArrayInterface.allowed_setindex!(df,imag(f(c1)) / epsilon_complex,i)
205-
ArrayInterface.allowed_setindex!(c1,c1_old,i)
202+
c1_old = ArrayInterfaceCore.allowed_getindex(c1,i)
203+
ArrayInterfaceCore.allowed_setindex!(c1,c1_old+im*epsilon_complex,i)
204+
ArrayInterfaceCore.allowed_setindex!(df,imag(f(c1)) / epsilon_complex,i)
205+
ArrayInterfaceCore.allowed_setindex!(c1,c1_old,i)
206206
end
207207
else
208208
fdtype_error(returntype)

‎src/hessians.jl‎

Lines changed: 22 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -68,27 +68,27 @@ function finite_difference_hessian!(H,f,x,
6868
end
6969

7070
for i = 1:n
71-
xi = ArrayInterface.allowed_getindex(x,i)
71+
xi = ArrayInterfaceCore.allowed_getindex(x,i)
7272
epsilon = compute_epsilon(Val(:hcentral), xi, relstep, absstep)
7373

7474
if inplace === Val(true)
75-
ArrayInterface.allowed_setindex!(xpp,xi + epsilon,i)
76-
ArrayInterface.allowed_setindex!(xmm,xi - epsilon,i)
75+
ArrayInterfaceCore.allowed_setindex!(xpp,xi + epsilon,i)
76+
ArrayInterfaceCore.allowed_setindex!(xmm,xi - epsilon,i)
7777
else
7878
_xpp = Base.setindex(xpp,xi + epsilon, i)
7979
_xmm = Base.setindex(xmm,xi - epsilon, i)
8080
end
8181

82-
ArrayInterface.allowed_setindex!(H,(f(_xpp) - 2*fx + f(_xmm)) / epsilon^2,i,i)
82+
ArrayInterfaceCore.allowed_setindex!(H,(f(_xpp) - 2*fx + f(_xmm)) / epsilon^2,i,i)
8383
epsiloni = compute_epsilon(Val(:central), xi, relstep, absstep)
8484
xp = xi + epsiloni
8585
xm = xi - epsiloni
8686

8787
if inplace === Val(true)
88-
ArrayInterface.allowed_setindex!(xpp,xp,i)
89-
ArrayInterface.allowed_setindex!(xpm,xp,i)
90-
ArrayInterface.allowed_setindex!(xmp,xm,i)
91-
ArrayInterface.allowed_setindex!(xmm,xm,i)
88+
ArrayInterfaceCore.allowed_setindex!(xpp,xp,i)
89+
ArrayInterfaceCore.allowed_setindex!(xpm,xp,i)
90+
ArrayInterfaceCore.allowed_setindex!(xmp,xm,i)
91+
ArrayInterfaceCore.allowed_setindex!(xmm,xm,i)
9292
else
9393
_xpp = Base.setindex(xpp,xp,i)
9494
_xpm = Base.setindex(xpm,xp,i)
@@ -97,30 +97,30 @@ function finite_difference_hessian!(H,f,x,
9797
end
9898

9999
for j = i+1:n
100-
xj = ArrayInterface.allowed_getindex(x,j)
100+
xj = ArrayInterfaceCore.allowed_getindex(x,j)
101101
epsilonj = compute_epsilon(Val(:central), xj, relstep, absstep)
102102
xp = xj + epsilonj
103103
xm = xj - epsilonj
104104

105105
if inplace === Val(true)
106-
ArrayInterface.allowed_setindex!(xpp,xp,j)
107-
ArrayInterface.allowed_setindex!(xpm,xm,j)
108-
ArrayInterface.allowed_setindex!(xmp,xp,j)
109-
ArrayInterface.allowed_setindex!(xmm,xm,j)
106+
ArrayInterfaceCore.allowed_setindex!(xpp,xp,j)
107+
ArrayInterfaceCore.allowed_setindex!(xpm,xm,j)
108+
ArrayInterfaceCore.allowed_setindex!(xmp,xp,j)
109+
ArrayInterfaceCore.allowed_setindex!(xmm,xm,j)
110110
else
111111
_xpp = Base.setindex(_xpp,xp,j)
112112
_xpm = Base.setindex(_xpm,xm,j)
113113
_xmp = Base.setindex(_xmp,xp,j)
114114
_xmm = Base.setindex(_xmm,xm,j)
115115
end
116116

117-
ArrayInterface.allowed_setindex!(H,(f(_xpp) - f(_xpm) - f(_xmp) + f(_xmm))/(4*epsiloni*epsilonj),i,j)
117+
ArrayInterfaceCore.allowed_setindex!(H,(f(_xpp) - f(_xpm) - f(_xmp) + f(_xmm))/(4*epsiloni*epsilonj),i,j)
118118

119119
if inplace === Val(true)
120-
ArrayInterface.allowed_setindex!(xpp,xj,j)
121-
ArrayInterface.allowed_setindex!(xpm,xj,j)
122-
ArrayInterface.allowed_setindex!(xmp,xj,j)
123-
ArrayInterface.allowed_setindex!(xmm,xj,j)
120+
ArrayInterfaceCore.allowed_setindex!(xpp,xj,j)
121+
ArrayInterfaceCore.allowed_setindex!(xpm,xj,j)
122+
ArrayInterfaceCore.allowed_setindex!(xmp,xj,j)
123+
ArrayInterfaceCore.allowed_setindex!(xmm,xj,j)
124124
else
125125
_xpp = Base.setindex(_xpp,xj,j)
126126
_xpm = Base.setindex(_xpm,xj,j)
@@ -130,10 +130,10 @@ function finite_difference_hessian!(H,f,x,
130130
end
131131

132132
if inplace === Val(true)
133-
ArrayInterface.allowed_setindex!(xpp,xi,i)
134-
ArrayInterface.allowed_setindex!(xpm,xi,i)
135-
ArrayInterface.allowed_setindex!(xmp,xi,i)
136-
ArrayInterface.allowed_setindex!(xmm,xi,i)
133+
ArrayInterfaceCore.allowed_setindex!(xpp,xi,i)
134+
ArrayInterfaceCore.allowed_setindex!(xpm,xi,i)
135+
ArrayInterfaceCore.allowed_setindex!(xmp,xi,i)
136+
ArrayInterfaceCore.allowed_setindex!(xmm,xi,i)
137137
end
138138
end
139139
LinearAlgebra.copytri!(H,'U')

‎src/iteration_utils.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ end
2929
end
3030

3131
#override default setting of using findstructralnz
32-
_use_findstructralnz(sparsity) = ArrayInterface.has_sparsestruct(sparsity)
32+
_use_findstructralnz(sparsity) = ArrayInterfaceCore.has_sparsestruct(sparsity)
3333
_use_findstructralnz(::SparseMatrixCSC) = false
3434

3535
# test if J, sparsity are both SparseMatrixCSC and have the same sparsity pattern of stored values

0 commit comments

Comments
 (0)