Help with using ModellingToolkit.jl for evolutionary analysis

Specific Domains Modelling & Simulations
,
1

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?

1 Like
2

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 

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

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)
3

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.
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)

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)

?

4

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

5

@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?

6

You’d use NLsolve.jl for that.

7

@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.

8

okay thanks will look into that! Thanks!