Catalyst.jl icon indicating copy to clipboard operation
Catalyst.jl copied to clipboard

DSL does not work for parameters/species/variables with complicated interdependencies

Open TorkelE opened this issue 10 months ago • 2 comments

Sometimes parameters/species/variables depend on each other. In the DSL, we create a single block each, in order:

@parameters begin

end
@variables begin

end
@species begin

end

Here, something declared, which depends on a subsequent thing, will throw an error. Right now, this is pretty much only a problem for when default values depend on other quantities. I don't think there is any relevant metadata which includes other parameters/variables/species, but if that gets added, that could be a problem.

I think the only "natural" case is a species/variable with a parametric initial condition, which is fine. But someone wanting to implement a model where something like this happen is possible, so would be nice to fix some time.

Here is the full list of cases which generate errors.

# Parameters depending on variables.
rn = @reaction_network begin
    @parameters a = A
    @variables A
end
# Parameters depending on species.
rn = @reaction_network begin
    @parameters a = X
    @species X
end
# Variables depending on species.
rn = @reaction_network begin
    @variables V = X
    @species X
end

# Parameters depending on subsequent parameters.
rn = @reaction_network begin
    @parameters a = b b
end

# Variables depending on subsequent variables.
rn = @reaction_network begin
    @variables A = B B
end

# Species depending on subsequent species.
rn = @reaction_network begin
    @species X = Y Y
end

TorkelE avatar Apr 15 '24 14:04 TorkelE

Was just about to make an issue related to this, it's a particular problem when using remake to update the initial conditions of an ODEProblem(rs, ...; remove_conserved = true).

  • This is because the Γ conservation law constants are defined as parameters, which have symbolic dependencies on the initial conditions.
  • Perhaps an overload is necessary for the internal SciMLBase.remake function here:
function _updated_u0_p_symmap(prob, u0, ::Val{true}, p, ::Val{true})
    isu0dep = any(symbolic_type(v) !== NotSymbolic() for (_, v) in u0)
    ispdep = any(symbolic_type(v) !== NotSymbolic() for (_, v) in p)

    if !isu0dep && !ispdep
        return remake_buffer(prob, state_values(prob), u0),
        remake_buffer(prob, parameter_values(prob), p)
    end
    if !isu0dep
        u0 = remake_buffer(prob, state_values(prob), u0)
        return _updated_u0_p_symmap(prob, u0, Val(false), p, Val(true))
    end
    if !ispdep
        p = remake_buffer(prob, parameter_values(prob), p)
        return _updated_u0_p_symmap(prob, u0, Val(true), p, Val(false))
    end

    varmap = merge(u0, p)
    u0 = anydict(k => symbolic_type(v) === NotSymbolic() ? v : symbolic_evaluate(v, varmap)
    for (k, v) in u0)
    p = anydict(k => symbolic_type(v) === NotSymbolic() ? v : symbolic_evaluate(v, varmap)
    for (k, v) in p)
    return remake_buffer(prob, state_values(prob), u0),
    remake_buffer(prob, parameter_values(prob), p)
end
  • Here the main issue I've found are that the ispdep check doesn't find that Γ in the parameters is a symbolic type, because it is evaluated to numerical values in the ODEProblem. The only place where the symbolic expressions for Γ are exposed are in defaults field of the ODESystem.
  • This obstructs the ostensible allowance for interdependencies between u0 and p, seen with the varmap = merge(u0, p) line.
  • Currently my workaround is calling symbolic_evaluate prior to calling remake using the symbolic expressions from the defaults map:
function compute_Γ(tunable_var_value_map, constraint_syms, sys_defs::T) where T

 #* Merge the default system definitions with the new symbolic map
 merged_defs = merge(sys_defs, Dict(tunable_var_value_map))

 NTuple{4, eltype(tunable_var_value_map)}(constraint => symbolic_evaluate(merged_defs[constraint], merged_defs) for constraint in constraint_syms)
end

I can make a separate issue for this in this repo (not sure if it pertains most here, or ModelingToolkit.jl or SciMLBase.jl) if this is truly a problem and I'm not missing anything.

jonathanfischer97 avatar Apr 16 '24 19:04 jonathanfischer97

Sounds good, yes, feel free to make a separate issue as well.

Right now there are a couple of bugs related to remake circulating within MTK (but I expect these to be sorted out soon). Once these are fixed, I will start having a look at your issue.

TorkelE avatar Apr 16 '24 20:04 TorkelE