JuliaMath / Interpolations.jl

Fast, continuous interpolation of discrete datasets in Julia
http://juliamath.github.io/Interpolations.jl/
Other
533 stars 109 forks source link

Type promotion #359

Open haampie opened 4 years ago

haampie commented 4 years ago

It seems interpolating static arrays of eltype Float32 causes a type promotion to Float64:

julia> interpolate([@SVector(rand(Float32, 3)) for i=1:5], BSpline(Linear()))
5-element interpolate(::Array{SArray{Tuple{3},Float32,1,3},1}, BSpline(Linear())) with element type SArray{Tuple{3},Float64,1,3}:
 Float32[0.40708745, 0.08114183, 0.64380157] 
 Float32[0.889678, 0.40728843, 0.3275386]    
 Float32[0.064697385, 0.27480924, 0.08783579]
 Float32[0.21843159, 0.546651, 0.7467506]    
 Float32[0.6662456, 0.8931141, 0.17646503]

It says: with element type SVector{3,Float64}

Is that expected?

It seems to be hitting tweight(A::AbstractArray) = Float64.

haampie commented 4 years ago

Similar issue with scale and without static arrays:

julia> itp = scale(interpolate([1f0, 2f0, 3f0], BSpline(Quadratic(Natural(OnCell())))), range(0.0f0, 10.0f0, length=3));

julia> itp(5.0f0)
1.9999999701976776

julia> typeof(ans)
Float64
haampie commented 4 years ago

There's multiple issues. Regarding scale:

diff --git a/src/scaling/scaling.jl b/src/scaling/scaling.jl
index b700824..5457814 100644
--- a/src/scaling/scaling.jl
+++ b/src/scaling/scaling.jl
@@ -104,7 +104,7 @@ coordlookup(r::AbstractUnitRange, x) = x - first(r) + oneunit(eltype(r))
 coordlookup(r::StepRange, x) = (x - r.start) / r.step + oneunit(eltype(r))

 coordlookup(r::AbstractRange, x) = (x - first(r)) / step(r) + oneunit(eltype(r))
-boundstep(r::AbstractRange) = 0.5*step(r)
+boundstep(r::AbstractRange) = step(r) / 2
 rescale_gradient(r::AbstractRange, g) = g / step(r)

 basetype(::Type{ScaledInterpolation{T,N,ITPT,IT,RT}}) where {T,N,ITPT,IT,RT} = ITPT

and

diff --git a/src/b-splines/b-splines.jl b/src/b-splines/b-splines.jl
index e066cce..4b88b54 100644
--- a/src/b-splines/b-splines.jl
+++ b/src/b-splines/b-splines.jl
@@ -116,8 +116,8 @@ ubound(ax::AbstractRange, bs::BSpline)   = ubound(ax, degree(bs))
 ubound(ax::AbstractRange, deg::Degree)   = last(ax)
 ubound(ax::AbstractRange, deg::DegreeBC) = ubound(ax, deg, deg.bc.gt)

-lbound(ax::AbstractUnitRange, ::DegreeBC, ::OnCell) = first(ax) - 0.5
-ubound(ax::AbstractUnitRange, ::DegreeBC, ::OnCell) = last(ax) + 0.5
+lbound(ax::AbstractUnitRange, ::DegreeBC, ::OnCell) = first(ax) - 0.5f0
+ubound(ax::AbstractUnitRange, ::DegreeBC, ::OnCell) = last(ax) + 0.5f0
 lbound(ax::AbstractUnitRange, ::DegreeBC, ::OnGrid) = first(ax)
 ubound(ax::AbstractUnitRange, ::DegreeBC, ::OnGrid) = last(ax)

The first is of course fine for a PR, but in the second one the actual element type is not propagated to the lbound and ubound functions, so I figure the proper solution would be to do that and use T(0.5).

timholy commented 4 years ago

Yeah, this could use a cleanup.

maxfreu commented 4 years ago

I also bumped into this issue while using Augmentor.jl. Unfortunately I dont see how I could fix it.