Plotting Catalyst systems no longer works with SciMLBase 2.11

[ranocha]

Describe the bug 🐞

I followed the instructions in the tutorial at https://docs.sciml.ai/Overview/stable/. Instead of getting the plot, I got an error.

Expected behavior

The code in the tutorial produces the same plot as shown on the website.

Minimal Reproducible Example 👇 Error & Stacktrace ⚠️

julia> using Plots, Catalyst, OrdinaryDiffEq

julia> rn = @reaction_network begin
           α, S + I --> 2I
           β, I --> R
       end
Model ##ReactionSystem#297
States (3):
  S(t)
  I(t)
  R(t)
Parameters (2):
  α
  β

julia> p     = [:α => .1/1000, :β => .01]
2-element Vector{Pair{Symbol, Float64}}:
 :α => 0.0001
 :β => 0.01

julia> tspan = (0.0,250.0)
(0.0, 250.0)

julia> u0    = [:S => 999.0, :I => 1.0, :R => 0.0]
3-element Vector{Pair{Symbol, Float64}}:
 :S => 999.0
 :I => 1.0
 :R => 0.0

julia> op    = ODEProblem(rn, u0, tspan, p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 250.0)
u0: 3-element Vector{Float64}:
 999.0
   1.0
   0.0

julia> sol   = solve(op, Tsit5())       # use Tsit5 ODE solver
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 25-element Vector{Float64}:
   0.0
   0.07042204340497453
   0.7746424774547197
   2.5720784698621917
   5.370475868328381
   8.892790675662736
  13.397314385835262
  18.77312249641433
  25.067606451778733
  32.21876025749868
  40.216010355336614
  49.06258942401493
  58.86717147955988
  69.82965358602912
  82.10838805404742
  95.46849615110189
 109.7545467360122
 122.01827800502546
 135.51407012931605
 149.4289751648907
 164.9581387700312
 182.1077023701446
 201.8902038451328
 225.42950830558135
 250.0
u: 25-element Vector{Vector{Float64}}:
 [999.0, 1.0, 0.0]
 [998.9929425461802, 1.006350999490624, 0.0007064543291988011]
 [998.9198582366232, 1.0721192430449888, 0.008022520331768075]
 [998.710958959454, 1.2601038170445216, 0.028937223501382614]
 [998.3106311383737, 1.6203391495804873, 0.06902971204575792]
 [997.6406610464425, 2.2231763292975812, 0.136162624259896]
 [996.4102628260856, 3.3301676362184667, 0.25956953769586744]
 [994.1215268974047, 5.388941226022392, 0.4895318765728288]
 [989.6067337013325, 9.44855116622177, 0.9447151324456879]
 [980.3144906488483, 17.797375650380932, 1.8881337007707701]
 [960.3226806148012, 35.72878549263647, 3.948533892562241]
 [916.1743022636258, 75.17093793996762, 8.654759796406582]
 [819.8113586901718, 160.42086871916172, 19.767772590666446]
 [633.230666598227, 321.17781962619586, 45.591513775577155]
 [374.89419038836917, 527.0968834174673, 98.00892619416355]
 [167.74733089130217, 653.8265939143656, 178.4260751943322]
 [64.78984789048314, 661.6683505558885, 273.5418015536284]
 [29.481085972467977, 618.2436612642134, 352.2752527633187]
 [13.346763344107229, 555.1339446439661, 431.5192920119268]
 [6.455508523840237, 489.39148441083967, 504.15300706532014]
 [3.184454885297364, 421.99665141068317, 574.8188937040195]
 [1.6352844539587053, 356.89846909139345, 641.466246454648]
 [0.8610792936956242, 293.53283900937345, 705.606081696931]
 [0.4650119511113042, 232.31441070419152, 767.2205773446973]
 [0.28026965071333754, 181.86671916839683, 817.8530111808899]

julia> plot(sol, lw=2)
ERROR: MethodError: no method matching cleansyms(::Vector{SymbolicUtils.BasicSymbolic{Real}})

Closest candidates are:
  cleansyms(::Tuple)
   @ SciMLBase ~/.julia/packages/SciMLBase/MMAmp/src/symbolic_utils.jl:74
  cleansyms(::LinearIndices)
   @ SciMLBase ~/.julia/packages/SciMLBase/MMAmp/src/symbolic_utils.jl:76
  cleansyms(::CartesianIndices)
   @ SciMLBase ~/.julia/packages/SciMLBase/MMAmp/src/symbolic_utils.jl:77
  ...

Stacktrace:
 [1] macro expansion
   @ ~/.julia/packages/SciMLBase/MMAmp/src/solutions/solution_interface.jl:180 [inlined]
 [2] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, sol::SciMLBase.AbstractTimeseriesSolution)
   @ SciMLBase ~/.julia/packages/RecipesBase/BRe07/src/RecipesBase.jl:300
 [3] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
   @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/user_recipe.jl:38
 [4] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
   @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/RecipesPipeline.jl:72
 [5] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
   @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:223
 [6] plot(args::Any; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})
   @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:102
 [7] top-level scope
   @ REPL[9]:1

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [479239e8] Catalyst v13.5.1
  [1dea7af3] OrdinaryDiffEq v6.63.0
  [91a5bcdd] Plots v1.39.0

[isaacsas]

Plotting still works fine on SciMLBase 2.10.

[ChrisRackauckas]

No don't overdramatize it, it's just a bug.

[isaacsas]

If it helps, I can also confirm the issues occurs if we directly use symbolic maps too, i.e.

...
rn = complete(rn)
p     = [rn.α => .1/1000, rn.β => .01]
u0    = [rn.S => 999.0, rn.I => 1.0, rn.R => 0.0]
...

[ChrisRackauckas]

I don't think that will make a difference. The real issue here is that plot recipes have some spotty coverage since plot functions are called during the tests but not plotting itself. That can be improved because that was setup due to some limitations in Julia v0.6 times, and I think we can run all of the plotting headless now. But reguardless, we need to fix the symbolic indexing in the plot recipe which seems to have regressed with the SymbolicIndexingInterface v0.3 update.

But there's a large issue in here which is really that the plot recipe has its own symbolic indexing feature which isn't the actual "symbolic indexing", so it's a bit janky right now. We need to throw out effectively 90% of the plot recipe and just make it use the now existent symbolic indexing features of SciML directly. If you look at the code, the vast majority of it is simply implementing a symbolic indexing feature because that all existed in the plot recipe about a few years before even ModelingToolkit existed, so it just doesn't work like the rest of the system. We should definitely get a hotfix in today ASAP however we can (@AayushSabharwal), but in the longer term the plot recipe needs an overhaul which would make it a lot more stable by using the features that actually exist in the library now,

i.e. idxs = rn.R should just use sol[rn.R] to get its values, but since that didn't exist it does not and does some weird side channel to do the same thing.

[AayushSabharwal]

I'll take a look. The plot recipes weren't really touched during the SII PR, just the obvious stuff like removing issymbollike

EDIT: Sheesh this is a mess

[isaacsas]

If plotting in tests now works via a headless mode, we can add that to Catalyst at least to catch such issues in the future.

[ChrisRackauckas]

Yes, we should just add some blank plotting to the tests. One thing at a time though, or accepting PRs 😅

[ChrisRackauckas]

EDIT: Sheesh this is a mess

Yes absolutely. In no way is the current plotting code any good, it's all technical debt and I'm pretty serious when I say you might be able to straight up delete 90% of the code just by using sol[x] 😅 . If you want to take this on then you're a hero.