Recommended way of lifting monotonic functions

[tpapp]

If I have a univariate monotonic function, what's the recommended way to define it for Interval? Is there a utility function for lifting? Eg

using ValidatedNumerics
oddsratio(x) = x==Inf ? one(x) : x/(1+x)
oddsratio(0.0)                  # 0.0
oddsratio(Inf)                  # 1.0
oddsratio(0..Inf)               # of course gives [0,∞], but we could do better

Naturally

oddsratio(i::Interval) = oddsratio(i.lo)..oddsratio(i.hi)

works but I am wondering if there is a shorthand (writing many functions).

[dpsanders]

There was actually an issue about this:
#198
but I closed it since it seemed unnecessary (!). I'll reopen it.

What do you think a good syntax would be? Just

@monotone oddsratio

to define the lifted Interval version?

[tpapp]

I would need to lift a callable object (currently done manually here) but I guess @monotone would also work with that, so it would be fine.

[dpsanders]

Hmmm, I'm not sure that this is as trivial as it seems. For example, what do you want to happen if someone calls your function with negative values? Maybe you need to specify the domain.

If you have many functions like this, you can use metaprogramming, something like

monotone_functions = (:oddsratio, :other_one)

for f in monotone_functions
  $f(i::Interval) = $f(i.lo)..$f(i.hi))
end 

However, note that the function as written will not give correctly rounded values.

[lbenet]

I guess that using try and catch is helpful here; for instance

julia> function monot_interval(f, a::Interval)
    res = try
        Interval(f(a.lo), f(a.hi))
    catch
        Interval(f(a.hi), f(a.lo))
    end
    return res
end

monot_interval (generic function with 1 method)

julia> monot_interval( x-> x^3, 0..2)
[0, 8]

julia> monot_interval( x-> -x^3, 0..2)
[-8, -0]

(Awful function name, I know, but you can change it.)

I guess this could be wrapped nicely in a macro.

[tpapp]

@dpsanders: of course I will add domain checks, but that is orthogonal.

How would I define a version with correct rounding, using a macro? I tried

"Apply a (monotone) transformation to an interval by endpoints."
macro monotone_interval(S)
    :((f::$S)(x::Interval) = ValidatedNumerics.@interval(f(x.lo), f(x.hi)))
end

but

julia>  macroexpand(quote @monotone_interval Logistic end)
quote  # REPL[17], line 1:
    (f::Logistic)(#13#x::Interval) = begin  # REPL[16], line 2:
            (ValidatedNumerics.Interval)(f(((ValidatedNumerics.convert)((ValidatedNumerics.esc
)(ValidatedNumerics.T)),#13#x.lo)).lo,f(((ValidatedNumerics.convert)((ValidatedNumerics.esc)(V
alidatedNumerics.T)),#13#x.hi)).hi)                                                          
        end
end

and those free T's are causing a problem. I am still lost when it comes to macro hygiene in Julia, but I suspect the problem is with @interval. Should I use a different interface?

(I would like to avoid using eval for now, and just have a macro that lifts).

[dpsanders]

Good question. Indeed it seems to be a macro hygiene issue in @interval. I'll try to fix this.

[dpsanders]

Does @interval(f(x.lo), f(x.hi)) work of it's not inside a macro?

[dpsanders]

Please also update to the latest version of the package that was released yesterday! (0.7)

[dpsanders]

With the latest version of the package, the following seems to solve your problem. Note the use of ::Type, which will clean up your code I think:

julia v0.5> immutable InvOddsRatio end

julia v0.5> (::Type{InvOddsRatio})(x::Number) = x == Inf ? one(x) : x/(1+x)

julia v0.5> InvOddsRatio(3)
0.75

julia v0.5> macro monotone_interval(S)
                :((::Type{$S})(x::Interval) = ValidatedNumerics.@interval($S(x.lo), $S(x.hi)))
            end

julia v0.5> @monotone_interval InvOddsRatio

julia v0.5> InvOddsRatio(3..4)
[0.75, 0.800001]

By the way, if the functions you need are simple one-liners, it would be possible to define a macro to do

@define_my_function InvOddsRatio x == Inf ? one(x) : x/(1+x)

and have it automatically generate everything!

[dpsanders]

Hmm, this is almost there but is giving a stack overflow. I'll have to think about this some more.

[tpapp]

My current solution is here. I wonder if you could look at it and tell me if I am doing the right thing.

Also, I imagine that rounding is overdone for expressions like x/(1+x), but I don't know how to deal with that, so I am planning to leave it as is for now.

[dpsanders]

I don't actually think that's completely right, since if you want to round down x/(1+x), you actually need (in principal) to round up y=1+x and then round down x/y.

The way to get around this problem is to use intervals. The best solution that I can think of is something like this:

f(x, dummy) = x == Inf ? one(x) : x/(1+x)
f(x::Number) = f(x, true)

f{T}(x::Interval{T}) = Interval(f(convert(Interval{T}, x.lo), true).lo, 
                                                f(convert(Interval{T}, x.hi), true).hi)

The idea is that the correct rounding is done by passing in an Interval to the function; the dummy variable is just to make a separate method to dispatch to:

julia> x = convert(Interval{Float64}, 3.1)
[3.09999, 3.10001]

julia> f(x, true)   # gives correct rounding
[0.756097, 0.756098]

We pass in separately the lo and hi components as intervals to round them correctly, and create the Interval from the result.

Compare the following: your method and my method give slightly different results:

julia> setrounding(Float64, RoundDown) do
           f(0.1, true)
       end
0.09090909090909091

julia> x = convert(Interval{Float64}, 0.1)
[0.0999999, 0.100001]

julia> f(x, true)
[0.090909, 0.0909091]

julia> setdisplay(:full)
6

julia> f(x, true)
Interval(0.09090909090909088, 0.09090909090909093)

I have a guaranteed lower and upper bound, whereas with your method I'm not sure if it's actually a lower bound. Indeed, it actually isn't:

julia> f(big"0.1", true)
9.09090909090909090909090909090909090909090909090909090909090909090909090909085e-02

[tpapp]

Thanks! I found that I can use invoke to call the right function:

using ValidatedNumerics

f(x::Real) = x == Inf ? one(x) : x/(1+x)
function f_singleton{T <: Real}(x::T)
    println("singleton interval around $x")
    invoke(f, Tuple{Real}, convert(Interval{T}, x))
end
f(x::Interval) = Interval(f_singleton(x.lo).lo, f_singleton(x.hi).hi)

and then

julia>  f(1.0)
0.5

julia>  f(3.1)
0.7560975609756099

julia>  f(Interval(3.1))
singleton interval around 3.1
singleton interval around 3.1
[0.756097, 0.756098]

I am still experimenting with this so I may have missed some corner case of dispatch, but this seems robust. What do you think?

[dpsanders]

I had thought of using invoke, but somehow I don't like it aesthetically. But it's a perfectly good solution, yes!

[dpsanders]

It's kind of annoying that it's necessary to write 3 methods for each function, but I don't see a way around that. And these functions can definitely be auto-generated with metaprogramming.

[dpsanders]

I'm not sure if there's anything to add to ValidatedNumerics.jl for this.

We could add a @monotone macro, but it's not clear to me how we could check if the function passed to it is actually monotone, since this requires bounding the derivative, which in general is a hard problem.

[tpapp]

I have realized that ValidatedNumerics.jl may not be the best match for what I am trying to do, so I will be implementing this differently. Thanks for all the help, closing the issue.

[dpsanders]

OK, no problem. What is it that you're trying to do? If you just need intervals as sets, and you don't need to worry about all the rounding issues, then you might want to look at
https://github.com/JuliaMath/IntervalSets.jl

Hopefully with traits it will be possible at some point to unify some of these ideas.

[lbenet]

The rigorous way of checking that a function is monotonic involves some tricks which involve the derivatives, with methods rather well established. I thought the idea was that the user takes responsibility, something like using @inbounds.

[dpsanders]

Yes of course if we delegate responsibility to the user, that is much easier ;) Do you have a reference for those techniques? I don't see that it would ever be possible to prove that eg. x/(1+x) is monotone over the real line.

[lbenet]

I guess one has to be careful around -1, but elsewhere not including it, I think it should be possible. If I have some time this afternoon, I work on this.

[dpsanders]

Actually that particular function I believe is only supposed to be for positive x.

[lbenet]

Then there should be no problem...

[tpapp]

@dpsanders answering the "what am I trying to do" question: initially I wanted "intervals as sets", then read about validated numerics and thought I could get lots of useful things like error analysis for free. Realized that I would like to dispatch on intervals being finite or infinite (one and two-sided), so I wrote yet another way of dealing with the problem. Hopefully with traits one day we can unify, but as things are at the moment working with either IntervalSets.jl or ValidatedNumerics.jl would require an outrageous amount of type piracy.