ModelingToolkit Example throwing errors

[RohanRajagopal]

Greetings,

Thanks very much for your package. I am just trying to run some of the examples in the link. The [MTK_example] is giving an error unfortunately. (https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/NME-MTK/#Neural-mass-equation-MTK)

using Revise, ModelingToolkit, LinearAlgebra
using DifferentialEquations, Plots
using BifurcationKit
const BK = BifurcationKit

indexof(sym, syms) = findfirst(isequal(sym),syms)

@variables t E(t) x(t) u(t) SS0(t) SS1(t) 	# independent and dependent variables
@parameters U0 τ J E0 τD U0 τF τS α    		# parameters
D = Differential(t) 				# define an operator for the differentiation w.r.t. time

# define the model
@named NMmodel = ODESystem([SS0 ~ J * u * x * E + E0,
	SS1 ~ α * log(1 + exp(SS0 / α)),
	D(E) ~ (-E + SS1) / τ,
	D(x) ~ (1.0 - x) / τD - u * x * E,
	D(u) ~ (U0 - u) / τF +  U0 * (1 - u) * E],
	defaults = Dict(E => 0.238616, x => 0.982747, u => 0.367876,
	α => 1.5, τ => 0.013, J => 3.07, E0 => -2.0, τD => 0.200, U0 => 0.3, τF => 1.5, τS => 0.007))

# get the vector field and jacobian
odeprob = ODEProblem(structural_simplify(NMmodel), [], (0.,10.), [], jac = true)
odefun = odeprob.f
F = (u,p) -> odefun(u,p,0)
J = (u,p) -> odefun.jac(u,p,0)

id_E0 = indexof(E0, parameters(NMmodel))
par_tm = odeprob.p

# we collect the differentials together in a problem
prob = BifurcationProblem(F, odeprob.u0, par_tm, (@lens _[id_E0]); J = J,
    record_from_solution = (x, p) -> (E = x[1], x = x[2], u = x[3]))

I get the following error.

ERROR: MethodError: no method matching isless(::Float64, ::Vector{Float64})

Closest candidates are:
  isless(::AbstractFloat, ::ForwardDiff.Dual{Ty}) where Ty
   @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/dual.jl:145
  isless(::T, ::T) where T<:Union{Float16, Float32, Float64}
   @ Base float.jl:548
  isless(::AbstractFloat, ::AbstractFloat)
   @ Base operators.jl:177
  ...

Stacktrace:
 [1] <(x::Float64, y::Vector{Float64})
   @ Base ./operators.jl:352
 [2] <=(x::Float64, y::Vector{Float64})
   @ Base ./operators.jl:401
 [3] iterate(it::ContIterable{…}; _verbosity::Int64)
   @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:262
 [4] iterate(it::ContIterable{…})
   @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:250
 [5] continuation(it::ContIterable{…})
   @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:474
 [6] continuation(prob::BifurcationProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…}; linear_algo::Nothing, bothside::Bool, kwargs::@Kwargs{…})
   @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:549

I tried to modify the model, creating only one ''parameter'' but I end up with a new error, isless(::Nothing, ::Int64). Am I missing something in the implementation?

ERROR: MethodError: no method matching isless(::Nothing, ::Int64)

Closest candidates are:
  isless(::ForwardDiff.Dual{Tx}, ::Integer) where Tx
   @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/dual.jl:144
  isless(::Num, ::Real)
   @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/num.jl:131
  isless(::ForwardDiff.Dual{Tx}, ::Real) where Tx
   @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/dual.jl:144
  ...

Stacktrace:
  [1] <(x::Nothing, y::Int64)
    @ Base ./operators.jl:352
  [2] <=(x::Nothing, y::Int64)
    @ Base ./operators.jl:401
  [3] getindex(buf::ModelingToolkit.MTKParameters{…}, i::Nothing)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/parameter_buffer.jl:425
  [4] get
    @ ~/.julia/packages/Setfield/PdKfV/src/lens.jl:202 [inlined]
  [5] getparam(pb::BifurcationProblem{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Problems.jl:225
  [6] iterate(it::ContIterable{…}; _verbosity::Int64)
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:254
  [7] iterate(it::ContIterable{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:250
  [8] continuation(it::ContIterable{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:474
  [9] continuation(prob::BifurcationProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…}; linear_algo::Nothing, bothside::Bool, kwargs::@Kwargs{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:549

Thank you very much.

[rveltz]

Thank you for reporting. The docs were built with this configuration: https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/#Reproducibility

You can perhaps try this interface https://docs.sciml.ai/ModelingToolkit/stable/tutorials/bifurcation_diagram_computation/

[rveltz]

maybe @TorkelE knows

[TorkelE]

Since its recent update, ModelingToolkit requires the system independent variable to be given to ODESystems. I.e. you have to add t after adding the equations:

# define the model
@named NMmodel = ODESystem([SS0 ~ J * u * x * E + E0,
	SS1 ~ α * log(1 + exp(SS0 / α)),
	D(E) ~ (-E + SS1) / τ,
	D(x) ~ (1.0 - x) / τD - u * x * E,
	D(u) ~ (U0 - u) / τF +  U0 * (1 - u) * E],
        t,
	defaults = Dict(E => 0.238616, x => 0.982747, u => 0.367876,
	α => 1.5, τ => 0.013, J => 3.07, E0 => -2.0, τD => 0.200, U0 => 0.3, τF => 1.5, τS => 0.007))

Also, at an early stage J is defined as a parameter, and later overwritten as the Jacobian. The code works perfectly fine, but might be a bit inadvisable, as if someone goes back and modifies the system, or runs it again or something, J can no longer be used to create equations (I messed this up myself here, hence I am mentioning it).

[rveltz]

Also, at an early stage J is defined as a parameter, and later overwritten as the Jacobian.

Oh yes!

[rveltz]

I guess related to https://github.com/bifurcationkit/BifurcationKit.jl/issues/154

[rveltz]

Should we close this?

[rveltz]

Should we close this?