Skip to content

WIP: memoization in the trees #9

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
ranocha wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

WIP: memoization in the trees #9

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 16 additions & 4 deletions src/BSeries.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,13 +29,21 @@ end
"""
substitute(b, a, t::RootedTree)

Compute the coefficient correspoding to the tree `t` of the B-series that is
Compute the coefficient corresponding to the tree `t` of the B-series that is
formed by substituting the B-series `b` into the B-series `a`.
"""
function substitute(b, a, t::RootedTree)
substitute(b, a, t::RootedTree) = substitute!((nothing, nothing), b, a, t)

"""
substitute!(cache, b, a, t::RootedTree)

Like [`substitute`](@ref) but with possible speed-up from memoization
(at the cost of more memory usage in the `cache`).
"""
function substitute!(cache, b, a, t::RootedTree)
result = zero(first(values(a)) * first(values(b)))

for (forest, skeleton) in PartitionIterator(t)
for (forest, skeleton) in PartitionIterator(t, cache)
result += reduce(*, b[tree] for tree in forest) * a[skeleton]
end

Expand Down Expand Up @@ -222,6 +230,10 @@ function modifying_integrator(A::AbstractMatrix, b::AbstractVector, c::AbstractV
t = rootedtree([1])
series[t] = series_rk[t]

cache = create_cache(PartitionIterator, t)
# cache = (first(cache), nothing) # TODO: Why is this faster?
# cache = (nothing, nothing)

# Recursively solve `substitute(series, series_rk, t) == series_ex[t]`.
# This works because
# substitute(series, series_rk, t) = series[t] + lower order terms
Expand All @@ -231,7 +243,7 @@ function modifying_integrator(A::AbstractMatrix, b::AbstractVector, c::AbstractV
# t = copy(_t)
# which are slightly less efficient due to additional computations and allocations.
for t in keys(series)
series[t] += series_ex[t] - substitute(series, series_rk, t)
series[t] += series_ex[t] - substitute!(cache, series, series_rk, t)
end

return series
Expand Down