JuliaIntervals / IntervalArithmetic.jl

Library for validated numerics using interval arithmetic
https://juliaintervals.github.io/IntervalArithmetic.jl/
Other
298 stars 71 forks source link

What to do with mod? #129

Open dpsanders opened 6 years ago

dpsanders commented 6 years ago

How should we deal with functions like mod acting on intervals?

E.g.

One solution is to define an interval union type; see e.g. https://link.springer.com/article/10.1007/s10543-016-0632-y

Really the correct solution is probably to use contractors.

lbenet commented 6 years ago

Speaking about Interval, I think it should behave as tan. For the specific example you give, I think it should return 0..1, which corresponds to the hull.

A complementary possibility is to have mod(a::Interval) returning the hull, and also implement extended_mod(a::Interval), among other examples (e.g. extended_inv(a::Interval)).

extended_mod(a) would return a Vector{Interval) (or perhaps NTuple{N,Interval}) containing the list of (connected) intervals, as a way of representing the union in the sense you use above. If the interval a does not contain the discontinuity, the result would correspond to [mod(a)].

Perhaps, this is what you mean by using contractors.

lbenet commented 6 years ago

In #145 I have implemented something that is equivalent to mod(0.8..2.2, 1) --> 0..1.

We could also implement some extended_mod which would return a tuple of intervals (two or three, this we should discuss) which somehow mimics extended_div.

dpsanders commented 6 years ago

Great, thanks. Why would there be three intervals?

Is this version of mod (which just returns 0..1) useful? I tried it with root finding and it didn't seem to be too useful.

On Thu, 3 May 2018, 3:01 pm Luis Benet, notifications@github.com wrote:

In #145 https://github.com/JuliaIntervals/IntervalArithmetic.jl/pull/145 I have implemented something that is equivalent to mod(0.8..2.2, 1) --> 0..1.

We could also implement some extended_mod which would return a tuple of intervals (two or three, this we should discuss) which somehow mimics extended_div.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/129#issuecomment-386402374, or mute the thread https://github.com/notifications/unsubscribe-auth/AALtTniH72rkhmw7slmA8eqx0WR-1EQJks5tu1QFgaJpZM4TXHeO .

lbenet commented 6 years ago

Why would there be three intervals?

Considering your example, mod(0.8..2.2, 1), I would say that extended_mod would return (0.8..1.0) ∪ (0.0..1.0) ∪(0.0..0.2), if we consider each segment.

Is this version of mod (which just returns 0..1) useful?

This may be useful for things like find_quadrants; for root-finding methods, we would need extended_mod. I'll push a new commit dealing with extended_mod later.

lbenet commented 6 years ago

Just pushed another commit to #145 with extended_mod implemented. Can you retest your example with root finding, using now extended_mod? Maybe you could post it so I can play with it.

dpsanders commented 6 years ago

The problem is that all of the root finding functions assume that a single box in gives a single box out. Basically extended mod etc lead to interval unions, which is not yet implemented.

On Fri, 4 May 2018, 11:49 am Luis Benet, notifications@github.com wrote:

Just pushed another commit to #145 https://github.com/JuliaIntervals/IntervalArithmetic.jl/pull/145 with extended_mod implemented. Can you retest your example with root finding, using now extended_mod? Maybe you could post it so I can play with it.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/129#issuecomment-386643737, or mute the thread https://github.com/notifications/unsubscribe-auth/AALtTnG9ekSMND3_12RR-cNS2dhdwRIDks5tvHihgaJpZM4TXHeO .

lbenet commented 6 years ago

I'm not sure I understand your point.

I agree that extended_mod returns an interval union (as a tuple), which was precisely the suggestion made above. It mimics what Newton method does by extending the division, thus returning a tuple of disjoint intervals, if the derivative contains a zero. Similarly, extended_mod returns a tuple; each interval corresponds to the range of mod when the domain has no discontinuities.

dpsanders commented 6 years ago

I just mean that the current versions of the algorithms can't deal with interval unions.

Datseris commented 1 year ago

Is there any update here, or at least instructions of how to handle mods? I have a simple dynamical system that involves the mod operation

function dissipative_standard_map_rule(u, p, n)
    x, y = u
    (; ν, f₀) = p
    s = x + y
    xn = mod2pi(s)
    yn = (1 - ν)*y + f₀*(sin(s))
    return SVector(xn, yn)
end

and I would like to find its fixed points using interval arithmetic (root finding) via utilizing Newton's method (I am using ChaosTools.fixedpoints). However, inevitably, the mod is called on an interval and I get

ERROR: MethodError: no method matching rem2pi(::IntervalArithmetic.Interval{Float64}, ::RoundingMode{:Down})

Since I am anyways only searching for roots for x bounded in [0, 2pi], would you say it is correct to simply remove the mod operation from the above function?

lbenet commented 1 year ago

I guess the work you are looking for is in #145, but somehow we got distracted with other stuff...

Datseris commented 1 year ago

Seems that https://github.com/JuliaIntervals/IntervalArithmetic.jl/pull/525 is much more recent and perhaps could be merged with relatively little effort? I would agree with the OP that the (interval, real) method is the most common.

Kolaru commented 1 year ago

Since I am anyways only searching for roots for x bounded in [0, 2pi], would you say it is correct to simply remove the mod operation from the above function?

To me it looks like you could use mod at the very end of the computation, or not at all, since for what matters it goes into sin anyway. Feels like the mod2pi is here for numerical stability.

I merged #525 for now. It will not remove the error directly, because it doesn't add mod2pi. It can be avoided with mod(x, 2pi). Also it may cause a new problem : if you get an interval that overlap with one of the bound, the mod will get you the full range (e.g. mod(-0.0001 .. 0.0002, 2pi) == (0..2pi)),

What you want may rather be something like xn = mod2pi(mid(s)) ± radius(s), which doesn't guarantee inclusion in the range [0, 2pi], but will keep xn close enough, without risking a huge increase in the interval radius. Disclaimer : I haven't check if the formula correctly include the numerical error, I would need a bit more time to come up with something that guarantee it.

Datseris commented 1 year ago

Thank you yes. In the meantime I indeed went ahead and removed the mod from my equations but unfortunately in the end of the day it didn't help me, turns out the Jacobian matrix isn't invertible so it was all for naught xD But thanks a lot for progressing with the mod!