diff --git a/doc/optimizationApplication.md b/doc/optimizationApplication.md index 5c8b3655..57f3a0c4 100644 --- a/doc/optimizationApplication.md +++ b/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 diff --git a/examples/parpeamici/steadystate/parpeExampleSteadystateBasic.ipynb b/examples/parpeamici/steadystate/parpeExampleSteadystateBasic.ipynb index cb7b6f80..f7217b50 100644 --- a/examples/parpeamici/steadystate/parpeExampleSteadystateBasic.ipynb +++ b/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 -} \ No newline at end of file +} diff --git a/examples/parpeamici/steadystate/steadyStateMultiConditionDataprovider.cpp b/examples/parpeamici/steadystate/steadyStateMultiConditionDataprovider.cpp index 5d33657d..e8fcd64c 100644 --- a/examples/parpeamici/steadystate/steadyStateMultiConditionDataprovider.cpp +++ b/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); } diff --git a/src/parpeamici/multiConditionProblem.cpp b/src/parpeamici/multiConditionProblem.cpp index 0ef2ec37..57c3bb10 100644 --- a/src/parpeamici/multiConditionProblem.cpp +++ b/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