Record parameter calibration values progress in output.

[TorkelE]

When calibrating a parameter set to data, and the optimizer iterates through parameter values: It would be useful to receive (if desired) the entire path in the output (and not just the best parameter set and loss function value). Either just all the loss function values or possibly the corresponding parameter sets as well. Is there a way to get this in the output? I feel that it is often useful to plot this to gain some idea of the performance of the optimising run.

[sebapersson]

I agree with you that this is useful information. Currently for the Optimizer's I wrap in PEtab.jl (Fides, Ipopt and Optim with box-constraints) for which I support multi-start optimization this information is not available. However, it is definitely possible to retrieve by writing a callback (how to write the callback function depends on the Optimization-package though so takes some time).

Overall, if we would like to carry out the parameter-estimation in PEtab (which is probably a good idea as we can write a method that takes PEtabODEProblem and then calibrates the model) a better output struct from a single calibration run is likely needed, kind of something like:

struct PEtabOptimisationResult
     trace::Vector{Vector{Float64}} # Parameter vectors (if user wants to save them)
     cost_trace::Vector{Float64} # Likelihood value (if user wants to save them)
     n_iterations # Number of iterations optimiser
     cost_best # Best optimised value 
     p_best # Last parameter value 
     convergence_criteria # If user wants to 
     run_time # Always fun :)
end

Is there anymore output you think currently is needed? Otherwise I could take a look at writing a better optimization struct.

[TorkelE]

Would it make sense to have a single interface for when multi-start approaches are used? In these cases it would be useful to have all information for e.g. parameters plots:

Then each start would have this kind of information, with the output structure containing these one a vector, as weöö as the best parameter and its cost value?

[sebapersson]

You are thinking that we have a function like perform_multistart_callibration(petab_problem, optimizer_alg) which then outputs a struct-like?

struct PEtabMultistartOptimisationResult
     best_vector::Vector{Vector{Float64}} # Parameter vectors (if user wants to save them)
     best_cost::Vector{Float64} # Likelihood value (if user wants to save them)
     n_multistarts::Int
     run_results::Vector{PEtabOptimisationResult} # See above
end

From a struct like this it would then be easy to write a function for generating the plot above. If something like this looks like a good format I can code a prototype.

[TorkelE]

Yes, that looks ideal :)

[sebapersson]

I have now added a struct and interface for multistart-parameter estimation:

struct PEtabMultistartOptimisationResult
    xMin::Vector{Float64} # Parameter vectors (if user wants to save them)
    fMin::Float64 # Likelihood value (if user wants to save them)
    nMultistarts::Int
    alg::Symbol
    multiStartMethod::String
    dirSave::Union{String, Nothing}
    run::Vector{PEtabOptimisationResult} # See above
end

And now it should, via package extensions, be possible to perform parameter estimation either single-start callibrateModel or multistart callibrateModelMultistart with Fides (requires PyCall), Ipopt or Optim; as in this example:

using Optim
using PyCall
using Ipopt
using QuasiMonteCarlo
using CSV
using Catalyst
using DataFrames
using Distributions
using PEtab

# Define reaction network model
rn = @reaction_network begin
    @parameters a0
    @species A(t)=a0
    (k1, k2), A <--> B
end
state_map = [:B => 1.0] # Constant initial value for B

# Setup PEtab problem 
measurements = DataFrame(simulation_id=["c0", "c0", "c1", "c1"],
                         obs_id=["obs_a", "obs_a", "obs_a", "obs_a"],
                         time=[0, 10.0, 0, 10.0],
                         measurement=[0.7, 0.1, 0.8, 0.2])
simulation_conditions = Dict("c0" => Dict(:a0 => 0.8),
                             "c1" => Dict(:a0 => 0.9))
petab_parameters = [PEtabParameter(:k1, value=0.8, scale=:lin),
                    PEtabParameter(:k2, value=0.6, scale=:lin)]
@unpack A = rn
observables = Dict(["obs_a" => PEtabObservable(A, :lin, 1.0)])
petab_model = readPEtabModel(rn, simulation_conditions, observables, measurements,
                            petab_parameters, verbose=false, stateMap=state_map)
petab_problem = createPEtabODEProblem(petab_model, verbose=false)

# Callibrate the model with different methods
saveTrace = true
p0 = petab_problem.θ_nominalT .* 0.4
# Fides = NewtonTrustRegion with box-constraints
res = callibrateModel(petab_problem, p0, Fides(verbose=true);
                      saveTrace=true)
res = callibrateModel(petab_problem, p0, Fides(:BFGS, verbose=true);
                      saveTrace=true)
# Optim Interior point Newton                      
res = callibrateModel(petab_problem, p0, Optim.IPNewton();
                      saveTrace=true)
# Ipopt Interior point Newton                       
res = callibrateModel(petab_problem, p0, IpoptOptimiser(approximateHessian=false);
                      saveTrace=true, options=IpoptOptions(print_level=5))
res = callibrateModel(petab_problem, p0, IpoptOptimiser(approximateHessian=true);
                      saveTrace=true, options=IpoptOptions(print_level=5))

# Mutlistart optimisation using LathinHypercube sampling, performing 100 multistarts 
dirSave = joinpath(@__DIR__, "Callibration")
res = callibrateModelMultistart(petab_problem, Optim.IPNewton(), 100, dirSave, saveTrace=true)

saveTrace allows the trace to be saved and works for all algorithms except Fides (I will file an issue on this). Any more output you could be useful for optimisation?

For optimisation it is also possible for the user to set options via the option-interface for the package, e.g the following is allowed:

# Mutlistart optimisation using LathinHypercube sampling, performing 100 multistarts 
dirSave = joinpath(@__DIR__, "Callibration")
res = callibrateModelMultistart(petab_problem, Optim.IPNewton(), 100, dirSave, saveTrace=true, 
                                                 options=Optim.Options(x_abstol=1e-8))

[TorkelE]

Awsome! I hope to start on Catalyst documentation for this this weekend.

[TorkelE]

I think this is more or less solved now, but then, maybe it is not:

When running calibrate_model_multistart, the output is saved to a file, and a PEtabMultistartOptimisationResult optimisation struct is returned. We are now writing plotting functions that take a PEtabMultistartOptimisationResult, and it seems possible that these will be used as the input for other types of analysis. Would it make sense for calibrate_model_multistart to either

1. When(/if) the runs finish, save the PEtabMultistartOptimisationResult to the output folder as well.
2. Create a constructor PEtabMultistartOptimisationResult(folderpath::String) which for a given folder (to which calibrate_model_multistart wrote the results) return a PEtabMultistartOptimisationResult.
   ?
   (I tried to see if this was already possible, but didn't find anything)

[sebapersson]

Sorry did not see the reply here.

I think 2 is a good solution PEtabMultistartOptimisationResult(folderpath::String) - I basically have functionality for this and will add it.

[TorkelE]

Sounds great :)