Fix CI errors from jcn variables being set with initial conditions

[agchesebro]

No description provided.

[agchesebro]

@helmutstrey this fixes the warning so we'll see if the test passes now. Locally it didn't fail originally just went through the solve with a warning, so unsure if this is the only fix needed.

[MasonProtter]

This failure looks like it might be from a downstream dependency on the MTK side. It looks like one of the differential equations got turned into a differential algebraic equation that it didn't know how to handle. Investigating now.

[agchesebro]

@MasonProtter so I changed one line in neural_masses.jl to fix the issue Helmut mentioned and now a GraphDynamics test is failing. This mass shouldn't be interacting with GraphDynamics at all though so I'm confused what happened.

[agchesebro]

Jinx lol ignore my comment

[MasonProtter]

I had hoped that capping the MTK compat would prevent things like this, but perhaps what happened is that mtk changed something and then a downstream package was changed in response to that change and then when that other package changed it broke something in our tests.

This with fully updated dependancies, this MWE now fails on the MTK backend (works fine with the GraphDynamics backend):

julia> let tspan=(0.0, 2.0)
           
           # First focus is on producing panels from Figure 1 of the PING network paper.

           # Setup parameters from the supplemental material
           μ_E = 1.5
           σ_E = 0.15
           μ_I = 0.8
           σ_I = 0.08

           # Define the PING network neuron numbers
           NE_driven = 2
           NE_other = 14
           NI_driven = 4
           N_total = NE_driven + NE_other + NI_driven

           # First, create the 20 driven excitatory neurons
           exci_driven = [PINGNeuronExci(name=Symbol("ED$i"), I_ext=rand(Normal(μ_I, σ_I))) for i in 1:NE_driven]
           exci_other  = [PINGNeuronExci(name=Symbol("EO$i")) for i in 1:NE_other]
           inhib       = [PINGNeuronInhib(name=Symbol("ID$i"), I_ext=rand(Normal(μ_I, σ_I))) for i in 1:NI_driven]

           # Create the network
           g = MetaDiGraph()
           add_blox!.(Ref(g), vcat(exci_driven, exci_other, inhib))

           # Extra parameters
           N=N_total
           g_II=0.2
           g_IE=0.6
           g_EI=0.6

           for i = 1:NE_driven+NE_other
               for j = NE_driven+NE_other+1:N_total
                   add_edge!(g, i, j, Dict(:weight => g_EI/N))
                   add_edge!(g, j, i, Dict(:weight => g_IE/N))
               end
           end

           for i = NE_driven+NE_other+1:N_total
               for j = NE_driven+NE_other+1:N_total
                   add_edge!(g, i, j, Dict(:weight => g_II/N))
               end
           end
           @named sys = system_from_graph(g)
           prob = ODEProblem(sys, [], tspan, [])
           solve(prob, Tsit5())
       end
┌ Warning: Initialization system is overdetermined. 20 equations for 0 unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/zfOUk/src/systems/diffeqs/abstractodesystem.jl:1291
ERROR: MethodError: no method matching (::ModelingToolkit.GetUpdatedMTKParameters{…})(::SciMLBase.NonlinearSolution{…})
The object of type `ModelingToolkit.GetUpdatedMTKParameters{SymbolicIndexingInterface.MultipleGetters{ContinuousTimeseries, Vector{Any}}, SymbolicIndexingInterface.ParameterHookWrapper{SymbolicIndexingInterface.MultipleSetters{Vector{Any}}, Vector{Any}}}` exists, but no method is defined for this combination of argument types when trying to treat it as a callable object.

Closest candidates are:
  (::ModelingToolkit.GetUpdatedMTKParameters)(::Any, ::Any)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/zfOUk/src/systems/problem_utils.jl:398

Stacktrace:
  [1] get_initial_values(prob::ODEProblem{…}, valp::OrdinaryDiffEqCore.ODEIntegrator{…}, f::Function, alg::SciMLBase.OverrideInit{…}, iip::Val{…}; nlsolve_alg::Nothing, abstol::Float64, reltol::Float64, kwargs::@Kwargs{})
    @ SciMLBase ~/.julia/packages/SciMLBase/NtgCQ/src/initialization.jl:213
  [2] _initialize_dae!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, prob::ODEProblem{…}, alg::SciMLBase.OverrideInit{…}, isinplace::Val{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/initialize_dae.jl:159
  [3] _initialize_dae!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, prob::ODEProblem{…}, alg::OrdinaryDiffEqCore.DefaultInit, x::Val{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/initialize_dae.jl:50
  [4] initialize_dae!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, initializealg::OrdinaryDiffEqCore.DefaultInit)
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/initialize_dae.jl:40
  [5] __init(prob::ODEProblem{…}, alg::Tsit5{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/solve.jl:514
  [6] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/solve.jl:11 [inlined]
  [7] #__solve#62
    @ ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/solve.jl:6 [inlined]
  [8] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/solve.jl:1 [inlined]
  [9] #solve_call#44
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:632 [inlined]
 [10] solve_call
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:589 [inlined]
 [11] #solve_up#53
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1120 [inlined]
 [12] solve_up
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1099 [inlined]
 [13] #solve#51
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1036 [inlined]
 [14] solve(prob::ODEProblem{…}, args::Tsit5{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1026
 [15] top-level scope
    @ REPL[4]:46
Some type information was truncated. Use `show(err)` to see complete types.

[MasonProtter]

oops, that was a badly named commit, I meant to call it "don't set PING jcn to zero".

This is related to the breaking changes that MTK made a while ago where they suddenly decided that overriding initial conditions isn't allowed. Something somewhere is now confused by the jcn=0 to the point where it causes the PING network system to be detected as an algebraic equation (I think?) and then some other bug I can't track down triggers causing the whole thing to blow up.

[agchesebro]

Gotta love the new MTK bugs. I think that's the last of the neural masses that has that notation. Interesting that it didn't fail the test, but that mass is an SDE so maybe it somehow evades the check?

[MasonProtter]

Ugh. I guess those =0s are important somehow for the discrete blox?

[agchesebro]

Did they just not receive inputs? Because then there's no default value for jcn. Maybe BloxConnector needs to detect all jcns and automatically set to zero if there's no rhs?

[MasonProtter]

yeah probably. that'll be kinda annoying to do though since it can only be done at the end of system construction once you know whether or not there's any inputs.

For now lets see what happens if I just leave the discrete blox as-is

[agchesebro]

What problem was this even supposed to fix for MTK to make initial conditions immutable?

[MasonProtter]

No idea.

[ChrisRackauckas]

This is related to the breaking changes that MTK made a while ago where they suddenly decided that overriding initial conditions isn't allowed. Something somewhere is now confused by the jcn=0 to the point where it causes the PING network system to be detected as an algebraic equation (I think?) and then some other bug I can't track down triggers causing the whole thing to blow up.

Wait, what's the issue?

[ChrisRackauckas]

Oh, you're just hitting SciML/SciMLBase.jl#866. It's a snag from moving the initialization system from being something for ODEs to being something for all systems, but there was an issue with the testing setup when moving/generalizing that bit of code as mentioned in SciML/Catalyst.jl#1127 (comment). Hopefully our CI is unclogged and it's patched within the hour and you shouldn't need any changes.

[ChrisRackauckas]

(works fine with the GraphDynamics backend)

Is it properly building an initialization system / checking if one is not required? Your system looks truly overdetermined, so it's a correctness issue if it's not hitting the initialize check.

[MasonProtter]

Wait, what's the issue?

SciML/ModelingToolkit.jl#3047

It used to be okay to set an initial value, and then change it with equations but that's now treated as separate equations so it causes it to be overdetermined. Change was considered a bugfix, but it's definitely behaviour we were relying on here (right or wrong).

Oh, you're just hitting SciML/SciMLBase.jl#866. It's a snag from moving the initialization system from being something for ODEs to being something for all systems, but there was an issue with the testing setup when moving/generalizing that bit of code as mentioned in SciML/Catalyst.jl#1127 (comment). Hopefully our CI is unclogged and it's patched within the hour and you shouldn't need any changes.

I think we're only hitting that because of the change to how overdetermined systems are handled, not sure though, that's my naive assumption from there being stuff about DAEs in the stacktrace and this isn't supposed to be a DAE

(works fine with the GraphDynamics backend)

Is it properly building an initialization system / checking if one is not required? Your system looks truly overdetermined, so it's a correctness issue if it's not hitting the initialize check.

The GraphDynamics version isn't hitting anything overdetermined since it handles junctions differently. That =0 value is actually the init argument to a reduction.

[MasonProtter]

Reading that old issue again, we should try and see if guess=0 helps with the issues encountered doing this to the discrete blox

[MasonProtter]

Update: guess=0 didn't help

[ChrisRackauckas]

It used to be okay to set an initial value, and then change it with equations

I'm not sure what you mean by this, MWE?

[MasonProtter]

MWE is in SciML/ModelingToolkit.jl#3047

In earlier MTK 0.9.x the linked example did not error

[harisorgn]

Update: guess=0 didn't help

I was expecting this to work, but to be honest we don't need an initial value or a guess for these input variables since we initialize add a jcn ~ 0 equation when we initialize the connection making.

[ChrisRackauckas]

In earlier MTK 0.9.x the linked example did not error

In the earlier version it wasn't properly validating the equations so there was a correctness bug that got through. Now there is a correctness bug only with GraphDynamics.jl which is why it's let through.

[MasonProtter]

Again, GraphDynamics.jl does something completely different here that's not analogous and is not incorrect (as far as I'm aware).

Rather, it uses that 0 as an init argument to a mapreduce over all connections.