I am working on a problem where I have a continuous real valued function, e.g. a solution to an ODE, and I want to be able to update it multiple times by adding other continuous functions, which are calculated by code running elsewhere (i.e. not user defined).
For example, the workflow looks something like: solve an ODE and save the solution as a function of time, f(t), do some calculations to find the “update” function, g(t), and then update f by f(t) += g(t).
I’ve put together a minimal working example to try and illustrate what I’ve tried. The idea I’ve had is to define a SolutionObject type which contains an at_time field, which is the function I want to update. The tspan field defines the range over which this function is defined but is not directly relevant to the problem. I saw this answer on StackOverflow which shows how to add together anonymous functions using the let keyword:
julia> mutable struct SolutionObject{T<:AbstractFloat}
at_time # Function (1-d continuous)
tspan::Tuple{T, T}
end
julia> sol_object = SolutionObject(t -> t^2, (0., 10.))
julia> update_function = t -> 2t
# add together the outputs of the two functions and save it
# back in the original function (not type stable)
julia> sol_object.at_time = let old_function = sol_object.at_time
t -> old_function(t) + update_function(t)
end
# this gives the expected result
julia> println(sol_object.at_time(3.0)) # 15.0 (= 9.0 + 6.0)
julia> using BenchmarkTools
# but is slow compared to the direct definition
julia> @btime sol_object.at_time(3.0) # 71.368 ns (4 allocations: 64 bytes)
julia> added = t -> t^2 + 2t
julia> @btime added(3.0) # 16.139 ns (1 allocation: 16 bytes)
This is pretty slow, and also not type stable. If I make a different constructor where the type of at_time is defined then the above code gives an error, which is the expected behaviour:
julia> mutable struct TypedSolutionObject{F<:Function, T<:AbstractFloat}
at_time::F
tspan::Tuple{T,T}
end
julia> typed_sol_object = TypedSolutionObject(t -> t^2, (0., 10.))
julia> typed_sol_object.at_time = let old_function = typed_sol_object.at_time
t -> old_function(t) + update_function(t)
end
ERROR: MethodError: Cannot `convert` an object of type var"#15#16"{var"#13#14"} to an object of type var"#13#14"
Closest candidates are:
convert(::Type{T}, ::T) where T at Base.jl:61
Stacktrace:
[1] setproperty!(x::TypedSolutionObject{var"#13#14", Float64}, f::Symbol, v::Function)
@ Base ./Base.jl:39
[2] top-level scope
@ ~/projects/mwe/anon_function_types.jl:31
One workaround here is to define a new TypedSolutionObject with the updated_f, but this gives the same performance as the original code (perhaps unsurprisingly?):
julia> updated_f = let old_function = typed_sol_object.at_time
t -> old_function(t) + update_function(t)
end
julia> updated_typed_solution_object = TypedSolutionObject(updated_f, typed_sol_object.tspan)
julia> @btime updated_typed_solution_object.at_time(3.0) # 74.553 ns (4 allocations: 64 bytes)
In my code I have currently implemented the workaround and it is way too slow and still not type stable, and this approach somehow seems wrong / overly complicated. Does anyone have suggestions of how I could solve this problem, perhaps with a different approach than what I have already tried? I’m happy to provide more info if needed.