Hi everyone, I am trying to use modelingtoolkitize to regenerate an ODEProblem with the analytical solution to the sparse Jacobian , as in this tutorial: Automated Sparse Analytical Jacobians · ModelingToolkit.jl
I am modelling gene regulation in a tissue of cells with possible diffusion of gene products. The gene regulatory network in each cell consists of three genes, which can inhibit or activate each others expression. The networks respond to a steady state morphogen gradient defined by m. Here is my code:
using ModelingToolkit, DifferentialEquations
const Nc = 100
const Ng = 3
const L = 1.
const θ = 5.
const c0 = 1.
const λm = 0.4
const tissue = range(0,L,length = Nc)
limit(a, N) = a == N+1 ? N : a == 0 ? 1 : a
m(x) = c0*exp(-x/λm)
σ(I) = 1/(1+exp(θ-θ*I))
# Method of lines: u_{j}(t) = u(x_j,t) where x_j = j*dx
function gene_regulation_1d(dg,g,p,t)
w, Dg, λg, dx = p
h = 1/dx^2
@inbounds for j in 1:Nc
x = tissue[j]
j_left, j_right = limit(j-1,Nc), limit(j+1,Nc)
dg[1,j] = σ(w[1,1]*g[1,j]+ w[1,2]*g[2,j] + w[1,3]*g[3,j] + w[1,4]*m(x)) + h*Dg[1]*(g[1,j_left] + g[1,j_right] - 2*g[1,j]) - λg[1]*g[1,j]
dg[2,j] = σ(w[2,1]*g[1,j]+ w[2,2]*g[2,j] + w[2,3]*g[3,j] + w[2,4]*m(x)) + h*Dg[2]*(g[2,j_left] + g[2,j_right] - 2*g[2,j]) - λg[2]*g[2,j]
dg[3,j] = σ(w[3,1]*g[1,j]+ w[3,2]*g[2,j] + w[3,3]*g[3,j] + w[3,4]*m(x)) + h*Dg[3]*(g[3,j_left] + g[3,j_right] - 2*g[3,j]) - λg[3]*g[3,j]
end
end
function init_gene_regulation_1d(start_conc)
g0 = zeros(Ng,Nc)
for j in 1:Nc
for i in 1:Ng
g0[i,j] = start_conc
end
end
return g0
end
const w_ex = [0.0 0.0 -0.0588219 2.08643;
0.388038 0.0 0.0 0.0;
-1.47171 0.0672493 0.0 0.0]
const Dg_ex = zeros(Ng)
const λg_ex = 0.05 .* ones(Ng)
Note I have set diffusion to zero in this example. I then create an ODE problem as per the tutorial. Since I am interested in finding the staeady state of the networks, I use an infinite time span with an appropriate call back function to terminate the steady state:
p = (w_ex,Dg_ex,λg_ex,step(tissue))
g0 = init_gene_regulation_1d(0.01);
tspan = (0,Inf)
prob_ode = ODEProblem(gene_regulation_1d,g0,tspan,p)
sol_ode = solve(prob_ode,AutoTsit5(Rosenbrock23()), isoutofdomain=(u,p,t) -> any(x -> x < 0, u), callback = TerminateSteadyState(1e-8,1e-6),maxiters = 1e5)
This works fine, and the solution is correctly terminated. I now want to improve the performance by generating the sparse jacobian. As per the tutorial, I use modeling toolkitize and generate the sparse problem:
sys = modelingtoolkitize(prob_ode);
sparseprob = ODEProblem(sys,Pair[],tspan,jac=true,sparse=true);
This works correctly. However now when I attempt to solve I get an error:
solve(sparseprob,AutoTsit5(Rosenbrock23()), isoutofdomain=(u,p,t) -> any(x -> x < 0, u), callback = TerminateSteadyState(1e-8,1e-6),maxiters = 1e5)
Error code:
Output exceeds the size limit. Open the full output data in a text editor.
MethodError: no method matching +(::Float64, ::Matrix{Float64}) For element-wise addition, use broadcasting with dot syntax: scalar .+ array
Can provide full output if useful. Any ideas what I am doing wrong? In addition, any pointers for performance upgrades in general for this model would be appreciated.