Add InfiniteInterpolate as an extension (#382)

* Update Fishing example docs

* Adjust unit test of fishing example to test output values within tolerance

* Update Hanging Chain example docs

* Update Jennings example docs

* Added unit tests to examples (excluding pandemic control)

* Update tolerance to 1E-06

* Add InfiniteInterpolate with basic value function extension based on linear interpolation + test case

* Add argument so that interpolation method is user defined

* Add logic that throws an error for unsupported interpolation methods and handling for variables that don't rely on infinite parameters

* Added unit tests for catching irregular grid method errors and checking interpolated values

* Include InfiniteInterpolate test in runtests.jl

* Add InfiniteInterpolate with basic value function extension based on linear interpolation + test case

* Add argument so that interpolation method is user defined

* Add logic that throws an error for unsupported interpolation methods and handling for variables that don't rely on infinite parameters

* Added unit tests for catching irregular grid method errors and checking interpolated values

* Include InfiniteInterpolate test in runtests.jl

* Add docstring to InfiniteInterpolate ext & unit test for unsupported interpolation method error

* Add InfiniteInterpolate with basic value function extension based on linear interpolation + test case

* Add argument so that interpolation method is user defined

* Add logic that throws an error for unsupported interpolation methods and handling for variables that don't rely on infinite parameters

* Added unit tests for catching irregular grid method errors and checking interpolated values

* Include InfiniteInterpolate test in runtests.jl

* Add InfiniteInterpolate with basic value function extension based on linear interpolation + test case

* Added unit tests for catching irregular grid method errors and checking interpolated values

* Include InfiniteInterpolate test in runtests.jl

* Add docstring to InfiniteInterpolate ext & unit test for unsupported interpolation method error

* Add InfiniteInterpolate to Project.toml + use InfiniteInteprolate value func in hovercraft example + shorten doc string

* Change Ipopt to mock optimizer in InfiniteInterpolate test case

* Add semi-infinite, point, finite & derivative variables to InfiniteInterpolate test case

* Add a small example on using InfiniteInterpolate in docs > src > guide > result.md

* Add Interpolations.jl to [extras] and [targets] in Project.toml

* Move Interpolations.jl import statement from runtests.jl to InfiniteInterpolate test file

* Update results.md doc about type piracy issue with InfiniteInterpolate

* Change Matrix datatype -> Array in InfiniteInterpolate wrapper function

* Adjust documentation for InfiniteInterpolate in guide/result.md and Interpolations version in project.toml files

* Add orthogonal collocation points to InfiniteInterpolate test case

* Add degree + fallback functions to InfiniteInterpolate extension; update docs as necessary

* Adjust InfiniteInterpolate documentation

* Adjust make.jl to load in InfiniteInterpolate extension

* Adjust hyperlinks to InfiniteInterpolate extension in technical manual

* Simplify InfiniteInterpolate case in hovercraft example + InfiniteInterpolate documentation wording

* Get rid of if statements in InfiniteInterpolate code and better utilize multiple dispatch + update docs accordingly

* Add finite parameter unit test for InfiniteInterpolate case

* Break up code lines to be around 80 char limit

* Adjust line breaks in make.jl

* Uncomment unit tests in runtests.jl

---------

Co-authored-by: Joshua Pulsipher <[email protected]>

Author: wenwen0231

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,12 @@ 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,24 @@ 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 are used as the interpolation method for both x1 and x2:
+using Interpolations
+xFunc = value.(x, cubic_spline_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 +88,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.6386021193421085, 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,57 @@ 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. 
 
+## Interpolation-Based Continuous Values
+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,11 @@ JuMP.reduced_cost(::GeneralVariableRef)
 JuMP.optimizer_index(::GeneralVariableRef)
 ```
 
+## [InfiniteInterpolate](@id infiniteInterpolate)
+```@docs
+JuMP.value(::GeneralVariableRef, ::Interpolations.InterpolationType)
+```
+
 ## Constraints
 ```@docs
 JuMP.has_duals(::InfiniteModel)

ext/InfiniteInterpolate.jl

@@ -0,0 +1,169 @@
+module InfiniteInterpolate
+
+import JuMP
+import InfiniteOpt
+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::Interpolations.InterpolationType);
+        [kwargs...])
+
+Extend `JuMP.value` to return `vref` as an interpolation object from
+Interpolations.jl, based on `method` which specifies the interpolation method. 
+Currently supported method(s) are:
+- `linear_interpolation`
+- `constant_interpolation`
+- `cubic_spline_interpolation`
+
+All methods support equidistant grid points. However, irregular grid points can 
+only be used with `linear_interpolation` and `constant_interpolation`.
+
+```julia
+JuMP.value(vref::GeneralVariableRef, degree::Interpolations.Degree; kwargs...)
+```
+Extend `JuMP.value` to return `vref` as an interpolation object from 
+Interpolations.jl, based on `degree` which specifies the degree of interpolation.
+The currently supported degrees are:
+- `Linear()`
+- `Constant()`
+- `Cubic()`
+
+All methods support equidistant grid points. However, irregular grid points can 
+only be used with `Linear()` and `Constant()`.
+
+**Examples**
+```julia-repl
+julia> zFunc = value(z, cubic_spline_interpolation)
+
+julia> zFunc(5.4);
+42.0
+
+julia> zFunc2 = value(z, Cubic())
+
+julia> zFunc2(5.4)
+42.0
+```
+"""
+function JuMP.value(vref::InfiniteOpt.GeneralVariableRef,
+                    method::IP.InterpolationType; kwargs...)
+    throw(ArgumentError("Unsupported interpolation type: $(method). 
+        Supported interpolation types are: linear_interpolation, 
+        constant_interpolation, and cubic_spline_interpolation."))
+end
+
+function JuMP.value(vref::InfiniteOpt.GeneralVariableRef,
+        interpMethod::typeof(IP.cubic_spline_interpolation);
+        kwargs...)
+    InfiniteOpt._check_result_is_current(JuMP.owner_model(vref), JuMP.value)
+
+    # Get infinite parameter references for which the variable depends on
+    prefs = InfiniteOpt.parameter_refs(vref)
+
+    if isempty(prefs)
+        # If no infinite parameters, return the value directly
+        return InfiniteOpt._get_value(vref,
+                                      InfiniteOpt._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
+            suppsLabel = first(InfiniteOpt.core_object(pref).supports)[2]
+            if !(InfiniteOpt.UniformGrid in suppsLabel)
+                throw(ArgumentError("Interpolation method $(interpMethod) does 
+                not support irregular grids for supports. Please specify a 
+                uniform grid or choose a different interpolation method."))
+            end
+
+            # Create an equidistant range for support values
+            paramVals = InfiniteOpt._get_value(pref,
+                                    InfiniteOpt._index_type(pref);
+                                    kwargs...)
+            numSupps = length(paramVals)
+            paramRange = LinRange(paramVals[1], paramVals[end], numSupps)
+            push!(Vparams, paramRange)
+        end
+
+        # Ensure Vparams is a tuple for interpolation
+        Vparams = Tuple(Vparams)
+
+        # Get the variable supports
+        Vsupps = InfiniteOpt._get_value(vref,
+                                        InfiniteOpt._index_type(vref);
+                                        kwargs...)
+
+        # Pass the parameter and variable values to interpolation function
+        varFunc = interpMethod(Vparams, Vsupps)
+        return varFunc
+    end
+end
+
+function JuMP.value(vref::InfiniteOpt.GeneralVariableRef,
+    interpMethod::_irregularGridMethods; kwargs...)
+    InfiniteOpt._check_result_is_current(JuMP.owner_model(vref), JuMP.value)
+
+    # Get infinite parameter references for which the variable depends on
+    prefs = InfiniteOpt.parameter_refs(vref)
+
+    if isempty(prefs)
+        # If no infinite parameters, return the value directly
+        return InfiniteOpt._get_value(vref, InfiniteOpt._index_type(vref); kwargs...)
+    else
+        # Get the parameter supports
+        Vparams = []
+        for pref in prefs
+            paramVals = InfiniteOpt._get_value(pref,
+                                               InfiniteOpt._index_type(pref);
+                                               kwargs...)
+            # Directly pass in a vector of support values, which may be irregularly spaced
+            push!(Vparams, paramVals)
+        end
+
+        # Ensure Vparams is a tuple for interpolation
+        Vparams = Tuple(Vparams)
+
+        # Get the variable supports
+        Vsupps = InfiniteOpt._get_value(vref, 
+                                        InfiniteOpt._index_type(vref); 
+                                        kwargs...)
+
+        # Pass the parameter and variable values to interpolation function
+        varFunc = interpMethod(Vparams, Vsupps)
+        return varFunc
+    end
+end
+
+# Fallback for unsupported interpolation degrees
+function JuMP.value(vref::InfiniteOpt.GeneralVariableRef,
+    degree::IP.Degree; kwargs...)
+    throw(ArgumentError("Unsupported interpolation degree: $(degree). Supported 
+    interpolation degrees are: Linear(), Constant(), and Cubic()."))
+end
+
+function JuMP.value(vref::InfiniteOpt.GeneralVariableRef,
+                    degree::IP.Linear;
+                    kwargs...)
+    return JuMP.value(vref, IP.linear_interpolation; kwargs...) 
+end
+
+function JuMP.value(vref::InfiniteOpt.GeneralVariableRef,
+                    degree::IP.Constant;
+                    kwargs...)
+    return JuMP.value(vref, IP.constant_interpolation; kwargs...) 
+end
+
+function JuMP.value(vref::InfiniteOpt.GeneralVariableRef,
+                    degree::IP.Cubic;
+                    kwargs...)
+    return JuMP.value(vref, IP.cubic_spline_interpolation; kwargs...) 
+end
+end