Add InfiniteInterpolate as an extension

Project.toml

@@ -12,6 +12,12 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
+[weakdeps]
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+
+[extensions]
+InfiniteInterpolate = "Interpolations"
+
 [compat]
 DataStructures = "0.14.2 - 0.18, 0.19"
 Distributions = "0.19 - 0.25"
@@ -20,12 +26,14 @@ JuMP = "1.18"
 MutableArithmetics = "1"
 Reexport = "0.2, 1"
 julia = "^1.6"
+Interpolations = "0.16"
 
 [extras]
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 
 [targets]
-test = ["Pkg", "Test", "Random", "Suppressor"]
+test = ["Pkg", "Test", "Random", "Suppressor", "Interpolations"]

docs/Project.toml

@@ -11,6 +11,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 
 [compat]
 Distributions = "0.25"
@@ -23,3 +24,4 @@ JuMP = "1.23"
 Literate = "2.18"
 Plots = "1"
 SpecialFunctions = "2"
+Interpolations = "0.16"

docs/make.jl

@@ -1,4 +1,4 @@
-using Documenter, InfiniteOpt, Distributions, Literate, Random
+using Documenter, InfiniteOpt, Distributions, Literate, Random, Interpolations
 
 if !@isdefined(EXAMPLE_DIR)
     const EXAMPLE_DIR = joinpath(@__DIR__, "src", "examples")
@@ -96,7 +96,8 @@ makedocs(;
         ],
     sitename = "InfiniteOpt.jl",
     authors = "Joshua Pulsipher and Weiqi Zhang",
-    modules = [InfiniteOpt],
+    modules = [InfiniteOpt, isdefined(Base, :get_extension) ? Base.get_extension(InfiniteOpt, :InfiniteInterpolate) :
+        InfiniteOpt.InfiniteInterpolate,],
     checkdocs = :none,
     linkcheck = true,
     linkcheck_ignore = [r"https://www.youtube.com/.*"],

docs/src/examples/Optimal Control/hovercraft.jl

@@ -61,15 +61,21 @@ end)
 # Optimize the model:
 optimize!(m)
 
-# Extract the results:
-x_opt = value.(x);
+# Extract the results. The [`InfiniteInterpolate`](@ref infiniteInterpolate) extension can be used to get a smooth interpolated function for x, which is invoked when both the `Interpolations` and `InfiniteOpt` packages are imported. Here, cubic splines and linear were chosen as the interpolation methods for x1 and x2, respectively:
+using Interpolations
+xFunc = value.(x, (cubic_spline_interpolation, linear_interpolation));
+
+# Query our interpolated function for the values of x1 and x2:
+tvals = LinRange(0, 60, 100)
+x1Vals = xFunc[1](tvals)
+x2Vals = xFunc[2](tvals);
 
 # Plot the results:
 using Plots
-scatter(xw[1,:], xw[2,:], label = "Waypoints", background_color = :transparent)
-plot!(x_opt[1], x_opt[2], label = "Trajectory")
+scatter(xw[1,:], xw[2,:], label = "Waypoints")
 xlabel!("x_1")
 ylabel!("x_2")
+plot!(x1Vals, x2Vals, label = "Trajectory")
 
 # That's it, now we have our optimal trajectory!
 
@@ -79,8 +85,10 @@ using Test
 tol = 1E-6
 @test termination_status(m) == MOI.LOCALLY_SOLVED
 @test has_values(m)
-@test x_opt isa Vector{<:Vector{<:Real}}
+@test x1Vals isa Vector{<:Real}
+@test x2Vals isa Vector{<:Real}
 @test isapprox(objective_value(m), 0.043685293177035435, atol=tol)
 @test isapprox(value(u[1])[end], -0.010503853944039986, atol=tol)
 @test isapprox(value(u[2])[end], 0.005456780217220367, atol=tol)
-
+@test isapprox(x1Vals[15], 1.3837274935883543, atol=tol)
+@test isapprox(x2Vals[15], 1.6392615199958938, atol=tol)

docs/src/guide/result.md

@@ -137,7 +137,7 @@ parameter dependencies.
     This also holds true for many other methods in `InfiniteOpt`. For example, 
     adding the dot also vectorizes `dual` and `set_binary`.
 
-We can also query the dual of a constraint via 
+We can query the dual of a constraint via 
 [`dual`](@ref JuMP.dual(::InfOptConstraintRef)) if a model has duals available
 as indicated by [`has_duals`](@ref):
 ```jldoctest results
@@ -184,6 +184,48 @@ These again all have a 1-to-1 correspondence.
     parameters, an n-dimensional array will typically be returned
     whose dimensions correspond to the supports of the infinite parameters. 
 
+## InfiniteInterpolate queries
+We can also get a continuous representation of a variable as an interpolations object from the `Interpolations.jl` package. This is based on the [`InfiniteInterpolate`](@ref infiniteInterpolate) extension, which is automatically loaded in when both `InfiniteOpt` and `Interpolations` are imported. The current supported methods are `linear_interpolation`, `constant_interpolation` and `cubic_spline_interpolation`.
+```jldoctest results
+julia> using Interpolations
+
+julia> yFunc = value(y, linear_interpolation)
+10-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
+ 42.0
+ 20.999999995627107
+ 20.999999995628606
+ 20.999999995628603
+ 20.999999995628592
+ 20.999999995628603
+ 20.999999995634035
+ 20.999999995620904
+ 20.99999999562204
+ 20.9999999956286
+
+julia> yFunc(5.12)
+20.9999999956286
+```
+Alternatively, we can specify the degree of interpolation with `Linear()`, `Constant()` or `Cubic()`. This will call the corresponding interpolation method.
+```jldoctest results
+julia> yFunc = value(y, Linear()) # equivalent to value(y, linear_interpolation)
+10-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
+ 42.0
+ 20.999999995627107
+ 20.999999995628606
+ 20.999999995628603
+ 20.999999995628592
+ 20.999999995628603
+ 20.999999995634035
+ 20.999999995620904
+ 20.99999999562204
+ 20.9999999956286
+
+julia> yFunc(5.12)
+20.9999999956286
+```
+!!! warning
+    There is a type piracy conflict between JuMP and OffsetArrays (a dependency of Interpolations.jl). As a result, type piracy issues may arise when Interpolations is loaded in. Please keep this in mind when using the InfiniteInterpolate extension.
+
 ## Termination Queries
 Termination queries are those that question about how the infinite model was 
 solved and what its optimized state entails. Programmatically, such queries on 

docs/src/manual/result.md

@@ -36,6 +36,12 @@ JuMP.reduced_cost(::GeneralVariableRef)
 JuMP.optimizer_index(::GeneralVariableRef)
 ```
 
+## [InfiniteInterpolate](@id infiniteInterpolate)
+```@docs
+JuMP.value(::GeneralVariableRef, ::Union{typeof(Interpolations.constant_interpolation), typeof(Interpolations.cubic_spline_interpolation), typeof(Interpolations.linear_interpolation)})
+JuMP.value(::GeneralVariableRef, ::Interpolations.Degree)
+```
+
 ## Constraints
 ```@docs
 JuMP.has_duals(::InfiniteModel)

ext/InfiniteInterpolate.jl

@@ -0,0 +1,100 @@
+module InfiniteInterpolate
+
+import JuMP
+import InfiniteOpt as infOpt
+import Interpolations as IP
+
+const _irregularGridMethods = Union{typeof(IP.linear_interpolation),
+                            typeof(IP.constant_interpolation)}
+
+const _convenienceConstructors = Union{typeof(IP.linear_interpolation),
+                            typeof(IP.constant_interpolation),
+                            typeof(IP.cubic_spline_interpolation)}
+
+"""
+    JuMP.value(vref::GeneralVariableRef, method::Union{typeof(Interpolations.constant_interpolation), typeof(Interpolations.cubic_spline_interpolation), typeof(Interpolations.linear_interpolation)}; [kwargs...])
+
+Extend `JuMP.value` to return `vref` as an interpolation object from Interpolations.jl, where the `method` argument specifies the interpolation method to be used.
+The currently supported methods are `linear_interpolation`, `constant_interpolation` and `cubic_spline_interpolation`.
+
+**Example**
+```julia-repl
+julia> zFunc = value(z, cubic_spline_interpolation)
+julia> zFunc(5.4)
+42.0
+```
+"""
+function JuMP.value(vref::infOpt.GeneralVariableRef, interpMethod::_convenienceConstructors; kwargs...)
+    infOpt._check_result_is_current(JuMP.owner_model(vref), JuMP.value)
+
+    # Get infinite parameter references for which the variable depends on
+    prefs = infOpt.parameter_refs(vref)
+
+    if isempty(prefs)
+        # If no infinite parameters, return the value directly
+        return infOpt._get_value(vref, infOpt._index_type(vref); kwargs...)
+    else
+        # Get the parameter supports
+        Vparams = []
+        for pref in prefs
+            # If user defined irregular grid points for supports, throw an error if the specified method doesn't support it
+            suppsLabel = first(infOpt.core_object(pref).supports)[2]
+            if !(infOpt.UniformGrid in suppsLabel) && !(interpMethod isa _irregularGridMethods)
+                throw(ArgumentError("Interpolation method $(interpMethod) does not support irregular grids for supports. Please specify a uniform grid or choose a different interpolation method."))
+            end
+
+            paramVals = infOpt._get_value(pref, infOpt._index_type(pref); kwargs...)
+            numSupps = length(paramVals)
+            if interpMethod isa _irregularGridMethods
+                # Can directly pass in a vector of support values, which may be irregularly spaced
+                push!(Vparams, paramVals)
+            else
+                # Create an equidistant range for support values
+                paramRange = LinRange(paramVals[1], paramVals[end], numSupps)
+                push!(Vparams, paramRange)
+            end
+        end
+
+        # Ensure Vparams is a tuple for interpolation
+        Vparams = Tuple(Vparams)
+
+        # Get the variable supports
+        Vsupps = infOpt._get_value(vref, infOpt._index_type(vref); kwargs...)
+
+        # Pass the parameter and variable values to interpolation function
+        varFunc = interpMethod(Vparams, Vsupps)
+        return varFunc
+    end
+end
+
+"""
+    JuMP.value(vref::GeneralVariableRef, degree::Interpolations.Degree; [kwargs...])
+
+Extend `JuMP.value` to return `vref` as an interpolation object from Interpolations.jl, where the `degree` argument specifies the degree of interpolation. The corresponding interpolation method is called depending on the degree.
+The currently supported degrees are `Linear()`, `Constant()` and `Cubic()`.
+
+**Example**
+```julia-repl
+julia> zFunc = value(z, Cubic())
+julia> zFunc(5.4)
+42.0
+```
+"""
+function JuMP.value(vref::infOpt.GeneralVariableRef, degree::IP.Degree; kwargs...)
+    if degree == IP.Linear()
+        return JuMP.value(vref, IP.linear_interpolation; kwargs...)
+    elseif degree == IP.Constant()
+        return JuMP.value(vref, IP.constant_interpolation; kwargs...)
+    elseif degree == IP.Cubic()
+        return JuMP.value(vref, IP.cubic_spline_interpolation; kwargs...)
+    else
+        throw(ArgumentError("Unsupported interpolation degree type: $(degree). Supported degree types are: Linear(), Constant(), and Cubic()."))
+    end
+end
+
+# Fallback for unsupported interpolation types
+function JuMP.value(vref::infOpt.GeneralVariableRef, method::IP.InterpolationType; kwargs...)
+    throw(ArgumentError("Unsupported interpolation type: $(method). Supported interpolation types are: Linear(), Constant(), and Cubic()."))
+end
+
+end