Support FixedPoint

[timholy]

Fixes #18

[codecov]

Codecov Report

Merging #19 (8b6d6fa) into master (3476a71) will increase coverage by 5.79%.
The diff coverage is 100.00%.

❗ Current head 8b6d6fa differs from pull request most recent head f29ff51. Consider uploading reports for the commit f29ff51 to get more accurate results

@@            Coverage Diff             @@
##           master      #19      +/-   ##
==========================================
+ Coverage   83.33%   89.13%   +5.79%     
==========================================
  Files           1        1              
  Lines          36       46      +10     
==========================================
+ Hits           30       41      +11     
+ Misses          6        5       -1     

Impacted Files	Coverage Δ
src/Ratios.jl	89.13% <100.00%> (+5.79%)	⬆️

---

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 3476a71...f29ff51. Read the comment docs.

[mkitti]

Conversion to floating point is practical, but it seems to defeat the purpose of using either SimpleRatio or Normed to begin with. If I understand correctly, I believe any fixed point number can be expressed exactly as a SimpleRatio.

0.035N0f8 is exactly equivalent to SimpleRatio{UInt8}(0x09, 0xff). That is we could write the following.

SimpleRatio(x::Normed{T,f}) where {T,f} = SimpleRatio{T}(x.i, 2^f-1)

From a mathematical perspective, promotion of FixedPoint to SimpleRatio seems better, I think.

I have not tested how this works for the original problem yet.

[timholy]

If I understand correctly, I believe any fixed point number can be expressed exactly as a SimpleRatio.

That is true, but:

julia> convert(N0f8, SimpleRatio(8, 3) * 0.8N0f8)
ERROR: ArgumentError: N0f8 is an 8-bit type representing 256 values from 0.0 to 1.0; cannot represent 2.1333334
Stacktrace:
 [1] throw_converterror(#unused#::Type{N0f8}, x::Float32)
   @ FixedPointNumbers ~/.julia/packages/FixedPointNumbers/HAGk2/src/FixedPointNumbers.jl:326
 [2] _convert
   @ ~/.julia/packages/FixedPointNumbers/HAGk2/src/normed.jl:76 [inlined]
 [3] FixedPoint
   @ ~/.julia/packages/FixedPointNumbers/HAGk2/src/FixedPointNumbers.jl:58 [inlined]
 [4] convert(#unused#::Type{N0f8}, x::Float32)
   @ Base ./number.jl:7
 [5] top-level scope
   @ REPL[4]:1

For Interpolations, I guess all the cases you care about are bounded by ±1? If so then this exact issue won't come up in practice. There's still the issue that, e.g., 7 is not a factor of the denominator for any of the FixedPoint types.

julia> Float32(convert(N0f8, SimpleRatio(6, 7) * 0.8N0f8)) == Float32((6*0.8N0f8)/7)
false

[timholy]

Oh, I see, you were interested in the other direction. Yes, that is actually kind of intriguing...

[timholy]

See if you like this better. One danger, though, is that

Ratios.jl/src/Ratios.jl

Lines 31 to 32 in 3476a71

	+(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)

will quickly lead to overflows.

[mkitti]

The overflow potential does seem like it could be problematic if misapplied but may be acceptable for something called a SimpleRatio. If someone wanted an OverflowError perhaps then Rational would be a better choice.

For the application in question, I would guess that the FixedPoint number being used might be of shallower bit-depth than Int64.

I wonder if there are some reasonable circumstances where a warning could be provided.

[timholy]

Perhaps the best defense would be to check whether x.den == y.den. That would slow arithmetic a bit but not catastrophically so. If Interpolations is careful to use a common denominator for all coefficients then this will go a long, long ways to preventing overflow in any realistic situation.

[timholy]

Let me know if you want other changes or are happy with this. We can address the arithmetic in a separate PR.

[mkitti]

Returning to @johnnychen94 's original issue in JuliaMath/Interpolations.jl#413, the following still does not work:

using FixedPointNumbers, Interpolations
interpolate(rand(N0f8, 4, 4),  BSpline(Quadratic(Flat(OnGrid()))))
ERROR: ArgumentError: N0f8 is an 8-bit type representing 256 values from 0.0 to 1.0; cannot represent 1.1411269974768714
Stacktrace:
  [1] throw_converterror(#unused#::Type{N0f8}, x::Float64)
    @ FixedPointNumbers ~/.julia/packages/FixedPointNumbers/HAGk2/src/FixedPointNumbers.jl:326
  [2] _convert
    @ ~/.julia/packages/FixedPointNumbers/HAGk2/src/normed.jl:76 [inlined]
  [3] FixedPoint
    @ ~/.julia/packages/FixedPointNumbers/HAGk2/src/FixedPointNumbers.jl:58 [inlined]
  [4] convert
    @ ./number.jl:7 [inlined]
  [5] setindex!
    @ ./array.jl:841 [inlined]
  [6] setindex!
    @ ./multidimensional.jl:639 [inlined]
  [7] _A_ldiv_B_md!(dest::Matrix{N0f8}, F::LinearAlgebra.LU{Float64, LinearAlgebra.Tridiagonal{Float64, Vector{Float64}}}, src::Matrix{N0f8}, R1::CartesianIndices{0, Tuple{}}, R2::CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, b::Vector{N0f8})
    @ Interpolations ~/.julia/dev/Interpolations/src/filter1d.jl:39
  [8] _A_ldiv_B_md!(dest::Matrix{N0f8}, W::WoodburyMatrices.Woodbury{Float64, LinearAlgebra.LU{Float64, LinearAlgebra.Tridiagonal{Float64, Vector{Float64}}}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Matrix{Float64}, Matrix{Float64}}, src::Matrix{N0f8}, R1::CartesianIndices{0, Tuple{}}, R2::CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, b::Vector{N0f8})
...

That's because quadratic interpolation is not guaranteed to be constrained between the values being interpolated. If I choose a carefully chosen set of values the it seems to work fine now as opposed to the master branch.

julia> r = reinterpret(N0f8,[0x73, 0x94, 0x66, 0x1e])
4-element reinterpret(N0f8, ::Vector{UInt8}):
 0.451N0f8
 0.58N0f8
 0.4N0f8
 0.118N0f8

julia> itp = interpolate(r,  BSpline(Quadratic(Flat(OnGrid()))))
4-element interpolate(OffsetArray(::Array{N0f8,1}, 0:5), BSpline(Quadratic(Flat(OnGrid())))) with element type Float64:
 Ratios.SimpleRatio{Int64}(1914336000, 4244832000)
 Ratios.SimpleRatio{Int64}(2465748000, 4244832000)
 Ratios.SimpleRatio{Int64}(1695852000, 4244832000)
 Ratios.SimpleRatio{Int64}(499392000, 4244832000)

julia> x = 1:0.1:4; y = itp(x);

julia> Plots.plot(x,y)

julia> Plots.scatter!(1:4, r)

To address the original question, we need to clamp the "image":

julia> r = rand(N0f8, 4, 4)
4×4 reinterpret(N0f8, ::Matrix{UInt8}):
 0.275  1.0    0.467  0.298
 0.953  0.635  0.773  0.706
 0.38   0.471  0.039  0.925
 0.243  0.278  0.227  0.694

julia> r = r * 0.3N0f16 .+ 0.35N0f16
4×4 Array{N0f16,2} with eltype N0f16:
 0.43235  0.64999  0.49     0.4394
 0.63587  0.54058  0.58175  0.56176
 0.46412  0.49117  0.36176  0.62763
 0.42293  0.43352  0.41823  0.55822

julia> itp = interpolate(r,  BSpline(Quadratic(Flat(OnGrid()))))
4×4 interpolate(OffsetArray(::Array{N0f16,2}, 0:5, 0:5), BSpline(Quadratic(Flat(OnGrid())))) with element type Float64:
 SimpleRatio{Int64}(5372212362956767232, 2533270495428608)   …  SimpleRatio{Int64}(-4891770677985542144, 2533270495428608)
 SimpleRatio{Int64}(5312948799365185536, 2533270495428608)      SimpleRatio{Int64}(7026829071363342336, 2533270495428608)
 SimpleRatio{Int64}(683411144106311680, 2533270495428608)       SimpleRatio{Int64}(4225819746047623168, 2533270495428608)
 SimpleRatio{Int64}(-7070532375658102784, 2533270495428608)     SimpleRatio{Int64}(633194344873984000, 2533270495428608)

julia> ritp = [itp(x,y) for x in 1:0.1:4, y in 1:0.1:4];

julia> Plots.plot(Plots.heatmap.((r,ritp))..., layout=(1,2))

In this regard, I'm pretty happy with this.

[mkitti]

How about inferring the parameter S from the Fixed type?

SimpleRatio(x::Normed{S}) where S = SimpleRatio{S}( x )
SimpleRatio(x::Fixed{S, f}) where {S,f} = x == -1 ? SimpleRatio{S}( 1, -1 ) : SimpleRatio{S}( -reinterpret(x), -1 << f )

[timholy]

We could add that. For the particular use here the promotion specifies the specific type and so of course you need to respect that, but for anyone who just calls SimpleRatio(x::Fixed) this would make sense.

[timholy]

That's because quadratic interpolation is not guaranteed to be constrained between the values being interpolated.

Interpolations either needs to pick a safe element type (e.g., Float32 or widen it while allowing negative values) or allow the user to control the storage type. Since there is no signed direct analog of N0f8 (we don't have a S7f8 which I think would be the ideal type, see JuliaMath/FixedPointNumbers.jl#143), really I think Float32 is your best option. I want to evaluate making N0f8 storage-only for the purpose of image processing, anyway, because the overflow issues with N0f8 are pretty troubling. The only issue is the 4x increase in memory, which would not be ideal.