Remove SmoothInterpolation as a dependency (#1649)

Fixes #1625. Some upstream
issues came out of this:

SciML/DataInterpolations.jl#301
SciML/DataInterpolations.jl#302
SciML/DataInterpolations.jl#303

---------

Co-authored-by: Martijn Visser <[email protected]>

Project.toml

@@ -41,11 +41,11 @@ OteraEngine = "b2d7f28f-acd6-4007-8b26-bc27716e5513"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
+ReadOnlyArrays = "988b38a3-91fc-5605-94a2-ee2116b3bd83"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Ribasim = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635"
 SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-SmoothInterpolation = "4412c578-8aba-4b04-8c8c-6ae573ccb792"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

core/Project.toml

@@ -30,9 +30,9 @@ LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
+ReadOnlyArrays = "988b38a3-91fc-5605-94a2-ee2116b3bd83"
 SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-SmoothInterpolation = "4412c578-8aba-4b04-8c8c-6ae573ccb792"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
@@ -51,7 +51,7 @@ ComponentArrays = "0.13, 0.14, 0.15"
 Configurations = "0.17"
 DBInterface = "2.4"
 DataFrames = "1.4"
-DataInterpolations = "=5.2.0"
+DataInterpolations = "=5.3.1"
 DataStructures = "0.18"
 Dates = "<0.0.1, 1"
 DiffEqCallbacks = "3.6"
@@ -71,9 +71,9 @@ MetaGraphsNext = "0.6, 0.7"
 OrdinaryDiffEq = "6.7"
 PreallocationTools = "0.4"
 ReTestItems = "1.20"
+ReadOnlyArrays = "0.2"
 SQLite = "1.5.1"
 SciMLBase = "2.36"
-SmoothInterpolation = "0.1.0"
 SparseArrays = "<0.0.1, 1"
 StructArrays = "0.6.13"
 TOML = "<0.0.1, 1"

core/src/Ribasim.jl

@@ -25,7 +25,8 @@ using Accessors: @set
 using Arrow: Arrow, Table
 using CodecZstd: ZstdCompressor
 using ComponentArrays: ComponentVector
-using DataInterpolations: LinearInterpolation, derivative, integral
+using DataInterpolations:
+    LinearInterpolation, LinearInterpolationIntInv, invert_integral, derivative, integral
 using Dates: Dates, DateTime, Millisecond, @dateformat_str
 using DBInterface: execute
 using DiffEqCallbacks:
@@ -50,6 +51,7 @@ using MetaGraphsNext:
     inneighbor_labels
 using OrdinaryDiffEq: OrdinaryDiffEq, OrdinaryDiffEqRosenbrockAdaptiveAlgorithm, get_du
 using PreallocationTools: DiffCache, get_tmp
+using ReadOnlyArrays: ReadOnlyVector
 using SciMLBase:
     init,
     solve!,
@@ -64,7 +66,6 @@ using SciMLBase:
     ODESolution,
     VectorContinuousCallback,
     get_proposed_dt
-using SmoothInterpolation
 using SparseArrays: SparseMatrixCSC, spzeros
 using SQLite: SQLite, DB, Query, esc_id
 using StructArrays: StructVector

core/src/parameter.jl

@@ -73,7 +73,13 @@ end
 
 Base.to_index(id::NodeID) = Int(id.value)
 
-const ScalarInterpolation = LinearInterpolation{Vector{Float64}, Vector{Float64}, Float64}
+const ScalarInterpolation = LinearInterpolation{
+    ReadOnlyVector{Float64, Vector{Float64}},
+    ReadOnlyVector{Float64, Vector{Float64}},
+    Vector{Float64},
+    Vector{Float64},
+    Float64,
+}
 
 """
 Store information for a subnetwork used for allocation.
@@ -261,7 +267,12 @@ end
     current_area::T = zeros(length(node_id))
     # Discrete values for interpolation
     storage_to_level::Vector{
-        LinearInterpolationIntInv{Vector{Float64}, Vector{Float64}, Float64},
+        LinearInterpolationIntInv{
+            ReadOnlyVector{Float64, Vector{Float64}},
+            Vector{Float64},
+            ScalarInterpolation,
+            Float64,
+        },
     }
     level_to_area::Vector{ScalarInterpolation}
     # Demands for allocation if applicable

core/src/read.jl

@@ -580,8 +580,7 @@ function Basin(db::DB, config::Config, graph::MetaGraph, chunk_sizes::Vector{Int
         error("Invalid Basin / profile table.")
     end
 
-    level_to_area = SmoothedLinearInterpolation.(area, level; extrapolate = true, λ = 0.01)
-    level_to_area = LinearInterpolation.(level_to_area)
+    level_to_area = LinearInterpolation.(area, level; extrapolate = true)
     storage_to_level = invert_integral.(level_to_area)
 
     return Basin(;
@@ -901,7 +900,12 @@ function user_demand_static!(
         for row in group
             priority_idx = findsorted(priorities, row.priority)
             demand_row = coalesce(row.demand, 0.0)
-            demand_itp[user_demand_idx][priority_idx].u .= demand_row
+            demand_itp_old = demand_itp[user_demand_idx][priority_idx]
+            demand_itp[user_demand_idx][priority_idx] = LinearInterpolation(
+                fill(demand_row, 2),
+                demand_itp_old.t;
+                extrapolate = true,
+            )
             demand[user_demand_idx, priority_idx] = demand_row
         end
     end
@@ -974,7 +978,9 @@ function UserDemand(db::DB, config::Config, graph::MetaGraph)::UserDemand
     demand_reduced = zeros(n_user, n_priority)
     trivial_timespan = [0.0, prevfloat(Inf)]
     demand_itp = [
-        [LinearInterpolation(zeros(2), trivial_timespan) for i in eachindex(priorities)] for j in eachindex(node_ids)
+        ScalarInterpolation[
+            LinearInterpolation(zeros(2), trivial_timespan) for i in eachindex(priorities)
+        ] for j in eachindex(node_ids)
     ]
     demand_from_timeseries = fill(false, n_user)
     allocated = fill(Inf, n_user, n_priority)

core/src/util.jl

@@ -104,8 +104,7 @@ function get_scalar_interpolation(
         push!(parameter, parameter[end])
     end
 
-    itp = SmoothedLinearInterpolation(parameter, times; extrapolate = true, λ = 0.1)
-    return LinearInterpolation(itp), allunique(times)
+    return LinearInterpolation(parameter, times; extrapolate = true), allunique(times)
 end
 
 """
@@ -121,8 +120,7 @@ function qh_interpolation(node_id::NodeID, table::StructVector)::ScalarInterpola
     pushfirst!(level, first(level) - 1)
     pushfirst!(flow_rate, first(flow_rate))
 
-    itp = SmoothedLinearInterpolation(flow_rate, level; extrapolate = true)
-    return LinearInterpolation(itp)
+    return LinearInterpolation(flow_rate, level; extrapolate = true)
 end
 
 """

core/test/allocation_test.jl

@@ -584,8 +584,8 @@ end
     (; allocation_models) = p.allocation
     (; basin, level_demand, graph) = p
 
-    fill!(level_demand.max_level[1].u, Inf)
-    fill!(level_demand.max_level[2].u, Inf)
+    fill!(level_demand.max_level[1].u.parent, Inf)
+    fill!(level_demand.max_level[2].u.parent, Inf)
 
     # Given a max_level of Inf, the basin capacity is 0.0 because it is not possible for the basin level to be > Inf
     @test Ribasim.get_basin_capacity(allocation_models[1], u, p, t, basin.node_id[1]) == 0.0

core/test/run_models_test.jl

@@ -296,8 +296,7 @@ end
 
 @testitem "Profile" begin
     import Tables
-    using DataInterpolations: LinearInterpolation, integral
-    using SmoothInterpolation: invert_integral
+    using DataInterpolations: LinearInterpolation, integral, invert_integral
 
     "Shorthand for Ribasim.get_area_and_level"
     function lookup(profile, S)
@@ -313,7 +312,7 @@ end
     storage = range(0.0, 1000.0, n_interpolations)
 
     # Covers interpolation for constant and non-constant area, extrapolation for constant area
-    A = [0.0, 100.0, 100.0]
+    A = [1e-9, 100.0, 100.0]
     h = [0.0, 10.0, 15.0]
     S = integral.(Ref(LinearInterpolation(A, h)), h)
     profile = (; S, A, h)
@@ -341,7 +340,7 @@ end
     end
 
     # Covers extrapolation for non-constant area
-    A = [0.0, 100.0]
+    A = [1e-9, 100.0]
     h = [0.0, 10.0]
     S = integral.(Ref(LinearInterpolation(A, h)), h)
 

core/test/utils_test.jl

@@ -12,8 +12,7 @@ end
 @testitem "bottom" begin
     using StructArrays: StructVector
     using Ribasim: NodeID
-    using DataInterpolations: LinearInterpolation, integral
-    using SmoothInterpolation: invert_integral
+    using DataInterpolations: LinearInterpolation, integral, invert_integral
 
     # create two basins with different bottoms/levels
     area = [[0.01, 1.0], [0.01, 1.0]]
@@ -46,8 +45,7 @@ end
     using StructArrays: StructVector
     using Logging
     using Ribasim: NodeID
-    using DataInterpolations: LinearInterpolation
-    using SmoothInterpolation: invert_integral
+    using DataInterpolations: LinearInterpolation, invert_integral
 
     level = [
         0.0,

core/test/validation_test.jl

@@ -384,7 +384,7 @@ end
     parameters = model.integrator.p
 
     (; graph, tabulated_rating_curve, basin) = parameters
-    tabulated_rating_curve.table[1].t[1] = invalid_level
+    tabulated_rating_curve.table[1].t.parent[1] = invalid_level
 
     logger = TestLogger()
     with_logger(logger) do
@@ -438,6 +438,9 @@ end
         Ribasim.main(toml_path)
     end
 
-    @test occursin("Warning: Convergence bottlenecks in descending order of severity:", output)
+    @test occursin(
+        "Warning: Convergence bottlenecks in descending order of severity:",
+        output,
+    )
     @test occursin("Basin #11 = ", output)
 end

pixi.toml

@@ -31,7 +31,7 @@ install = { depends_on = [
 # Julia
 update-registry-julia = "julia --eval='using Pkg; Registry.update()'"
 update-manifest-julia = "julia --project --eval='using Pkg; Pkg.update()'"
-instantiate-julia = "julia --project --eval='using Pkg; Pkg.instantiate()'"
+instantiate-julia = { cmd = "julia --project --eval='using Pkg; Pkg.instantiate()'", env = { JULIA_SSL_CA_ROOTS_PATH = "" } }
 initialize-julia = { depends_on = [
     "update-registry-julia",
     "instantiate-julia",