Symbolic Jacobian sparsity

[visr]

This doesn't work yet, but I quickly tried this out yesterday, thought I might as well put it up to look if we should pursue this.
I think the main advantage is that we don't have to maintain our own code for determining the sparsity, which can more easily get out of sync. It does add some dependencies.

The error I get is with DataInterpolations not supporting ::SymbolicUtils.BasicSymbolic{Real}, which looks like this closed issue: SciML/DataInterpolations.jl#168. DataInterpolations has an extension on Symbolics for this purpose, so perhaps we just need to create a minimal example and file an issue.

ERROR: MethodError: no method matching (::DataInterpolations.LinearInterpolation{Vector{…}, Vector{…}, Float64})(::SymbolicUtils.BasicSymbolic{Real})

Closest candidates are:
  (::DataInterpolations.AbstractInterpolation)(::Symbolics.Num)
   @ DataInterpolationsSymbolicsExt C:\Users\visser_mn\.julia\packages\DataInterpolations\NoIUa\ext\DataInterpolationsSymbolicsExt.jl:15
  (::DataInterpolations.AbstractInterpolation)(::Number)
   @ DataInterpolations C:\Users\visser_mn\.julia\packages\DataInterpolations\NoIUa\src\DataInterpolations.jl:22
  (::DataInterpolations.AbstractInterpolation)(::Number, ::Integer)
   @ DataInterpolations C:\Users\visser_mn\.julia\packages\DataInterpolations\NoIUa\src\DataInterpolations.jl:23
  ...

Stacktrace:
  [1] get_area_and_level(basin::Ribasim.Basin{…}, state_idx::Int64, storage::Symbolics.Num)
    @ Ribasim d:\repo\ribasim\Ribasim\core\src\util.jl:57
  [2] set_current_basin_properties!(basin::Ribasim.Basin{…}, storage::SubArray{…})
    @ Ribasim d:\repo\ribasim\Ribasim\core\src\solve.jl:43
  [3] water_balance!(du::ComponentArrays.ComponentVector{…}, u::ComponentArrays.ComponentVector{…}, p::Ribasim.Parameters{…}, t::Float64)

It could be that there are a bunch more issues of this kind once this is resolved.

If it starts to work but we don't want it, we should perhaps use it to test our manual Jacobian sparsity.

[SouthEndMusic]

SciML/DataInterpolations.jl#310

[SouthEndMusic]

I did some digging. I found out these things:

• The problem you state above is likely not caused by DataInterpolations.jl, but by SmoothInterpolation.jl, because the problematic type is the output of a function call of that package and the problem goes away when I remove that call. That will be solved by Remove SmoothInterpolation as a dependency #1649
• Code that branches out based on sparsity-relevant parameters is incompatible, for instance in reduction_factor. I made a non-branching version of it
• PreallocationTools is not directly compatible with Symbolics, I found a workaround here: https://discourse.julialang.org/t/type-issue-to-get-preallocationtools-diffcache-work-with-symbolics-jacobian-sparsity/88542/3?u=bdekoning
• To trace trough the code to obtain the jacobian sparsity, the type parameter T in Parameters and its fields must (probably) be Vector{Symbolics.Num}, so the previous point might be irrelevant. I have not yet thought of a non-intrusive way to achieve this

[SouthEndMusic]

There seems to be a problem with DataInterpolations after all:

SciML/DataInterpolations.jl#317

[SouthEndMusic]

One thing to keep in mind here is that some basin dependencies are temporary and not 'active' at initialization. To get a proper sparsity pattern, all nodes should be active

[visr]

Left some comments. Let's have a look at the impact performance and binary size.

[visr]

Calculating the Jacobian seems fast, but runtime performance went down a lot. The wrap_forcing function leads to GC and runtime dispatch. One possible way to resolve this that we discussed is to switch the forcing ComponentVector to four separate Vectors. The issue is probably that the ComponentVector does not play well with similar_type used in LazyBufferCache, since the Axis of the components are part of the type.

[visr]

We should also have a look at https://github.com/adrhill/SparseConnectivityTracer.jl as an alternative, lighter dependency.

[SouthEndMusic]

SparseConnectivityTracer.jl needs some work: adrhill/SparseConnectivityTracer.jl#145

[SouthEndMusic]

SparseConnectivityTracer (the unsafe version) now works, only Aqua complains about type piracy on a custom overload, but that should be solved upstream

[visr]

Looks good to me! You'll have to approve yourself ;)

[SouthEndMusic]

https://nl.wikipedia.org/wiki/Wij_van_Wc-eend_adviseren_Wc-eend