Catalyst integration cont.

docs/Project.toml

@@ -4,6 +4,6 @@ PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
 SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
 
 [compat]
-Documenter = "0.27"
+Documenter = "1"
 PEtab = "2"
 SBML = "1"

docs/make.jl

@@ -5,9 +5,10 @@ DocMeta.setdocmeta!(PEtab, :DocTestSetup, :(using PEtab); recursive=true)
 
 makedocs(;
     modules=[PEtab],
-    strict=false,
     authors="Viktor Hasselgren, Sebastian Persson, Damiano Ognissanti, Rafael Arutjunjan",
     repo="https://github.com/sebapersson/PEtab.jl/blob/{commit}{path}#{line}",
+    checkdocs=:exports,
+    warnonly=false,
     sitename="PEtab.jl",
     format=Documenter.HTML(;
         prettyurls=get(ENV, "CI", "false") == "true",

docs/src/API_choosen.md

@@ -1,9 +1,7 @@
 # API
 
-```@docs
+```@docs; canonical=true
 PEtabModel
-PEtabModel
-PEtabODEProblem
 PEtabODEProblem
 PEtabObservable
 PEtabParameter
@@ -17,4 +15,5 @@ calibrate_model
 calibrate_model_multistart
 run_PEtab_select
 solve_SBML
+generate_startguesses
 ```

docs/src/Define_in_julia.md

@@ -7,7 +7,7 @@ Here we demonstrate how to define a parameter estimation problem using a simple
 1. **Dynamic Model**: You can use either a `ReactionSystem` defined in [Catalyst](https://petab.readthedocs.io/en/latest/) or an `ODESystem` defined in [ModellingToolkit](https://github.com/SciML/ModelingToolkit.jl).
 2. **Observable Formula**: To link the model to the measurement data, you need an observable formula. Since real-world data often comes with measurement noise, you also must specify a noise formula and noise distribution. This is specified as a `PEtabObservable`.
 3. **Parameters to Estimate**: Typically, you do not want to estimate all model parameters. Moreover, sometimes you might want to incorporate prior beliefs by assigning priors to certain parameters. Parameter information is provided as a vector of `PEtabParameter`.
-4. **Simulation Conditions**: Measurements are often taken under various experimental conditions, such as different substrate concentrations. These experimental conditions typically correspond to model control parameters, like the initial value of a model species. You specify these conditions as a `Dict` (see below).
+4. **Simulation Conditions**: Measurements are often taken under various experimental conditions, such as different substrate concentrations. These experimental conditions typically correspond to model control parameters, like the initial value of a model species. You specify these conditions as a `Dict` (see below). In case the model only has a single simulation conditions, `simulation_conditions` can be omitted when building the `PEtabModel`.
 5. **Measurement Data**: To calibrate the model, you need measurement data, which should be provided as a `DataFrame`. The data format is explained below.
 
 ## Defining the Dynamic Model
@@ -187,6 +187,9 @@ petab_problem = PEtabODEProblem(petab_model, verbose=true)
 
 The `PEtabODEProblem` contains all the necessary information to work with most available optimizers (see [here](@ref import_petab_problem)). Alternatively, if you want to perform parameter estimation using a multi-start approach, you can use the `calibrate_model_multistart` function (see see [Parameter estimation](@ref parameter_estimation)).
 
+!!! note
+    If the model does not have multiple simulation conditions (e.g., data is collected under a single condition), you can omit the `simulation_conditions` argument when constructing the `PEtabModel` and the `simulation_id` columns from the measurement data. Simply use the following format: `PEtabModel(sys, observables, measurements, parameters, <keyword arguments>)`.
+
 ## Where to Go Next
 
 This example has covered the fundamental aspects of setting up a parameter estimation problem directly Julia, but there are additional options:

docs/src/Parameter_estimation.md

@@ -105,11 +105,12 @@ res.runs[1].ftrace
 
 If we want to perform single-start parameter estimation instead of multistart, we can use the `calibrate_model` function. This function runs a single optimization from a given initial guess.
 
-Given a starting point `p0`, and that we want to use Ipopt for optimization the model can be parameter estimated via:
+Given a starting point `p0` which can be generated by the `generate_startguesses` function, and that we want to use Ipopt for optimization the model can be parameter estimated via:
 
 ```julia
+p0 = generate_startguesses(petab_problem, 1)
 res = calibrate_model(petab_problem, p0, IpoptOptimiser(false),
-                     options=IpoptOptions(max_iter = 1000))
+                      options=IpoptOptions(max_iter = 1000))
 print(res)
 ```
 ```
@@ -122,4 +123,12 @@ Run time              = 1.9e+00s
 Optimiser algorithm   = Ipopt_user_Hessian
 ```
 
-The results are returned as a `PEtabOptimisationResult`, which includes the following information: minimum parameter values found (`xmin`), smallest objective value (`fmin`), number of iterations, runtime, whether the optimizer converged, and optionally, the trace if `save_trace=true`.
+The results are returned as a `PEtabOptimisationResult`, which includes the following information: minimum parameter values found (`xmin`), smallest objective value (`fmin`), number of iterations, runtime, whether the optimizer converged, and optionally, the trace if `save_trace=true`.
+
+Note that `generate_startguesses` can generate several start-guesses within the bounds of the `PEtabODEProblem` with any sampling method from [QuasiMonteCarlo](https://github.com/SciML/QuasiMonteCarlo.jl). For example, to generate 10 start guesses with Sobol sampling do:
+
+```julia
+using QuasiMonteCarlo
+p0 = generate_startguesses(petab_problem, 10, 
+                           sampling_method=SobolSample())
+```

examples/0002/PEtab_select_forward_AIC.yaml

@@ -1,9 +1,9 @@
 criteria:
-  AIC: -4.705325994421834
-  NLLH: -4.352662997210917
+  AIC: -4.705325993742942
+  NLLH: -4.352662996871471
 estimated_parameters:
-  k1: 0.2016087805174873
-  sigma_x2: 0.11714062782969972
+  k1: 0.20160878221098633
+  sigma_x2: 0.11714134130776435
 model_hash: 4f0ae55358f04ad6ab79c4dff4c6a969b83846171dd05a97dae872b3264a0eb3b2e6bddd0a71ac602786c75202b8ad033d9e42762d792c27d4a00dc3b2f4c8d2
 model_id: 4f0ae55358f04ad6ab79c4dff4c6a969b83846171dd05a97dae872b3264a0eb3b2e6bddd0a71ac602786c75202b8ad033d9e42762d792c27d4a00dc3b2f4c8d2
 model_subspace_id: M1_3

examples/Parameter_estimation.jl

@@ -6,8 +6,7 @@ dir_save = nothing
 
 using Optim
 res1 = calibrate_model_multistart(petab_problem, IPNewton(), 10, dir_save,
-                                  options=Optim.Options(iterations = 200,
-                                  show_trace=true))
+                                  options=Optim.Options(iterations = 200, show_trace=true))
 
 using PyCall
 using QuasiMonteCarlo
@@ -16,13 +15,13 @@ res2 = calibrate_model_multistart(petab_problem, Fides(nothing; verbose=true), 1
 
 using Ipopt
 res3 = calibrate_model_multistart(petab_problem, IpoptOptimiser(false), 10, dir_save,
-                               save_trace=true,
-                               seed=123)
+                                  save_trace=true,
+                                  seed=123)
 println("x-trace = ", res3.runs[1].xtrace)
 println("f-trace = ", res3.runs[1].ftrace)
 res3.runs[1].xtrace
 res3.runs[1].ftrace
 
-p0 = petab_problem.θ_nominalT .* 0.5
+p0 = generate_startguesses(petab_problem, 1)
 res = calibrate_model(petab_problem, p0, IpoptOptimiser(false),
-                     options=IpoptOptions(max_iter = 1000, print_level=5))
+                      options=IpoptOptions(max_iter = 1000, print_level=5))

ext/CatalystExtension/Create_PEtab_model.jl

@@ -11,6 +11,19 @@ function PEtab.PEtabModel(system::ReactionSystem,
     return PEtab._PEtabModel(system, model_name, simulation_conditions, observables, measurements, 
                              petab_parameters, state_map, parameter_map, verbose)
 end
+function PEtab.PEtabModel(system::ReactionSystem,
+                          observables::Dict{String, PEtab.PEtabObservable},
+                          measurements::DataFrame,
+                          petab_parameters::Vector{PEtab.PEtabParameter};
+                          state_map::Union{Nothing, Vector{Pair{T1, Float64}}}=nothing,
+                          parameter_map::Union{Nothing, Vector{Pair{T2, Float64}}}=nothing,
+                          verbose::Bool=false)::PEtab.PEtabModel where {T1<:Union{Symbol, Num}, T2<:Union{Symbol, Num}}
+
+    simulation_conditions = Dict("__c0__" => Dict())                        
+    model_name = "ReactionSystemModel"                          
+    return PEtab._PEtabModel(system, model_name, simulation_conditions, observables, measurements, 
+                             petab_parameters, state_map, parameter_map, verbose)
+end
 
 
 function PEtab.get_default_values(system::ReactionSystem)

src/Model_callibration/Common.jl

@@ -137,44 +137,59 @@ gradient computation. If left blank, we automatically select appropriate options
 function run_PEtab_select end
 
 
-function save_partial_results(path_save_res::String,
-                              path_save_parameters::String,
-                              path_save_trace::Union{String, Nothing},
-                              res::PEtabOptimisationResult,
-                              θ_names::Vector{Symbol},
-                              i::Int64)::Nothing
+"""
+    generate_startguesses(petab_problem::PEtabODEProblem,
+                          n_multistarts::Int64;
+                          sampling_method::T=QuasiMonteCarlo.LatinHypercubeSample(),
+                          allow_inf_for_startguess::Bool=false,
+                          verbose::Bool=false)::Array{Float64} where T <: QuasiMonteCarlo.SamplingAlgorithm
 
-    df_save_res = DataFrame(fmin=res.fmin,
-                          alg=String(res.alg),
-                          n_iterations = res.n_iterations,
-                          run_time = res.runtime,
-                          converged=string(res.converged),
-                          Start_guess=i)
-    df_save_parameters = DataFrame(Matrix(res.xmin'), θ_names)
-    df_save_parameters[!, "Start_guess"] = [i]
-    CSV.write(path_save_res, df_save_res, append=isfile(path_save_res))
-    CSV.write(path_save_parameters, df_save_parameters, append=isfile(path_save_parameters))
+Generate `n_multistarts` initial parameter guesses within the parameter bounds in the `petab_problem` with `sampling_method`
 
-    if !isnothing(path_save_trace)
-        df_save_trace = DataFrame(Matrix(reduce(vcat, res.xtrace')), θ_names)
-        df_save_trace[!, "f_trace"] = res.ftrace
-        df_save_trace[!, "Start_guess"] = repeat([i], length(res.ftrace))
-        CSV.write(path_save_trace, df_save_trace, append=isfile(path_save_trace))
-    end
-    return nothing
-end
+Any sampling algorithm from QuasiMonteCarlo is supported, but `LatinHypercubeSample` is recomended as it usually 
+performs well.
+
+If `n_multistarts` is set to 1, a single random vector within the parameter bounds is returned. For 
+`n_multistarts > 1`, a matrix is returned, with each column representing a different initial guess.
 
+By default `allow_inf_startguess=false` - only initial guesses that result in finite cost evaluations are returned. 
+If `allow_inf_startguess=true`, initial guesses that result in `Inf` are allowed.
 
+## Example
+```julia
+# Generate a single initial guess within the parameter bounds
+start_guess = generate_startguesses(petab_problem, 1)
+```
+
+```julia
+# Generate 10 initial guesses using Sobol sampling
+start_guess = generate_startguesses(petab_problem, 10, 
+                                    sampling_method=QuasiMonteCarlo.SobolSample())
+```
+"""
 function generate_startguesses(petab_problem::PEtabODEProblem,
-                               sampling_method::T,
                                n_multistarts::Int64;
-                               verbose::Bool=false)::Matrix{Float64} where T <: QuasiMonteCarlo.SamplingAlgorithm
+                               sampling_method::T=QuasiMonteCarlo.LatinHypercubeSample(),
+                               allow_inf_for_startguess::Bool=false,
+                               verbose::Bool=false)::Array{Float64} where T <: QuasiMonteCarlo.SamplingAlgorithm
 
     verbose == true && @info "Generating start-guesses"
 
     # Nothing prevents the user from sending in a parameter vector with zero parameters
     if length(petab_problem.lower_bounds) == 0
-        return nothing
+        return Vector{Float64}(undef, 0)
+    end
+
+    if n_multistarts == 1
+        while true
+            _p::Vector{Float64} = [rand() * (petab_problem.upper_bounds[j] - petab_problem.lower_bounds[j]) + petab_problem.lower_bounds[j] for j in eachindex(petab_problem.lower_bounds)]
+            _cost = petab_problem.compute_cost(_p)
+            if allow_inf_for_startguess == true
+                return _p
+            elseif !isinf(_cost)
+                return _p
+            end
+        end
     end
 
     startguesses = Matrix{Float64}(undef, length(petab_problem.lower_bounds), n_multistarts)
@@ -194,7 +209,10 @@ function generate_startguesses(petab_problem::PEtabODEProblem,
         for i in 1:size(_samples)[2]
             _p = _samples[:, i]
             _cost = petab_problem.compute_cost(_p)
-            if !isinf(_cost)
+            if allow_inf_for_startguess == true
+                found_starts += 1
+                startguesses[:, found_starts] .= _p
+            elseif !isinf(_cost)
                 found_starts += 1
                 startguesses[:, found_starts] .= _p
             end
@@ -209,6 +227,34 @@ function generate_startguesses(petab_problem::PEtabODEProblem,
 end
 
 
+function save_partial_results(path_save_res::String,
+                              path_save_parameters::String,
+                              path_save_trace::Union{String, Nothing},
+                              res::PEtabOptimisationResult,
+                              θ_names::Vector{Symbol},
+                              i::Int64)::Nothing
+
+    df_save_res = DataFrame(fmin=res.fmin,
+                          alg=String(res.alg),
+                          n_iterations = res.n_iterations,
+                          run_time = res.runtime,
+                          converged=string(res.converged),
+                          Start_guess=i)
+    df_save_parameters = DataFrame(Matrix(res.xmin'), θ_names)
+    df_save_parameters[!, "Start_guess"] = [i]
+    CSV.write(path_save_res, df_save_res, append=isfile(path_save_res))
+    CSV.write(path_save_parameters, df_save_parameters, append=isfile(path_save_parameters))
+
+    if !isnothing(path_save_trace)
+        df_save_trace = DataFrame(Matrix(reduce(vcat, res.xtrace')), θ_names)
+        df_save_trace[!, "f_trace"] = res.ftrace
+        df_save_trace[!, "Start_guess"] = repeat([i], length(res.ftrace))
+        CSV.write(path_save_trace, df_save_trace, append=isfile(path_save_trace))
+    end
+    return nothing
+end
+
+
 function _calibrate_model_multistart(petab_problem::PEtabODEProblem,
                                      alg,
                                      n_multistarts,
@@ -239,7 +285,7 @@ function _calibrate_model_multistart(petab_problem::PEtabODEProblem,
         end
     end
 
-    startguesses = generate_startguesses(petab_problem, sampling_method, n_multistarts)
+    startguesses = generate_startguesses(petab_problem, n_multistarts; sampling_method=sampling_method)
     if !isnothing(path_save_x0)
         startguessesDf = DataFrame(Matrix(startguesses)', petab_problem.θ_names)
         startguessesDf[!, "Start_guess"] = 1:size(startguessesDf)[1]

src/PEtab.jl

@@ -85,7 +85,6 @@ include(joinpath("SBML", "Solve_SBML_model.jl"))
 
 # For correct struct printing
 include(joinpath("Show.jl"))
-
 # Reduce time for reading a PEtabModel and for building a PEtabODEProblem
 @setup_workload begin
     path_yaml = joinpath(@__DIR__, "..", "test", "Test_model3", "Test_model3.yaml")
@@ -96,7 +95,7 @@ include(joinpath("Show.jl"))
     end
 end
 
-export PEtabModel, PEtabODEProblem, ODESolver, SteadyStateSolver, PEtabModel, PEtabODEProblem, remake_PEtab_problem, Fides, solve_SBML, PEtabOptimisationResult, IpoptOptions, IpoptOptimiser, PEtabParameter, PEtabObservable, PEtabMultistartOptimisationResult
+export PEtabModel, PEtabODEProblem, ODESolver, SteadyStateSolver, PEtabModel, PEtabODEProblem, remake_PEtab_problem, Fides, solve_SBML, PEtabOptimisationResult, IpoptOptions, IpoptOptimiser, PEtabParameter, PEtabObservable, PEtabMultistartOptimisationResult, generate_startguesses
 
 # These are given as extensions, but their docstrings are availble in the 
 # general documentation 

src/PEtab_structs.jl

@@ -3,7 +3,7 @@
 
 A Julia-compatible representation of a PEtab-specified problem.
 
-For how to construct see below.
+For constructor see below.
 
 !!! note
     Most of the functions in `PEtabModel` are not intended to be accessed by the user. For example, `compute_h` 
@@ -240,9 +240,11 @@ end
 
 
 """
+    PEtabODEProblem
+
 Everything needed to setup an optimization problem (compute cost, gradient, hessian and  parameter bounds) for a PEtab model.
 
-For how to construct, see below.
+For constructor, see below.
 
 !!! note
     The parameter vector θ is always assumed to be on the parameter scale specified in the PEtab parameters file. If needed, θ is transformed to the linear scale inside the function call.
@@ -402,7 +404,6 @@ If `hessian_method=nothing`, the Hessian method from the `PEtabODEProblem` is us
 
 For more information on each method, see the Fides [documentation](https://fides-optimizer.readthedocs.io/en/latest/generated/fides.hessian_approximation.html).
 
-See also [`callibrateModel`](@ref) and [`callibrateModelMultistart`](@ref)
 
 ## Examples
 ```julia

src/Process_PEtab_files/Julia_tables_provided/Create_PEtab_model.jl

@@ -85,6 +85,32 @@ function PEtabModel(system::ODESystem,
     return _PEtabModel(system, model_name, simulation_conditions, observables, measurements, 
                        petab_parameters, state_map, parameter_map, verbose)
 end
+"""
+    PEtabModel(system::Union{ReactionSystem, ODESystem},
+               observables::Dict{String, PEtabObservable},
+               measurements::DataFrame,
+               petab_parameters::Vector{PEtabParameter};
+               state_map::Union{Nothing, Vector{Pair}=nothing,
+               parameter_map::Union{Nothing, Vector{Pair}=nothing,
+               verbose::Bool=false)::PEtabModel
+
+Create a PEtabModel directly in Julia from a Catalyst ReactionSystem or MTK ODESystem without simulation conditions.
+
+In case of simulation conditions, see above.
+"""
+function PEtabModel(system::ODESystem,
+                    observables::Dict{String, PEtab.PEtabObservable},
+                    measurements::DataFrame,
+                    petab_parameters::Vector{PEtab.PEtabParameter};
+                    state_map::Union{Nothing, Vector{Pair{T1, Float64}}}=nothing,
+                    parameter_map::Union{Nothing, Vector{Pair{T2, Float64}}}=nothing,
+                    verbose::Bool=false)::PEtab.PEtabModel where {T1<:Union{Symbol, Num}, T2<:Union{Symbol, Num}}
+
+    simulation_conditions = Dict("__c0__" => Dict())                        
+    model_name = "ODESystemModel"                          
+    return _PEtabModel(system, model_name, simulation_conditions, observables, measurements, 
+                       petab_parameters, state_map, parameter_map, verbose)
+end
 
 
 function _PEtabModel(system,