Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Integer overflow in innocent looking linear interpolation. Should some warning be issued? #457

Closed
hurak opened this issue Aug 28, 2021 · 24 comments · Fixed by timholy/Ratios.jl#24

Comments

@hurak
Copy link
Contributor

hurak commented Aug 28, 2021

Overflow is encountered during and innocent looking linear interpolation/extrapolation:

Interpolate linearly two points (1,1) and (10000,10000) and then find the value at 10001 by linear extrapolation.
The following code does it

julia> itpi = LinearInterpolation([1,10000],[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Int64}, Gridded(Linear())), Line()) with element type Float64:
 Ratios.SimpleRatio{Int64}(99980001, 99980001)
 Ratios.SimpleRatio{Int64}(999800010000, 99980001)

julia> convert(Float64,itpi(10001))
773.9376918087789

The result is clearly incorrect because the correct answer is 10001. Looking at the displayed ratios, I guess that an integer overflow is what happened here. There are no problems with smaller numbers:

julia> itp = LinearInterpolation([1,2],[1,2]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Int64}, Gridded(Linear())), Line()) with element type Float64:
 Ratios.SimpleRatio{Int64}(1, 1)
 Ratios.SimpleRatio{Int64}(2, 1)

julia> convert(Float64,itp(3))
3.0

or with the floating point data:

julia> itpf = LinearInterpolation([1.0,10000.0],[1.0,10000.0]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Line()) with element type Float64:
     1.0
 10000.0

julia> convert(Float64,itpf(10001.0))
10001.0

I am hesitating if this should be regarded (and filed) as an Issue. In general the responsibility for the overflow problems is on the application programmer, right? But shouldn’t the interpolation function give some warning here? The parameters of the problem look innocent, don’t they?

@mkitti
Copy link
Collaborator

mkitti commented Aug 30, 2021

If there is a bug here, it is that overflow occurs for such a simple case. If the user supplies an integer type, we stick with that integer type for the computation to take advantage of any performance gains along with any potential issues such as overflow.

Let's distill the bug to be a failure to simplify the SimpleRatio which makes us prematurely vulnerable to overflow in this case.

@timholy
Copy link
Member

timholy commented Aug 30, 2021

Can someone check whether it would be fixed with timholy/Ratios.jl#19 (comment)?

@mkitti
Copy link
Collaborator

mkitti commented Aug 30, 2021

Can someone check whether it would be fixed with timholy/Ratios.jl#19 (comment)?

julia> using Interpolations, Ratios

julia> itpi = LinearInterpolation([1,10000],[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Int64}, Gridded(Linear())), Line()) with element type Float64:
 SimpleRatio{Int64}(99980001, 99980001)
 SimpleRatio{Int64}(999800010000, 99980001)

julia> itpi(10001)
SimpleRatio{Int64}(7736281631652211921, 9996000599960001)

julia> convert(Float64, itpi(10001))
773.9376918087789

julia> Base.:+(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num + y.num, x.den) : SimpleRatio(x.num*y.den + x.den*y.num, x.den*y.den)

julia> Base.:-(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num - y.num, x.den) : SimpleRatio(x.num*y.den - x.den*y.num, x.den*y.den)

julia> itpi = LinearInterpolation([1,10000],[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Int64}, Gridded(Linear())), Line()) with element type Float64:
 SimpleRatio{Int64}(9999, 9999)
 SimpleRatio{Int64}(99990000, 9999)

julia> convert(Float64, itpi(10001))
10001.0

@mkitti
Copy link
Collaborator

mkitti commented Aug 30, 2021

Within Interpolations.jl, my thought would be to simplify the ratio before construction:

julia> using Interpolations, Ratios

julia> Interpolations.ratio(num::Integer, denom::Integer) = SimpleRatio(Base.divgcd(promote(num, denom)...)...)

julia> itpi = LinearInterpolation([1,10000],[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Int64}, Gridded(Linear())), Line()) with element type Float64:
 SimpleRatio{Int64}(1, 1)
 SimpleRatio{Int64}(10000, 1)

julia> itpi(10001)
SimpleRatio{Int64}(999899990001, 99980001)

julia> convert(Float64, itpi(10001))
10001.0

This also combines well with the above fix for Ratios.jl:

julia> Base.:+(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num + y.num, x.den) : SimpleRatio(x.num*y.den + x.den*y.num, x.den*y.den)

julia> Base.:-(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num - y.num, x.den) : SimpleRatio(x.num*y.den - x.den*y.num, x.den*y.den)

julia> itpi(10001)
SimpleRatio{Int64}(99999999, 9999)

Perhaps a simplify method would be generally useful?

julia> simplify(r::SimpleRatio) = (d = gcd(r.num, r.den); SimpleRatio(r.num ÷ d, r.den ÷ d))
simplify (generic function with 1 method)

julia> itpi(10001)
SimpleRatio{Int64}(99999999, 9999)

julia> simplify(itpi(10001))
SimpleRatio{Int64}(10001, 1)

At some point, it seems we are patching SimpleRatio until it becomes Rational.

@timholy
Copy link
Member

timholy commented Aug 31, 2021

At some point, it seems we are patching SimpleRatio until it becomes Rational.

Yes, and that would be a huge problem.

julia> using BenchmarkTools

julia> @btime Base.divgcd(7736281631652211921, 9996000599960001)
  127.862 ns (0 allocations: 0 bytes)
(7736281631652211921, 9996000599960001)

vs

julia> using Interpolations

julia> using Interpolations.Ratios

julia> Base.:+(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num + y.num, x.den) : SimpleRatio(x.num*y.den + x.den*y.num, x.den*y.den)

julia> Base.:-(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num - y.num, x.den) : SimpleRatio(x.num*y.den - x.den*y.num, x.den*y.den)

julia> itpi = LinearInterpolation([1,10000],[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Int64}, Gridded(Linear())), Line()) with element type Float64:
 SimpleRatio{Int64}(9999, 9999)
 SimpleRatio{Int64}(99990000, 9999)

julia> @btime $itpi(10001)
  7.331 ns (0 allocations: 0 bytes)
SimpleRatio{Int64}(99999999, 9999)

Do you see why I suggested that strategy first?

@mkitti
Copy link
Collaborator

mkitti commented Sep 1, 2021

Do you see why I suggested that strategy first?

It seems the denominator check seems to address this specific issue and is not too onerous. Do you want me to create an issue over at Ratios.jl? Would ifelse be better than a ternary operator to avoid branching logic?

@mkitti
Copy link
Collaborator

mkitti commented Sep 1, 2021

Regarding a warning or error that overflow occurred, perhaps working on the composition of this package, Ratios.jl, and SaferIntegers.jl would be useful.

Here's my attempt which requires some hot patching at the moment:

julia> using SaferIntegers, Interpolations, Ratios

julia> import SaferIntegers: ndigits0z, baseint, SaferIntegers

julia> for S in (:SafeSigned, :SafeUnsigned)
         @eval begin
               ndigits0z(x::T) where T<:$S = ndigits0z(baseint(x)) # do not reconvert
           ndigits0z(x::T1, b::T2) where T1<:$S where T2<:SafeInteger = ndigits0z(baseint(x), baseint(b)) # do not reconvert
           ndigits0z(x::T1, b::T2) where T1<:$S where T2<:Integer = ndigits0z(baseint(x), b) # do not reconvert
         end
       end

julia> Base.float(x::SafeInteger) = Base.float(SaferIntegers.baseint(x))

julia> Base.:-(x::SimpleRatio{T}) where {T<:SafeSigned} = SimpleRatio(-x.num, x.den)

julia> itpi = LinearInterpolation(SafeInt[1,10000],SafeInt[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{SafeInt64},), ::Vector{SafeInt64}, Gridded(Linear())), Line()) with element type Float64:
 SimpleRatio{SafeInt64}(99980001, 99980001)
 SimpleRatio{SafeInt64}(999800010000, 99980001)

julia> itpi(10001)
ERROR: OverflowError: 999800010000 * 99980001 overflowed for type Int64
Stacktrace:
 [1] throw_overflowerr_binaryop(op::Symbol, x::Int64, y::Int64)
   @ Base.Checked .\checked.jl:154
 [2] checked_mul
   @ .\checked.jl:288 [inlined]
 [3] *
   @ ~\.julia\packages\SaferIntegers\nOUgY\src\arith_ops.jl:21 [inlined]
 [4] +
   @ ~\.julia\packages\Ratios\uRs4y\src\Ratios.jl:31 [inlined]
 [5] extrapolate_axis
   @ ~\.julia\dev\Interpolations\src\extrapolation\extrapolation.jl:168 [inlined]
 [6] extrapolate_value
   @ ~\.julia\dev\Interpolations\src\extrapolation\extrapolation.jl:162 [inlined]
 [7] (::Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, SafeInt64, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{SafeInt64}}}, Gridded{Linear{Throw{OnGrid}}}, Line{Nothing}})(x::Int64)
   @ Interpolations ~\.julia\dev\Interpolations\src\extrapolation\extrapolation.jl:52
 [8] top-level scope
   @ REPL[7]:1

@timholy
Copy link
Member

timholy commented Sep 1, 2021

I think before making any decisions, all the options should be benchmarked, including against the (buggy) implementation we have now. SaferIntegers is appealing in principle, but what are the performance implications? Interpolations is an extraordinarily performance-sensitive package.

@mkitti
Copy link
Collaborator

mkitti commented Sep 1, 2021

The SaferIntegers approach is orthogonal to the other approaches. Once we provide support for it, the user can make a choice whether to use a SafeInteger to detect an overflow or not.

The gcd approach is probably better handled by Rational. Perhaps we should have a Preferences.jl mechanism to use Rational instead of SimpleRatio when performance is not the primary need of the user?

To summarize, I think we have three orthogonal approaches:

  1. Modify Ratios.jl to check if SimpleRatios has the same denominator during addition and subtraction. The added overhead of that check should be benchmarked carefully. Based on my informal experiments and the above times, the overhead looks reasonable to me.

  2. Modify SaferIntegers.jl and Ratios.jl to work together with Interpolations.jl. While benchmarks would be informative, I do not think this influences whether we should do this or not since the choice to use SaferIntegers.jl lies with the user.

  3. Modify Interpolations.jl to use Preferences.jl to modify the implementation of Interpolations.ratio between SimpleRatio or Rational.

@JeffreySarnoff
Copy link
Member

How would SaferIntegers.jl need to be modified to work with Interpolations.jl?

@mkitti
Copy link
Collaborator

mkitti commented Sep 2, 2021

How would SaferIntegers.jl need to be modified to work with Interpolations.jl?

Hi @JeffreySarnoff , adding Base.float for SafeInteger would help.

Base.float(x::SafeInteger) = Base.float(SaferIntegers.baseint(x))

See JeffreySarnoff/SaferIntegers.jl#19

@JeffreySarnoff
Copy link
Member

I incorporated this

Base.float(x::S) where S<:SafeSigned = Base.float(baseint(x))
Base.float(x::U) where U<:SafeUnsigned = Base.float(baseint(x))

into SaferIntegers/src/promote.jl, bumped the version and reregistered it
(I still need to add tests)

@mkitti
Copy link
Collaborator

mkitti commented Sep 3, 2021

On approach #2, we are then here:

julia> using Interpolations, SaferIntegers

julia> itpi = LinearInterpolation(SafeInt[1,10000],SafeInt[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{SafeInt64},), ::Vector{SafeInt64}, Gridded(Linear())), Line()) with element type Float64:
 Ratios.SimpleRatio{SafeInt64}(99980001, 99980001)
 Ratios.SimpleRatio{SafeInt64}(999800010000, 99980001)

julia> itpi(10001)
ERROR: MethodError: no method matching -(::Ratios.SimpleRatio{SafeInt64})
...

We just need to add the missing method to Ratios now since SafeInt64 is not Signed, but SafeSigned.

julia> using Interpolations.Ratios

julia> Base.:-(x::SimpleRatio{T}) where {T<:SafeSigned} = SimpleRatio(-x.num, x.den)

julia> itpi(10001)
ERROR: OverflowError: 999800010000 * 99980001 overflowed for type Int64
Stacktrace:
 [1] throw_overflowerr_binaryop(op::Symbol, x::Int64, y::Int64)
   @ Base.Checked .\checked.jl:154
 [2] checked_mul
   @ .\checked.jl:288 [inlined]
 [3] *
   @ ~\.julia\packages\SaferIntegers\zbJSp\src\arith_ops.jl:21 [inlined]
 [4] +
   @ ~\.julia\packages\Ratios\fgO34\src\Ratios.jl:35 [inlined]
 [5] extrapolate_axis
   @ ~\.julia\dev\Interpolations\src\extrapolation\extrapolation.jl:168 [inlined]
 [6] extrapolate_value
   @ ~\.julia\dev\Interpolations\src\extrapolation\extrapolation.jl:162 [inlined]
 [7] (::Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, SafeInt64, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{SafeInt64}}}, Gridded{Linear{Throw{OnGrid}}}, Line{Nothing}})(x::Int64)
   @ Interpolations ~\.julia\dev\Interpolations\src\extrapolation\extrapolation.jl:52
 [8] top-level scope
   @ REPL[7]:1

@JeffreySarnoff
Copy link
Member

We just need to add the missing method to Ratios
since SafeInt64 is not Signed, but SafeSigned

This structuring

SafeInt64 <: Integer
SafeInt64 <: SafeSigned
!(SafeInt64 <: Signed)

Was borrowed from an early version of RoundingIntegers.jl by @timholy .
Improvements + amendations inside Julia have allowed that to be reworked

export RSigned, RUnsigned, RInteger
abstract type RSigned   <: Signed end
abstract type RUnsigned <: Unsigned end

const RInteger = Union{RSigned,RUnsigned}

which recovers the desired, dispatch simplifying, subtyping relationships.

I should apply this approach to SaferIntegers -- unless there is something about the development of SaferIntegers that precluded doing so years ago. Honestly, I do recall considering it ... and do not recall precisely why it did not serve. There is a structural distinction in the implementation of RoundingIntegers.jl when compared to SaferIntegers.jl: less complexity enures to its overarching logical schema.

@mkitti
Copy link
Collaborator

mkitti commented Sep 3, 2021

I should apply this approach to SaferIntegers

I gave it a shot: JeffreySarnoff/SaferIntegers.jl#20

@mkitti
Copy link
Collaborator

mkitti commented Sep 11, 2021

You should be able to detect overflow now without modification to any of the packages.

julia> using SaferIntegers, Interpolations

julia> itpi = LinearInterpolation(SafeInt[1,10000],SafeInt[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{SafeInt64},), ::Vector{SafeInt64}, Gridded(Linear())), Line()) with element type Float64:
 SimpleRatio{SafeInt64}(99980001, 99980001)
 SimpleRatio{SafeInt64}(999800010000, 99980001)

julia> itpi(10001)
ERROR: OverflowError: 999800010000 * 99980001 overflowed for type Int64
Stacktrace:

@JeffreySarnoff
Copy link
Member

yes, that works for me using Julia 17ish

@stevengj
Copy link
Member

stevengj commented Sep 11, 2021

If the user supplies an integer type, we stick with that integer type for the computation to take advantage of any performance gains along with any potential issues such as overflow.

This doesn't make sense to me.

It seems like the interpolation inputs should essentially be promoted to the type of the output, which is a floating-point type for integer inputs.

Significant performance gains from keeping in Int64 form seem unlikely; floating-point arithmetic is fast (we no longer live in the 1970s–80s). And the danger of catastrophic overflow is real — correctness beats small performance gains.

@mkitti
Copy link
Collaborator

mkitti commented Sep 11, 2021

At the moment the interpolations in this package return a SimpleRatio from Ratios.jl rather than a float when given an integer. That is a design decision which preceeds me. Tim would be better able to comment on this than I.

Integer based output from integer input seems reasonable to me, although I can see how some may not expect it or specifically SimpleRatio. The user could select float output by giving float input, just as they will get unitful output from unitful input.

I agree there may be distinct expectations of this package. For this reason, I am considering using Preferences.jl to make Interpolations.ratio configurable at compile time.

@timholy
Copy link
Member

timholy commented Sep 12, 2021

Not just integers, either. The idea of using a ratio comes from

julia> 3//4 * 1.0
0.75

julia> 3//4 * 1.0f0
0.75f0

This allows the correct type to be returned from itp(x, y). It's much harder to do that if you pick a particular floating-point type.

That said, we should still do the den trick if the performance implications are not too bad. I forget whether this has been carefully benchmarked (I have not done it seriously myself). Detecting overflow is all fine and good, but avoiding it is even better.

@mkitti
Copy link
Collaborator

mkitti commented Sep 12, 2021

I was thinking that perhaps we should create a new ratio type for the common denominator detection.

@timholy
Copy link
Member

timholy commented Sep 15, 2021

It looks like Interpolations is the only consumer:

julia> using RegistryCompatTools

help?> held_back_by
search: held_back_by held_back_packages print_held_back

  held_back_by(pkgname::AbstractString, version::VersionNumber)

  Return a list of packages would hold back the given not-yet-registered version of pkgname.

  ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

  held_back_by(pkgname::AbstractString, d=held_back_packages())
  
  Return a list of packages that are holding back `pkgname`. `d` can be obtained from [`held_back_packages`](@ref).

julia> held_back_by("Ratios", v"1.0")
1-element Vector{String}:
 "Interpolations"

And it looks like it has no performance impact (the following must be run while using Revise):

julia> using Interpolations, BenchmarkTools
[ Info: Precompiling Interpolations [a98d9a8b-a2ab-59e6-89dd-64a1c18fca59]

julia> itp1 = interpolate(rand(15), BSpline(Quadratic(Flat(OnGrid()))));

julia> itp2 = interpolate(rand(15, 12), BSpline(Quadratic(Periodic(OnGrid()))));

julia> itp2l = interpolate(rand(15, 12), BSpline(Linear()));

julia> using Ratios

julia> SimpleRatio(2, 5) + SimpleRatio(1, 5)
SimpleRatio{Int64}(15, 25)

julia> @benchmark itp1(x) setup=(x=3+3*rand())
BenchmarkTools.Trial: 10000 samples with 995 evaluations.
 Range (min  max):  27.977 ns  705.963 ns  ┊ GC (min  max): 0.00%  95.72%
 Time  (median):     28.592 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   30.399 ns ±  17.046 ns  ┊ GC (mean ± σ):  1.45% ±  2.52%

   ▄█▆▂                             ▂▂ ▁▃▁                     ▁
  █████▇▆▅▄▅▃▃▄▄▅▄▅▅▆▅▂▃▃▄▃▃▃▃▂▅██▆████████▆▅▆▆▆▄▆▆▆▅▅▅▅▅▄▄▅▄▄ █
  28 ns         Histogram: log(frequency) by time      41.7 ns <

 Memory estimate: 32 bytes, allocs estimate: 2.

julia> @benchmark itp2(x, y) setup=(x=3+3*rand(); y=4+2*rand())
BenchmarkTools.Trial: 10000 samples with 979 evaluations.
 Range (min  max):  64.200 ns  932.400 ns  ┊ GC (min  max): 0.00%  91.49%
 Time  (median):     65.199 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   67.493 ns ±  26.095 ns  ┊ GC (mean ± σ):  1.20% ±  2.91%

  ▂▃▇█▆▃▁                                  ▁▂▂▁▁▂▁             ▂
  ███████▇▇▅▆███▆▇▇█▇▇▅▅▆▃▄▄▁▃▄▁▁▁▄▃▁▁▁▁▆▇█████████▆▇▇▇▆▇▆▅▇▇▇ █
  64.2 ns       Histogram: log(frequency) by time      81.9 ns <

 Memory estimate: 48 bytes, allocs estimate: 3.

julia> @benchmark itp2l(x, y) setup=(x=3+3*rand(); y=4+2*rand())
BenchmarkTools.Trial: 10000 samples with 994 evaluations.
 Range (min  max):  33.037 ns  896.178 ns  ┊ GC (min  max): 0.00%  94.73%
 Time  (median):     34.229 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   36.485 ns ±  25.571 ns  ┊ GC (mean ± σ):  2.18% ±  3.02%

  ▂▃▂█▇▅▃▆▄                                ▁▁▁▂▂▁              ▂
  █████████▇▆▆▆▇▇▅▆▅▆▆▆▆▄▆▄▃▁▄▃▁▁▁▁▄▃▁▃▅▇▇█████████▇▇▇▆▇▇▇▇▆▇▆ █
  33 ns         Histogram: log(frequency) by time      50.4 ns <

 Memory estimate: 48 bytes, allocs estimate: 3.

shell> git stash pop
On branch master
Your branch is up to date with 'origin/master'.

Changes not staged for commit:
  (use "git add <file>..." to update what will be committed)
  (use "git restore <file>..." to discard changes in working directory)
	modified:   src/Ratios.jl

no changes added to commit (use "git add" and/or "git commit -a")
Dropped refs/stash@{0} (da13fd30491a55a8f70425dbfb94122e543f35e7)

julia> SimpleRatio(2, 5) + SimpleRatio(1, 5)
SimpleRatio{Int64}(3, 5)

julia> @benchmark itp1(x) setup=(x=3+3*rand())
BenchmarkTools.Trial: 10000 samples with 995 evaluations.
 Range (min  max):  27.946 ns   1.055 μs  ┊ GC (min  max): 0.00%  96.76%
 Time  (median):     28.552 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   30.614 ns ± 26.329 ns  ┊ GC (mean ± σ):  2.25% ±  2.56%

   ▅█▆▂                            ▁▂▂▁▂▂▁                    ▁
  █████▆▅▅▃▅▅▄▄▄▄▄▄▃▄▅▆▆▆▅▃▄▄▄▂▇██▆███████▇▆▅▅▄▄▄▄▅▅▅▅▆▆▅▆▆▅▅ █
  27.9 ns      Histogram: log(frequency) by time      41.6 ns <

 Memory estimate: 32 bytes, allocs estimate: 2.

julia> @benchmark itp2(x, y) setup=(x=3+3*rand(); y=4+2*rand())
BenchmarkTools.Trial: 10000 samples with 979 evaluations.
 Range (min  max):  63.468 ns   1.133 μs  ┊ GC (min  max): 0.00%  92.65%
 Time  (median):     64.884 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   67.562 ns ± 32.278 ns  ┊ GC (mean ± σ):  1.48% ±  2.96%

  ▁▆█▇▃                  ▁▂▃▂▂                                ▂
  ██████▇▇▇██▇▇▇▇▆▆▅▅▄▄▄▇██████▇▇▇▇██▇▅▅▄▁▃▅▄▄▁▄▁▃▁▁▄▃▁▁▁▄▄▁▃ █
  63.5 ns      Histogram: log(frequency) by time      94.2 ns <

 Memory estimate: 48 bytes, allocs estimate: 3.

julia> @benchmark itp2l(x, y) setup=(x=3+3*rand(); y=4+2*rand())
BenchmarkTools.Trial: 10000 samples with 994 evaluations.
 Range (min  max):  33.057 ns   1.052 μs  ┊ GC (min  max): 0.00%  96.43%
 Time  (median):     34.306 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   36.839 ns ± 30.796 ns  ┊ GC (mean ± σ):  2.61% ±  3.04%

  ▂▃▂█▆▃▄▅▁                             ▁ ▁▂▁                 ▁
  █████████▆▅▆▆▆▅▅▆▇▅▅▆▆▅▅▅▃▄▅▄▃▄▃▃▃▅▆▇████████▆▅▆▄▅▆▆▆▆▇▆▆▄▅ █
  33.1 ns      Histogram: log(frequency) by time      51.9 ns <

 Memory estimate: 48 bytes, allocs estimate: 3.

shell> git diff
diff --git a/src/Ratios.jl b/src/Ratios.jl
index c7a4a98..75ffb4f 100644
--- a/src/Ratios.jl
+++ b/src/Ratios.jl
@@ -32,8 +32,10 @@ Rational{T}(r::SimpleRatio{S}) where {T<:Integer, S<:Integer} = convert(T, r.num
 /(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den, y.num)
 +(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den + y.num, y.den)
 -(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den - y.num, y.den)
-+(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.den + x.den*y.num, x.den*y.den)
--(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.den - x.den*y.num, x.den*y.den)
++(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num + y.num, x.den) :
+                                                     SimpleRatio(x.num*y.den + x.den*y.num, x.den*y.den)
+-(x::SimpleRatio, y::SimpleRatio) = x.den == y.den ? SimpleRatio(x.num - y.num, x.den) :
+                                                     SimpleRatio(x.num*y.den - x.den*y.num, x.den*y.den)
 ^(x::SimpleRatio, y::Integer) = SimpleRatio(x.num^y, x.den^y)
 
 -(x::SimpleRatio{T}) where {T<:Signed} = SimpleRatio(-x.num, x.den)

timholy added a commit to timholy/Ratios.jl that referenced this issue Sep 15, 2021
When implementing `r1 + r2`, check whether they have the same
denominator, and if so skip the generic multiplication-of-denominators.
This is not as comprehensive as using `gcd`, but it is *much* faster,
and that is a key design goal for this package.

For consumers, one can exploit this by ensuring that all coefficients
for a single axis have the same denominator.  For example, if you
wanted to construct a weighted average of `a`, `b`, and `c`, rather
than

    cedge, cmid = SimpleRatio(1, 4), SimpleRatio(1, 2)
    cedge*a + cmid*b + cedge*c

it would be better to use

    cedge, cmid = SimpleRatio(1, 4), SimpleRatio(2, 4)

Fixes JuliaMath/Interpolations.jl#457
@timholy
Copy link
Member

timholy commented Sep 15, 2021

timholy/Ratios.jl#24

timholy added a commit to timholy/Ratios.jl that referenced this issue Sep 16, 2021
When implementing `r1 + r2`, check whether they have the same
denominator, and if so skip the generic multiplication-of-denominators.
This is not as comprehensive as using `gcd`, but it is *much* faster,
and that is a key design goal for this package.

For consumers, one can exploit this by ensuring that all coefficients
for a single axis have the same denominator.  For example, if you
wanted to construct a weighted average of `a`, `b`, and `c`, rather
than

    cedge, cmid = SimpleRatio(1, 4), SimpleRatio(1, 2)
    cedge*a + cmid*b + cedge*c

it would be better to use

    cedge, cmid = SimpleRatio(1, 4), SimpleRatio(2, 4)

Fixes JuliaMath/Interpolations.jl#457
@mkitti
Copy link
Collaborator

mkitti commented Sep 16, 2021

Just to summarize, with Ratios.jl 0.4.2 we now get the correct answer without the need for overflow detection:

julia> using Interpolations

julia> itpi = LinearInterpolation(Int[1,10000],Int[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Int64}, Gridded(Linear())), Line()) with element type Float64:
 Ratios.SimpleRatio{Int64}(9999, 9999)
 Ratios.SimpleRatio{Int64}(99990000, 9999)

julia> convert(Float64, itpi(10001))
10001.0

julia> itpi(10001)
Ratios.SimpleRatio{Int64}(99999999, 9999)

With SaferIntegers 3.x, you can also still watch for overflow:

julia> using SaferIntegers, Interpolations

julia> itpi = LinearInterpolation(SafeInt[1,10000],SafeInt[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{SafeInt64},), ::Vector{SafeInt64}, Gridded(Linear())), Line()) with element type Float64:
 Ratios.SimpleRatio{SafeInt64}(9999, 9999)
 Ratios.SimpleRatio{SafeInt64}(99990000, 9999)

julia> itpi(10001)
Ratios.SimpleRatio{SafeInt64}(99999999, 9999)

The benchmark indicates the performance impact is either negligible or tolerable.

Also, if you prefer floating point, you can always still use floating point input:

julia> itpi = LinearInterpolation(Float64[1,10000], Float64[1,10000]; extrapolation_bc=Line())
2-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Line()) with element type Float64:
     1.0
 10000.0

julia> itpi(10001)
10001.0

I'm still interested in providing a mechanism to configure Interpolations.ratio so that the end user has some choice of the backend implementation. Alternatives to SimpleRatio from Ratios.jl include Float64 from Base, Rationals from Base, and FastRationals.jl.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

5 participants