Support for forward sensitivities (#393)

* Minor fix for enabling running FSA. Tested with Boehm_JProteomeRes2014. Works.
* Updated example
* Slightly changed output to indicate which sensitivity method was used

Closes #343

doc/optimizationApplication.md

@@ -40,7 +40,9 @@ Run the created executable with the `-h`/`--help` argument.
 - llh: Log-likelihood as provided by AMICI
 - status: AMICI status
 - time: simulation wall-time
-- gradientFlag: '+' indicates simulation with gradient, '-' indicates simulation without gradient   
+- gradientFlag: '-' indicates simulation without gradient computation;
+  'A' indicates gradient computation using adjoint sensitivity analysis;
+  'F' indicates gradient computation using forward sensitivity analysis
 
 ## Environment variables
 

examples/parpeamici/steadystate/parpeExampleSteadystateBasic.ipynb

@@ -34,6 +34,7 @@
     "import h5py\n",
     "from importlib import reload\n",
     "import parpe\n",
+    "from pprint import pprint\n",
     "\n",
     "# set paths\n",
     "parpe_source_root = os.path.abspath('../../../')\n",
@@ -1065,10 +1066,97 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 43,
    "metadata": {},
+   "source": [
+    "# let's create a copy of the previous input file\n",
+    "input_file_fsa = f'{example_data_dir}/example_data_fsa.h5'\n",
+    "!cp {input_file} {input_file_fsa}\n",
+    "\n",
+    "# first, we just look at the current settings\n",
+    "with h5py.File(input_file_fsa, 'r') as f:\n",
+    "    pprint(dict(f['/amiciOptions'].attrs))\n"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": [
+    "There is currently no detailed documentation of the attributes available, but they can be looked up in `amici::writeSolverSettingsToHDF5` in `deps/AMICI/src/hdf5.cpp`. There are additional settings available to the ones shown above.\n",
+    "The relevant attributes for sensitivity analysis are `sensi_meth` (and `sensi_meth_preeq`, if preequilibration is used). This is how you can change them:"
+   ]
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "with h5py.File(input_file_fsa, 'r+') as f:\n",
+    "    print(f'Before: sensi_meth={amici.SensitivityMethod(f[\"/amiciOptions\"].attrs[\"sensi_meth\"])!r}')\n",
+    "    print(f'Before: sensi_meth_preeq={amici.SensitivityMethod(f[\"/amiciOptions\"].attrs[\"sensi_meth_preeq\"])!r}')\n",
+    "    f[\"/amiciOptions\"].attrs[\"sensi_meth\"] = amici.SensitivityMethod.forward\n",
+    "    f[\"/amiciOptions\"].attrs[\"sensi_meth_preeq\"] = amici.SensitivityMethod.forward\n",
+    "    print(f'After: sensi_meth={amici.SensitivityMethod(f[\"/amiciOptions\"].attrs[\"sensi_meth\"])!r}')\n",
+    "    print(f'After: sensi_meth_preeq={amici.SensitivityMethod(f[\"/amiciOptions\"].attrs[\"sensi_meth_preeq\"])!r}')\n"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": "Other settings can be changed in the same way. The important ones will probably be the integration tolerances and the maximum number of integration steps."
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": [
+    "### FSA vs. ASA\n",
+    "\n",
+    "We can repeat the optimization with the updated input file and compare the results from forward sensitivity analysis to adjoint sensitivity analysis:"
+   ]
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": "!PARPE_NO_DEBUG=1 {example_binary_dir}/example_steadystate_multi -o deleteme-fsa/ {input_file_fsa}\n",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "filename = 'deleteme-fsa/_rank00000.h5'\n",
+    "trajectories_ipopt_asa = trajectories_ipopt.copy()\n",
+    "trajectories_ipopt_fsa = parpe.getCostTrajectories(filename)\n",
+    "ax = parpe.plotting.plotCostTrajectory(trajectories_ipopt_asa, log=False)\n",
+    "parpe.plotting.plotCostTrajectory(trajectories_ipopt_fsa, log=False, ax=ax);"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# pad trajectories to the same length\n",
+    "tmp = np.full((np.fmax(len(trajectories_ipopt_asa), len(trajectories_ipopt_fsa)), 2), np.nan)\n",
+    "tmp[:len(trajectories_ipopt_asa), 0] = trajectories_ipopt_asa.squeeze()\n",
+    "tmp[:len(trajectories_ipopt_fsa), 1] = trajectories_ipopt_fsa.squeeze()\n",
+    "\n",
+    "# create dataframe\n",
+    "df = pd.DataFrame(tmp, columns=['fval Ipopt ASA', 'fval Ipopt FSA'])\n",
+    "df[\"diff\"] = df['fval Ipopt ASA'] - df['fval Ipopt FSA']\n",
+    "df"
+   ],
    "outputs": [],
-   "source": []
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": "As expected, the differences are tiny."
   }
  ],
  "metadata": {
@@ -1110,4 +1198,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}

examples/parpeamici/steadystate/steadyStateMultiConditionDataprovider.cpp

@@ -27,8 +27,8 @@ void SteadyStateMultiConditionDataProvider::setupModelAndSolver() const {
     //model.requireSensitivitiesForAllParameters();
 
     solver_->setSensitivityOrder(amici::SensitivityOrder::first);
-    solver_->setSensitivityMethod(amici::SensitivityMethod::adjoint);
-    solver_->setMaxSteps(10000);
-    solver_->setNewtonMaxSteps(40);
+    // solver_->setSensitivityMethod(amici::SensitivityMethod::adjoint);
+    // solver_->setMaxSteps(10000);
+    // solver_->setNewtonMaxSteps(40);
 }
 

src/parpeamici/multiConditionProblem.cpp

@@ -213,14 +213,22 @@ void printSimulationResult(Logger *logger, int jobId,
         return;
     }
 
-    bool with_sensi = rdata->sensi >= amici::SensitivityOrder::first;
+    char sensi_mode = '-';
+    if(rdata->sensi >= amici::SensitivityOrder::first
+        && rdata->sensi_meth == amici::SensitivityMethod::adjoint) {
+        sensi_mode = 'A';
+    } else if(rdata->sensi >= amici::SensitivityOrder::first
+               && rdata->sensi_meth == amici::SensitivityMethod::forward) {
+        sensi_mode = 'F';
+    }
+    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)",
                        jobId, rdata->llh, rdata->status,
                        rdata->numsteps.empty()?-1:rdata->numsteps[rdata->numsteps.size() - 1],
                        rdata->numstepsB.empty()?-1:rdata->numstepsB[0],
                        timeSeconds,
-                       with_sensi?'+':'-');
+                       sensi_mode);
 
 
     // check for NaNs, only report first
@@ -368,7 +376,7 @@ AmiciSimulationRunner::AmiciResultPackageSimple runAndLogSimulation(
 
             if(solver->getSensitivityOrder() >= amici::SensitivityOrder::first
                 && solver->getSensitivityMethod()
-                       == amici::SensitivityMethod::adjoint) {
+                       != amici::SensitivityMethod::none) {
                 solver->setReturnDataReportingMode(amici::RDataReporting::likelihood);
             } else {
                 // unset sensitivity method, because `residuals` is not allowed