Improve the SparseIR.refine_grid function

[terasakisatoshi]

This PR improves the implementation of SparseIR.refine_grid and adds tests for the function.

[terasakisatoshi]

This change is required because the newn gets a negative value for empty input, i.e., grid=Float64[]; otherwise, refine_grid([]) will crash a Julia process.

julia> using SparseIR

julia> SparseIR.refine_grid(Float64[], Val(2))
ERROR: ArgumentError: invalid GenericMemory size: too large for system address width
Stacktrace:
 [1] GenericMemory
   @ ./boot.jl:516 [inlined]
 [2] Array
   @ ./boot.jl:578 [inlined]
 [3] refine_grid(grid::Vector{Float64}, ::Val{2})
   @ SparseIR ~/.julia/packages/SparseIR/iqwia/src/_roots.jl:45
 [4] top-level scope
   @ REPL[2]:1

[terasakisatoshi]

The current implementation (on the main branch) causes an error.

julia> SparseIR.refine_grid([1,2,3], Val(2))
ERROR: InexactError: Int64(1.5)
Stacktrace:
 [1] Int64
   @ ./float.jl:994 [inlined]
 [2] convert
   @ ./number.jl:7 [inlined]
 [3] setindex!
   @ ./array.jl:987 [inlined]
 [4] refine_grid(grid::Vector{Int64}, ::Val{2})
   @ SparseIR ~/.julia/packages/SparseIR/iqwia/src/_roots.jl:53
 [5] top-level scope
   @ REPL[3]:1

To avoid this issue, newgrid = Vector{float(T)}(undef, newn) is required.

[SamuelBadr]

Thanks, I didn't think of this edge case! However, the function is now type unstable for integer inputs. Also, in the empty case, the output aliases the input.
I'd propose using newn = max(..., 0) instead.

(on my phone, so haven't actually tested)

[terasakisatoshi]

The following implementation causes e.g., max(α * (n - 1) + 1, 0) can't treat for the case isempty(grid) is true

function refine_grid(grid::Vector{T}, ::Val{α}) where {T, α}
    n = length(grid)
    newn = max(α * (n - 1) + 1, 0)
    newgrid = Vector{float(T)}(undef, newn)

    @inbounds for i in 1:(n - 1)
        xb = grid[i]
        xe = grid[i + 1]
        Δx = (xe - xb) / α
        newi = α * (i - 1)
        for j in 1:α
            newgrid[newi + j] = xb + Δx * (j - 1)
        end
    end
    newgrid[end] = last(grid) ## This line can't treat for empty case
    return newgrid
end

refine_grid(Float64[], Val(2))

#=
julia> [][end]
ERROR: BoundsError: attempt to access 0-element Vector{Any} at index [0]
=#

[terasakisatoshi]

What about this?

From

function refine_grid(grid::Vector{T}, ::Val{α}) where {T, α}
    isempty(grid) && return grid

To

function refine_grid(grid::Vector{T}, ::Val{α}) where {T, α}
    isempty(grid) && return T[]

[SamuelBadr]

True, didn't even look at the last line. I like your solution, but suggest changing it to float(T)[] for type stability.

[terasakisatoshi]

but suggest changing it to float(T)[] for type stability.

Yes!!! That makes sense.

[terasakisatoshi]

I've inserted @inferred macro:

refined = @inferred SparseIR.refine_grid(int_grid, Val(2))

[SamuelBadr]

Thank you!