Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support for forward sensitivities #393

Merged
merged 1 commit into from
Oct 18, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion doc/optimizationApplication.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
94 changes: 91 additions & 3 deletions examples/parpeamici/steadystate/parpeExampleSteadystateBasic.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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": {
Expand Down Expand Up @@ -1110,4 +1198,4 @@
},
"nbformat": 4,
"nbformat_minor": 2
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

16 changes: 12 additions & 4 deletions src/parpeamici/multiConditionProblem.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading