Skip to content

Commit f1bb430

Browse files
Merge pull request #821 from ChrisRackauckas-Claude/narrow-unseed
Unseed only the chunk, and split unseed! out of seed!
2 parents 37bab0a + a337ee6 commit f1bb430

9 files changed

Lines changed: 170 additions & 71 deletions

File tree

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "ForwardDiff"
22
uuid = "f6369f11-7733-5829-9624-2563aa707210"
3-
version = "1.4.3"
3+
version = "1.4.4"
44

55
[deps]
66
CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"

src/apiutils.jl

Lines changed: 21 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ end
2727
function vector_mode_dual_eval!(f!::F, cfg::JacobianConfig, y, x) where {F}
2828
ydual, xdual = cfg.duals
2929
seed!(xdual, x, cfg.seeds)
30-
seed!(ydual, y)
30+
seed_zero_partials!(ydual, y)
3131
f!(ydual, xdual)
3232
return ydual
3333
end
@@ -70,14 +70,30 @@ function structural_eachindex(x::Diagonal, y::AbstractArray)
7070
return diagind(x)
7171
end
7272

73-
function seed!(duals::AbstractArray{Dual{T,V,N}}, x,
74-
seed::Partials{N,V} = zero(Partials{N,V})) where {T,V,N}
73+
# Copies the values of `x` into `duals` with zero partials. Used both to remove seeds `duals` is
74+
# currently carrying and to initialize a freshly allocated work buffer, whose elements must all be
75+
# written before the target function reads them.
76+
seed_zero_partials!(duals::AbstractArray{Dual{T,V,N}}, x) where {T,V,N} =
77+
_seed_zero_partials!(duals, x, structural_eachindex(duals, x))
78+
79+
# Zeroes the partials of `count` elements starting at structural position `index`. Chunk mode only
80+
# needs to clear the chunk it just seeded, so writing through to the end of the array would be O(n)
81+
# redundant work per chunk, i.e. O(n^2/N) per sweep. `count` mirrors the `chunksize` argument of
82+
# `seed!(duals, x, index, seeds, chunksize)`.
83+
function seed_zero_partials!(duals::AbstractArray{Dual{T,V,N}}, x, index,
84+
count = N) where {T,V,N}
85+
idxs = Iterators.take(Iterators.drop(structural_eachindex(duals, x), index - 1), count)
86+
return _seed_zero_partials!(duals, x, idxs)
87+
end
88+
89+
function _seed_zero_partials!(duals::AbstractArray{Dual{T,V,N}}, x, idxs) where {T,V,N}
90+
seed = zero(Partials{N,V})
7591
if isbitstype(V)
76-
for idx in structural_eachindex(duals, x)
92+
for idx in idxs
7793
duals[idx] = Dual{T,V,N}(x[idx], seed)
7894
end
7995
else
80-
for idx in structural_eachindex(duals, x)
96+
for idx in idxs
8197
if isassigned(x, idx)
8298
duals[idx] = Dual{T,V,N}(x[idx], seed)
8399
else
@@ -106,26 +122,6 @@ function seed!(duals::AbstractArray{Dual{T,V,N}}, x,
106122
return duals
107123
end
108124

109-
function seed!(duals::AbstractArray{Dual{T,V,N}}, x, index,
110-
seed::Partials{N,V} = zero(Partials{N,V})) where {T,V,N}
111-
offset = index - 1
112-
idxs = Iterators.drop(structural_eachindex(duals, x), offset)
113-
if isbitstype(V)
114-
for idx in idxs
115-
duals[idx] = Dual{T,V,N}(x[idx], seed)
116-
end
117-
else
118-
for idx in idxs
119-
if isassigned(x, idx)
120-
duals[idx] = Dual{T,V,N}(x[idx], seed)
121-
else
122-
Base._unsetindex!(duals, idx)
123-
end
124-
end
125-
end
126-
return duals
127-
end
128-
129125
function seed!(duals::AbstractArray{Dual{T,V,N}}, x, index,
130126
seeds::NTuple{N,Partials{N,V}}, chunksize = N) where {T,V,N}
131127
offset = index - 1

src/derivative.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ Set `check` to `Val{false}()` to disable tag checking. This can lead to perturba
2727
require_one_based_indexing(y)
2828
CHK && checktag(T, f!, x)
2929
ydual = cfg.duals
30-
seed!(ydual, y)
30+
seed_zero_partials!(ydual, y)
3131
f!(ydual, Dual{T}(x, one(x)))
3232
map!(value, y, ydual)
3333
return extract_derivative(T, ydual)
@@ -65,7 +65,7 @@ Set `check` to `Val{false}()` to disable tag checking. This can lead to perturba
6565
result isa DiffResult ? require_one_based_indexing(y) : require_one_based_indexing(result, y)
6666
CHK && checktag(T, f!, x)
6767
ydual = cfg.duals
68-
seed!(ydual, y)
68+
seed_zero_partials!(ydual, y)
6969
f!(ydual, Dual{T}(x, one(x)))
7070
result = extract_value!(T, result, y, ydual)
7171
result = extract_derivative!(T, result, ydual)

src/gradient.jl

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -127,22 +127,23 @@ function chunk_mode_gradient_expr(result_definition::Expr)
127127
# seed work vectors
128128
xdual = cfg.duals
129129
seeds = cfg.seeds
130-
seed!(xdual, x)
131130

132-
# do first chunk manually to calculate output type
131+
# do first chunk manually to calculate output type. Seeding the first chunk and zeroing the
132+
# remaining elements partitions `xdual`, so every element is initialized exactly once.
133133
seed!(xdual, x, 1, seeds)
134+
seed_zero_partials!(xdual, x, N + 1, xlen - N)
134135
ydual = f(xdual)
135136
$(result_definition)
136137
extract_gradient_chunk!(T, result, ydual, 1, N)
137-
seed!(xdual, x, 1)
138+
seed_zero_partials!(xdual, x, 1)
138139

139140
# do middle chunks
140141
for c in middlechunks
141142
i = ((c - 1) * N + 1)
142143
seed!(xdual, x, i, seeds)
143144
ydual = f(xdual)
144145
extract_gradient_chunk!(T, result, ydual, i, N)
145-
seed!(xdual, x, i)
146+
seed_zero_partials!(xdual, x, i)
146147
end
147148

148149
# do final chunk

src/jacobian.jl

Lines changed: 11 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -184,22 +184,24 @@ function jacobian_chunk_mode_expr(work_array_definition::Expr, compute_ydual::Ex
184184
$(work_array_definition)
185185
seeds = cfg.seeds
186186

187-
# do first chunk manually to calculate output type
187+
# do first chunk manually to calculate output type. Seeding the first chunk and zeroing the
188+
# remaining elements partitions `xdual`, so every element is initialized exactly once.
188189
seed!(xdual, x, 1, seeds)
190+
seed_zero_partials!(xdual, x, N + 1, xlen - N)
189191
$(compute_ydual)
190192
ydual isa AbstractArray || throw(JACOBIAN_ERROR)
191193
$(result_definition)
192194
out_reshaped = reshape_jacobian(result, ydual, xdual)
193195
extract_jacobian_chunk!(T, out_reshaped, ydual, 1, N)
194-
seed!(xdual, x, 1)
196+
seed_zero_partials!(xdual, x, 1)
195197

196198
# do middle chunks
197199
for c in middlechunks
198200
i = ((c - 1) * N + 1)
199201
seed!(xdual, x, i, seeds)
200202
$(compute_ydual)
201203
extract_jacobian_chunk!(T, out_reshaped, ydual, i, N)
202-
seed!(xdual, x, i)
204+
seed_zero_partials!(xdual, x, i)
203205
end
204206

205207
# do final chunk
@@ -214,41 +216,29 @@ function jacobian_chunk_mode_expr(work_array_definition::Expr, compute_ydual::Ex
214216
end
215217

216218
@eval function chunk_mode_jacobian(f::F, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
217-
$(jacobian_chunk_mode_expr(quote
218-
xdual = cfg.duals
219-
seed!(xdual, x)
220-
end,
219+
$(jacobian_chunk_mode_expr(:(xdual = cfg.duals),
221220
:(ydual = f(xdual)),
222221
:(result = similar(ydual, valtype(T, eltype(ydual)), length(ydual), xlen)),
223222
:()))
224223
end
225224

226225
@eval function chunk_mode_jacobian(f!::F, y, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
227-
$(jacobian_chunk_mode_expr(quote
228-
ydual, xdual = cfg.duals
229-
seed!(xdual, x)
230-
end,
231-
:(f!(seed!(ydual, y), xdual)),
226+
$(jacobian_chunk_mode_expr(:((ydual, xdual) = cfg.duals),
227+
:(f!(seed_zero_partials!(ydual, y), xdual)),
232228
:(result = similar(y, length(y), xlen)),
233229
:(map!(d -> value(T,d), y, ydual))))
234230
end
235231

236232
@eval function chunk_mode_jacobian!(result, f::F, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
237-
$(jacobian_chunk_mode_expr(quote
238-
xdual = cfg.duals
239-
seed!(xdual, x)
240-
end,
233+
$(jacobian_chunk_mode_expr(:(xdual = cfg.duals),
241234
:(ydual = f(xdual)),
242235
:(),
243236
:(extract_value!(T, result, ydual))))
244237
end
245238

246239
@eval function chunk_mode_jacobian!(result, f!::F, y, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
247-
$(jacobian_chunk_mode_expr(quote
248-
ydual, xdual = cfg.duals
249-
seed!(xdual, x)
250-
end,
251-
:(f!(seed!(ydual, y), xdual)),
240+
$(jacobian_chunk_mode_expr(:((ydual, xdual) = cfg.duals),
241+
:(f!(seed_zero_partials!(ydual, y), xdual)),
252242
:(),
253243
:(extract_value!(T, result, y, ydual))))
254244
end

test/AllocationsTest.jl

Lines changed: 11 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -7,22 +7,27 @@ include(joinpath(dirname(@__FILE__), "utils.jl"))
77

88
convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,Float64,8},4},2}, 1.3)
99

10-
@testset "Test seed! allocations" begin
10+
@testset "Test seed!/seed_zero_partials! allocations" begin
1111
x = rand(1000)
1212
cfg = ForwardDiff.GradientConfig(nothing, x)
1313
duals = cfg.duals
1414
seeds = cfg.seeds
15-
seed = cfg.seeds[1]
1615

1716
allocs_seed!(args...) = @allocated ForwardDiff.seed!(args...)
1817
allocs_seed!(duals, x, seeds)
1918
@test iszero(allocs_seed!(duals, x, seeds))
20-
allocs_seed!(duals, x, seed)
21-
@test iszero(allocs_seed!(duals, x, seed))
2219
allocs_seed!(duals, x, 1, seeds)
2320
@test iszero(allocs_seed!(duals, x, 1, seeds))
24-
allocs_seed!(duals, x, 1, seed)
25-
@test iszero(allocs_seed!(duals, x, 1, seed))
21+
22+
# the 4-arg form passes `count` as a runtime value, so it catches an inference regression at the
23+
# `_seed_zero_partials!` boundary that the forms defaulting `count` to `N` could hide
24+
allocs_szp!(args...) = @allocated ForwardDiff.seed_zero_partials!(args...)
25+
allocs_szp!(duals, x)
26+
@test iszero(allocs_szp!(duals, x))
27+
allocs_szp!(duals, x, 1)
28+
@test iszero(allocs_szp!(duals, x, 1))
29+
allocs_szp!(duals, x, 1, 4)
30+
@test iszero(allocs_szp!(duals, x, 1, 4))
2631

2732
allocs_convert_test_574() = @allocated convert_test_574()
2833
allocs_convert_test_574()

test/JacobianTest.jl

Lines changed: 21 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -298,18 +298,27 @@ end
298298
@test res == I
299299
end
300300

301-
# Unassigned (but unused) entry in the input and unassigned entries in the output
302-
resize!(x, 10)
303-
f = (y, x) -> copyto!(y, 1, x, 1, 9)
304-
for chunksize in (1, 2, 10)
305-
y = similar(x, 9)
306-
@test all(i -> !isassigned(y, i), eachindex(y))
307-
cfg = ForwardDiff.JacobianConfig(f, y, x, ForwardDiff.Chunk{chunksize}())
308-
res = ForwardDiff.jacobian(f, y, x, cfg)
309-
@test y == x[1:(end-1)]
310-
@test res isa Matrix{BigFloat}
311-
@test res[:, 1:(end-1)] == I
312-
@test all(iszero, res[:, end])
301+
# Unassigned (but unused) entry in the input and unassigned entries in the output. `hole` is
302+
# varied so the unassigned entry lands in a middle chunk as well as in the last one: only the
303+
# former reaches the `Base._unsetindex!` branch of the windowed seeding path, since the last
304+
# chunk is never cleared.
305+
@testset "unassigned input entry at $hole" for hole in (5, 10)
306+
x = Vector{BigFloat}(undef, 10)
307+
for i in eachindex(x)
308+
i == hole || (x[i] = BigFloat(i))
309+
end
310+
used = [i for i in eachindex(x) if i != hole]
311+
f = (y, x) -> (for (k, i) in enumerate(used); y[k] = x[i]; end; y)
312+
for chunksize in (1, 2, 10)
313+
y = similar(x, 9)
314+
@test all(i -> !isassigned(y, i), eachindex(y))
315+
cfg = ForwardDiff.JacobianConfig(f, y, x, ForwardDiff.Chunk{chunksize}())
316+
res = ForwardDiff.jacobian(f, y, x, cfg)
317+
@test y == x[used]
318+
@test res isa Matrix{BigFloat}
319+
@test res[:, used] == I
320+
@test all(iszero, res[:, hole])
321+
end
313322
end
314323
end
315324

test/SeedTest.jl

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
module SeedTest
2+
3+
import ForwardDiff
4+
using ForwardDiff: Partials
5+
using LinearAlgebra
6+
using Test
7+
8+
include("utils.jl")
9+
10+
# The windowed `seed_zero_partials!` is only ever called to clear a chunk that was just seeded, so
11+
# clearing *too much* is harmless and no test written against the public API can distinguish a
12+
# correctly bounded implementation from an unbounded one. These tests pin the window down directly:
13+
# they seed every structural position with a marker whose partials are all nonzero, clear a window,
14+
# and check exactly which positions lost their marker.
15+
#
16+
# The expected structural index sets are written out by hand rather than obtained from
17+
# `structural_eachindex`, so a bug in that iterator cannot hide inside the assertions depending on
18+
# it; one test ties the two together. Order is significant: `index` and `count` are positions along
19+
# the sequence, not array indices. The sets are heterogeneous by design — `Vector` and `Diagonal`
20+
# enumerate linear indices (the latter via `diagind`), `UpperTriangular` enumerates `CartesianIndex`
21+
# in column-major order.
22+
const SEED_CASES = (
23+
(rand(10), collect(1:10)),
24+
(UpperTriangular(rand(5, 5)), [CartesianIndex(i, j) for j in 1:5 for i in 1:j]),
25+
(Diagonal(rand(6, 6)), collect(1:7:36)),
26+
)
27+
28+
# Positions within `sidx` whose partials are zero.
29+
zeroed_positions(duals, sidx) =
30+
[i for (i, idx) in enumerate(sidx) if iszero(ForwardDiff.partials(duals[idx]))]
31+
32+
# Compares over *every* index of `x`, not just the structural ones, so a bug misplacing values
33+
# outside the structural set is visible. Off-structure reads are safe: the wrapper types return
34+
# `zero(Dual)` without touching the (uninitialized) parent storage.
35+
values_match(duals, x) = all(idx -> ForwardDiff.value(duals[idx]) == x[idx], eachindex(x))
36+
37+
function fill_marker!(duals, x, sidx, marker)
38+
D = eltype(duals)
39+
for idx in sidx
40+
duals[idx] = D(x[idx], marker)
41+
end
42+
return duals
43+
end
44+
45+
@testset "seed_zero_partials!: $(nameof(typeof(x)))" for (x, sidx) in SEED_CASES
46+
cfg = ForwardDiff.GradientConfig(nothing, x, ForwardDiff.Chunk{3}())
47+
duals, seeds = cfg.duals, cfg.seeds
48+
N = ForwardDiff.npartials(eltype(duals))
49+
marker = Partials(ntuple(i -> Float64(i), N))
50+
nstruct = length(sidx)
51+
52+
# everything below counts positions along `sidx`, so pin it to the implementation once
53+
@test collect(ForwardDiff.structural_eachindex(duals, x)) == sidx
54+
@test ForwardDiff.structural_length(x) == nstruct
55+
56+
# `count` defaults to N
57+
fill_marker!(duals, x, sidx, marker)
58+
ForwardDiff.seed_zero_partials!(duals, x, 4)
59+
@test zeroed_positions(duals, sidx) == collect(4:(4 + N - 1))
60+
@test values_match(duals, x)
61+
62+
# an explicit `count` narrows the window; a `count` overrunning the end is clamped by
63+
# `Iterators.take` rather than throwing; a zero-width window is a no-op, which is what makes
64+
# `xlen - N` safe as the `count` of chunk mode's tail clear
65+
@testset "index=$index count=$count" for (index, count, expected) in
66+
((4, 2, 4:5),
67+
(nstruct - 1, N, (nstruct - 1):nstruct),
68+
(1, 0, 1:0))
69+
fill_marker!(duals, x, sidx, marker)
70+
ForwardDiff.seed_zero_partials!(duals, x, index, count)
71+
@test zeroed_positions(duals, sidx) == collect(expected)
72+
@test values_match(duals, x)
73+
end
74+
75+
# the 2-arg form clears every structural position
76+
fill_marker!(duals, x, sidx, marker)
77+
ForwardDiff.seed_zero_partials!(duals, x)
78+
@test zeroed_positions(duals, sidx) == collect(1:nstruct)
79+
@test values_match(duals, x)
80+
81+
# `seed!` and `seed_zero_partials!` must agree on what "the chunk at `index`" is, or chunk mode
82+
# would leave stale seeds behind. `duals` enters each iteration fully cleared.
83+
@testset "round-trips seed! at index=$index" for index in unique((1, 4, nstruct - N + 1))
84+
ForwardDiff.seed!(duals, x, index, seeds)
85+
@test zeroed_positions(duals, sidx) ==
86+
[i for i in 1:nstruct if !(index <= i <= index + N - 1)]
87+
ForwardDiff.seed_zero_partials!(duals, x, index)
88+
@test zeroed_positions(duals, sidx) == collect(1:nstruct)
89+
@test values_match(duals, x)
90+
end
91+
end
92+
93+
end # module

test/runtests.jl

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,11 @@ Random.seed!(SEED)
4848
t = @elapsed include("MiscTest.jl")
4949
println("##### done (took $t seconds).")
5050
end
51+
@testset "Seeding" begin
52+
println("##### Testing seeding...")
53+
t = @elapsed include("SeedTest.jl")
54+
println("##### done (took $t seconds).")
55+
end
5156
@testset "Allocations" begin
5257
println("##### Testing allocations...")
5358
t = @elapsed include("AllocationsTest.jl")

0 commit comments

Comments
 (0)