"extended" vector with heterogeneous "padding"

General Usage
1

Hello,

I seem to frequently encounter situations where accounting for boundary/initial conditions requires “extending” vectors in a way similar to padding, but where padded values are not necessarily constant and might not even be of the same type as the values stored in the parent vector.

Latest case in point:

(Minimal code to make the example runnable) 
struct ExtendedVector{N, T, V}
    vec :: V
    ext :: NTuple{N, Pair{Int, T}}
end

extend_vector(vec, args...) = ExtendedVector(vec, args)

function Base.getindex(ev::ExtendedVector, i::Int)
    for (j, v) in ev.ext
        i == j && return v
    end
    return ev.vec[i]
end

julia> using JuMP
julia> m = Model();
julia> N = 5;
julia> @variable(m, x[1:N]);
       
julia> # Everything should behave as if
       #   x̄[0] == 0, and
       #   x̄[N+1] == 1
       x̄ = extend_vector(x, 0=>0, N+1=>1);

julia> # No need for explicit code to take care of boundary conditions
       @constraint(m, [i∈1:N],
                   x̄[i+1] - 2x̄[i] + x̄[i-1] <= 1)
5-element Vector{ConstraintRef{...}}:
 -2 x[1] + x[2] ≤ 1
 x[1] - 2 x[2] + x[3] ≤ 1
 x[2] - 2 x[3] + x[4] ≤ 1
 x[3] - 2 x[4] + x[5] ≤ 1
 x[4] - 2 x[5] ≤ 0

I’m not sure how to call such an “extended vector”. Does there exist anything like this in the ecosystem? Under what name?

I can use the minimal implementation above in my real use case, but I expect things might get more complicated very quickly if I ever need more features (like a full AbstractVector-like behavior). And if something more fancy already exists, I might as well use it instead.

1 Like
2

That kinda sounds like Interpolations.jl extrapolations?

3

Thanks!

To be clear, are you suggesting to refer to this as an “extrapolated” vector instead of an “extended” vector? Or are you implying that Interpolations.jl’s extrapolate function could actually be used for this as it is?

4

Both! But they don’t just pad one or two indices, it’s everything:

julia> using Interpolations

julia> etp = extrapolate(interpolate(1:7, NoInterp()), missing)
7-element extrapolate(interpolate(::UnitRange{Int64}, NoInterp()), missing) with element type Union{Missing, Int64}:
 1
 2
 3
 4
 5
 6
 7

julia> etp[begin]
1

julia> etp[0]
missing

julia> etp[end]
7

julia> etp[end+1]
missing

julia> etp[-1000]
missing

They also permit floating point indexing in ways you may not want:

julia> etp[1.0]
1

julia> etp[0.9]
missing

julia> etp[1.1]
ERROR: InexactError: Int64(1.1)

You may still want a specialize type that explicitly includes the padding as its axes and such.

5

Maybe you can use functionality from Stencils.jl although I have never tried to use it with JUMP.
But there are different ways to implement boundaries (periodic, always some fixed value, or also the option to embed it into a larger array (ghost cells) so the stencil never indexes out of bounds.

6

PS: actually i dont know if its possible to use asymmetrical padding but maybe you can work around it by using a halo and manually setting the value at the halo indices?