Comments (12)
JuMP turns parameters into tuples, so you probably don't want to use JuMP if the number of parameters gets large. That said, here's a working version of that code:
using ModelingToolkit
using OrdinaryDiffEq, DiffEqParamEstim, JuMP, NLopt
t_max = 10.0
dt = 0.1
t_out = 0:dt:t_max
tspan = (0,t_max)
@ModelingToolkit.parameters t, k1, k2
@ModelingToolkit.variables A(t) B(t)
@ModelingToolkit.derivatives D'~t
parameters = [k1 => 0.8, k2 => 0.6]
initial_conditions = [A => 1,
B => 0]
eqs = [D(A) ~ -k1*A + k2*B,
D(B) ~ k1*A - k2*B]
sys = ODESystem(eqs)
model_ode(p_) = OrdinaryDiffEq.ODEProblem(sys,initial_conditions,tspan,p_,jac=true)
solve_model(mp_) = OrdinaryDiffEq.solve(model_ode(mp_), Tsit5(),saveat=dt)
mock_data = Array(solve_model([0.8, 0.6]))
loss_objective(mp_, dat) = build_loss_objective(model_ode(collect(mp_)), Tsit5(), L2Loss(t_out,dat))
juobj(args...) = loss_objective(args, mock_data)(args)
jumodel = Model()
JuMP.register(jumodel, :juobj, 2, juobj, autodiff=true)
@JuMP.variable(jumodel, k1, start=1)
@JuMP.variable(jumodel, k2, start=1)
@JuMP.NLobjective(jumodel, Min, juobj(k1, k2))
set_optimizer(jumodel, NLopt.Optimizer)
set_optimizer_attribute(jumodel, "algorithm", :LD_MMA)
JuMP.optimize!(jumodel)
The other change was to make time floating point valued, since integer-valued time usually isn't a great idea. Hopefully that helps!
from sbml2julia.
build_loss_objective
isn't the issue here (and you could / probably should use DiffEqFlux.jl's more expressive form anyways). The issue is that is an ill-formed JuMP expression, so you couldn't write that with one experimental condition. As to why it says it's ill-formed, that's a JuMP question and you might want to ask a JuMP dev.
from sbml2julia.
That is fantastic!
from sbml2julia.
Thanks for the great suggestion @ChrisRackauckas
It would be nice to have an option to formulate the problem as a simulation-based problem using adjoint evaluation capability. And I think DifferentialEquations.jl would be the most convenient package for doing this 🙂 We can't promise anything at this point, but this seems like a great long-term plan.
Do you have an example code of hooking up DifferentialEquations.jl with NLP solvers via JuMP?
from sbml2julia.
https://diffeq.sciml.ai/stable/analysis/parameter_estimation/#Using-JuMP-with-DiffEqParamEstim is an example, but that will use forward-mode AD so it won't scale to huge problems as well. Mixing it with https://diffeq.sciml.ai/stable/analysis/sensitivity/#solve-Differentiation-Examples for the gradient definition is what you'd want.
from sbml2julia.
Thanks a lot for your interest, @ChrisRackauckas. I think your suggestion on Twitter to use a two step approach could make SBML2Julia much more powerful.
The current Julia code output of SBML2Julia looks like the tests/fixtures/jl_code_gold.jl
file. Do you think you would be able to contribute fixtures tests/fixtures/mtk_code_gold.jl
(or tests/fixtures/catalyst_code_gold.jl
), tests/fixtures/JuMP_code_gold.jl
, tests/fixtures/optim_code_gold.jl
, etc. which solve the optimization problem specified in tests/fixtures/jl_code_gold.jl
to the feature_mtk branch? You will probably see from the code structure that I know very little Julia, so please feel free to make improvements. I could then take it from there and modify the Python part to auto generate MTK, JuMP, optim etc. models.
from sbml2julia.
I am just trying to optimize a simple ODE system specified in MTK with JuMP:
using ModelingToolkit
using OrdinaryDiffEq, DiffEqParamEstim, JuMP, NLopt
t_max = 10
dt = 0.1
t_out = 0:dt:t_max
tspan = (0,t_max)
@ModelingToolkit.parameters t, k1, k2
@ModelingToolkit.variables A(t) B(t)
@ModelingToolkit.derivatives D'~t
parameters = [k1 => 0.8, k2 => 0.6]
initial_conditions = [A => 1,
B => 0]
eqs = [D(A) ~ -k1*A + k2*B,
D(B) ~ k1*A - k2*B]
sys = ODESystem(eqs)
model_ode(p_) = OrdinaryDiffEq.ODEProblem(sys,initial_conditions,tspan,p_,jac=true)
solve_model(mp_) = OrdinaryDiffEq.solve(model_ode(mp_), Tsit5(),saveat=dt)
mock_data = Array(solve_model([0.8, 0.6]))
loss_objective(mp_, dat) = build_loss_objective(model_ode(mp_), Tsit5(), L2Loss(t_out,dat))
juobj(args...) = loss_objective(args, mock_data)(args)
jumodel = Model()
JuMP.register(jumodel, :juobj, 2, juobj, autodiff=true)
@JuMP.variable(jumodel, k1, start=1)
@JuMP.variable(jumodel, k2, start=1)
@JuMP.NLobjective(jumodel, Min, juobj(k1, k2))
set_optimizer(jumodel, NLopt.Optimizer)
set_optimizer_attribute(jumodel, "algorithm", :LD_MMA)
JuMP.optimize!(jumodel)
Unfortunately, I get the following error:
ERROR: MethodError: no method matching lower_mapnames(::Tuple{Float64,Float64})
I am using Julia v1.5.3 and just added the used packages yesterday. Do you have ideas how to resolve this?
from sbml2julia.
Thanks! That helped a lot.
from sbml2julia.
I cannot find a way to get multiple experimental condition to work. What I tried for instance is:
mock_data_c1 = Array(solve_model([0.8, 0.6]))
mock_data_c2 = Array(solve_model([0.8, 0.5]))
loss_objective(mp_, dat) = build_loss_objective(model_ode(collect(mp_)), Tsit5(), L2Loss(t_out,dat))
juobj(args1, args2) = loss_objective(args1, mock_data1)(args1) + loss_objective(args2, mock_data2)(args2)
jumodel = Model()
JuMP.register(jumodel, :juobj, 2, juobj, autodiff=true)
@JuMP.variable(jumodel, k1_1, start=1)
@JuMP.variable(jumodel, k2_1, start=1)
@JuMP.variable(jumodel, k1_2, start=1)
@JuMP.variable(jumodel, k2_2, start=1)
@JuMP.NLobjective(jumodel, Min, juobj((k1_1, k2_1), (k1_2, k2_2)))
But I get the following error:
ERROR: Unexpected object (k1, k2) (of type Tuple{VariableRef,VariableRef} in nonlinear expression.
Stacktrace:
[1] error(::String) at ./error.jl:33
[2] _parse_NL_expr_runtime(::Model, ::Tuple{VariableRef,VariableRef}, ::Array{JuMP._Derivatives.NodeData,1}, ::Int64, ::Array{Float64,1}) at /home/paul/.julia/packages/JuMP/qhoVb/src/parse_nlp.jl:223
[3] top-level scope at /home/paul/.julia/packages/JuMP/qhoVb/src/parse_nlp.jl:115
[4] top-level scope at /home/paul/.julia/packages/JuMP/qhoVb/src/macros.jl:1368
[5] top-level scope at REPL[23]:1
If anyone knows how to solve this, please let me know.
from sbml2julia.
I would ask on #jump-dev as it's more of a JuMP question.
from sbml2julia.
Oh, I see. Thanks. I actually wanted to write an objective function that includes multiple experimental conditions and could be fed to any optimizer, not only JuMP. Is there really no way this can be done using build_loss_objective
?
from sbml2julia.
Little update from my side. I tried to do the first step of the 2-step approach and wrote SbmlInterface.jl, which converts SBML models to MTK and simulates them. Thanks a lot for your help on slack, @ChrisRackauckas . Next plan is to interface with optimizers.
from sbml2julia.
Related Issues (9)
- Restoration fails when using objectivePriors HOT 2
- Cannot increase weighting of likelihood in `examples/Shin_PLOS2019` HOT 6
- Reading measurementData into Julia HOT 2
- MUMPS requires more memory HOT 3
- Invalid number in NLP function or derivative detected HOT 2
- Cannot install linear solver HOT 2
- Failing Benchmark-Models-PEtab cases HOT 6
- Different time vectors for different exp. conditions HOT 7
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from sbml2julia.