Help with using ModellingToolkit.jl for evolutionary analysis

Spawned from earlier question here: SymPy vs Sympy.jl and/or symengine.jl

Hello all. I am working on a converting a model written in python using sympy into julia. I posted a previous thread on packages I should use but the discussion strayed to actual troubleshooting and it was recommended I make a new thread for specific help.

As part of my workflow I am building a system of equation to describe a load of state based transitions. For example F is a matrix where each entry is a frequency of that type. Kind of like a compartment based SIR model. But the compartments are abstract organism states F[1,1] could be a those individuals with one unit of territory and no group members, and F[2,2] would be those individuals with two units of territory and 1 group member. In sympy I was then constructing large equations for how each of those frequencies evolves over time (mortality, birth, migration,various interactions, etc.). I was then taking that matrix Fâ€™ (which had the equations for how frequency changed) and creating a lambdified function that took in some matrix variables and returned the frequency changes. Then using scipy.optimize I looked for roots (using some biologically informed initial guesses).

Here is my actual initial code showing the addition of random state swaps due to environmental variables:

``````using ModelingToolkit

# possible states for an individual
world = Dict(
:q => 2,
:n => 3
)

# F the frequencies of each state
# W the natural transition probability to a new state
@variables begin
F[1:world[:q], 1:world[:n]]
W[1:world[:q], 1:world[:n]]
end

# F prime the change in frequency matrix
Fp = Array{Expression, 2}(undef, world[:q], world[:n])
for q in 1:world[:q]
for n in 1:world[:n]
Fp[q,n] = 0
end
end

# function to add transitions. I know this is clunky but other
# functions won't add simple things like this the logic code
# for transitions can be much more complicated.
for q in 1:world[:q]
for n in 1:world[:n]
for newQ in 1:world[:q]
Fp[q, n] = Fp[q, n] - Tr[q, newQ] * F[q, n]
Fp[q, n] = Fp[q, n] + Tr[newQ, q] * F[newQ, n]
end
end
end
return Fp
end

# output how frequency changes
Fp1 = addTransF(Fp, F, W, world)

# in sympy I would just generate function by:
# funF = build_function(Fp1, [F,W])
# that doesn't work
# instead you need to flatten and concat all elements or all arguments
funF = build_function(Fp1, vcat([x for x in reshape(F, (6,1))], [x for x in reshape(W, (6,1))]))
# or splat
funF2 = build_function(simplify.(Fp1), [F..., W...])

# making a callable fun
callFun = eval(funF2[1])

tW = [0 0.2 0.1; 0.4 0 0.2]

callFun([F..., tW...])
``````

The next bit would be solving for callFun==0 when W is fixed at some value. Havenâ€™t got that far because I wanted to check I hadnâ€™t done anything dumb.

Hopefully this make sense. Please let me know if there is a better way to assemble the model. The reason I havenâ€™t just substituted in values for W (which I canâ€™t seem to work out how to anyway) is that some values are evolving so each timestep I need to solve with different constants.

tldr: Is splatting the best way to pass matrix variables? How do you substitute in matrix variables to a matrix of expressions? How do you find a root for a multivariate function in Julia using ModelingToolkit or otherwise?

1 Like

Played around a bit more today not much progress. Pretty sure there must be an easier way to do this.

Not sure how to use IntervalRootFinding.jl properly and canâ€™t find any good examples of how to deal with variables and constants in the function you want to root.

``````using ModelingToolkit

@parameters g, l, d, epsilon, k, b, B, deltaX, deltaY, dummy, t

world = Dict(
:q => 2,
:n => 3
)

@variables begin
F[1:world[:q], 1:world[:n]]
W[1:world[:q], 1:world[:n]]
R[1:world[:q], 1:world[:n]]
M[1:world[:q], 1:world[:n]]
Mf[1:world[:q], 1:world[:n]]
Ml[1:world[:q], 1:world[:n]]
P[1:world[:q], 1:world[:n]]
Pf[1:world[:q], 1:world[:n]]
Pl[1:world[:q], 1:world[:n]]
C[1:world[:q], 1:world[:n]]
Cf[1:world[:q], 1:world[:n]]
Cl[1:world[:q], 1:world[:n]]
Tr[1:world[:q], 1:world[:n]]
X[1:world[:q], 1:world[:n]]
Y[1:world[:q], 1:world[:n]]
Xf[1:world[:q], 1:world[:n]]
Yf[1:world[:q], 1:world[:n]]
Xl[1:world[:q], 1:world[:n]]
Yl[1:world[:q], 1:world[:n]]
end

Fp = Array{Expression, 2}(undef, world[:q], world[:n])
for q in 1:world[:q]
for n in 1:world[:n]
Fp[q,n] = 0
end
end

for q in 1:world[:q]
for n in 1:world[:n]
for newQ in 1:world[:q]
Fp[q, n] = Fp[q, n] - Tr[q, newQ] * F[q, n]
Fp[q, n] = Fp[q, n] + Tr[newQ, q] * F[newQ, n]
end
end
end
return Fp
end

Fp1 = addTransF(Fp, F, W, world)

# funF = build_function(Fp1, [F,W])
funF = build_function(Fp1, vcat([x for x in reshape(F, (6,1))], [x for x in reshape(W, (6,1))]))
funF2 = build_function(Fp1, [F..., W...])

callFun = eval(funF[1])
callFun2 = eval(funF2[1])

tF = [0 0.3 0.2; 0 0.2 0.3]
tW = [0 0.2 0; 0.4 0 0]

callFun([tF..., tW...])

using StaticArrays
using IntervalArithmetic
using IntervalRootFinding

function functionF(tup)
vec = [i for i in tup]
arr = reshape(vec, (2,3))
ans = callFun([tF..., [0 0.2 0; 0.4 0 0]...])
return(SVector(SMatrix{1, 6}(ans)))
end

X = 0..1

rts = roots(functionF, XĂ—XĂ—XĂ—XĂ—XĂ—X)
``````

Hey, sorry for the late response. Are you just trying to make a function of two inputs like:

``````using ModelingToolkit

# possible states for an individual
world = Dict(
:q => 2,
:n => 3
)

# F the frequencies of each state
# W the natural transition probability to a new state
@variables begin
F[1:world[:q], 1:world[:n]]
W[1:world[:q], 1:world[:n]]
end

# F prime the change in frequency matrix
Fp = Array{Expression, 2}(undef, world[:q], world[:n])
for q in 1:world[:q]
for n in 1:world[:n]
Fp[q,n] = 0
end
end

# function to add transitions. I know this is clunky but other
# functions won't add simple things like this the logic code
# for transitions can be much more complicated.
for q in 1:world[:q]
for n in 1:world[:n]
for newQ in 1:world[:q]
Fp[q, n] = Fp[q, n] - Tr[q, newQ] * F[q, n]
Fp[q, n] = Fp[q, n] + Tr[newQ, q] * F[newQ, n]
end
end
end
return Fp
end

# output how frequency changes
Fp1 = addTransF(Fp, F, W, world)

funF = build_function(Fp1, F, W)
callFun = eval(funF[1])

Fv = [0 0.2 0.1; 0.4 0 0.2]
tW = [0 0.2 0.1; 0.4 0 0.2]

callFun(Fv,tW)
``````

?

Can you give a simple example of what youâ€™re trying to do using normal mathematical notation?

@ChrisRackauckas Thanks for the reply sorry I only just had time to come back to this problem. Maybe I was on an older version or just didnâ€™t read the docs properly. But your code as is does exactly what I wanted.

Is there then a similar function to `scipy.optimize.root` to get the value of `Fv` that zeros the `callFun`?

Youâ€™d use NLsolve.jl for that.

@dpsanders Hey Iâ€™m really sorry I am not really familiar with proper notation so not sure what exactly you mean. From a programming point of view now that I have a function G that accepts a matrix tF and returns the change in frequency. I want to find zeros of that function (stationary points).

I found one of your previous posts on IntervalRootFinding and have been trying to use that but not to much success.

okay thanks will look into that! Thanks!