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

Kolaru · 2018-05-30T23:08:18Z

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 · 2018-05-31T20:03:56Z

src/roots.jl

-    outputs::Vector{Root{T}}
+struct RootSearchState{V, VR}
+    working::V  # Should be a container of the form  CONT{T}
+    outputs::VR  # Should ba a container of root of the form CONT{Root{T}}


I don't think it's necessary for outputs to be a fancy container -- isn't a standard list (vector) good enough, since there won't be complicated manipulations here?

dpsanders · 2018-05-31T20:05:16Z

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 · 2018-05-31T20:08:46Z

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 · 2018-06-02T13:44:07Z

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 · 2018-06-04T19:26:01Z

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 · 2018-06-07T11:55:35Z

src/roots.jl

-    sizehint!(state.outputs, 100)
-    sizehint!(state.working, 1000)
-    return state
+    return iter.state


This implementation does not conform to the Base.Iterator specification: start(iter) should return the initial state of the iterator; not the last state that iterator was in. I can imagine this leading to problems down the line when someone tries to restart a RootSearch after it has run for a few iterations. Ideally, start() should create a fresh state. (Resetting state would probably also lead to some confusion because of side effects.)

gwater · 2018-06-07T12:43:44Z

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 · 2018-06-07T14:30:13Z

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 · 2018-06-08T08:36:38Z

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

codecov-io · 2018-07-05T14:50:10Z

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% <ø> (ø)`	⬆️
src/roots.jl	`75% <84.21%> (-13.14%)`	⬇️
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 · 2018-07-05T14:58:42Z

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 · 2018-07-09T09:08:12Z

@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 · 2018-07-09T12:35:45Z

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 · 2018-07-09T13:02:56Z

@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 · 2018-07-12T13:00:26Z

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 · 2018-07-17T12:31:19Z

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 · 2018-07-17T12:35:53Z

src/roots.jl

@@ -119,100 +54,162 @@ IntervalLike{T} = Union{Interval{T}, IntervalBox{T}}
 NewtonLike = Union{Type{Newton}, Type{Krawczyk}}

 """
-    roots(f, X, contractor, tol=1e-3)
-    roots(f, deriv, X, contractor, tol=1e-3)
+    roots(f, X, contractor=Newton, strategy=BreadthFirstSearch, tol=1e-15)


In line L80 the same value (for tolerance) is defined as a constant. It might be worth moving the constant up to reuse it here (and on the following line).

gwater · 2018-07-17T12:40:27Z

src/rootsearch_iterator.jl

+is always modified. If a root is found, it is added to `state.outputs`.
+"""
+function step!(state::RootSearchState, contractor, searchstrat, tolerance)
+    X = searchstrat.retrieve!(state.working)


This line introduces a type instability because julia cannot infer the return type of retrieve!() for some reason. Three options:

put the remaining code into a separate function which is then fully inferred

ignore the issue for now because the heavy lifting is done in contractor() which should be fully inferrible

somehow fix type inferrence for retrieve!() (any ideas?)

For your convenience; this is how I tested option 1:

function step_inner!(X, working, outputs, store!, contractor, tolerance) status, output = contractor(X, tolerance) if status == :empty return nothing elseif status == :unique store!(outputs, Root(output, :unique)) elseif diam(output) < tolerance store!(outputs, Root(output, :unknown)) else # branch X1, X2 = bisect(X) store!(working, X1, X2) end return nothing end """ step!(state::RootSearchState, contractor, tolerance) Progress `state` by treating one of its `working` regions. Note: `state.working` is always modified. If a root is found, it is added to `state.outputs`. """ function step!(state::RootSearchState, contractor, searchstrat, tolerance) X = searchstrat.retrieve!(state.working) step_inner!(X, state.working, state.outputs, searchstrat.store!, contractor, tolerance) return nothing end

Here, step!() is still unstable but step_inner!() is fully inferred.

Nice catch! Thinking about it, this is probably why each function has its own type: parametrization on function seems to allow to infer returned type in struct body. Consider the following

struct A func::Function end struct B{F} func::F end a = A(x -> x^2) b = B(x -> x^2) test(cont, x) = cont.func(x)

In the case of test(a, 3.2) the return type can not be inferred, but for test(b, 3.2) it can. So this is the way to go.

Yes, functions inside structs must "always" be parametrized.

gwater · 2018-07-17T12:44:32Z

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)

Allow to have type predictability for the full iteration cycle.

- Introduce default strategies - Rework interface to allow using them - Changes break tests for complex Krawczyk method which expectation has been lowered

dpsanders · 2018-08-19T16:22:59Z

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

dpsanders · 2018-08-21T06:55:35Z

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

Two general comments:

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
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 · 2018-08-21T15:11:54Z

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 · 2018-08-21T22:35:53Z

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 · 2019-03-08T01:31:58Z

Superseeded by #97.

dpsanders reviewed May 31, 2018

View reviewed changes

gwater reviewed Jun 7, 2018

View reviewed changes

Kolaru force-pushed the simple_custom_strat branch from 343f7de to 5fccaac Compare July 12, 2018 12:38

Kolaru changed the title ~~[WIP] Allow control over the way the intervals are selected to be processed~~ Allow control over the way the intervals are selected to be processed Jul 12, 2018

gwater reviewed Jul 17, 2018

View reviewed changes

Kolaru added 9 commits July 19, 2018 14:52

Implement SearchStrategy for root finding

16a0a92

Fix bug for structures without sizehint!

dd333cb

Decouple the iterator and its state

cf0570e

Parametrize RootSearchState on container type

709c83a

Allow to have type predictability for the full iteration cycle.

Redefine RootSearch with parametric SearchStrategy

8a1f1d0

Create separate file for the RootSearch iterator

0a9218a

Rewrite documentation for RootSearch related objects

38f0982

Rework roots interface to allow choosing a strategy

dcc2bc1

- Introduce default strategies - Rework interface to allow using them - Changes break tests for complex Krawczyk method which expectation has been lowered

Add simple test for SearchStrategy

f6b959c

Fix type instability and some other minor issues

3eadb00

Kolaru force-pushed the simple_custom_strat branch from 5fccaac to 3eadb00 Compare July 19, 2018 13:31

Kolaru mentioned this pull request Sep 29, 2018

Separated the logic of branch and bound search #97

Merged

Kolaru closed this Mar 8, 2019

Kolaru deleted the simple_custom_strat branch March 8, 2019 01:32

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

Kolaru commented May 30, 2018

dpsanders May 31, 2018

dpsanders commented May 31, 2018

dpsanders commented May 31, 2018 •

edited

Loading

Kolaru commented Jun 2, 2018

Kolaru commented Jun 4, 2018

gwater Jun 7, 2018 •

edited

Loading

gwater commented Jun 7, 2018

Kolaru commented Jun 7, 2018

gwater commented Jun 8, 2018

codecov-io commented Jul 5, 2018 •

edited

Loading

Kolaru commented Jul 5, 2018 •

edited

Loading

gwater commented Jul 9, 2018

gwater commented Jul 9, 2018

Kolaru commented Jul 9, 2018

Kolaru commented Jul 12, 2018 •

edited

Loading

gwater commented Jul 17, 2018

gwater Jul 17, 2018 •

edited

Loading

gwater Jul 17, 2018 •

edited

Loading

gwater Jul 17, 2018

Kolaru Jul 17, 2018

dpsanders Aug 19, 2018

gwater commented Jul 17, 2018

dpsanders commented Aug 19, 2018

dpsanders commented Aug 21, 2018 •

edited

Loading

dpsanders commented Aug 21, 2018

Kolaru commented Aug 21, 2018

Kolaru commented Mar 8, 2019

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

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

Conversation

Kolaru commented May 30, 2018

dpsanders May 31, 2018

Choose a reason for hiding this comment

dpsanders commented May 31, 2018

dpsanders commented May 31, 2018 • edited Loading

Kolaru commented Jun 2, 2018

Kolaru commented Jun 4, 2018

gwater Jun 7, 2018 • edited Loading

Choose a reason for hiding this comment

gwater commented Jun 7, 2018

Kolaru commented Jun 7, 2018

gwater commented Jun 8, 2018

codecov-io commented Jul 5, 2018 • edited Loading

Codecov Report

Kolaru commented Jul 5, 2018 • edited Loading

gwater commented Jul 9, 2018

gwater commented Jul 9, 2018

Kolaru commented Jul 9, 2018

Kolaru commented Jul 12, 2018 • edited Loading

gwater commented Jul 17, 2018

gwater Jul 17, 2018 • edited Loading

Choose a reason for hiding this comment

gwater Jul 17, 2018 • edited Loading

Choose a reason for hiding this comment

gwater Jul 17, 2018

Choose a reason for hiding this comment

Kolaru Jul 17, 2018

Choose a reason for hiding this comment

dpsanders Aug 19, 2018

Choose a reason for hiding this comment

gwater commented Jul 17, 2018

dpsanders commented Aug 19, 2018

dpsanders commented Aug 21, 2018 • edited Loading

dpsanders commented Aug 21, 2018

Kolaru commented Aug 21, 2018

Kolaru commented Mar 8, 2019

dpsanders commented May 31, 2018 •

edited

Loading

gwater Jun 7, 2018 •

edited

Loading

codecov-io commented Jul 5, 2018 •

edited

Loading

Kolaru commented Jul 5, 2018 •

edited

Loading

Kolaru commented Jul 12, 2018 •

edited

Loading

gwater Jul 17, 2018 •

edited

Loading

gwater Jul 17, 2018 •

edited

Loading

dpsanders commented Aug 21, 2018 •

edited

Loading