Merge pull request #2671 from hersle/offset_indices

Fix structural simplification of non-1-indexed variable arrays

Project.toml

@@ -36,6 +36,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"

src/structural_transformation/symbolics_tearing.jl

@@ -1,3 +1,5 @@
+using OffsetArrays: Origin
+
 # N.B. assumes `slist` and `dlist` are unique
 function substitution_graph(graph, slist, dlist, var_eq_matching)
     ns = length(slist)
@@ -585,7 +587,12 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
         Symbolics.shape(lhs) !== Symbolics.Unknown() || continue
         arg1 = arguments(lhs)[1]
         haskey(obs_arr_subs, arg1) && continue
-        obs_arr_subs[arg1] = [arg1[i] for i in eachindex(arg1)]
+        obs_arr_subs[arg1] = [arg1[i] for i in eachindex(arg1)] # e.g. p => [p[1], p[2]]
+        index_first = eachindex(arg1)[1]
+
+        # respect non-1-indexed arrays
+        # TODO: get rid of this hack together with the above hack, then remove OffsetArrays dependency
+        obs_arr_subs[arg1] = Origin(index_first)(obs_arr_subs[arg1])
     end
     for i in eachindex(neweqs)
         neweqs[i] = fast_substitute(neweqs[i], obs_arr_subs; operator = Symbolics.Operator)

test/odesystem.jl

@@ -1164,6 +1164,13 @@ for sys in [sys1, sys2]
     end
 end
 
+@testset "Non-1-indexed variable array (issue #2670)" begin
+    @variables x(t)[0:1] # 0-indexed variable array
+    @named sys = ODESystem([x[0] ~ 0.0, D(x[1]) ~ x[0]], t, [x], [])
+    @test_nowarn sys = structural_simplify(sys)
+    @test equations(sys) == [D(x[1]) ~ 0.0]
+end
+
 # Namespacing of array variables
 @variables x(t)[1:2]
 @named sys = ODESystem(Equation[], t)