Skip to content

Change initialization of tolerances of fixed point iterations #95

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.

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

Merged
merged 4 commits into from
Jan 20, 2019

Conversation

devmotion
Copy link
Member

This PR changes how the tolerances of fixed point iterations are initialized. If unspecified, the error tolerances of the integrator are used without any conversions - the correct types are supposed to be set by the logic in OrdinaryDiffEq (https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve.jl#L112-L134). If they are specified, then it is only ensured that they are real numbers (similar to https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve.jl#L121 and https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve.jl#L133). Hence if units are used and tolerances are specified they always have to be given with correct units. Currently unitless numbers were converted to the correct units but in my opinion we shouldn't do that since it is not transparent to the user.

This PR, however, has bigger implications: when used with Flux currently fixedpoint_reltol_internal is initialized as TrackedReal{TrackedReal{Float64}} which is incorrect and causes a stack overflow. With this PR the following code (modification of JuliaLang/www.julialang.org#274) runs without any issues:

using DelayDiffEq
using DiffEqFlux
using Flux
using Plots

gr()

function delay_lotka_volterra(du,u,h,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx =- β*y)*h(p,t-0.1)[1]
  du[2] = dy =*x - γ)*y
end

h(p,t) = ones(eltype(p),2)

prob = DDEProblem(delay_lotka_volterra,[1.0,1.0],h,(0.0,10.0),constant_lags=[0.1])

p = param([2.2, 1.0, 2.0, 0.4])
params = Flux.Params([p])

reduction(sol) = sol[1,:]

function predict_rd_dde()
  diffeq_rd(p,reduction,prob,MethodOfSteps(Tsit5()),saveat=0.1)
end

loss_rd_dde() = sum(abs2,x-1 for x in predict_rd_dde())

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function ()
  display(loss_rd_dde())
  display(plot(solve(remake(prob,p=Flux.data(p)),MethodOfSteps(Tsit5()),saveat=0.1),
          ylim=(0,6)))
end

cb()

Flux.train!(loss_rd_dde, params, data, opt, cb = cb)

abstol = integrator.opts.abstol
if typeof(abstol) <: Number
fixedpoint_abstol_internal = abstol
else
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure whether this specialization is required. To me it seems recursivecopy(::Number) falls back to deepcopy which both give a more complex output with e.g. @code_native recursivecopy(3) and @code_native deepcopy(3) than @code_native copy(3)? Should we just define recursivecopy(x::Number) = copy(x) in RecursiveArrayTools which by default just returns x (https://github.com/JuliaLang/julia/blob/master/base/number.jl#L88)?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. There shouldn't be a need to deepcopy a number anymore. This could improve performance in some places as well.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should make sure Measurements.jl still works after changing to copy. It should.

@devmotion devmotion merged commit 7b62c3f into master Jan 20, 2019
@devmotion devmotion deleted the tolerances branch January 20, 2019 20:51
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants