Hello,
I am working with a small non-linear ODE problem to model some protein dynamics data. I am fitting my model to the data and estimating the parameters using Optimization.jl. So far everything worked fine, my problem occurs when I need to retrieve a desired parameter from the optimisation solution.
I built my ODE model with @mtkmodel and defined the ODE problem with ODEProblem. Here is a simpler example of my problem:
@mtkmodel protein_dynamics begin
@parameters begin
x1_act
x1_inact
x2_act
x2_inact
end
@variables begin
x1(t)
x2(t)
end
@equations begin
D(x1) ~ x1_act - x1_inact * x1
D(x2) ~ x2_act - x2_inact * x2
end
end
@mtkbuild protein_model = protein_dynamics()
true_param = [protein_model.x1_act => 0.02, protein_model.x1_inact => 0.2, protein_model.x2_act => 0.04, protein_model.x2_inact => 0.4]
simul_data_prob = ODEProblem(protein_model, [protein_model.x1 => 0.5, protein_model.x2 => 1], (0.0, 200.0), true_param, jac = true)
simul_data_sol = solve(simul_data_prob, saveat=1, dtmax=0.5)
# Import the data as array
simul_data_x1 = simul_data_sol[protein_model.x1]
simul_data_x2 = simul_data_sol[protein_model.x2]
When I pass the array of parameters to ODEProblem, it automatically builds an MTKParameters object which shuffles the order of the parameters.
println(simul_data_prob.p)
> ModelingToolkit.MTKParameters{Tuple{Vector{Float64}}, Tuple{}, Tuple{}, Tuple{}, Tuple{}, Nothing, Nothing}(([0.4, 0.02, 0.2, 0.04],), (), (), (), (), nothing, nothing)
This order is kept during the optimisation process and the parameters in the solution array are also in this order.
# Define MSE loss function to estimate parameters
function loss(new_parameters)
# Update the ODE problem with the new set of parameters
new_ODE_prob = remake(simul_data_prob, p=new_parameters)
# Solve the new ODE problem
new_sol = solve(new_ODE_prob, saveat = 1, dtmax=0.5)
new_sol_array_x1 = new_sol[protein_model.x1]
new_sol_array_x2 = new_sol[protein_model.x2]
# Determine the loss by computing the sum of squares
if length(simul_data_x1) == length(new_sol_array_x1)
loss = sum(abs2, new_sol_array_x1 .- simul_data_x1) +
sum(abs2, new_sol_array_x2 .- simul_data_x2)
else
loss = 1000
end
return loss, new_sol
end
# Define optimization function and differentiation method
optim_f = Optimization.OptimizationFunction((x, p) -> loss(x), Optimization.AutoForwardDiff())
# Initial parameter guess
param_guess = [0.1, 0.1, 0.1, 0.1]
# Define the optimisation problem
optim_prob = Optimization.OptimizationProblem(optim_f, param_guess)
# Run the optimisation
results = Optimization.solve(optim_prob, PolyOpt(), maxiters = 100)
println(results.u)
> [0.39999999267413, 0.020000000029211443, 0.19999999991817688, 0.03999999927883234]
I could not figure out how to obtain the index mapping between the initial array and the MTKParameters indexing order although I tried using SymbolicIndexingInterface.jl or SciMLStructures.jl as explained in the FAQ of ModelingToolkit.jl.
Does anyone know how to handle that indexing?