Allow control over the way the intervals are selected to be processed

[Kolaru]

This PR introduces the SearchStrategy type that store the methods used to store and retrieve the intervals in roots and modifies the RootSearchState type to allow for any type of containers. This should in principle allow to have any strategy and thus lay the foundation of what is needed to fix #22 .

Here is a small working example using a Deque to implement a bread-first strategy:

using ValidatedNumerics
using DataStructures

import IntervalRootFinding: SearchStrategy, RootSearchState

f(x) = x^2 - 2
region = -4..4
contractor = Newton(f, x -> ForwardDiff.derivative(f, x))
tol = 1e-3
R = typeof(region)

strat = SearchStrategy(push!, shift!)
working = Deque{R}() 
push!(working, region)
outputs = Deque{Root{R}}()

rss = RootSearchState(working, outputs)
rsearch = RootSearch(rss, contractor, strat, tol)

for state in rsearch end

println(rsearch.state.outputs)

One problem is that due to the current type system (it will change in 0.7 if I get it correctly) it is not easy to write something like

struct RootSearchState{CONTAINER, T}
    working::CONTAINER{T}
    outputs::CONTAINER{Root{T}}
end

For now the RootSearchState is therefore more general than needed and no check is performed to see if the two fields are compatible.

Another problem is that the syntax is yet very verbose. One of the reason for that is that the outputs set must be provided. In principle, its type can be inferred from the type of the working set, however DataStructures do not provide a nice way of doing as far as I know (Base provides the similar function which does exactly that). I do not know what is the best strategy for us, requesting both fields or implementing the similar function (or something similar... ^^') for supported container types.

I post this here as a work in progress as I would appreciate feedback and be sure that we can agree on this design before going further to implement what is missing:

• Cleaner syntax
• Several default strategies
• Integration with the roots function
• Examples
• Tests

[dpsanders]

This is interesting, thanks!

Note that as far as I know, standard Julia arrays can actually be thought of as deques, since pushing and popping (shift! and unshift!, renamed to pushfirst! and popfirst! in Julia 0.7) are efficient -- the data is stored in the middle of the space reserved for the array.

[dpsanders]

One problem is that due to the current type system (it will change in 0.7 if I get it correctly) it is not easy to write something like

struct RootSearchState{CONTAINER, T} working::CONTAINER{T} outputs::CONTAINER{Root{T}} end

I don't think you can actually do this in 0.7 either?

[Kolaru]

Note that as far as I know, standard Julia arrays can actually be thought of as deques, since pushing and popping (shift! and unshift!, renamed to pushfirst! and popfirst! in Julia 0.7) are efficient -- the data is stored in the middle of the space reserved for the array.

Despite the claim of better performance by the DataStructure documentation (which did make me go for it in the first place), my quick tests do not show very noticeable differences so it is probably not worth adding the dependency here for defaults strategies.

I don't think you can actually do this in 0.7 either?

Well I think I saw something along that line, but I can't find it again so I may have just dreamed of it.

[Kolaru]

I don't think you can actually do this in 0.7 either?

Well I think I saw something along that line, but I can't find it again so I may have just dreamed of it.

Okay after checking it may rather be that it can be implemented using SimpleTraits.jl which will get a more usable syntax in 0.7. A simple custom implementation could also do it and avoid a new dependency.

[gwater]

More generally, let's discuss why state is proposed as a field of RootSearch: As I understand it, that becomes necessary in order to implement multiple strategies while implementing only a single Iterator. Semantically, the strategy is proposed as a property of the state, not the search Iterator. I find that conceptually flawed: Different search strategies produce different paths to (hopefully) the same final state. So why not implement multiple Iterators; one for each strategy?

[Kolaru]

The reason I put stateas a field of RootSearch is for the start function: it allows to define the working and outputs containers outside of it. Your other comment made me rethink it and I come to the conclusion that all this could be avoided by modifying SearchStrategy to include a constructor function for the working container:

struct SearchStrategy
    constructor::Function
    retrieve!::Function
    store!::Function
end

constructor takes an element type as argument and returning an empty container (in principle it could take the type T of container as argument and use the syntax T{E}() to produce a the empty container with element type E but I fear some data structures may be incompatible with that syntax).

If I get it correctly, this should solve your concern about the state of the iterator as it is totally determined and produced by the strategy used. Therefore the strategy becomes a property of the search iterator while the state is restored as a simple bookkeeping object.

Finally implementing a RootSearch type for each strategy seems like an overkill to me as what change from one strategy to one other is a well defined subset of the RootSearch properties and is independent of the other properties.

[gwater]

@Kolaru That sounds like an elegant solution to my concerns, thanks!

[codecov-io]

Codecov Report

Merging #72 into master will decrease coverage by 0.54%. The diff coverage is 84.61%.

@@            Coverage Diff             @@
##           master      #72      +/-   ##
==========================================
- Coverage   64.97%   64.42%   -0.55%     
==========================================
  Files          11       12       +1     
  Lines         551      565      +14     
==========================================
+ Hits          358      364       +6     
- Misses        193      201       +8

Impacted Files	Coverage Δ
src/IntervalRootFinding.jl	4.65% <ø> (ø)	:arrow_up:
src/roots.jl	75% <84.21%> (-13.14%)	:arrow_down:
src/rootsearch_iterator.jl	84.84% <84.84%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 5a0027a...3eadb00. Read the comment docs.

[Kolaru]

I finally got the time to implement this. The example now looks much simpler:

using IntervalArithmetic
using IntervalRootFinding
using DataStructures

f(x) = x^2 - 2
region = -4..4
contractor = Newton(f, x -> ForwardDiff.derivative(f, x))
tol = 1e-10

strat = SearchStrategy(R -> Deque{R}(), push!, shift!)
rootsearch = RootSearch(region, contractor, strat, tol)

global state
for state in rootsearch end

println(state.outputs)

@gwater In your implementation you defined the eltype of the RootSearch iterator. Could you explain me the rationale behind it ? I wasn't able to find much information about eltype of Iterator. As for now, I had to change its definition and it does not return a concrete type anymore. I wonder if it could impact performance significantly.

[gwater]

@Kolaru I think it could because we lose type stability in the loop. Can you check what @code_warntype next(rootsearch, state) prints? I imagine there will be a few red indicators without eltype().

I think the eltype problem arises from the constructor field in SearchStrategy. I think it may be possible to make SearchStrategy parametric on the container type rather than just storing a reference to the constructor. I think we can instruct users to only use containers with constructors like Container{ElementType}(). Vector and Deque both have that structure.

[gwater]

I'm thinking of something like

struct SearchStrategy{CONT}
    store!::Function
    retrieve!::Function
end

function SearchStrategy(CONT::Type, push!::Function, retrieve!::Function)
    SearchStrategy{CONT}(push!, retrieve!)
end

struct RootSearch{R <: Union{Interval,IntervalBox}, CONT, C <: Contractor, T <: Real}
    region::R
    contractor::C
    strategy::SearchStrategy{CONT}
    tolerance::T
end

struct RootSearchState{T, CONT}
    working::CONT{T}  # Should be a container of the form  
    outputs::CONT{Root{T}}  # Should ba a container of root of the form 
end

eltype(::Type{RS}) where {T, CONT, RS <: RootSearch{T, CONT}} = RootSearchState{T, CONT}

function RootSearchState(region::R, strat::SearchStrategy{C}) where {R <: Union{Interval,IntervalBox}, CONT}
    working = CONT{R}()
    outputs = CONT{Root{R}}()
    strat.store!(working, region)
    return RootSearchState(working, outputs)
end

(untested)

[Kolaru]

@gwater Looks like your last message and my commit just crossed, but we arrived at (mostly) the same conclusion.

Except from the nomenclature, the main difference is that I defined the RootSearch struct more poorly (which make the definition of several other functions a bit akward).

Concerning eltype, I think that it is not used at all during the iteration. However, I still had a problem, because the type of the state returned by the start function could not be inferred. This is now fixed.

[Kolaru]

So here is the changes brought by the last commit:

• PR rebased into master.
• Better definition of the RootSearch object and better use of the fact that SearchStrategy is parametrized by the container type used.
• Separation of the RootSearch iterator and roots functions in two files, so the roots.jl file only contains the (numerous) roots functions and branch_and_prune.
• Definition of two simple SearchStrategy: BreadfirstSearch and DepthfirstSearch.
• Addition of the possibility to give a SearchStrategy to the roots function (which lead to the next bullet).
• Rework of the roots interface to keep all optional parameters in the same place and definition of an internal _roots function which handle all concrete cases. Definition of constants for defaults parameters. It may be better to have contractor, strategyand tol as keywords argument of the roots function, but this would be a breaking change and is beyond the scope of this PR.
• Addition or rewriting a several docstrings to keep them up to date.
• Addition of basic tests for SearchStrategy. They are not very interesting yet, but implementing a limit on the number of working regions should provide more vibrant differences between search strategies.
• Lowering of the requirement for Complex tests. For some (currently mysterious) reasons the changes make computing roots with autodiff or with explicit derivative not exactly match for Krawczyk method.

Note that the ability to give a sizehint! to the container has been lost in the process, but I do not know how much impact this may have.

Quick benchmark using the Smiley 2.2 example show no difference in performance between this PR and current master.

[gwater]

I finally got a round to play with this version and I think it works fine (as the tests demonstrate). It's very cool to see how depth-first vs breadth-first affects the number of working intervals.

I noticed a few issues, which I will annotate in the code.

[gwater]

The changes also affect one example in examples/roots.jl. On line 28 one now needs to declare the search strategy, i.e.:

search = RootSearch(-10..10, contractor, DepthFirstSearch, 1e-3)

[dpsanders]

Sorry for the delay in getting to this. I will have a detailed look.

[dpsanders]

This looks great, thanks a lot! Is it ready to merge?

Two general comments:

1. The possible options are now getting so complicated that it may be term to introduce some kind of RootProblem object to hold the various options, analogous to ODEProblem in DifferentialEquations.jl. This should simplify a lot of function signatures as well. cc @ChrisRackauckas

2. The possibility of using different search strategies seems to be even more crucial for IntervalOptimisation.jl. Maybe this should be refactored out into some kind of separate package? Something like IntervalBranchAndBound.jl.

[dpsanders]

It seems to me that the outputs field in RootSearchState does not need to be a container of any particular type -- a simple Vector should be enough for that?

Or, rather, it should be some kind of tree object that reflects how a subinterval(box) is related to a bigger interval(box).

[Kolaru]

This looks great, thanks a lot! Is it ready to merge?

Thanks ! It is ready to merged if we are still targeting Julia 0.6, but it is still untested with the recent revision of IntervalArithmetic or Julia 0.7.

The possible options are now getting so complicated that it may be term to introduce some kind of RootProblem object to hold the various options, analogous to ODEProblem in DifferentialEquations.jl. This should simplify a lot of function signatures as well. cc @ChrisRackauckas

As far as I can see, this would come down to promoting the roots function to a type. From an implementation point of view it doesn't make much difference, but I am not sure how it would change usability for users.

The possibility of using different search strategies seems to be even more crucial for IntervalOptimisation.jl. Maybe this should be refactored out into some kind of separate package? Something like IntervalBranchAndBound.jl.

Actually, the SearchStrategy object has nothing to do with intervals, it is just a container for two methods and the RootSearch object itself could be generalized to any data. So it could be something shorter like BranchAndBound.jl. I think that such generalization indeed make it a convincing candidate to become a separate package.

It seems to me that the outputs field in RootSearchState does not need to be a container of any particular type -- a simple Vector should be enough for that?

Or, rather, it should be some kind of tree object that reflects how a subinterval(box) is related to a bigger interval(box).

That's interesting, the output container could indeed convey more informations. Since intervals and intervals boxes are spatially arrange relative to each other, representing them as a network of neighbors could be more useful than representing them as a tree (for example to help merge interval boxes with no solution when plotting them).

Do you have any particular use case in mind for a tree-like structure ?

[Kolaru]

Superseeded by #97.