Any examples utilizing PEtab.jl for parameter fitting with MethodOfLines.jl?
If not, then could an ODESystem be (re)constructed from an ODEProblem or a PDESystem?
Use symbolic_discretize instead of discretize to get the ODESystem from MethodOfLines.jl and it should just compose.
and it does, thanks!!
with some further guessing I managed to move to almost the end but this last error I don’t really understand. It happens when executing p0 = generate_startguesses(petab_problem, 1)
this is an example, as simple as possible - the values and conditions are arbitrary:
using ModelingToolkit, MethodOfLines, Plots, DomainSets, OrdinaryDiffEq, PEtab, DataFrames
pythonplot()
# Parameters, variables, and derivatives
@parameters k a To
@variables t x T(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
Tamb = 273 #K
To_test = 1000 #K
h = 30 #W/m^2K arbritrary value
L = 0.1 #m
ρ = 7500 #kg/m^3
cp = 420 #J/kgK
kin = 15 #W/m K
ain = kin/ρ/cp #m^2/s
param = [k => kin, a => ain, To => To_test]
# 1D PDE and boundary conditions
eq = Dt(T(t, x)) ~ a * Dxx(T(t, x))
bcs = [T(0, x) ~ Tamb, #initial condition
Dx(T(t, 0)) ~ -h/k *(T(t, 0) - Tamb), #cold end
T(t, L) ~ To] #hot end
# Space and time domains
domains = [t ∈ Interval(0.0, 300.0),
x ∈ Interval(0.0, L)]
# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [T(t, x)], param)
# Method of lines discretization
dx = L/100
order = 2
discretization = MOLFiniteDifference([x => dx], t)
# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)
# Solve ODE problem
sol = solve(prob, Tsit5(), saveat=1.)
# Plot the solution
px = collect(sol[x])
pt = sol[t]
pT = sol[T(t,x)]
plt = plot(sol[t], sol[T(t,x)][:, 1], xlabel="Time", ylabel="Temperature", title="1D Heat Equation")
cnt = contourf(px, pt, pT, xlabel="Position", ylabel="Time", title="1D Heat Equation")
display(cnt)
#extract simulated experiments to create measurements for PEtab
idx =[3, 50, 97]
mx = px[idx]
mT = pT[:, idx]
#Create PEtab problem
odesys = symbolic_discretize(pdesys, discretization)
simpsys = structural_simplify(odesys[1])
@unpack T, k, a, To = simpsys
obs_T3 = PEtabObservable(T[3], 0.5)
observables = Dict("obs_T3" => obs_T3)
_k = PEtabParameter(:k, lb=0.1, ub=100., scale=:lin)
_a = PEtabParameter(:a, lb=1E-7, ub=1E-4, scale=:log10)
params = [_k , _a]
E0 = Dict(:To => 1000.)
E1 = Dict(:To => 900.)
sim_cond = Dict("c0" => E0, "c1" => E1)
measurements = DataFrame(
simulation_id = ["c0", "c1"],
obs_id=["obs_T3", "obs_T3"],
time = [300., 300.],
measurement=[400., 420.])
petab_model = PEtabModel(simpsys, sim_cond, observables, measurements, params, verbose=false)
petab_problem = PEtabODEProblem(petab_model, verbose=false)
p0 = generate_startguesses(petab_problem, 1)
res = calibrate_model(petab_problem, p0, IpoptOptimiser(false),
options=IpoptOptions(max_iter = 1000))
println(res)
and this is the error:
ERROR: CanonicalIndexError: setindex! not defined for Symbolics.Arr{Num, 1}
Stacktrace:
[1] error_if_canonical_setindex(::IndexCartesian, A::Symbolics.Arr{Num, 1}, ::Int64)
@ Base .\abstractarray.jl:1403
[2] setindex!(A::Symbolics.Arr{Num, 1}, v::Float64, I::Int64)
@ Base .\abstractarray.jl:1392
[3] macro expansion
@ .\none:7 [inlined]
[4] macro expansion
@ C:\Users\kkakosim\.julia\packages\RuntimeGeneratedFunctions\Yo8zx\src\RuntimeGeneratedFunctions.jl:163 [inlined]
[5] macro expansion
@ .\none:0 [inlined]
[6] generated_callfunc
@ .\none:0 [inlined]
[7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…})
@ RuntimeGeneratedFunctions C:\Users\kkakosim\.julia\packages\RuntimeGeneratedFunctions\Yo8zx\src\RuntimeGeneratedFunctions.jl:150
[8] change_ode_parameters!(p_ode_problem::Vector{…}, u0::Vector{…}, θ::Vector{…}, θ_indices::PEtab.ParameterIndices, petab_model::PEtabModel{…})
@ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Common.jl:252
[9] compute_cost_solve_ODE(θ_dynamic::SubArray{…}, θ_sd::SubArray{…}, θ_observable::SubArray{…}, θ_non_dynamic::SubArray{…}, ode_problem::ODEProblem{…}, ode_solver::ODESolver, ss_solver::SteadyStateSolver{…}, petab_model::PEtabModel{…}, simulation_info::PEtab.SimulationInfo{…}, θ_indices::PEtab.ParameterIndices, measurement_info::PEtab.MeasurementsInfo{…}, parameter_info::PEtab.ParametersInfo, petab_ODE_cache::PEtab.PEtabODEProblemCache{…}, petab_ODESolver_cache::PEtab.PEtabODESolverCache; compute_cost::Bool, compute_hessian::Bool, compute_gradient_θ_dynamic::Bool, compute_residuals::Bool, exp_id_solve::Vector{…})
@ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Objective\Objective.jl:78
[10] compute_cost_solve_ODE
@ C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Objective\Objective.jl:36 [inlined]
[11] compute_cost(θ_est::Vector{…}, ode_problem::ODEProblem{…}, ode_solver::ODESolver, ss_solver::SteadyStateSolver{…}, petab_model::PEtabModel{…}, simulation_info::PEtab.SimulationInfo{…}, θ_indices::PEtab.ParameterIndices, measurement_info::PEtab.MeasurementsInfo{…}, parameter_info::PEtab.ParametersInfo, prior_info::PEtab.PriorInfo, petab_ODE_cache::PEtab.PEtabODEProblemCache{…}, petab_ODESolver_cache::PEtab.PEtabODESolverCache, exp_id_solve::Vector{…}, compute_cost::Bool, compute_hessian::Bool, compute_residuals::Bool)
@ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Objective\Objective.jl:19
[12] #451
@ C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\PEtabODEProblem\Create.jl:358 [inlined]
[13] generate_startguesses(petab_problem::PEtabODEProblem{…}, n_multistarts::Int64; sampling_method::QuasiMonteCarlo.LatinHypercubeSample, sample_from_prior::Bool, allow_inf_for_startguess::Bool, verbose::Bool)
@ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Calibrate\Common.jl:288
[14] generate_startguesses(petab_problem::PEtabODEProblem{…}, n_multistarts::Int64)
@ PEtab C:\Users\kkakosim\.julia\packages\PEtab\xvijy\src\Calibrate\Common.jl:255
[15] top-level scope
@ p:\Chemical Engineering\Air_Quality_Eng_R_Team\!GITHUB\tamuq-chen-secarelab-receiver\aysha\MWE PEtab.jl:82
Some type information was truncated. Use `show(err)` to see complete types.
Looks like PETab.jl doesn’t support array variables. That’s an issue to open there.