Roots.jl icon indicating copy to clipboard operation
Roots.jl copied to clipboard

RFC: storing intermediate values

Open tpapp opened this issue 6 years ago • 5 comments

Sometimes I solve nonlinear systems of the form

f(g(x)) ≈ 0 # could be a scalar or a vector

where g(x) is a large, costly, and ultimately interesting object, so one would call it on the root, but that is a bit wasteful. It would be nice to have an interface that returns it. I am thinking of something like

(root = x, value = g(x)) = find_zero(f, arg_or_bracket, method; by = identity)

where by would be g, and a NamedTuple is returned. But this of course would be breaking, so a new name would perhaps be better.

tpapp avatar Jan 03 '19 14:01 tpapp

Why can't you just use a closure or something and store it yourself?

KristofferC avatar Jan 03 '19 14:01 KristofferC

local val
find_zero(x -> (val = g(x); f(val)))

KristofferC avatar Jan 03 '19 14:01 KristofferC

I don't know how to get the type inferred correctly. Eg

using Roots
function find_zero_ext(f, g, bracket)
    local val
    x = find_zero(x -> (val = g(x); f(val)), bracket)
    (x, val)
end

results in

julia> @code_warntype find_zero_ext(x -> x^3, cbrt, (-1, 1))
Body::Tuple{Float64,Any}
1 ─ %1  = %new(Core.Box)::Core.Box
│   %2  = %new(getfield(Main, Symbol("##9#10")){getfield(Main, Symbol("##15#16")),typeof(cbrt)}, f, g, %1)::getfield(Main, Symbol("##9#10")){getfield(Main, Symbol("##15#16")),typeof(cbrt)}
│   %3  = Roots.find_zero::typeof(find_zero)
│   %4  = invoke Roots.:(#find_zero#8)($(QuoteNode(Roots.NullTracks()))::Roots.NullTracks, false::Bool, $(QuoteNode(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}()))::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, %3::Function, %2::Function, _4::Tuple{Int64,Int64}, $(QuoteNode(Bisection()))::Bisection)::Float64
│   %5  = (Core.isdefined)(%1, :contents)::Bool
└──       goto #3 if not %5
2 ─       goto #4
3 ─       $(Expr(:throw_undef_if_not, :val, false))
4 ┄ %9  = (Core.getfield)(%1, :contents)::Any
│   %10 = (Core.tuple)(%4, %9)::Tuple{Float64,Any}
└──       return %10

tpapp avatar Jan 03 '19 15:01 tpapp

That's the old "boxes in closures"-issue. Can be worked around with something like:

julia> function find_zero_ext(f, g, bracket)
           val = Ref{Float64}()
           x = find_zero(x -> (val[] = g(x); f(val[])), bracket)
           (x, val[])
       end

KristofferC avatar Jan 03 '19 15:01 KristofferC

Assume I don't know the result type of g(x) (or that it's a nested NamedTuple of epic horribleness).

tpapp avatar Jan 03 '19 15:01 tpapp