PEtab pipeline: Only compute sensitivities with respect to optimized …

…parameters * Make use of amici::ExpData::plist to compute only relevant sensitivities for each simulation condition. * Update logging for nplist * Include Bruno_JExpBot2016 on GHA, which is a great test case for plist Closes #380
ICB-DCM · Oct 19, 2024 · b3a6845 · b3a6845
1 parent c76c7d3
commit b3a6845
Show file tree

Hide file tree

Showing 6 changed files with 54 additions and 16 deletions.
diff --git a/benchmark_collection/all.sh b/benchmark_collection/all.sh
@@ -10,6 +10,7 @@ set -eou pipefail
 expected_to_work="
 Boehm_JProteomeRes2014
 Borghans_BiophysChem1997
+Bruno_JExpBio2016
 Crauste_CellSystems2017
 Elowitz_Nature2000
 Fujita_SciSignal2010
@@ -25,7 +26,6 @@ Zheng_PNAS2012
 #
 # Alkan_SciSignal2018
 # Blasi_CellSystems2016
-# Bruno_JExpBio2016
 # Giordano_Nature2020
 # Okuonghae_ChaosSolitonsFractals2020
 # Perelson_Science1996

diff --git a/doc/optimizationApplication.md b/doc/optimizationApplication.md
@@ -42,7 +42,9 @@ Run the created executable with the `-h`/`--help` argument.
 - time: simulation wall-time
 - gradientFlag: '-' indicates simulation without gradient computation;
   'A' indicates gradient computation using adjoint sensitivity analysis;
-  'F' indicates gradient computation using forward sensitivity analysis
+  'F' indicates gradient computation using forward sensitivity analysis;
+  this is followed by a number indicating w.r.t. how many parameters the
+  gradient was computed.
 
 ## Environment variables
 

diff --git a/include/parpeamici/multiConditionDataProvider.h b/include/parpeamici/multiConditionDataProvider.h
@@ -311,6 +311,13 @@ class MultiConditionDataProviderHDF5 : public MultiConditionDataProvider
     std::unique_ptr<amici::ExpData> getExperimentalDataForCondition(
       int conditionIdx) const override;
 
+    /**
+     * @brief Get list of parameters w.r.t. which we need sensitivities
+     * @param mapping Mapping from model simulation to objective parameters
+     * @return AMICI's 'plist'
+     */
+    std::vector<int> getSensitivityParameterList(std::vector<int> const& mapping) const;
+
     std::vector<std::vector<double>> getAllMeasurements() const override;
     std::vector<std::vector<double>> getAllSigmas() const override;
 

diff --git a/python/parpe/hdf5_pe_input.py b/python/parpe/hdf5_pe_input.py
@@ -277,6 +277,7 @@ def _generate_simulation_to_optimization_parameter_mapping(self) -> None:
         self.parameter_mapping = self.petab_problem \
             .get_optimization_to_simulation_parameter_mapping(
             warn_unmapped=False, scaled_parameters=False,
+            fill_fixed_parameters=True,
             allow_timepoint_specific_numeric_noise_parameters=True)
 
         variable_par_ids = self.amici_model.getParameterIds()

diff --git a/src/parpeamici/multiConditionDataProvider.cpp b/src/parpeamici/multiConditionDataProvider.cpp
@@ -110,18 +110,16 @@ MultiConditionDataProviderHDF5::mapSimulationToOptimizationGradientAddMultiply(
   double coefficient) const
 {
     auto mapping = getSimulationToOptimizationParameterMapping(conditionIdx);
-
+    auto plist = getSensitivityParameterList(mapping);
     // Need to consider varying scaling
     auto scaleOpt = getParameterScaleOpt();
     auto scaleSim = getParameterScaleSim(conditionIdx);
 
-    for (int i = 0; i < model_->np(); ++i) {
-        // some model parameter are not mapped if there is no respective data
-        if (mapping[i] >= 0) {
-            double newGrad = applyChainRule(
-              simulation[i], parameters[i], scaleSim[i], scaleOpt[mapping[i]]);
-            optimization[mapping[i]] += coefficient * newGrad;
-        }
+    for(std::vector<int>::size_type i_plist = 0; i_plist < plist.size(); ++i_plist) {
+        auto i_p = plist[i_plist];
+        double newGrad = applyChainRule(simulation[i_plist], parameters[i_p],
+                                        scaleSim[i_p], scaleOpt[mapping[i_p]]);
+        optimization[mapping[i_p]] += coefficient * newGrad;
     }
 }
 
@@ -153,6 +151,8 @@ MultiConditionDataProviderHDF5::mapAndSetOptimizationToSimulationVariables(
     }
 
     for (int i = 0; i < model_->np(); ++i) {
+        // a negative index means there is no mapping to a problem parameter,
+        //  but that there is a fixed value
         if (mapping[i] >= 0) {
             // map from optimization parameters
             simulation[i] = getScaledParameter(
@@ -310,11 +310,33 @@ std::unique_ptr<amici::ExpData> MultiConditionDataProviderHDF5::getExperimentalD
                     + std::to_string(simulationIdx)));
     edata->setObservedData(getMeasurementForSimulationIndex(simulationIdx));
     edata->setObservedDataStdDev(getSigmaForSimulationIndex(simulationIdx));
+    // fixed parameters and state variable indices for reinitialization
     updateFixedSimulationParameters(simulationIdx, *edata);
 
+    auto mapping = getSimulationToOptimizationParameterMapping(simulationIdx);
+    edata->plist = getSensitivityParameterList(mapping);
+
     return edata;
 }
 
+std::vector<int> MultiConditionDataProviderHDF5::getSensitivityParameterList(
+    const std::vector<int> &mapping) const {
+    // We can construct plist simply from the parameter mapping:
+    //  If a model parameter does not map to a problem parameter (w.r.t. all
+    //  of which we need to compute sensitivities), there is a negative index
+    //  in the mapping. That means, plist is the list of indices for which
+    //  the mapping is non-negative.
+    int nplist = std::count_if(mapping.begin(), mapping.end(), [](int i) { return i >= 0; });
+    int i_plist = 0;
+    std::vector<int> plist(nplist);
+    for (int i_p = 0; i_p < model_->np(); ++i_p) {
+        if (mapping[i_p] >= 0) {
+            plist[i_plist++] = i_p;
+        }
+    }
+    return plist;
+}
+
 std::vector<std::vector<double>>
 MultiConditionDataProviderHDF5::getAllMeasurements() const
 {
@@ -618,6 +640,7 @@ MultiConditionDataProviderDefault ::
     // TODO redundant
     auto mapping = getSimulationToOptimizationParameterMapping(conditionIdx);
 
+    // TODO respect plist
     for (int i = 0; i < model_->np(); ++i) {
         optimization[mapping[i]] = coefficient * simulation[i];
     }

diff --git a/src/parpeamici/multiConditionProblem.cpp b/src/parpeamici/multiConditionProblem.cpp
@@ -223,17 +223,17 @@ void printSimulationResult(Logger *logger, int jobId,
     }
     bool with_sensi = rdata->sensi >= amici::SensitivityOrder::first && rdata->sensi_meth != amici::SensitivityMethod::none;
 
-    logger->logmessage(loglevel::debug, "Result for %d: %g (%d) (%d/%d/%.4fs/%c)",
+    logger->logmessage(loglevel::debug, "Result for %d: %g (%d) (%d/%d/%.4fs/%c%d)",
                        jobId, rdata->llh, rdata->status,
                        rdata->numsteps.empty()?-1:rdata->numsteps[rdata->numsteps.size() - 1],
                        rdata->numstepsB.empty()?-1:rdata->numstepsB[0],
                        timeSeconds,
-                       sensi_mode);
-
+                       sensi_mode,
+                       rdata->nplist);
 
     // check for NaNs, only report first
     if (with_sensi) {
-        for (int i = 0; i < rdata->np; ++i) {
+        for (int i = 0; i < rdata->sllh.size(); ++i) {
             if (std::isnan(rdata->sllh[i])) {
                 logger->logmessage(loglevel::debug,
                                    "Gradient contains NaN at %d", i);
@@ -272,13 +272,17 @@ void saveSimulation(const H5::H5File &file, std::string const& pathStr,
                 file, fullGroupPath, "jobId",
                 gsl::make_span<const int>(&jobId, 1));
 
-    if (!gradient.empty()) {
+    // save sllh; but sllh may vary depending on the condition-specific plist
+    if (gradient.size() == parameters.size()) {
         hdf5CreateOrExtendAndWriteToDouble2DArray(
             file, fullGroupPath, "simulationLogLikelihoodGradient",
                     gradient);
     } else if(!parameters.empty()) {
         double dummyGradient[parameters.size()];
         std::fill_n(dummyGradient, parameters.size(), NAN);
+        if(!gradient.empty()) {
+            std::copy_n(gradient.begin(), std::min(gradient.size(), parameters.size()), dummyGradient);
+        }
         hdf5CreateOrExtendAndWriteToDouble2DArray(
             file, fullGroupPath, "simulationLogLikelihoodGradient",
             gsl::make_span<const double>(dummyGradient, parameters.size()));
@@ -474,7 +478,8 @@ AmiciSimulationRunner::AmiciResultPackageSimple runAndLogSimulation(
             rdata->llh,
             timeSeconds,
             (solverTemplate.getSensitivityOrder()
-             > amici::SensitivityOrder::none)
+             > amici::SensitivityOrder::none
+             && solverTemplate.getSensitivityMethod() != amici::SensitivityMethod::none)
                 ? rdata->sllh : std::vector<double>(),
             rdata->y, rdata->sigmay,
             sendStates ? rdata->x : std::vector<double>(),