Merge pull request #374 from JuliaDynamics/hw/bump_sciml

bump sciml versions

Author: hexaeder

NEWS.md

@@ -1,6 +1,6 @@
 # NetworkDynamics Release Notes
 
-## v0.10.18 Changelog
+## v1.0.0 Changelog
 
 - **ModelingToolkit v11 compatibility** ([#344](https://github.com/JuliaDynamics/NetworkDynamics.jl/pull/344)):
   - Raises compat from `ModelingToolkit` v10 to v11. MTK v11 split into `ModelingToolkitBase` (MIT-licensed) and `ModelingToolkit` (AGPL-licensed); NetworkDynamics.jl now only requires the MIT-licensed `ModelingToolkitBase` as a direct dependency, dropping the AGPL dependency.
@@ -15,6 +15,7 @@
   - The default initial state for fixpoint search now applies guesses and formulas automatically (`guess=true, apply_formulas=true`).
 - **`NWState` constructor** gains `guess`, `apply_formulas`, and `verbose` keyword arguments. With `default=true` (the default), values are filled in order: defaults/inits → `InitFormula`s → guesses (if `guess=true`) → `GuessFormula`s. `apply_formulas=true` by default; `guess=false` by default.
 - `doctor` check: added smoketest for the observable function (`obsf`) of each component.
+- `SciMLBase@3` and `OrdinaryDfifEq@7` compat. Most notably, upstream `DiffEq` changed its default `initializealg` to *check* rather than *reinit* inconsistent initial conditions. To preserve the previous behavior, the `ODEProblem(nw::Network, ...)` constructor now defaults to `initializealg=BrownFullBasicInit()` (deliberately differing from the OrdinaryDiffEq default). Pass `initializealg=...` to the constructor to override. Also, the interface of the `VectorContinuousComponentCallback` changed to fit that of `VectorContinuousCallback`: the affect now receives a buffer with per-element up or down crossing information rather than a single event idx.
 
 ## v0.10.17 Changelog
 - **Open-loop linearization** ([#341](https://github.com/JuliaDynamics/NetworkDynamics.jl/pull/341)):

NetworkDynamicsInspector/Project.toml

@@ -36,8 +36,8 @@ GraphMakie = "0.6"
 Graphs = "1.12.0"
 NetworkDynamics = "1, 0.10"
 Observables = "0.5.5"
-OrderedCollections = "1.8.0"
+OrderedCollections = "1.8.0, 2"
 PrecompileTools = "1.2.1"
-SciMLBase = "2.75.1"
+SciMLBase = "3"
 WGLMakie = "0.13"
 julia = "1.10"

NetworkDynamicsInspector/test/Project.toml

@@ -28,7 +28,7 @@ Graphs = "≥0.0.1"
 ModelingToolkitBase = "≥0.0.1"
 NetworkDynamics = "≥0.0.1"
 NetworkDynamicsInspector = "≥0.0.1"
-OrdinaryDiffEqTsit5 = "1.1.0"
+OrdinaryDiffEqTsit5 = "1.1.0, 2"
 SafeTestsets = "0.1.0"
 SciCompDSL = "≥0.0.1"
 Test = "1"

Project.toml

@@ -12,6 +12,7 @@ Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -67,9 +68,10 @@ Accessors = "0.1.43"
 Adapt = "4.0.4"
 ArgCheck = "2.3.0"
 Atomix = "1"
-CUDA = "5.5.2"
+CUDA = "5.5.2, 6"
 ConstructionBase = "1.6.0"
 DataFrames = "1"
+DiffEqBase = "7"
 DiffEqCallbacks = "4.2.2"
 DifferentiationInterface = "0.7"
 DocStringExtensions = "0.9.3"
@@ -80,29 +82,29 @@ Hungarian = "0.7.0"
 InteractiveUtils = "1"
 KernelAbstractions = "0.9.18"
 LinearAlgebra = "1"
-LinearSolve = "3.46.1"
+LinearSolve = "3,4"
 MacroTools = "0.5.15"
 Mixers = "0.1.2"
 ModelingToolkitBase = "1.21"
 Moshi = "0.3"
 NNlib = "0.9.13"
 NonlinearSolve = "4"
-OrderedCollections = "1.8.0"
+OrderedCollections = "1.8.0, 2"
 Polyester = "0.7.12"
-PreallocationTools = "0.4.23, 1"
+PreallocationTools = "1"
 PrecompileTools = "1.2.1"
 Printf = "1.10.0"
 Random = "1"
-RecursiveArrayTools = "3.27.0"
+RecursiveArrayTools = "4"
 RuntimeGeneratedFunctions = "0.5.15"
-SciMLBase = "2"
+SciMLBase = "3"
 SparseArrays = "1"
 SparseConnectivityTracer = "0.6, 1"
 Static = "1.1.1"
 StaticArrays = "1.9.4"
 SteadyStateDiffEq = "2.2.0"
 StyledStrings = "1.0.3"
-SymbolicIndexingInterface = "0.3.27"
+SymbolicIndexingInterface = "0.3.36"
 SymbolicUtils = "4"
 Symbolics = "7"
 TimerOutputs = "0.5.23"

benchmark/Project.toml

@@ -31,7 +31,7 @@ Chairmarks = "1.3.1"
 Graphs = "≥0.0.1"
 Makie = "0.24"
 NetworkDynamics = "≥0.0.1"
-OrderedCollections = "1.7.0"
+OrderedCollections = "1.7.0, 2"
 PkgBenchmark = "0.2.12"
 PrettyTables = "3"
 ProgressMeter = "1.10.2"

docs/Project.toml

@@ -45,7 +45,7 @@ Bonito = "≥0.0.1"
 CairoMakie = "0.15"
 Chairmarks = "1.3.1"
 DataFrames = "≥0.0.1"
-DataInterpolations = "7, 8"
+DataInterpolations = "9"
 DiffEqCallbacks = "4.2.2"
 Distributions = "0.25.117"
 Documenter = "1.8.0"
@@ -60,10 +60,10 @@ ModelingToolkitBase = "≥0.0.1"
 NetworkDynamics = "≥0.0.1"
 NetworkDynamicsInspector = "≥0.0.1"
 NonlinearSolve = "≥0.0.1"
-OrdinaryDiffEqNonlinearSolve = "1.3.0"
-OrdinaryDiffEqRosenbrock = "1.4.0"
-OrdinaryDiffEqSDIRK = "1"
-OrdinaryDiffEqTsit5 = "1.1.0"
+OrdinaryDiffEqNonlinearSolve = "2"
+OrdinaryDiffEqRosenbrock = "2"
+OrdinaryDiffEqSDIRK = "2"
+OrdinaryDiffEqTsit5 = "2"
 Plots = "1.40.9"
 Printf = "≥0.0.1"
 Random = "≥0.0.1"

docs/examples/getting_started_with_network_dynamics.jl

@@ -224,7 +224,7 @@ nd_2 = Network(g, nd_diffusion_vertex_2, nd_diffusion_edge_2)
 We define the first diffusion as x ~ $N(0,1)^\mathrm{2}$ and the second diffusion as ϕ ~ $N(0,1)$. So the propagation
 of the diffusion from node 0 to node 2 is given by creating a vector:
 =#
-x0_2 = vec(transpose([randn(rng, N) .^ 2 randn(rng, N)]))
+x0_2 = collect(vec(transpose([randn(rng, N) .^ 2 randn(rng, N)])))
 
 # We then define the `ODEProblem`:
 ode_prob_2 = ODEProblem(nd_2, x0_2, (0.0, 3.0))
@@ -290,7 +290,7 @@ g = barabasi_albert(N, k)
 nd_diffusion_vertex_2 = VertexModel(; f=diffusionvertex_f!, g=1:2, dim=2, sym=[:x, :ϕ])
 nd_diffusion_edge_2 = EdgeModel(; g=AntiSymmetric(diffusionedge_g!), outsym=[:flow_x, :flow_ϕ])
 nd_2 = Network(g, nd_diffusion_vertex_2, nd_diffusion_edge_2)
-x0_2 = vec(transpose([randn(rng, N) .^ 2 randn(rng, N)]))
+x0_2 = collect(vec(transpose([randn(rng, N) .^ 2 randn(rng, N)])))
 ode_prob_2 = ODEProblem(nd_2, x0_2, (0.0, 3.0))
 Main.test_execution_styles(ode_prob_2) #src
 sol_2 = solve(ode_prob_2, Tsit5());

docs/src/callbacks.md

@@ -28,7 +28,7 @@ The main entry points are the types [`ContinuousComponentCallback`](@ref),
 [`VectorContinuousComponentCallback`](@ref) and [`DiscreteComponentCallback`](@ref). All of those objects combine a [`ComponentCondition`](@ref) with an [`ComponentAffect`](@ref).
 
 The "normal" [`ContinuousComponentCallback`](@ref) and [`DiscreteComponentCallback`](@ref) have a condition which returns a single value. The corresponding affect is triggered when the return value hits zero.
-In contrast, the "vector" version has an in-place condition which writes `len` outputs. When any of those outputs hits zero, the affect is triggered with an additional argument `event_idx` which tells the effect which dimension encountered the zerocrossing.
+In contrast, the "vector" version has an in-place condition which writes `len` outputs. When any of those outputs hits zero, the affect is triggered once with an additional argument `event_signs`: a length-`len` vector where entry `i` is `0` (no crossing), `+1` (upcrossing) or `-1` (downcrossing) for the `i`-th output. The affect resolves the crossing direction (and any simultaneous crossings) itself. This mirrors the `affect!(integrator, event_signs)` interface of the underlying [`VectorContinuousCallback`](@extref SciMLBase.VectorContinuousCallback).
 
 There is a special type [`PresetTimeComponentCallback`](@ref) which has no explicit condition and triggers the affect at given times.
 This internally generates a [`PresetTimeCallback`](@extref DiffEqCallbacks.PresetTimeCallback) object from `DiffEqCallbacks.jl`.
@@ -64,11 +64,15 @@ affect = ComponentAffect([:u], [:p]) do u, p, ctx
    obs = NWState(ctx.integrator)[VIndex(ctx.vidx, :obs)] # extract some observed state from context
    println("Trigger affect at t=$t")
 end
-vectoraffect = ComponentAffect([:u], [:p]) do u, p, event_idx, ctx
-    if event_idx == 1
-        u[:u] = 0 # change state
-    else
-        u[:p] = 0 # change parameter
+vectoraffect = ComponentAffect([:u], [:p]) do u, p, event_signs, ctx
+    for i in eachindex(event_signs)
+        event_signs[i] == 0 && continue # skip outputs that did not cross
+        if i == 1
+            u[:u] = 0 # first output crossed: change state
+        else
+            u[:p] = 0 # second output crossed: change parameter
+        end
+        # event_signs[i] is +1 for an upcrossing and -1 for a downcrossing
     end
 end
 ```
@@ -78,7 +82,7 @@ However the affect gets passed a `ctx` "context" object, which is a named tuple
 Lastly we need to define the actual callback object using [`ContinuousComponentCallback`](@ref)/[`VectorContinuousComponentCallback`](@ref):
 ```julia
 ccb  = ContinuousComponentCallback(condition, affect; kwargs...)
-vccb = VectorContinuousComponentCallback(condition, affect; kwargs...)
+vccb = VectorContinuousComponentCallback(condition, affect, len; kwargs...)
 ```
 where the `kwargs` are passed to the underlying [`SciMLBase.VectorContinuousCallback`](@extref) to finetune the zerocrossing-detection.
 

docs/src/mtk_integration.md

@@ -397,41 +397,38 @@ function attach_limint_callback!(cf, namespace)
         end
     end
 
-    upcrossing_affect = ComponentAffect([], [satmax, satmin]) do u, p, eventidx, ctx
-        if eventidx == 1
+    # A single affect handles all three conditions. `event_signs[i]` tells us for each
+    # condition output `i` whether it crossed and in which direction: `0` no crossing,
+    # `+1` upcrossing, `-1` downcrossing.
+    affect = ComponentAffect([], [satmax, satmin]) do u, p, event_signs, ctx
+        if event_signs[1] > 0 # upper saturation upcrossing
             println("$namespace: /⎺ reached upper saturation at $(round(ctx.t, digits=4))s")
             p[satmax] = 1.0
             p[satmin] = 0.0
-        elseif eventidx == 2
+        end
+        if event_signs[2] > 0 # lower saturation upcrossing
             println("$namespace: \\_ reached lower saturation at $(round(ctx.t, digits=4))s")
             p[satmax] = 0.0
             p[satmin] = 1.0
-        elseif eventidx == 3
-            # upcrossing means, forcing went from negative to positive, i.e. we leave lower saturation
+        end
+        if event_signs[3] > 0
+            # forcing went from negative to positive, i.e. we leave lower saturation
             insatmin = !iszero(p[satmin])
             if insatmin
                 println("$namespace: _/ left lower saturation at $(round(ctx.t, digits=4))s")
                 p[satmin] = 0.0
             end
-        else
-            error("Unknown event index $eventidx")
-        end
-    end
-
-    downcrossing_affect = ComponentAffect([],[satmax]) do u, p, eventidx, ctx
-        if eventidx == 3 # downcrossing means nothing for saturation affects
-            # downcrossing means, forcing went from positive to negative, i.e. we leave upper saturation
+        elseif event_signs[3] < 0
+            # forcing went from positive to negative, i.e. we leave upper saturation
             insatmax = !iszero(p[satmax])
             if insatmax
                 println("$namespace: ⎺\\ left upper saturation at $(round(ctx.t, digits=4))s")
                 p[satmax] = 0.0
             end
-        else
-            error("Unknown event index $eventidx")
         end
     end
 
-    cb = VectorContinuousComponentCallback(condition, upcrossing_affect, 3; affect_neg! = downcrossing_affect)
+    cb = VectorContinuousComponentCallback(condition, affect, 3)
 
     # finally add callback to component
     NetworkDynamics.add_callback!(cf, cb)

ext/NetworkDynamicsCUDAExt.jl

@@ -155,9 +155,9 @@ _adapt_rgf(f) = f
 function _adapt_diffcache(to, c::DiffCache)
     du = adapt(to, c.du)
     dual_du = adapt(to, c.dual_du)
-    DiffCache(du, dual_du, c.any_du)
-    # N = length(c.dual_du) ÷ length(c.du) - 1
-    # DiffCache(du, N)
+    # PreallocationTools 1.x added the `warn_on_resize` field to `DiffCache`; pass all
+    # fields positionally to hit the struct constructor (not the `(u, siz, chunks)` one).
+    DiffCache(du, dual_du, c.any_du, c.warn_on_resize)
 end
 
 end