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.
function addTransF(Fp, F, Tr, world)
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?