sebapersson / PEtab.jl

Create parameter estimation problems for ODE models
https://sebapersson.github.io/PEtab.jl/stable/
MIT License
33 stars 5 forks source link

Error when declaring (noise) parameters with `σ` in name #156

Open TorkelE opened 7 months ago

TorkelE commented 7 months ago

I have a simple CRN which I wish to fit to data. There is noise in the observations. I want to incooperate this by adding some parameters to the observation formulas. In the example below I want to set noise which scales with the concentration of the species. Depending on the name of my parameter, I might get errors. does not work, ns does.

This works

using Catalyst, PEtab, DataFrames

# Model.
rs = @reaction_network begin
    (k1, k2), X1 <--> X2 
end

# Known initial conditions.
u0 = [:X1 => 1.0]

# Observables.
@parameters ns
@unpack X1 = rs
obs_X1 = PEtabObservable(X1, X1*ns)
observables = Dict("obs_X1" => obs_X1)

# Parameters.
par_k1 = PEtabParameter(:k1)
par_k2 = PEtabParameter(:k2)
par_noise = PEtabParameter(:ns)
params = [par_k1, par_k2, par_noise]

# Conditions.
c1 = Dict(:X2 => 1.0)
c2 = Dict(:X2 => 2.0)
simulation_conditions = Dict("c1" => c1, "c2" => c2)

# Measurements
using DataFrames
m_c1 = DataFrame(simulation_id = "c1", obs_id="obs_X1", time=[1.0, 2.0, 3.0], measurement=[1.1, 1.2, 1.3])
m_c2 = DataFrame(simulation_id = "c2", obs_id="obs_X1", time=[1.0, 2.0, 3.0], measurement=[1.2, 1.4, 1.6])
measurements = vcat(m_c1, m_c2)

# Make optimisation problem.
petab_model = PEtabModel(rs, simulation_conditions, observables, measurements, params; state_map=u0)

But this:

using Catalyst, PEtab, DataFrames

# Model.
rs = @reaction_network begin
    (k1, k2), X1 <--> X2 
end

# Known initial conditions.
u0 = [:X1 => 1.0]

# Observables.
@parameters nσ
@unpack X1 = rs
obs_X1 = PEtabObservable(X1, X1*nσ)
observables = Dict("obs_X1" => obs_X1)

# Parameters.
par_k1 = PEtabParameter(:k1)
par_k2 = PEtabParameter(:k2)
par_noise = PEtabParameter(:nσ)
params = [par_k1, par_k2, par_noise]

# Conditions.
c1 = Dict(:X2 => 1.0)
c2 = Dict(:X2 => 2.0)
simulation_conditions = Dict("c1" => c1, "c2" => c2)

# Measurements
using DataFrames
m_c1 = DataFrame(simulation_id = "c1", obs_id="obs_X1", time=[1.0, 2.0, 3.0], measurement=[1.1, 1.2, 1.3])
m_c2 = DataFrame(simulation_id = "c2", obs_id="obs_X1", time=[1.0, 2.0, 3.0], measurement=[1.2, 1.4, 1.6])
measurements = vcat(m_c1, m_c2)

# Make optimisation problem.
petab_model = PEtabModel(rs, simulation_conditions, observables, measurements, params; state_map=u0)

gives a

ERROR: StringIndexError: invalid index [24], valid nearby indices [23]=>'σ', [25]=>','
Stacktrace:
  [1] string_index_err(s::String, i::Int64)
    @ Base ./strings/string.jl:12
  [2] getindex
    @ Base ./strings/string.jl:470 [inlined]
  [3] create_top_function_h(model_state_names::Vector{…}, parameter_info::PEtab.ParametersInfo, p_ode_problem_names::Vector{…}, θ_non_dynamic_names::Vector{…})
    @ PEtab ~/.julia/packages/PEtab/fBBay/src/Process_PEtab_files/Observables/Create_u0_h_sigma.jl:191
  [4] create_h_function(model_name::String, dir_model::String, model_state_names::Vector{…}, parameter_info::PEtab.ParametersInfo, p_ode_problem_names::Vector{…}, θ_non_dynamic_names::Vector{…}, observables_data::CSV.File, model_SBML::SBMLImporter.ModelSBML, write_to_file::Bool)
    @ PEtab ~/.julia/packages/PEtab/fBBay/src/Process_PEtab_files/Observables/Create_u0_h_sigma.jl:111
  [5] create_σ_h_u0_file(model_name::String, system::ReactionSystem{…}, experimental_conditions::CSV.File, measurements_data::CSV.File, parameters_data::CSV.File, observables_data::CSV.File, state_map::Vector{…})
    @ PEtab ~/.julia/packages/PEtab/fBBay/src/Process_PEtab_files/Observables/Create_u0_h_sigma.jl:73
  [6] macro expansion
    @ ~/.julia/packages/PEtab/fBBay/src/Process_PEtab_files/Julia_tables_provided/Create_PEtab_model.jl:192 [inlined]
  [7] macro expansion
    @ ./timing.jl:395 [inlined]
  [8] _PEtabModel(system::ReactionSystem{…}, model_name::String, simulation_conditions::Dict{…}, observables::Dict{…}, measurements::DataFrame, petab_parameters::Vector{…}, state_map::Vector{…}, parameter_map::Nothing, events::Nothing, verbose::Bool)
    @ PEtab ~/.julia/packages/PEtab/fBBay/src/Process_PEtab_files/Julia_tables_provided/Create_PEtab_model.jl:191
  [9] #PEtabModel#493
    @ PEtab ~/.julia/packages/PEtab/fBBay/src/Process_PEtab_files/Julia_tables_provided/Create_PEtab_model.jl:130 [inlined]
 [10] top-level scope
    @ ~/Desktop/Research Projects/KineticNetworkIdentifiabilityLevels/playground.jl:123
Some type information was truncated. Use `show(err)` to see complete types.

error.

sigma in rates is fine. E.g. this:

rs = @reaction_network begin
    (k1, kσ), X1 <--> X2 
end

works.

sebapersson commented 7 months ago

Thanks for noticing this. This seems to be due to using String processing when internally building the compute_σ as σ seems to not be a valid Char.

At the moment use normal characters in noise formula like s or sigma, and hopefully I should patch this before next weeks end.