Skip to content
This repository has been archived by the owner on Feb 9, 2020. It is now read-only.

Commit

Permalink
Merge pull request #35 from mbauman/siunits
Browse files Browse the repository at this point in the history
Allow interpolation with units!
  • Loading branch information
timholy committed Jul 16, 2014
2 parents ad2f1d7 + 06c4ff9 commit a099b87
Showing 1 changed file with 12 additions and 11 deletions.
23 changes: 12 additions & 11 deletions src/interp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -263,16 +263,16 @@ function valgradhess{T}(g::Vector{T}, h::Matrix{T}, G::InterpGrid{T}, x::Real...
end

## Vectorized evaluation at multiple points
function getindex{T,R<:Real}(G::AbstractInterpGrid{T,1}, x::AbstractVector{R})
function getindex{T,R<:Number}(G::AbstractInterpGrid{T,1}, x::AbstractVector{R})
n = length(x)
v = Array(T, n)
for i = 1:n
v[i] = getindex(G,x[i])
end
v
end
getindex{T,N,R<:Real}(G::AbstractInterpGrid{T,N}, x::AbstractVector{R}) = error("Linear indexing not supported")
function getindex{T,R<:Real}(G::AbstractInterpGrid{T,2}, x::AbstractVector{R}, y::AbstractVector{R})
getindex{T,N,R<:Number}(G::AbstractInterpGrid{T,N}, x::AbstractVector{R}) = error("Linear indexing not supported")
function getindex{T,R<:Number}(G::AbstractInterpGrid{T,2}, x::AbstractVector{R}, y::AbstractVector{R})
nx, ny = length(x), length(y)
v = Array(T, nx, ny)
for i = 1:nx
Expand All @@ -282,7 +282,7 @@ function getindex{T,R<:Real}(G::AbstractInterpGrid{T,2}, x::AbstractVector{R}, y
end
v
end
function getindex{T,N,R<:Real}(G::AbstractInterpGrid{T,N}, x::AbstractVector{R}, xrest::AbstractVector{R}...)
function getindex{T,N,R<:Number}(G::AbstractInterpGrid{T,N}, x::AbstractVector{R}, xrest::AbstractVector{R}...)
if length(xrest) != N-1
error("Dimensionality mismatch")
end
Expand All @@ -308,17 +308,17 @@ end
# Currently supports only 1d, nearest-neighbor or linear
# Consequently, the internal representation may change in the future
# BCperiodic and BCreflect not supported
type InterpIrregular{T<:FloatingPoint, S, N, BC<:BoundaryCondition, IT<:InterpType} <: AbstractInterpGrid{S,N,BC,IT}
type InterpIrregular{T<:Number, S, N, BC<:BoundaryCondition, IT<:InterpType} <: AbstractInterpGrid{S,N,BC,IT}
grid::Vector{Vector{T}}
coefs::Array{S,N}
x::Vector{T}
fillval::S # used only for BCfill (if ever)
end
InterpIrregular{T<:FloatingPoint, BC<:BoundaryCondition, IT<:Union(InterpNearest,InterpLinear)}(grid::Vector{T}, A::AbstractVector, ::Type{BC}, ::Type{IT}) =
InterpIrregular{T<:Number, BC<:BoundaryCondition, IT<:Union(InterpNearest,InterpLinear)}(grid::Vector{T}, A::AbstractArray, ::Type{BC}, ::Type{IT}) =
InterpIrregular(Vector{T}[grid], A, BC, IT) # special 1d syntax
InterpIrregular{T<:FloatingPoint, BC<:BoundaryCondition, IT<:Union(InterpNearest,InterpLinear)}(grid::(Vector{T}...), A::AbstractVector, ::Type{BC}, ::Type{IT}) =
InterpIrregular{T<:Number, BC<:BoundaryCondition, IT<:Union(InterpNearest,InterpLinear)}(grid::(Vector{T}...), A::AbstractArray, ::Type{BC}, ::Type{IT}) =
InterpIrregular(Vector{T}[grid...], A, BC, IT)
function InterpIrregular{T<:FloatingPoint, S, N, BC<:BoundaryCondition, IT<:Union(InterpNearest,InterpLinear)}(grid::Vector{Vector{T}}, A::AbstractArray{S, N}, ::Type{BC}, ::Type{IT})
function InterpIrregular{T<:Number, S, N, BC<:BoundaryCondition, IT<:Union(InterpNearest,InterpLinear)}(grid::Vector{Vector{T}}, A::AbstractArray{S, N}, ::Type{BC}, ::Type{IT})
if length(grid) != 1
error("Sorry, for now only 1d is supported")
end
Expand All @@ -342,23 +342,24 @@ function InterpIrregular{IT<:InterpType}(grid, A::Array, f::Number, ::Type{IT})
iu
end

function _getindexii{T,S,BC<:Union(BCfill,BCna,BCnan)}(G::InterpIrregular{T,S,1,BC}, x::Real)
function _getindexii{T,S,BC<:Union(BCfill,BCna,BCnan)}(G::InterpIrregular{T,S,1,BC}, x::Number)
g = G.grid[1]
i = (x == g[1]) ? 2 : searchsortedfirst(g, x)
(i == 1 || i == length(g)+1) ? G.fillval : _interpu(x, g, i, G.coefs, interptype(G))
end
function _getindexii{T,S}(G::InterpIrregular{T,S,1,BCnil}, x::Real)
function _getindexii{T,S}(G::InterpIrregular{T,S,1,BCnil}, x::Number)
g = G.grid[1]
i = (x == g[1]) ? 2 : searchsortedfirst(g, x)
(i == 1 || i == length(g)+1) ? error(BoundsError) : _interpu(x, g, i, G.coefs, interptype(G))
end
function _getindexii{T,S}(G::InterpIrregular{T,S,1,BCnearest}, x::Real)
function _getindexii{T,S}(G::InterpIrregular{T,S,1,BCnearest}, x::Number)
g = G.grid[1]
i = (x == g[1]) ? 2 : searchsortedfirst(g, x)
i == 1 ? G.coefs[1] : i == length(g)+1 ? G.coefs[end] : _interpu(x, g, i, G.coefs, interptype(G))
end
# This next is necessary for precedence
getindex(G::InterpIrregular, x::Real) = _getindexii(G, x)
getindex(G::InterpIrregular, x::Number) = _getindexii(G, x)

_interpu(x, g, i, coefs, ::Type{InterpNearest}) = (x-g[i-1] < g[i]-x) ? coefs[i-1] : coefs[i]
function _interpu(x, g, i, coefs, ::Type{InterpLinear})
Expand Down

0 comments on commit a099b87

Please # to comment.