Merge pull request #172 from sebapersson/Bugfix_prior

Fix bug in gradient for prior, test gradient, and test block Hessian
sebapersson · Feb 16, 2024 · 3222155 · 3222155
2 parents 536198d + 0a25f81
commit 3222155
Show file tree

Hide file tree

Showing 23 changed files with 1,208 additions and 129 deletions.
diff --git a/ext/PEtabPlotsExtension/optimisation_trajectory_recipes.jl b/ext/PEtabPlotsExtension/optimisation_trajectory_recipes.jl
@@ -1,10 +1,12 @@
 # The possible types of plots avaiable for PEtabOptimisationResult.
 const plot_types = [:objective, :best_objective]
-const plot_types_ms = [:objective,
-                       :best_objective,
-                       :waterfall,
-                       :runtime_eval,
-                       :parallel_coordinates]
+const plot_types_ms = [
+    :objective,
+    :best_objective,
+    :waterfall,
+    :runtime_eval,
+    :parallel_coordinates
+]
 
 # Plots the objective function progression for a PEtabOptimisationResult.
 @recipe function f(res::PEtabOptimisationResult{T}; plot_type = :best_objective) where {T}

diff --git a/src/Derivatives/Gradient.jl b/src/Derivatives/Gradient.jl
@@ -22,6 +22,9 @@ function compute_gradient_autodiff!(gradient::Vector{Float64},
     # We need to track a variable if ODE system could be solve as checking retcode on solution array it not enough.
     # This is because for ForwardDiff some chunks can solve the ODE, but other fail, and thus if we check the final
     # retcode we cannot catch these cases
+    # We need to track a variable if ODE system could be solve as checking retcode on solution array it not enough.
+    # This is because for ForwardDiff some chunks can solve the ODE, but other fail, and thus if we check the final
+    # retcode we cannot catch these cases
     simulation_info.could_solve[1] = true
 
     # Case where based on the original PEtab file read into Julia we do not have any parameter vectors fixated.
@@ -213,10 +216,10 @@ function compute_gradient_prior!(gradient::Vector{Float64},
                                  θ::Vector{<:Real},
                                  θ_indices::ParameterIndices,
                                  prior_info::PriorInfo)::Nothing
-    _evalPriors = (θ_est) -> begin
+    _eval_priors = (θ_est) -> begin
         θ_estT = transformθ(θ_est, θ_indices.θ_names, θ_indices)
         return -1.0 * compute_priors(θ_est, θ_estT, θ_indices.θ_names, prior_info) # We work with -loglik
     end
-    gradient .+= ForwardDiff.gradient(_evalPriors, θ)
+    gradient .+= ForwardDiff.gradient(_eval_priors, θ)
     return nothing
 end
diff --git a/src/Objective/Priors.jl b/src/Objective/Priors.jl
@@ -4,10 +4,6 @@ function compute_priors(θ_parameter_scale::Vector{T},
                         θ_linear_scale::Vector{T},
                         θ_names::Vector{Symbol},
                         prior_info::PriorInfo)::T where {T <: Real}
-    if prior_info.has_priors == false
-        return 0.0
-    end
-
     prior_value = 0.0
     for (θ_name, logpdf) in prior_info.logpdf
         iθ = findfirst(x -> x == θ_name, θ_names)
@@ -18,6 +14,5 @@ function compute_priors(θ_parameter_scale::Vector{T},
             prior_value += logpdf(θ_linear_scale[iθ])
         end
     end
-
     return prior_value
 end
diff --git a/src/PEtabODEProblem/Create.jl b/src/PEtabODEProblem/Create.jl
@@ -265,6 +265,8 @@ function PEtabODEProblem(petab_model::PEtabModel;
         return measurement_info.simulated_values
     end
 
+    # Computing nllh along with the gradient is needed for efficient Bayesian
+    # inference, as for example AdvancedHMC.jl needs both nllh and gradient
     # Computing nllh along with the gradient is needed for efficient Bayesian
     # inference, as for example AdvancedHMC.jl needs both nllh and gradient
     # in its evaluations.
@@ -631,7 +633,9 @@ function create_gradient_function(which_method::Symbol,
                                                                             prior_info,
                                                                             cfg,
                                                                             petab_ODE_cache,
-                                                                            exp_id_solve = [:all],
+                                                                            exp_id_solve = [
+                                                                                :all
+                                                                            ],
                                                                             split_over_conditions = split_over_conditions,
                                                                             isremade = isremade)
         end
@@ -810,8 +814,7 @@ function create_hessian_function(which_method::Symbol,
 
             _chunksize = isnothing(chunksize) ? ForwardDiff.Chunk(θ_dynamic) :
                          ForwardDiff.Chunk(chunksize)
-            cfg = ForwardDiff.HessianConfig(compute_cost_θ_dynamic, θ_dynamic,
-                                            ForwardDiff.Chunk(chunksize))
+            cfg = ForwardDiff.HessianConfig(compute_cost_θ_dynamic, θ_dynamic, _chunksize)
 
             _compute_hessian! = let compute_cost_θ_not_ODE = compute_cost_θ_not_ODE,
                 compute_cost_θ_dynamic = compute_cost_θ_dynamic,
@@ -990,7 +993,9 @@ function create_hessian_function(which_method::Symbol,
                                                                            cfg,
                                                                            cfg_not_solve_ode,
                                                                            petab_ODE_cache,
-                                                                           exp_id_solve = [:all],
+                                                                           exp_id_solve = [
+                                                                               :all
+                                                                           ],
                                                                            reuse_sensitivities = reuse_sensitivities,
                                                                            return_jacobian = return_jacobian,
                                                                            split_over_conditions = split_over_conditions,

diff --git a/src/Process_input/Julia_input.jl b/src/Process_input/Julia_input.jl
@@ -257,8 +257,8 @@ function parse_petab_measurements(petab_measurements::DataFrame,
                                   simulation_conditions,
                                   petab_parameters::Vector{PEtabParameter})::DataFrame
     allowed_names = ["time", "obs_id", "observable_id", "noise_parameters", "measurement",
-                     "simulation_id", "pre_equilibration_id", "pre_eq_id",
-                     "observable_parameters"]
+        "simulation_id", "pre_equilibration_id", "pre_eq_id",
+        "observable_parameters"]
     column_names = names(petab_measurements)
     for name in column_names
         if name ∈ allowed_names

diff --git a/src/Process_input/Observables/h_sigma_derivatives.jl b/src/Process_input/Observables/h_sigma_derivatives.jl
@@ -64,7 +64,7 @@ function create_derivative_σ_h_file(model_name::String,
                                         parameter_map, state_map, experimental_conditions)
 
     # Dummary variables to keep PEtab importer happy even as we are not providing any PEtab files
-    model_SBML = SBMLImporter.ModelSBML()
+    model_SBML = SBMLImporter.ModelSBML("")
 
     ∂h∂u_str, ∂h∂p_str = PEtab.create∂h∂_function(model_name, @__DIR__, model_state_names,
                                                   parameter_info, p_ode_problem_names,

diff --git a/src/Process_input/Observables/u0_h_sigma.jl b/src/Process_input/Observables/u0_h_sigma.jl
@@ -75,7 +75,7 @@ function create_σ_h_u0_file(model_name::String,
                                         parameter_map, state_map, experimental_conditions)
 
     # Dummary variables to keep PEtab importer happy even as we are not providing any PEtab files
-    model_SBML = SBMLImporter.ModelSBML()
+    model_SBML = SBMLImporter.ModelSBML("")
 
     h_str = PEtab.create_h_function(model_name, @__DIR__, model_state_names, parameter_info,
                                     p_ode_problem_names,
@@ -227,7 +227,7 @@ end
                        dir_model::String,
                        parameter_info::ParametersInfo,
                        p_ode_problem_names::Vector{String},
-                       state_map, 
+                       state_map,
                        model_SBML;
                        inplace::Bool=true)
 
@@ -270,7 +270,7 @@ function create_u0_function(model_name::String,
     # Write the formula for each initial condition to file
     _model_state_names = [replace.(string.(state_map[i].first), "(t)" => "")
                           for i in eachindex(state_map)]
-    # If we create from Julia model (e.g.) Catalyst this is not applicable and step is skipped 
+    # If we create from Julia model (e.g.) Catalyst this is not applicable and step is skipped
     if !isempty(model_SBML.parameters)
         model_state_names1 = filter(x -> x ∈ keys(model_SBML.species) &&
                                         model_SBML.species[x].assignment_rule == false,

diff --git a/src/Process_input/Table_input/Read_tables.jl b/src/Process_input/Table_input/Read_tables.jl
@@ -75,14 +75,16 @@ function check_petab_files(experimental_conditions::CSV.File, measurements_data,
     transformation_types = ["lin", "log", "log10"]
     distribution_types = ["laplace", "normal"]
     estimateTypes = [0, 1]
-    priorParameterTypes = ["uniform",
-                           "normal",
-                           "laplace",
-                           "logNormal",
-                           "logLaplace",
-                           "parameterScaleUniform",
-                           "parameterScaleNormal",
-                           "parameterScaleLaplace"]
+    priorParameterTypes = [
+        "uniform",
+        "normal",
+        "laplace",
+        "logNormal",
+        "logLaplace",
+        "parameterScaleUniform",
+        "parameterScaleNormal",
+        "parameterScaleLaplace"
+    ]
 
     # Check experimental_conditions
     columns_check = ["conditionId", "conditionName"]
@@ -106,74 +108,88 @@ function check_petab_files(experimental_conditions::CSV.File, measurements_data,
     end
 
     # Check measurements_data
-    columns_check = ["observableId",
-                     "simulationConditionId",
-                     "measurement",
-                     "time",
-                     "preequilibrationConditionId",
-                     "observableParameters",
-                     "noiseParameters",
-                     "datasetId",
-                     "replicateId"]
-    allowed_types = [string_types,
-                     string_types,
-                     number_types,
-                     number_types,
-                     string_types,
-                     string_number_types,
-                     string_number_types,
-                     string_number_types,
-                     string_number_types]
+    columns_check = [
+        "observableId",
+        "simulationConditionId",
+        "measurement",
+        "time",
+        "preequilibrationConditionId",
+        "observableParameters",
+        "noiseParameters",
+        "datasetId",
+        "replicateId"
+    ]
+    allowed_types = [
+        string_types,
+        string_types,
+        number_types,
+        number_types,
+        string_types,
+        string_number_types,
+        string_number_types,
+        string_number_types,
+        string_number_types
+    ]
     required_columns = ["observableId", "simulationConditionId", "measurement", "time"]
     check_df_columns(measurements_data, "measurements_data", columns_check, allowed_types,
                      required_columns)
 
     # Check parameters_data
-    columns_check = ["parameterId",
-                     "parameterScale",
-                     "lowerBound",
-                     "upperBound",
-                     "nominalValue",
-                     "estimate",
-                     "parameterName",
-                     "initializationPriorType",
-                     "initializationPriorParameters",
-                     "objectivePriorType",
-                     "objectivePriorParameters"]
+    columns_check = [
+        "parameterId",
+        "parameterScale",
+        "lowerBound",
+        "upperBound",
+        "nominalValue",
+        "estimate",
+        "parameterName",
+        "initializationPriorType",
+        "initializationPriorParameters",
+        "objectivePriorType",
+        "objectivePriorParameters"
+    ]
     # Some models have missing parameters in their bounds even though it's mandatory, so we add Missing as an allowed DataType for these columns.
-    allowed_types = [string_types,
-                     transformation_types,
-                     number_types,
-                     number_types,
-                     number_types,
-                     estimateTypes,
-                     string_types,
-                     priorParameterTypes,
-                     string_types,
-                     priorParameterTypes,
-                     string_types]
-    required_columns = ["parameterId",
-                        "parameterScale",
-                        "lowerBound",
-                        "upperBound",
-                        "nominalValue",
-                        "estimate"]
+    allowed_types = [
+        string_types,
+        transformation_types,
+        number_types,
+        number_types,
+        number_types,
+        estimateTypes,
+        string_types,
+        priorParameterTypes,
+        string_types,
+        priorParameterTypes,
+        string_types
+    ]
+    required_columns = [
+        "parameterId",
+        "parameterScale",
+        "lowerBound",
+        "upperBound",
+        "nominalValue",
+        "estimate"
+    ]
     check_df_columns(parameters_data, "parameters_data", columns_check, allowed_types,
                      required_columns)
 
     # Check observables_data
-    columns_check = ["observableId",
-                     "observableFormula",
-                     "noiseFormula",
-                     "observableName",
-                     "observableTransformation",
-                     "noiseDistribution"]
-    allowed_types = [string_types,
-                     string_types,
-                     string_number_types,
-                     string_types,
-                     transformation_types,
-                     distribution_types]
+    columns_check = [
+        "observableId",
+        "observableFormula",
+        "noiseFormula",
+        "observableName",
+        "observableTransformation",
+        "noiseDistribution"
+    ]
+    allowed_types = [
+        string_types,
+        string_types,
+        string_number_types,
+        string_types,
+        transformation_types,
+        distribution_types
+    ]
     required_columns = ["observableId", "observableFormula", "noiseFormula"]
     check_df_columns(observables_data, "observables_data", columns_check, allowed_types,
                      required_columns)

diff --git a/test/Julia_import/Example.jl b/test/Julia_import/Example.jl