diff --git a/Project.toml b/Project.toml index 482e4aa..3f87d6c 100644 --- a/Project.toml +++ b/Project.toml @@ -8,8 +8,8 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] -julia = "1" StaticArrays = "0.9, 0.10, 0.11, 0.12, 1" +julia = "1" [extras] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" diff --git a/src/GridInterpolations.jl b/src/GridInterpolations.jl index a40f1c6..bf516aa 100644 --- a/src/GridInterpolations.jl +++ b/src/GridInterpolations.jl @@ -13,7 +13,7 @@ mutable struct RectangleGrid{D} <: AbstractGrid{D} cut_counts::Vector{Int} cuts::Vector{Float64} - function RectangleGrid{D}(cutPoints...) where D + function RectangleGrid{D}(cutPoints...) where {D} cut_counts = Int[length(cutPoints[i]) for i = 1:length(cutPoints)] cuts = vcat(cutPoints...) myCutPoints = Array{Vector{Float64}}(undef, length(cutPoints)) @@ -21,7 +21,7 @@ mutable struct RectangleGrid{D} <: AbstractGrid{D} @assert numDims == D for i = 1:numDims if length(Set(cutPoints[i])) != length(cutPoints[i]) - error(@sprintf("Duplicates cutpoints are not allowed (duplicates observed in dimension %d)",i)) + error(@sprintf("Duplicates cutpoints are not allowed (duplicates observed in dimension %d)", i)) end if !issorted(cutPoints[i]) error("Cut points must be sorted") @@ -43,7 +43,7 @@ mutable struct SimplexGrid{D} <: AbstractGrid{D} ilo::Vector{Int} # indices of cuts below point n_ind::Vector{Int} - function SimplexGrid{D}(cutPoints...) where D + function SimplexGrid{D}(cutPoints...) where {D} cut_counts = Int[length(cutPoints[i]) for i = 1:length(cutPoints)] cuts = vcat(cutPoints...) myCutPoints = Array{Vector{Float64}}(undef, length(cutPoints)) @@ -51,7 +51,7 @@ mutable struct SimplexGrid{D} <: AbstractGrid{D} @assert numDims == D for i = 1:numDims if length(Set(cutPoints[i])) != length(cutPoints[i]) - error(@sprintf("Duplicates cutpoints are not allowed (duplicates observed in dimension %d)",i)) + error(@sprintf("Duplicates cutpoints are not allowed (duplicates observed in dimension %d)", i)) end if !issorted(cutPoints[i]) error("Cut points must be sorted") @@ -72,8 +72,8 @@ Base.length(grid::RectangleGrid) = prod(grid.cut_counts) Base.size(grid::RectangleGrid) = Tuple(grid.cut_counts) Base.length(grid::SimplexGrid) = prod(grid.cut_counts) -dimensions(grid::AbstractGrid{D}) where D = D -Base.ndims(grid::AbstractGrid{D}) where D = D +dimensions(grid::AbstractGrid{D}) where {D} = D +Base.ndims(grid::AbstractGrid{D}) where {D} = D label(grid::RectangleGrid) = "multilinear interpolation grid" label(grid::SimplexGrid) = "simplex interpolation grid" @@ -101,21 +101,21 @@ end function ind2x!(grid::AbstractGrid, ind::Int, x::AbstractArray) # Populates x with the value at ind. - # In-place version of ind2x. - # Example: - # rgrid = RectangleGrid([2,5],[20,50]) - # x = [0,0] - # ind2x!(rgrid,4,x) # x now contains [5,50] + # In-place version of ind2x. + # Example: + # rgrid = RectangleGrid([2,5],[20,50]) + # x = [0,0] + # ind2x!(rgrid,4,x) # x now contains [5,50] # @show x # displays [5,50] ndims = dimensions(grid) stride = grid.cut_counts[1] - for i=2:ndims-1 + for i = 2:ndims-1 stride *= grid.cut_counts[i] end - for i=(ndims-1):-1:1 - rest = rem(ind-1, stride) + 1 - x[i + 1] = grid.cutPoints[i + 1][div(ind - rest, stride) + 1] + for i = (ndims-1):-1:1 + rest = rem(ind - 1, stride) + 1 + x[i+1] = grid.cutPoints[i+1][div(ind - rest, stride)+1] ind = rest stride = div(stride, grid.cut_counts[i]) end @@ -127,8 +127,8 @@ end # masked interpolation ignores points that are masked function maskedInterpolate(grid::AbstractGrid, data::DenseArray, x::AbstractVector, mask::BitArray{1}) index, weight = interpolants(grid, x) - val = 0 - totalWeight = 0 + val = zero(eltype(data)) + totalWeight = zero(eltype(data)) for i = 1:length(index) if mask[index[i]] continue @@ -139,13 +139,13 @@ function maskedInterpolate(grid::AbstractGrid, data::DenseArray, x::AbstractVect return val / totalWeight end -interpolate(grid::AbstractGrid, data::Matrix, x::AbstractVector) = interpolate(grid, map(Float64, data[:]), x) +interpolate(grid::AbstractGrid, data::Matrix, x::AbstractVector) = interpolate(grid, map(eltype(data), data[:]), x) function interpolate(grid::AbstractGrid, data::DenseArray, x::AbstractVector) index, weight = interpolants(grid, x) - v = 0.0 - for (i,data_ind) in enumerate(index) - v += data[data_ind]*weight[i] + v = zero(eltype(data)) + for (i, data_ind) in enumerate(index) + v += data[data_ind] * weight[i] end return v end @@ -157,14 +157,19 @@ function interpolants(grid::RectangleGrid, x::AbstractVector) cut_counts = grid.cut_counts cuts = grid.cuts - # Reset the values in index and weight: - index = @MVector(ones(Int, 2^dimensions(grid))) - index2 = @MVector(ones(Int, 2^dimensions(grid))) - weight = @MVector(zeros(eltype(x), 2^dimensions(grid))) - weight2 = @MVector(zeros(eltype(x), 2^dimensions(grid))) - index[1] = 1 - index2[1] = 1 + num_points = 2^dimensions(grid) + index = MVector{num_points, Int}(undef) + index2 = MVector{num_points, Int}(undef) + weight = MVector{num_points, eltype(x)}(undef) + weight2 = MVector{num_points, eltype(x)}(undef) + + # Note: these values are set explicitly because we have not verified that the logic below is independent of the initial values. See discussion in PR #47. These can be removed if it can be proved that the logic is independent of the initial values. + index .= 1 + index2 .= 1 + weight .= zero(eltype(weight)) + weight2 .= zero(eltype(weight2)) + weight[1] = one(eltype(weight)) weight2[1] = one(eltype(weight2)) @@ -174,7 +179,7 @@ function interpolants(grid::RectangleGrid, x::AbstractVector) n = 1 for d = 1:length(x) coord = x[d] - lasti = cut_counts[d]+cut_i-1 + lasti = cut_counts[d] + cut_i - 1 ii = cut_i if coord <= cuts[ii] @@ -188,42 +193,42 @@ function interpolants(grid::RectangleGrid, x::AbstractVector) if cuts[ii] == coord i_lo, i_hi = ii, ii else - i_lo, i_hi = (ii-1), ii + i_lo, i_hi = (ii - 1), ii end end # the @inbounds are needed below to prevent allocation if i_lo == i_hi for i = 1:l - @inbounds index[i] += (i_lo - cut_i)*subblock_size + @inbounds index[i] += (i_lo - cut_i) * subblock_size end else - low = (1 - (coord - cuts[i_lo])/(cuts[i_hi]-cuts[i_lo])) + low = (1 - (coord - cuts[i_lo]) / (cuts[i_hi] - cuts[i_lo])) for i = 1:l - @inbounds index2[i ] = index[i] + (i_lo-cut_i)*subblock_size - @inbounds index2[i+l] = index[i] + (i_hi-cut_i)*subblock_size + @inbounds index2[i] = index[i] + (i_lo - cut_i) * subblock_size + @inbounds index2[i+l] = index[i] + (i_hi - cut_i) * subblock_size end @inbounds index[:] = index2 for i = 1:l - @inbounds weight2[i ] = weight[i]*low - @inbounds weight2[i+l] = weight[i]*(1-low) + @inbounds weight2[i] = weight[i] * low + @inbounds weight2[i+l] = weight[i] * (1 - low) end @inbounds weight[:] = weight2 - l = l*2 - n = n*2 + l = l * 2 + n = n * 2 end cut_i = cut_i + cut_counts[d] - subblock_size = subblock_size*(cut_counts[d]) + subblock_size = subblock_size * (cut_counts[d]) end - v = min(l,length(index)) - return view(SVector(index),1:v), view(SVector(weight),1:v) + v = min(l, length(index)) + return view(SVector(index), 1:v), view(SVector(weight), 1:v) end function interpolants(grid::SimplexGrid, x::AbstractVector) - weight = MVector{dimensions(grid)+1, eltype(x)}(undef) - index = MVector{dimensions(grid)+1, Int}(undef) + weight = MVector{dimensions(grid) + 1,eltype(x)}(undef) + index = MVector{dimensions(grid) + 1,Int}(undef) x_p = grid.x_p # residuals ihi = grid.ihi # indicies of cuts above point @@ -238,17 +243,17 @@ function interpolants(grid::SimplexGrid, x::AbstractVector) for i = 1:dimensions(grid) # find indicies of coords if match coord = x[i] - lasti = cut_counts[i]+cut_i-1 + lasti = cut_counts[i] + cut_i - 1 ii = cut_i # check bounds, snap to closest if out if coord <= cuts[ii] ihi[i] = ii ilo[i] = ii - x_p[i] = 0.0 + x_p[i] = zero(eltype(x)) elseif coord >= cuts[lasti] ihi[i] = lasti ilo[i] = lasti - x_p[i] = 0.0 + x_p[i] = zero(eltype(x)) else # increment through cut points if in bounds while cuts[ii] < coord @@ -258,10 +263,10 @@ function interpolants(grid::SimplexGrid, x::AbstractVector) if cuts[ii] == coord ilo[i] = ii ihi[i] = ii - x_p[i] = 0.0 + x_p[i] = zero(eltype(x)) else # if between cuts assign lo and high indecies and translate - ilo[i] = ii-1 + ilo[i] = ii - 1 ihi[i] = ii lo = cuts[ilo[i]] hi = cuts[ihi[i]] @@ -272,7 +277,9 @@ function interpolants(grid::SimplexGrid, x::AbstractVector) end # initialize sort indecies - for i = 1:length(n_ind); n_ind[i] = i; end + for i = 1:length(n_ind) + n_ind[i] = i + end # sort translated and scaled x values sortperm!(n_ind, x_p, rev=true) ############################################# killer of speed x_p = x_p[n_ind] @@ -282,7 +289,7 @@ function interpolants(grid::SimplexGrid, x::AbstractVector) for i = 1:(length(x_p)+1) if i == 1 weight[i] = 1 - x_p[i] - elseif i == length(x_p)+1 + elseif i == length(x_p) + 1 weight[i] = x_p[i-1] else weight[i] = x_p[i-1] - x_p[i] @@ -307,7 +314,7 @@ function interpolants(grid::SimplexGrid, x::AbstractVector) else index[i] += (ilo[k] - 1 - ct) * siz end - siz = siz*cut_counts[k] + siz = siz * cut_counts[k] ct += cut_counts[k] end index[i] += 1 @@ -323,14 +330,14 @@ function vertices(grid::AbstractGrid) n_dims = dimensions(grid) mem = Array{Float64,2}(undef, n_dims, length(grid)) - for idx = 1 : length(grid) - this_idx::Int = idx-1 + for idx = 1:length(grid) + this_idx::Int = idx - 1 # Get the correct index into each dimension # and populate vertex index with corresponding cut point - for j = 1 : n_dims + for j = 1:n_dims cut_idx::Int = this_idx % grid.cut_counts[j] - this_idx = div(this_idx,grid.cut_counts[j]) + this_idx = div(this_idx, grid.cut_counts[j]) mem[j, idx] = grid.cutPoints[j][cut_idx+1] end end @@ -341,7 +348,7 @@ function vertices(grid::AbstractGrid) (http://juliaarrays.github.io/StaticArrays.jl/stable/pages/ api.html#Arrays-of-static-arrays-1), and tests should catch these errors. =# - return reshape(reinterpret(SVector{n_dims, Float64}, mem), (length(grid),)) + return reshape(reinterpret(SVector{n_dims,Float64}, mem), (length(grid),)) end @@ -353,12 +360,12 @@ using Base.Sort # for sortperm! const DEFAULT_UNSTABLE = QuickSort function sortperm!(x::Vector{I}, v::AbstractVector; alg::Algorithm=DEFAULT_UNSTABLE, - lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward) where {I<:Integer} - sort!(x, alg, Perm(ord(lt,by,rev,order),v)) + lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward) where {I<:Integer} + sort!(x, alg, Perm(ord(lt, by, rev, order), v)) end function Base.iterate(iter::RectangleGrid, state::Int64=1) - return state<=length(iter) ? (ind2x(iter, state), state+1) : nothing + return state <= length(iter) ? (ind2x(iter, state), state + 1) : nothing end Base.getindex(grid::RectangleGrid, key::CartesianIndex) = ind2x(grid, LinearIndices(Dims((grid.cut_counts...,)))[key]) diff --git a/test/runtests.jl b/test/runtests.jl index 85a6cc0..5cfc075 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,36 +13,36 @@ function compareToGrid(testType::Symbol=:random, numDims::Int=3, pointsPerDim::I # Compare interpolation results from GridInterpolations with the Grid Module: # create a multi-dimensional test Grid with data, cut points, etc. # (dataM, gridM, gridI) = compareToGrid(2, 3, 5); - cutsDim = pointsPerDim*ones(Int,numDims) + cutsDim = pointsPerDim * ones(Int, numDims) gridData = rand(cutsDim...) # Set up the data structures: gridI = Interpolations.interpolate(gridData, BSpline(Linear())) # This is for the Interpolations package gridI = Interpolations.extrapolate(gridI, Flat()) cutPointMat = 0:(1/(pointsPerDim-1)):1 - cutPoints = [cutPointMat for i=1:numDims] + cutPoints = [cutPointMat for i = 1:numDims] gridM = RectangleGrid(tuple(cutPoints...)...) # Mykel's interpolation package # Package the data in a way that can be read by the interpolation package: dataM = Array{Float64}(undef, length(gridM)) - for i=1:length(gridM) + for i = 1:length(gridM) # x = ind2x(gridM, i) x = ind2x(gridM, i) - dataM[i] = gridData[getFractionalIndexes(gridM,x)...] + dataM[i] = gridData[getFractionalIndexes(gridM, x)...] end # Check a finer grid of coordinates, for comprehensiveness? - if testType==:random + if testType == :random # Interpolate in GridInterpolations with the same points, report any discrepancy # Select a number of random points and compare the interpolation results. testPoint = zeros(numDims) - for i=1:numRandomTests + for i = 1:numRandomTests rand!(testPoint) - testI = gridI(((pointsPerDim.-1)*testPoint.+1)...) - testM = GridInterpolations.interpolate(gridM,dataM,testPoint) + testI = gridI(((pointsPerDim .- 1) * testPoint .+ 1)...) + testM = GridInterpolations.interpolate(gridM, dataM, testPoint) - if (abs(testI-testM)>eps) + if (abs(testI - testM) > eps) display("Failed Random Point Interpolation Test") display("Non-matching interpolation points:") display("Test point = $testPoint") @@ -55,15 +55,15 @@ function compareToGrid(testType::Symbol=:random, numDims::Int=3, pointsPerDim::I return true end - if testType==:extrapNeg + if testType == :extrapNeg # Test behavior outside the grid, i.e. extrapolation: testPoint = zeros(numDims) - for i=1:numRandomTests + for i = 1:numRandomTests testPoint = -rand!(testPoint) - testI = gridI(((pointsPerDim-1)*testPoint+1)...) - testM = interpolate(gridM,dataM,testPoint) + testI = gridI(((pointsPerDim - 1) * testPoint + 1)...) + testM = interpolate(gridM, dataM, testPoint) - if (abs(testI-testM)>eps) + if (abs(testI - testM) > eps) display("Failed Random Negative Point Extrapolation Test") display("Non-matching interpolation points:") display("Test point = $testPoint") @@ -76,13 +76,13 @@ function compareToGrid(testType::Symbol=:random, numDims::Int=3, pointsPerDim::I return true end - if testType==:extrapPos - for i=1:numRandomTests - testPoint = rand(numDims).+1 - testI = gridI(((pointsPerDim-1)*testPoint.+1)...) - testM = interpolate(gridM,dataM,testPoint) + if testType == :extrapPos + for i = 1:numRandomTests + testPoint = rand(numDims) .+ 1 + testI = gridI(((pointsPerDim - 1) * testPoint .+ 1)...) + testM = interpolate(gridM, dataM, testPoint) - if (abs(testI-testM)>eps) + if (abs(testI - testM) > eps) display("Failed Random Positive Point Extrapolation Test") display("Non-matching interpolation points:") display("Test point = $testPoint") @@ -104,7 +104,7 @@ function getFractionalIndexes(g::AbstractGrid, s::Array) fracInd = Array{Int64}(undef, length(g.cutPoints)) - for i=1:length(g.cutPoints) + for i = 1:length(g.cutPoints) gridDisc = g.cutPoints[i] fracInd[i] = getFracIndex(gridDisc, s[i]) end @@ -126,12 +126,12 @@ function getFracIndex(vararray::Array, value::Float64) return argmax(vararray) end - i=1 + i = 1 while value > vararray[i+1] - i+=1 + i += 1 end - return i+(value-vararray[i])÷(vararray[i+1] - vararray[i]) + return i + (value - vararray[i]) ÷ (vararray[i+1] - vararray[i]) end @@ -140,9 +140,9 @@ function simplexMagic(NDISC::Int=20, NPOINTS::Int=3, checkFileName::AbstractStri checkFileName = joinpath(@__DIR__, "simplexMagicTest20.txt") end - val = transpose([ [8.,1,6] [3,5,7] [4,9,2] ]) # transposed magic(3) from matlab + val = transpose([[8.0, 1, 6] [3, 5, 7] [4, 9, 2]]) # transposed magic(3) from matlab - sGrid = SimplexGrid(1.:NPOINTS, 1.:NPOINTS) + sGrid = SimplexGrid(1.0:NPOINTS, 1.0:NPOINTS) sInterpVal = zeros(NDISC, NDISC) for i = 1:NDISC, j = 1:NDISC @@ -153,14 +153,14 @@ function simplexMagic(NDISC::Int=20, NPOINTS::Int=3, checkFileName::AbstractStri sInterpValTest = readdlm(checkFileName) - testErr = sum(abs, sInterpVal-sInterpValTest) + testErr = sum(abs, sInterpVal - sInterpValTest) if (testErr > eps) display("Failed Simplex Comparison Test") display("Test error = $testErr") end - return testErr x == 1, weights) - @test length(full_weight_indices) == 1 - @test weights[full_weight_indices[1]] == 1.0 + @test length(full_weight_indices) == 1 + @test weights[full_weight_indices[1]] == 1.0 - indices, weights = interpolants(grid, [1.5,3]) + indices, weights = interpolants(grid, [1.5, 3]) full_weight_indices = findall(x -> isapprox(x, 0.666667, rtol=1e-5), weights) - @test length(full_weight_indices) == 1 - @test isapprox(weights[full_weight_indices[1]], 0.666667, rtol=1e-5) == true + @test length(full_weight_indices) == 1 + @test isapprox(weights[full_weight_indices[1]], 0.666667, rtol=1e-5) == true full_weight_indices = findall(x -> isapprox(x, 0.333333, rtol=1e-5), weights) - @test length(full_weight_indices) == 1 - @test isapprox(weights[full_weight_indices[1]], 0.333333, rtol=1e-5) == true + @test length(full_weight_indices) == 1 + @test isapprox(weights[full_weight_indices[1]], 0.333333, rtol=1e-5) == true - @test interpolate(grid, [1,2,3,4], [1,1]) == 1.0 - @test maskedInterpolate(grid, [1,2,3,4], [1,1], BitArray([false, false, false, false])) == 1.0 - @test isnan(maskedInterpolate(grid, [1,2,3,4], [1,1], BitArray([true, false, false, false]))) + @test interpolate(grid, [1, 2, 3, 4], [1, 1]) == 1.0 + @test maskedInterpolate(grid, [1, 2, 3, 4], [1, 1], BitArray([false, false, false, false])) == 1.0 + @test isnan(maskedInterpolate(grid, [1, 2, 3, 4], [1, 1], BitArray([true, false, false, false]))) - @test isapprox(interpolate(grid, [1,2,3,4], [1.5,3]), 1.66666666, rtol=1e-5) - @test maskedInterpolate(grid, [1,2,3,4], [1.5,3], BitArray([true, false, false, false])) == 3 - @test maskedInterpolate(grid, [1,2,3,4], [1.5,3], BitArray([false, false, true, false])) == 1 + @test isapprox(interpolate(grid, [1, 2, 3, 4], [1.5, 3]), 1.66666666, rtol=1e-5) + @test maskedInterpolate(grid, [1, 2, 3, 4], [1.5, 3], BitArray([true, false, false, false])) == 3 + @test maskedInterpolate(grid, [1, 2, 3, 4], [1.5, 3], BitArray([false, false, true, false])) == 1 end @testset "Simplex Impl" begin test_simplex_implemented() @@ -322,14 +322,14 @@ end # check whether rectangle & simplex expect sorted order of cutpoints on dimensions function test_ordering(grid) - @test ind2x(grid,1) == [2,18] # 1 - @test ind2x(grid,2) == [5,18] # 2 - @test ind2x(grid,3) == [2,15] # 3 - @test ind2x(grid,4) == [5,15] # 4 - @test ind2x(grid,5) == [2,12] # 5 - @test ind2x(grid,6) == [5,12] # 6 - - return interpolate(grid, [1,2,3,4,5,6], [6, 20]) == 2.0 + @test ind2x(grid, 1) == [2, 18] # 1 + @test ind2x(grid, 2) == [5, 18] # 2 + @test ind2x(grid, 3) == [2, 15] # 3 + @test ind2x(grid, 4) == [5, 15] # 4 + @test ind2x(grid, 5) == [2, 12] # 5 + @test ind2x(grid, 6) == [5, 12] # 6 + + return interpolate(grid, [1, 2, 3, 4, 5, 6], [6, 20]) == 2.0 end # tests the order of vertices returned by vertices() @@ -340,18 +340,18 @@ function test_vertices_ordering(grid) @test length(grid_verts) == length(grid) - for i = 1 : length(grid) - @test grid_verts[i] == ind2x(grid,i) + for i = 1:length(grid) + @test grid_verts[i] == ind2x(grid, i) end return true end @testset "Ordering" begin - @test_throws ErrorException test_ordering( RectangleGrid([2,5], [18,15,12]) ) - @test_throws ErrorException test_ordering( SimplexGrid([2,5], [18,15,12]) ) + @test_throws ErrorException test_ordering(RectangleGrid([2, 5], [18, 15, 12])) + @test_throws ErrorException test_ordering(SimplexGrid([2, 5], [18, 15, 12])) - @test test_vertices_ordering( RectangleGrid([2,5], [12,15,18]) ) == true - @test test_vertices_ordering( SimplexGrid([1,4,6,10], [8,12,17]) ) == true + @test test_vertices_ordering(RectangleGrid([2, 5], [12, 15, 18])) == true + @test test_vertices_ordering(SimplexGrid([1, 4, 6, 10], [8, 12, 17])) == true end @testset "NaN" begin @@ -361,46 +361,55 @@ end @testset "Iteration and Indexing" begin @testset "2d" begin - g2 = GridInterpolations.RectangleGrid(1:3,4:5) - pts = ([1.,4.],[2.,4.],[3.,4.],[1.,5.],[2.,5.],[3.,5.]) + g2 = GridInterpolations.RectangleGrid(1:3, 4:5) + pts = ([1.0, 4.0], [2.0, 4.0], [3.0, 4.0], [1.0, 5.0], [2.0, 5.0], [3.0, 5.0]) @test length(g2) == length(pts) - @test size(g2) == (3,2) + @test size(g2) == (3, 2) @test all(pts .== [x for x in g2]) - @test iterate(g2,7) === nothing - @test all(pts .== [g2[ci] for ci in reshape(CartesianIndices((3,2)), length(g2))]) - @test all(pts .== reshape([g2[i,j] for i=1:3, j=1:2], length(g2))) - @test_throws BoundsError g2[CartesianIndex(0,0)] - @test_throws BoundsError g2[CartesianIndex(4,2)] - @test_throws BoundsError g2[CartesianIndex(3,3)] - @test_throws BoundsError g2[0,0] - @test_throws BoundsError g2[4,2] - @test_throws BoundsError g2[3,3] - @test_throws BoundsError g2[1,1,2] + @test iterate(g2, 7) === nothing + @test all(pts .== [g2[ci] for ci in reshape(CartesianIndices((3, 2)), length(g2))]) + @test all(pts .== reshape([g2[i, j] for i = 1:3, j = 1:2], length(g2))) + @test_throws BoundsError g2[CartesianIndex(0, 0)] + @test_throws BoundsError g2[CartesianIndex(4, 2)] + @test_throws BoundsError g2[CartesianIndex(3, 3)] + @test_throws BoundsError g2[0, 0] + @test_throws BoundsError g2[4, 2] + @test_throws BoundsError g2[3, 3] + @test_throws BoundsError g2[1, 1, 2] end @testset "3d" begin - g3 = GridInterpolations.RectangleGrid(1:3,4:5,6:10) - pts = reshape([[x,y,z] for x=1.:3., y=4.:5., z=6.:10.], length(g3)) + g3 = GridInterpolations.RectangleGrid(1:3, 4:5, 6:10) + pts = reshape([[x, y, z] for x = 1.0:3.0, y = 4.0:5.0, z = 6.0:10.0], length(g3)) @test length(g3) == length(pts) - @test size(g3) == (3,2,5) + @test size(g3) == (3, 2, 5) @test all(pts .== [x for x in g3]) - @test iterate(g3,31) === nothing - @test all(pts .== [g3[ci] for ci in reshape(CartesianIndices((3,2,5)), length(g3))]) - @test all(pts .== reshape([g3[i,j,k] for i=1:3, j=1:2, k=1:5], length(g3))) - @test_throws BoundsError g3[CartesianIndex(0,0,0)] - @test_throws BoundsError g3[CartesianIndex(4,2,5)] - @test_throws BoundsError g3[CartesianIndex(3,3,5)] - @test_throws BoundsError g3[CartesianIndex(3,2,6)] - @test_throws BoundsError g3[0,0,0] - @test_throws BoundsError g3[4,2,5] - @test_throws BoundsError g3[3,3,5] - @test_throws BoundsError g3[3,2,6] + @test iterate(g3, 31) === nothing + @test all(pts .== [g3[ci] for ci in reshape(CartesianIndices((3, 2, 5)), length(g3))]) + @test all(pts .== reshape([g3[i, j, k] for i = 1:3, j = 1:2, k = 1:5], length(g3))) + @test_throws BoundsError g3[CartesianIndex(0, 0, 0)] + @test_throws BoundsError g3[CartesianIndex(4, 2, 5)] + @test_throws BoundsError g3[CartesianIndex(3, 3, 5)] + @test_throws BoundsError g3[CartesianIndex(3, 2, 6)] + @test_throws BoundsError g3[0, 0, 0] + @test_throws BoundsError g3[4, 2, 5] + @test_throws BoundsError g3[3, 3, 5] + @test_throws BoundsError g3[3, 2, 6] if VERSION.major > 0 # this test is known to fail for julia v0.7 - @test_throws BoundsError g3[1,1] + @test_throws BoundsError g3[1, 1] end - @test_throws BoundsError g3[1,1,1,2] + @test_throws BoundsError g3[1, 1, 1, 2] end end +@testset "Data types" begin + grid = RectangleGrid(1.0:10, 1.0:10) + data = rand(Complex{Float64}, 10, 10) + @test typeof(interpolate(grid, data, [-1.5, 1.5])) == Complex{Float64} + + data = rand(Float32, 10, 10) + @test typeof(interpolate(grid, data, [-1.5f0, 1.5f0])) == Float32 +end + include(joinpath(@__DIR__, "..", "bench", "interpBenchmarks.jl")) compareBenchmarks(4, 10, 100, quiet=false)