I am trying to program a function that tries different root-finding algorithms depending on whether they succeed or not. How can I accomplish that? I programmed a prototype below, but it does not seem to work:
using Roots
algos = (Roots.Order0(), Roots.Order1())
function solveroots(f::F, init, algos) where F
for alg in algos
println("trying algorithm $alg")
try
root = Roots.find_zero(f, init, alg)
catch
@warn "Algorithm $alg errored, continue with another one"
end
end
return root
end
f(x) = x^5 - x + 1/2
but this errors:
julia> solveroots(f, algos, 0.1)
trying algorithm Order0()
trying algorithm Secant()
ERROR: UndefVarError: root not defined
This seems to work, I think there was a problem with the scope of the for loop.
function solveroots(f::F, init, algos) where F
for alg in algos
println("trying algorithm $alg")
root = try
Roots.find_zero(f, init, alg)
catch
@warn "Algorithm $alg errored, continue with another one"
end
return root
end
end
Define root outside of the for loop to a guard value (e.g. nothing), assign the value of the try block to root (note that @warn returns nothing), if root is not nothing break out of the loop, return root at the end whatever it is (may be nothing if there are no working algorithms)
julia> function solveroots(f, init, algos)
root = nothing
for alg in algos
println("trying algorithm $alg")
root = try
Roots.find_zero(f, init, alg)
catch
@warn "Algorithm $alg errored, continue with another one"
end
isnothing(root) || break
end
return root
end
solveroots (generic function with 1 method)
julia> solveroots(f, 0.1, algos)
trying algorithm Order0()
0.5506065793341349
The try block is concerning performance-wise to start with And yes, the fact that root is changing value affects performance, but what would you like to do otherwise?
Just wondering if there is a better way to do the same thing, ie try different algorithms and return the first one that works. If the problem was say in Optim I can use the returned struct to check for convergence, for example; but here find_root returns a float or throws an error and I am not familiar with handling errors.
function solveroots(f, init, algos)
for alg in algos
println("trying algorithm $alg")
root = try
Roots.find_zero(f, init, alg)
catch
@warn "Algorithm $alg errored, continue with another one"
end
isnothing(root) || return root
end
end
If you want to try several different algo in parallel and only return the first one that succeed. Then use the find treasure parallel algorithmn below.
Proof of concept below
using Base.Threads
struct WorkDetail
tid::Int64
data::Vector{Float64}
end
@show rawdata = rand(10)
NumOfThreads = 4
chunksize = length(rawdata) ÷ NumOfThreads
WorkDetailArray = WorkDetail[]
for t = 1:NumOfThreads
if t < NumOfThreads
push!( WorkDetailArray, WorkDetail(t, rawdata[chunksize*(t-1)+1:chunksize*(t)+1-1]) )
else
push!( WorkDetailArray, WorkDetail(t, rawdata[chunksize*(t-1)+1:length(rawdata)]) )
end
end
@show WorkDetailArray
Completed = Bool[ false for _ = 1:NumOfThreads ]
Accepted = Bool[ false for _ = 1:NumOfThreads ]
Success = Bool[ false for _ = 1:NumOfThreads ]
WorkTask(t) = begin
result = nothing
for c = 1:length(WorkDetailArray[t].data)
println("Thread $(t), Working with data = $(WorkDetailArray[t].data[c])")
sleep((NumOfThreads - t + 1) * rand())
FoundTreasureBool = rand() < 0.20
if FoundTreasureBool
Success[t] = true
result = WorkDetailArray[t].data[c]
break
else
# We have not found the Treasure
# now check if others have
# and stop working if they did
if any(Success)
break
end
end
end
"Signal that this thread has completed"
Completed[t] = true
println("Thread $(t) completion has been signaled")
return result
end
Workers = Task[]
for t = 1:NumOfThreads
push!(Workers, Threads.@spawn WorkTask(t))
println("Thread $(t) has been spawned")
end
# while Accepted contains at least one element which is false
treasurefound = false
while ! all(Accepted) && ! treasurefound
for t = 1:length(Workers)
if xor(Completed[t],Accepted[t])
println("Thread $(t) completion has been accepted")
"Mark that controller has accepted the thread completion"
Accepted[t] = true
"Check if the worker is successful in finding the treasure"
if Success[t]
treasure = Threads.fetch(Workers[t])
println("Thread $(t) has found treasure! Its value is ",treasure)
global treasurefound = true
break
else
Threads.wait(Workers[t])
end
end
end
sleep(0.05) # Sleep for 50 milisecond
end
# Now wait for any remainging worker
while ! all(Accepted)
for t = 1:length(Workers)
if xor(Completed[t],Accepted[t])
println("Thread $(t) completion has been accepted")
"Mark that controller has accepted the thread completion"
Accepted[t] = true
Threads.wait(Workers[t])
end
end
sleep(0.05) # Sleep for 50 milisecond
end
println("Done")
exit()
And yes, the fact that