From cdbcddf80554d124323141ec7f3af3d36875f4df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Thu, 16 May 2024 10:55:06 +0200 Subject: [PATCH 1/4] add initialization identifier --- include/amici/edata.h | 8 +++++++- include/amici/rdata.h | 7 ++++++- src/edata.cpp | 1 + 3 files changed, 14 insertions(+), 2 deletions(-) diff --git a/include/amici/edata.h b/include/amici/edata.h index 5d613c6262..e135bedd9e 100644 --- a/include/amici/edata.h +++ b/include/amici/edata.h @@ -398,10 +398,15 @@ class ExpData : public SimulationParameters { void clear_observations(); /** - * @brief Arbitrary (not necessarily unique) identifier. + * @brief Unique identifier. */ std::string id; + /** + * @brief Unique Identifier of the ExpData to use for initialization. + */ + std::string initialization_id; + protected: /** * @brief resizes observedData, observedDataStdDev, observedEvents and @@ -480,6 +485,7 @@ inline bool operator==(ExpData const& lhs, ExpData const& rhs) { return *dynamic_cast(&lhs) == *dynamic_cast(&rhs) && lhs.id == rhs.id && lhs.nytrue_ == rhs.nytrue_ + && lhs.initialization_id == rhs.initialization_id && lhs.nztrue_ == rhs.nztrue_ && lhs.nmaxevent_ == rhs.nmaxevent_ && is_equal(lhs.observed_data_, rhs.observed_data_) && is_equal(lhs.observed_data_std_dev_, rhs.observed_data_std_dev_) diff --git a/include/amici/rdata.h b/include/amici/rdata.h index b5b4c269b4..fe7186888a 100644 --- a/include/amici/rdata.h +++ b/include/amici/rdata.h @@ -95,10 +95,15 @@ class ReturnData : public ModelDimensions { Model& model, Solver const& solver, ExpData const* edata ); /** - * @brief Arbitrary (not necessarily unique) identifier. + * @brief Unique identifier. */ std::string id; + /** + * @brief Unique Identifier of the ExpData to use for initialization. + */ + std::string initialization_id; + /** * timepoints (shape `nt`) */ diff --git a/src/edata.cpp b/src/edata.cpp index 2408e4becd..fee49d5717 100644 --- a/src/edata.cpp +++ b/src/edata.cpp @@ -126,6 +126,7 @@ ExpData::ExpData( } id = rdata.id; + initialization_id = rdata.initialization_id; } void ExpData::setTimepoints(std::vector const& ts) { From cc620e15061572ea3bccb910bf7eccf7df9635eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Thu, 16 May 2024 15:59:32 +0200 Subject: [PATCH 2/4] add timeseries class --- python/sdist/amici/swig_wrappers.py | 52 ++++++++++++++++++++++++----- 1 file changed, 43 insertions(+), 9 deletions(-) diff --git a/python/sdist/amici/swig_wrappers.py b/python/sdist/amici/swig_wrappers.py index 2e90891df3..95d1235d38 100644 --- a/python/sdist/amici/swig_wrappers.py +++ b/python/sdist/amici/swig_wrappers.py @@ -73,10 +73,35 @@ def runAmiciSimulation( return numpy.ReturnDataView(rdata) +class ExpDataTimeseries: + """ + Prototype for timeseries data storage. + """ + + _exp_datas: list[amici.ExpData] + + def __len__(self): + return len(self._exp_datas) + + def __init__(self, exp_datas: list[amici.ExpData]): + self._exp_datas = exp_datas + + def __getitem__(self, key): + return next(edata for edata in self._exp_datas if edata.id == key) + + def __setitem__(self, key, value): + self._exp_datas = [ + edata if edata.id != key else value for edata in self._exp_datas + ] + + def get_ptr_vector(self): + return amici_swig.ExpDataPtrVector(self._exp_datas) + + def runAmiciSimulations( model: AmiciModel, solver: AmiciSolver, - edata_list: AmiciExpDataVector, + edata_list: AmiciExpDataVector | ExpDataTimeseries, failfast: bool = True, num_threads: int = 1, ) -> list["numpy.ReturnDataView"]: @@ -105,14 +130,23 @@ def runAmiciSimulations( stacklevel=1, ) - edata_ptr_vector = amici_swig.ExpDataPtrVector(edata_list) - rdata_ptr_list = amici_swig.runAmiciSimulations( - _get_ptr(solver), - edata_ptr_vector, - _get_ptr(model), - failfast, - num_threads, - ) + if isinstance(edata_list, ExpDataTimeseries): + edata_ptr_vector = edata_list.get_ptr_vector() + rdata_ptr_list = amici_swig.runAmiciSimulationsTimeseries( + _get_ptr(solver), + edata_ptr_vector, + _get_ptr(model), + failfast, + ) + else: + edata_ptr_vector = amici_swig.ExpDataPtrVector(edata_list) + rdata_ptr_list = amici_swig.runAmiciSimulations( + _get_ptr(solver), + edata_ptr_vector, + _get_ptr(model), + failfast, + num_threads, + ) for rdata in rdata_ptr_list: _log_simulation(rdata) if solver.getReturnDataReportingMode() == amici.RDataReporting.full: From 601aa755a4f07a8a08397e79cc24fce70ee06205 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Fri, 17 May 2024 14:01:57 +0100 Subject: [PATCH 3/4] prototype implementation --- include/amici/backwardproblem.h | 6 ++++++ include/amici/edata.h | 6 ++++++ src/edata.cpp | 5 +++++ 3 files changed, 17 insertions(+) diff --git a/include/amici/backwardproblem.h b/include/amici/backwardproblem.h index a5d32e3c47..64395e00fe 100644 --- a/include/amici/backwardproblem.h +++ b/include/amici/backwardproblem.h @@ -70,6 +70,12 @@ class BackwardProblem { */ AmiVector const& getAdjointState() const { return xB_; } + /** + * @brief Accessor for dxB + * @return xB + */ + AmiVector const& getAdjointStateDerivative() const { return dxB_; } + /** * @brief Accessor for xQB * @return xQB diff --git a/include/amici/edata.h b/include/amici/edata.h index e135bedd9e..3d7433b1ec 100644 --- a/include/amici/edata.h +++ b/include/amici/edata.h @@ -445,6 +445,12 @@ class ExpData : public SimulationParameters { std::vector const& input, char const* fieldname ) const; + /** + * @brief check whether this instance implements only equilibration + * without simulation + */ + bool isEquilibration() const; + /** @brief number of observables */ int nytrue_{0}; diff --git a/src/edata.cpp b/src/edata.cpp index fee49d5717..3bb0a062fc 100644 --- a/src/edata.cpp +++ b/src/edata.cpp @@ -386,6 +386,11 @@ void ExpData::checkEventsDimension( ); } +bool ExpData::isEquilibration() const { + return ts_.size() == 1 && std::isinf(ts_[0]); +} + + void checkSigmaPositivity( std::vector const& sigmaVector, char const* vectorName ) { From cc035a6daa405d2c5b877ba502bf4f6e8a3de227 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Mon, 20 May 2024 14:36:56 +0100 Subject: [PATCH 4/4] prototype implementation (unfinished) --- include/amici/amici.h | 14 ++ src/amici.cpp | 376 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 390 insertions(+) diff --git a/include/amici/amici.h b/include/amici/amici.h index 93b8daed9b..6d7576c9bc 100644 --- a/include/amici/amici.h +++ b/include/amici/amici.h @@ -38,6 +38,20 @@ std::vector> runAmiciSimulations( Model const& model, bool failfast, int num_threads ); +/** + * @brief Same as runAmiciSimulation, but for a timeseries of ExpData instances. + * + * @param solver Solver instance + * @param edatas experimental data objects + * @param model model specification object + * @param failfast flag to allow early termination + * @return vector of pointers to return data objects + */ +std::vector> runAmiciSimulationsTimeseries( + Solver const& solver, std::vector const& edatas, + Model const& model, bool failfast +); + /** * @brief Get the string representation of the given simulation status code * (see ReturnData::status). diff --git a/src/amici.cpp b/src/amici.cpp index c74c88e3a5..ea476c2a59 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -301,6 +301,382 @@ std::vector> runAmiciSimulations( return results; } + +std::unique_ptr runAmiciSimulationPeriod( + Solver& solver, ExpData& const edata, Model& model, + &std::map(std::string, ForwardProblem) fwds, + &std::map(std::string, SteadystateProblem) peqs, + bool rethrow +) { + + // create a temporary logger instance for Solver and Model to capture + // messages from only this simulation + Logger logger; + solver.logger = &logger; + model.logger = &logger; + // prevent dangling pointer + auto _ = gsl::finally([&solver, &model] { + solver.logger = model.logger = nullptr; + }); + + CpuTimer cpu_timer; + solver.startTimer(); + + /* Applies condition-specific model settings and restores them when going + * out of scope */ + ConditionContext cc(&model, edata, FixedParameterContext::simulation); + + std::unique_ptr rdata + = std::make_unique(solver, model); + rdata->id = edata->id; + + try { + auto peq = peqs.find(edata.id); + auto fwd = fwds.find(edata.id); + // simulate + if (peq) { + // steadystate + peq.workSteadyStateProblem(solver, model, -1); + } else { + // forward + fwd.workForwardProblem(); + } + + rdata->status = AMICI_SUCCESS; + } catch (amici::IntegrationFailure const& ex) { + if (ex.error_code == AMICI_RHSFUNC_FAIL && solver.timeExceeded()) { + rdata->status = AMICI_MAX_TIME_EXCEEDED; + if (rethrow) + throw; + logger.log( + LogSeverity::error, "MAXTIME_EXCEEDED", + "AMICI forward simulation failed at t = %g: " + "Maximum time exceeded in forward solve.", + ex.time + ); + } else { + rdata->status = ex.error_code; + if (rethrow) + throw; + logger.log( + LogSeverity::error, "FORWARD_FAILURE", + "AMICI forward simulation failed at t = %g: %s", ex.time, + ex.what() + ); + } + } catch (amici::AmiException const& ex) { + rdata->status = AMICI_ERROR; + if (rethrow) + throw; + logger.log( + LogSeverity::error, "OTHER", "AMICI simulation failed: %s", + ex.what() + ); +#ifndef NDEBUG + logger.log( + LogSeverity::debug, "BACKTRACE", + "The previous error occurred at:\n%s", ex.getBacktrace() + ); +#endif + } catch (std::exception const& ex) { + rdata->status = AMICI_ERROR; + if (rethrow) + throw; + logger.log( + LogSeverity::error, "OTHER", "AMICI simulation failed: %s", + ex.what() + ); + } + + rdata->cpu_time = cpu_timer.elapsed_milliseconds(); + rdata->cpu_time_total = rdata->cpu_time; + + rdata->messages = logger.items; + return rdata; +} + +std::unique_ptr runAmiciSimulationPeriodB( + Solver& solver, ExpData& const edata, Model& model, + const &std::map(std::string, ForwardProblem) fwds, + const &std::map(std::string, SteadystateProblem) peqs, + &std::map(std::string, BackwardProblem) bwds, + bool rethrow +) { + Logger logger; + solver.logger = &logger; + model.logger = &logger; + // prevent dangling pointer + auto _ = gsl::finally([&solver, &model] { + solver.logger = model.logger = nullptr; + }); + + CpuTimer cpu_timerB; + solver.startTimer(); + + /* Applies condition-specific model settings and restores them when going + * out of scope */ + ConditionContext cc(&model, edata, FixedParameterContext::simulation); + + bool bwd_success = false; + + try { + if (solver.computingASA()) { + bwd = std::make_unique(*fwd, posteq.get()); + auto peq = peqs.find(edata.id); + auto fwd = fwds.find(edata.id); + + if (peq != peqs.end()) { + // initialize from preequilibration + peq->getAdjointUpdates(model, *edata); + peq->workSteadyStateBackwardProblem( + solver, model, bwd.get() + ); + + } else (fwd != fwds.end()) { + // initialize from simulation + fwd->getAdjointUpdates(model, *edata); + bwd->workBackwardProblem(); + } + } + bwd_success = true; + } catch (amici::IntegrationFailureB const& ex) { + if (ex.error_code == AMICI_RHSFUNC_FAIL && solver.timeExceeded()) { + rdata->status = AMICI_MAX_TIME_EXCEEDED; + if (rethrow) + throw; + logger.log( + LogSeverity::error, "MAXTIME_EXCEEDED", + "AMICI backward simulation failed when trying to solve until " + "t = %g: Maximum time exceeded in backward solve.", + ex.time + ); + + } else { + rdata->status = ex.error_code; + if (rethrow) + throw; + logger.log( + LogSeverity::error, "BACKWARD_FAILURE", + "AMICI backward simulation failed when trying to solve until t " + "= %g" + " (check debug logs for details): %s", + ex.time, ex.what() + ); + } + } catch (amici::AmiException const& ex) { + rdata->status = AMICI_ERROR; + if (rethrow) + throw; + logger.log( + LogSeverity::error, "OTHER", "AMICI simulation failed: %s", + ex.what() + ); +#ifndef NDEBUG + logger.log( + LogSeverity::debug, "BACKTRACE", + "The previous error occurred at:\n%s", ex.getBacktrace() + ); +#endif + } catch (std::exception const& ex) { + rdata->status = AMICI_ERROR; + if (rethrow) + throw; + logger.log( + LogSeverity::error, "OTHER", "AMICI simulation failed: %s", + ex.what() + ); + } + + try { + rdata->processSimulationObjects( + nullptr, fwd.get(), bwd_success ? bwd.get() : nullptr, + nullptr, model, solver, edata + ); + } catch (std::exception const& ex) { + rdata->status = AMICI_ERROR; + if (rethrow) + throw; + logger.log( + LogSeverity::error, "OTHER", "AMICI simulation failed: %s", + ex.what() + ); + } + + rdata->cpu_timeB = cpu_timerB.elapsed_milliseconds(); + rdata->cpu_time_total += rdata->cpu_timeB; + + gsl_EnsuresDebug(rdata->cpu_timeB <= rdata->cpu_time_total); + + gsl_EnsuresDebug( + std::is_sorted(rdata->numstepsB.begin(), rdata->numstepsB.end()) + || rdata->status != AMICI_SUCCESS + ); + + // TODO(ff): append rather than set + rdata->messages = logger.items; + return rdata; +} + +void initializePeriodFromPreviousPeriod( + Model& model, std::string initialization_id, + std::map(std::string, ForwardProblem) fwds, + std::map(std::string, SteadystateProblem) peqs, +){ + realtype t0; + amiVector x0; + amiVector dx0; + amiVector sx0; + amiVector sdx0; + // initialize from initialization_id + if (edata.initialization_id == "") + return + + auto peq = peqs.find(edata.initialization_id); + auto fwd = fwds.find(edata.initialization_id); + if (peq != peqs.end()) { + // initialize from preequilibration + t0 = edata.t_start_; + x0 = peq.getState(); + dx0 = peq.getStateDerivative(); + } else if (fwd != fwds.end()) { + // initialize from simulation + t0 = fwd.getTime(); + x0 = fwd.getState(); + dx0 = fwd.getStateDerivative(); + if solver.computingFSA() { + sx0 = fwd.getStateSensitivity(); + sdx0 = fwd.getStateDerivativeSensitivity(); + } + } else { + throw amici::AmiException("Invalid initialization condition %s in simulation condition %s", edata.initialization_id, edata.id); + } + // TODO(ff) do we need a context manager here? + if (edata.isEquilibration()) { + // steadystate + peqs[edata.id] = std::make_unique(solver, model); + peqs.initialize(x0, dx0, sx0, sdx0); + } else { + // forward + fwds[edata.id] = std::make_unique( + edata, &model, &solver, nullptr + ); + fwd.initialize(x0, dx0, sx0, sdx0); + } + // TODO(ff) necessary? + model.setT0(t0); +} + +void initializePeriodFromNextPeriods( + Model& model, std::string initialization_id, + const std::vector edatas, + const std::map(std::string, ForwardProblem) fwds, + const std::map(std::string, SteadystateProblem) peqs, +) { + realtype t0; + amiVector xB0; + amiVector dxB0; + amiVector xQB0; + // initialize from initialization_id + for (int i = 0; i < (int)edatas.size(); ++i) { + auto edata = edatas[i]; + if (edata.initialization_id != initialization_id) + continue; + + auto bwd = bwds.find(edata.initialization_id); + if (peq != peqs.end()) { + // initialize from preequilibration + xB0 += peq.getAdjointState(); + dxB0 += peq.getAdjointStateDerivative(); + xQB0 += peq.getAdjointStateQuadrature(); + } else if (fwd != fwds.end()) { + // initialize from simulation + t0 = fwd.getTime(); + x0 = fwd.getState(); + dx0 = fwd.getStateDerivative(); + if solver + .computingFSA() { + sx0 = fwd.getStateSensitivity(); + sdx0 = fwd.getStateDerivativeSensitivity(); + } + } else { + throw amici::AmiException( + "Invalid initialization condition %s in simulation condition %s", + edata.initialization_id, edata.id + ); + } + } + // TODO(ff) do we need a context manager here? + if (edata.isEquilibration()) { + // steadystate + peqs[edata.id] = std::make_unique(solver, model); + peqs.initialize(x0, dx0, sx0, sdx0); + } else { + // forward + fwds[edata.id] = std::make_unique( + edata, &model, &solver, nullptr + ); + fwd.initialize(x0, dx0, sx0, sdx0); + } + // TODO(ff) necessary? + model.setT0(t0); +} + +std::vector> runAmiciSimulationsTimeseries( + Solver const& solver, std::vector const& edatas, + Model const& model, bool failfast +) { + std::vector> results(edatas.size()); + // is set to true if one simulation fails and we should skip the rest. + // shared across threads. + + std::map(std::string, ForwardProblem) fwds(edatas.size()); + std::map(std::string, SteadystateProblem) peqs(edatas.size()); + std::map(std::string, BackwardProblem) bwds(edatas.size()); + + bool skipThrough = false; + + // forward pass + for (int i = 0; i < (int)edatas.size(); ++i) { + auto edata = edatas[i]; + auto mySolver = std::unique_ptr(solver.clone()); + auto myModel = std::unique_ptr(model.clone()); + + if (skipThrough) { + ConditionContext conditionContext(myModel.get(), edata); + results[i] + = std::unique_ptr(new ReturnData(solver, model)); + continue + } + + initializeModelFromPreviousPeriod(myModel, edata, fwds, peqs); + results[i] = runAmiciSimulationPeriod(mySolver, edata, myModel); + + skipThrough |= failfast && results[i]->status < 0; + } + + // backward pass + for (int i = ((int)edatas.size()) -1; i >= 0; --i) { + auto edata = edatas[i]; + auto mySolver = std::unique_ptr(solver.clone()); + auto myModel = std::unique_ptr(model.clone()); + + if (skipThrough) { + ConditionContext conditionContext(myModel.get(), edata); + results[i] + = std::unique_ptr(new ReturnData(solver, model)); + continue + } + + initializeModelFromPreviousPeriod(myModel, edata, fwds, peqs); + results[i] = runAmiciSimulationPeriod(mySolver, edata, myModel); + + skipThrough |= failfast && results[i]->status < 0; + + } + + return results; +} + std::string simulation_status_to_str(int status) { try { return simulation_status_to_str_map.at(status);