dpsanders commented 6 years ago

@Kolaru: I'm having trouble with the complex functions...

dpsanders commented 6 years ago

@Kolaru: I think I've fixed complex functions. Would you please be able to take a look / review? Don't worry if not.

dpsanders commented 6 years ago

Great comments, thanks a lot! I may not be able to get to them today.

lbenet commented 6 years ago

Do you understand the following (current) behavior?

julia> using IntervalRootFinding, IntervalArithmetic

julia> find_roots(sin, interval(-4,4), IntervalRootFinding.newton)
3-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([-3.1416, -3.14159], :unique)        
 Root([-3.23118e-27, 3.23118e-27], :unique)
 Root([3.14159, 3.1416], :unique)          

julia> roots(sin, interval(-4,4), Newton)
4-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([3.14112, 3.14176], :unique)  
 Root([0, 0.000325267], :unknown)   
 Root([-0.000325267, 0], :unknown)  
 Root([-3.14176, -3.14112], :unique)

Note that using find_roots, three :unique Roots are obtained; using roots four are given, among which two are :unknown (and correspond to the root at 0).

Kolaru commented 6 years ago

It is because roots use bisect which does not use guarded_mid (instead it uses mid). The first iteration bisects exactly on zero, causing the problem. This can be seen by giving a non symmetric interval:

julia> roots(sin, -4..4.0001, Newton, eps(1.0))
3-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([3.14159, 3.1416], :unique)          
 Root([-7.11885e-17, 5.26982e-18], :unique)
 Root([-3.1416, -3.14159], :unique)

Note that shifting the center when bisecting does not totally solve the problem since we also have (same problem with find_roots):

julia> roots(sin, 0..4, Newton, eps(1.0))
2-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([3.14159, 3.1416], :unique)
 Root([0, 7.40155e-17], :unknown)

Also the default tolerance for roots and find_roots are not the same.

Speaking about tolerance, it is a bit odd that the default tolerance is given per function, especially since the value is always the same. Defining it globally would probably be better (either via a constant or with some mechanics which allows the user to change it globally).

lbenet commented 6 years ago

Also the default tolerance for roots and find_roots are not the same.

Good point!

It is because roots use bisect which does not use guarded_mid (instead it uses mid). The first iteration bisects exactly on zero, causing the problem.

Does it make sense to propose/suggest to use guarded_mid in roots?

Kolaru commented 6 years ago

Does it make sense to propose/suggest to use guarded_mid in roots?

I do not think so, as it would introduce a lot of complications. In order to work guarded_mid need to know the function used to test if the center is a root, which will be tricky to implement for multi-dimensional functions. In fact testing the center of an interval will not help at all for multi-dimensional function, unless the presence of a solution is tested on the whole border.

Some possible general fixes are:

Always shift the middle. It makes the convergence a bit worse but solve the problem for "obvious" roots (like 0 or 1).
Makes the documentation clear that roots (with Newton method) will not find roots on the border of the given Interval/IntervalBox and that roots of the form [Root([a, b], :unkown), Root([b, c], :unkown)] may indicate that b is a root.
Check explicitly for the latter and try to automate it.
Solve this at the theoretical level (I may have a shot at this if I have time to read the references in detail) either by providing other method or finding a way to characterize roots on interval borders.

dpsanders commented 6 years ago

Should the Root type have its own file ?

Sure, why not. I'm also not sure if having a Root type is a good solution. Since the only two possibilities are :unique and :unknown, it's more common to just return two arrays, one with the intervals which give :unique and another with those that are :unknown.

Is the name deriv better than f_prime for the derivative/jacobian ? (I prefer the latter but I may only be a matter of taste)

I don't mind too much. It's just a bit strange to write f_prime if the function is not called f.

dpsanders commented 6 years ago

There is no way that interval methods can locate roots that fall precisely on the boundary of the interval, intuitively since it can "never know what happens on the either side of that boundary". I think that #43 is the simplest solution for this problem.

lbenet commented 6 years ago

I'm also not sure if having a Root type is a good solution. Since the only two possibilities are :unique and :unknown, it's more common to just return two arrays, one with the intervals which give :unique and another with those that are :unknown.

I am actually in favor of having Root as a type that encodes the returned roots. I think it is a nice way of having in the same object the actual interval having the (plausible) root and the info about it (:unique or :unknown). Maybe we should consider to have, instead of returning a symbol, a closure. And obviously we should add lots of methods to deal comfortably with Roots.

lbenet commented 6 years ago

Related to the case of having a root precisely at the border of the interval, which is one of the cases we label with :unknown, we could document that the user should check whether this is the case by its own, or provide some methods to do so.

EDIT: Changed the label of a Root at the border

dpsanders commented 6 years ago

A root at the border is labelled :unknown.

lbenet commented 6 years ago

Ups... I meant so; I just edited my comment

Kolaru commented 6 years ago

I don't mind too much. It's just a bit strange to write f_prime if the function is not called f.

You are right, f_prime would only look nice in the documentation. Another alternative would be to use jacobian or J.

I'm also not sure if having a Root type is a good solution.

Keeping a Root type is probably required to (elegantly) solve #29, since the information that the root is unique need to be given one way or another.

I also reviewed the last commits and except for the failing newton1d test I think that everything is fine.

dpsanders commented 6 years ago

Regarding the name of the keyword, I have left it as deriv for the moment, since that seems to be the most generic. We can change it later if we find a better alternative.

Will merge once green.

codecov-io commented 6 years ago

Codecov Report

Merging #55 into master will decrease coverage by 36.18%. The diff coverage is 73.4%.

@@             Coverage Diff             @@
##           master      #55       +/-   ##
===========================================
- Coverage   80.53%   44.34%   -36.19%     
===========================================
  Files           8        9        +1     
  Lines         339      327       -12     
===========================================
- Hits          273      145      -128     
- Misses         66      182      +116

Impacted Files	Coverage Δ
src/IntervalRootFinding.jl	`2.04% <ø> (-86.85%)`	:arrow_down:
src/bisect.jl	`100% <ø> (ø)`	:arrow_up:
src/newton.jl	`0% <0%> (-68.5%)`	:arrow_down:
src/krawczyk.jl	`0% <0%> (-95.46%)`	:arrow_down:
src/root_object.jl	`12.5% <12.5%> (ø)`
src/roots.jl	`81.63% <81.63%> (ø)`
src/contractors.jl	`93.33% <93.33%> (ø)`
... and 3 more

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 14d97bf...5f2aebc. Read the comment docs.

lbenet commented 6 years ago

Thanks for all the work in this PR!

JuliaIntervals / IntervalRootFinding.jl

Refactor and cleanup. Removes NewtonContractor #55

Codecov Report