Merge remote-tracking branch 'upstream/master'

.github/workflows/CompatHelper.yml

@@ -13,7 +13,7 @@ jobs:
   CompatHelper:
     runs-on: ubuntu-latest
     steps:
-      - uses: julia-actions/setup-julia@3645a07f58c7f83b9f82ac8e0bb95583e69149e6
+      - uses: julia-actions/setup-julia@780022b48dfc0c2c6b94cfee6a9284850107d037
         with:
           version: 1.3
       - name: Pkg.add("CompatHelper")

.github/workflows/Downgrade.yml

@@ -28,7 +28,7 @@ jobs:
           - windows-latest
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v2.2.0
+      - uses: julia-actions/setup-julia@v2.3.0
         with:
           version: ${{ matrix.version }}
       - uses: julia-actions/julia-downgrade-compat@v1

Project.toml

@@ -1,13 +1,12 @@
 name = "DataInterpolations"
 uuid = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
-version = "5.3.0"
+version = "6.0.0"
 
 [deps]
 FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
-ReadOnlyArrays = "988b38a3-91fc-5605-94a2-ee2116b3bd83"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
@@ -16,6 +15,7 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 RegularizationTools = "29dad682-9a27-4bc3-9c72-016788665182"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extensions]
 DataInterpolationsChainRulesCoreExt = "ChainRulesCore"
@@ -25,22 +25,22 @@ DataInterpolationsSymbolicsExt = "Symbolics"
 
 [compat]
 Aqua = "0.8"
-ChainRulesCore = "1.18"
+ChainRulesCore = "1.24"
 FindFirstFunctions = "1.1"
 FiniteDifferences = "0.12.31"
 ForwardDiff = "0.10.36"
 LinearAlgebra = "1.10"
 Optim = "1.6"
 PrettyTables = "2"
 QuadGK = "2.9.1"
-ReadOnlyArrays = "0.2.0"
 RecipesBase = "1.3"
 Reexport = "1"
 RegularizationTools = "0.6"
 SafeTestsets = "0.1"
 StableRNGs = "1"
 Symbolics = "5.29"
 Test = "1"
+Zygote = "0.6.70"
 julia = "1.10"
 
 [extras]
@@ -55,6 +55,7 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics"]
+test = ["Aqua", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Zygote"]

docs/Project.toml

@@ -1,15 +1,23 @@
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RegularizationTools = "29dad682-9a27-4bc3-9c72-016788665182"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-DataInterpolations = "5"
+DataInterpolations = "6"
 Documenter = "1"
+ModelingToolkit = "9"
+ModelingToolkitStandardLibrary = "2"
 Optim = "1"
+OrdinaryDiffEq = "6"
 Plots = "1"
 RegularizationTools = "0.6"
-StableRNGs = "1"
+StableRNGs = "1"
+Symbolics = "5.29"

docs/make.jl

@@ -13,6 +13,7 @@ makedocs(modules = [DataInterpolations],
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/DataInterpolations/stable/"),
     pages = ["index.md", "Methods" => "methods.md",
-        "Interface" => "interface.md", "Manual" => "manual.md", "Inverting Integrals" => "inverting_integrals.md"])
+        "Interface" => "interface.md", "Using with Symbolics/ModelingToolkit" => "symbolics.md",
+        "Manual" => "manual.md", "Inverting Integrals" => "inverting_integrals.md"])
 
 deploydocs(repo = "github.com/SciML/DataInterpolations.jl"; push_preview = true)

docs/src/index.md

@@ -1,10 +1,6 @@
 # DataInterpolations.jl
 
-DataInterpolations.jl is a library for performing interpolations of one-dimensional data. By
-"data interpolations" we mean techniques for interpolating possibly noisy data, and thus
-some methods are mixtures of regressions with interpolations (i.e. do not hit the data
-points exactly, smoothing out the lines). This library can be used to fill in intermediate
-data points in applications like timeseries data.
+DataInterpolations.jl is a library for performing interpolations of one-dimensional data. Interpolations are a very important component of many modeling workflows. Often, sampled or measured inputs need to be transformed into continuous functions or smooth curves for simulation purposes. In many scientific machine learning workflows, interpolating data is essential to learn continuous models. DataInterpolations.jl can be used for facilitating these types of workflows. By "data interpolations" we mean techniques for interpolating possibly noisy data, and thus some methods are mixtures of regressions with interpolations (i.e. do not hit the data points exactly, smoothing out the lines).
 
 ## Installation
 

docs/src/interface.md

@@ -35,22 +35,7 @@ A2(300.0)
 
     The values computed beyond the range of the time points provided during interpolation will not be reliable, as these methods only perform well within the range and the first/last piece polynomial fit is extrapolated on either side which might not reflect the true nature of the data.
 
-The keyword `safetycopy = false` can be passed to make sure no copies of `u` and `t` are made when initializing the interpolation object.
-
-```@example interface
-A3 = QuadraticInterpolation(u, t; safetycopy = false)
-
-# Check for same memory
-u === A3.u.parent
-```
-
-Note that this does not prevent allocation in every interpolation constructor call, because parameter values are cached for all interpolation types except [`ConstantInterpolation`](@ref).
-
-Because of the caching of parameters which depend on `u` and `t`, this data should not be mutated. Therefore `u` and `t` are wrapped in a `ReadOnlyArray` from [ReadOnlyArrays.jl](https://github.com/JuliaArrays/ReadOnlyArrays.jl).
-
-```@repl interface
-A3.t[2] = 3.14
-```
+The keyword `cache_parameters = true` can be passed to precalculate parameters at initialization, making evalations cheaper to compute. This is not compatible with modifying `u` and `t`. The default `cache_parameters = false` does however not prevent allocation in every interpolation constructor call.
 
 ## Derivatives
 

docs/src/symbolics.md

@@ -0,0 +1,65 @@
+# Using DataInterpolations.jl with Symbolics.jl and ModelingToolkit.jl
+
+All interpolation methods can be integrated with [Symbolics.jl](https://symbolics.juliasymbolics.org/stable/) and [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) seamlessly.
+
+## Using with Symbolics.jl
+
+### Expressions
+
+```@example symbolics
+using DataInterpolations, Symbolics
+using Test
+
+u = [0.0, 1.5, 0.0]
+t = [0.0, 0.5, 1.0]
+A = LinearInterpolation(u, t)
+
+@variables τ
+
+# Simple Expression
+ex = cos(τ) * A(τ)
+@test substitute(ex, Dict(τ => 0.5)) == cos(0.5) * A(0.5) # true
+```
+
+### Symbolic Derivatives
+
+```@example symbolics
+D = Differential(τ)
+
+ex1 = A(τ)
+
+# Derivative of interpolation
+ex2 = expand_derivatives(D(ex1))
+
+@test substitute(ex2, Dict(τ => 0.5)) == DataInterpolations.derivative(A, 0.5) # true
+
+# Higher Order Derivatives
+ex3 = expand_derivatives(D(D(A(τ))))
+
+@test substitute(ex3, Dict(τ => 0.5)) == DataInterpolations.derivative(A, 0.5, 2) # true
+```
+
+## Using with ModelingToolkit.jl
+
+Most common use case with [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) is to plug in interpolation objects as input functions. This can be done using `TimeVaryingFunction` component of [ModelingToolkitStandardLibrary.jl](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/).
+
+```@example mtk
+using DataInterpolations
+using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using OrdinaryDiffEq
+
+us = [0.0, 1.5, 0.0]
+times = [0.0, 0.5, 1.0]
+A = LinearInterpolation(us, times)
+
+@named src = TimeVaryingFunction(A)
+vars = @variables x(t) out(t)
+eqs = [out ~ src.output.u, D(x) ~ 1 + out]
+@named sys = ODESystem(eqs, t, vars, []; systems = [src])
+
+sys = structural_simplify(sys)
+prob = ODEProblem(sys, [x => 0.0], (times[1], times[end]))
+sol = solve(prob)
+```

ext/DataInterpolationsChainRulesCoreExt.jl

@@ -1,27 +1,97 @@
 module DataInterpolationsChainRulesCoreExt
-
 if isdefined(Base, :get_extension)
     using DataInterpolations: _interpolate, derivative, AbstractInterpolation,
+                              LinearInterpolation, QuadraticInterpolation,
                               LagrangeInterpolation, AkimaInterpolation,
-                              BSplineInterpolation, BSplineApprox
+                              BSplineInterpolation, BSplineApprox, get_idx, get_parameters,
+                              _quad_interp_indices
     using ChainRulesCore
 else
     using ..DataInterpolations: _interpolate, derivative, AbstractInterpolation,
+                                LinearInterpolation, QuadraticInterpolation,
                                 LagrangeInterpolation, AkimaInterpolation,
-                                BSplineInterpolation, BSplineApprox
+                                BSplineInterpolation, BSplineApprox, get_parameters,
+                                _quad_interp_indices
     using ..ChainRulesCore
 end
 
+function ChainRulesCore.rrule(
+        ::Type{LinearInterpolation}, u, t, I, p, extrapolate, cache_parameters)
+    A = LinearInterpolation(u, t, I, p, extrapolate, cache_parameters)
+    function LinearInterpolation_pullback(ΔA)
+        df = NoTangent()
+        du = ΔA.u
+        dt = NoTangent()
+        dI = NoTangent()
+        dp = NoTangent()
+        dextrapolate = NoTangent()
+        dcache_parameters = NoTangent()
+        df, du, dt, dI, dp, dextrapolate, dcache_parameters
+    end
+
+    A, LinearInterpolation_pullback
+end
+
+function ChainRulesCore.rrule(
+        ::Type{QuadraticInterpolation}, u, t, I, p, mode, extrapolate, cache_parameters)
+    A = QuadraticInterpolation(u, t, I, p, mode, extrapolate, cache_parameters)
+    function LinearInterpolation_pullback(ΔA)
+        df = NoTangent()
+        du = ΔA.u
+        dt = NoTangent()
+        dI = NoTangent()
+        dp = NoTangent()
+        dmode = NoTangent()
+        dextrapolate = NoTangent()
+        dcache_parameters = NoTangent()
+        df, du, dt, dI, dp, dmode, dextrapolate, dcache_parameters
+    end
+
+    A, LinearInterpolation_pullback
+end
+
+function u_tangent(A::LinearInterpolation, t, Δ)
+    out = zero(A.u)
+    idx = get_idx(A.t, t, A.idx_prev[])
+    t_factor = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx])
+    out[idx] = Δ * (one(eltype(out)) - t_factor)
+    out[idx + 1] = Δ * t_factor
+    out
+end
+
+function u_tangent(A::QuadraticInterpolation, t, Δ)
+    out = zero(A.u)
+    i₀, i₁, i₂ = _quad_interp_indices(A, t, A.idx_prev[])
+    t₀ = A.t[i₀]
+    t₁ = A.t[i₁]
+    t₂ = A.t[i₂]
+    Δt₀ = t₁ - t₀
+    Δt₁ = t₂ - t₁
+    Δt₂ = t₂ - t₀
+    out[i₀] = Δ * (t - A.t[i₁]) * (t - A.t[i₂]) / (Δt₀ * Δt₂)
+    out[i₁] = -Δ * (t - A.t[i₀]) * (t - A.t[i₂]) / (Δt₀ * Δt₁)
+    out[i₂] = Δ * (t - A.t[i₀]) * (t - A.t[i₁]) / (Δt₂ * Δt₁)
+    out
+end
+
+function u_tangent(A, t, Δ)
+    NoTangent()
+end
+
 function ChainRulesCore.rrule(::typeof(_interpolate),
         A::Union{
+            LinearInterpolation,
+            QuadraticInterpolation,
             LagrangeInterpolation,
             AkimaInterpolation,
             BSplineInterpolation,
             BSplineApprox
         },
         t::Number)
     deriv = derivative(A, t)
-    interpolate_pullback(Δ) = (NoTangent(), NoTangent(), deriv * Δ)
+    function interpolate_pullback(Δ)
+        (NoTangent(), Tangent{typeof(A)}(; u = u_tangent(A, t, Δ)), deriv * Δ)
+    end
     return _interpolate(A, t), interpolate_pullback
 end
 

ext/DataInterpolationsOptimExt.jl

@@ -18,9 +18,8 @@ function Curvefit(u,
         box = false,
         lb = nothing,
         ub = nothing;
-        extrapolate = false,
-        safetycopy = false)
-    u, t = munge_data(u, t, safetycopy)
+        extrapolate = false)
+    u, t = munge_data(u, t)
     errfun(t, u, p) = sum(abs2.(u .- model(t, p)))
     if box == false
         mfit = optimize(p -> errfun(t, u, p), p0, alg)