diff --git a/include/amici/rdata.h b/include/amici/rdata.h index 9945771dd7..145d9b378a 100644 --- a/include/amici/rdata.h +++ b/include/amici/rdata.h @@ -468,19 +468,6 @@ class ReturnData : public ModelDimensions { */ int status = AMICI_NOT_RUN; - /** - * @brief Number of state variables. - * - * (alias `nx_rdata`, kept for backward compatibility) - */ - int nx{0}; - - /** - * Number of state variables in the unaugmented system - * (alias nxtrue_rdata, kept for backward compatibility) - */ - int nxtrue{0}; - /** Number of parameters w.r.t. which sensitivities were requested */ int nplist{0}; diff --git a/include/amici/serialization.h b/include/amici/serialization.h index d2240266fc..6b322c53fe 100644 --- a/include/amici/serialization.h +++ b/include/amici/serialization.h @@ -208,8 +208,8 @@ void serialize( ) { ar& dynamic_cast(r); ar & r.id; - ar & r.nx; - ar & r.nxtrue; + ar & r.nx_rdata; + ar & r.nxtrue_rdata; ar & r.nplist; ar & r.nmaxevent; ar & r.nt; diff --git a/models/model_calvetti_py/main.cpp b/models/model_calvetti_py/main.cpp index ecdff85a46..28b10015d2 100644 --- a/models/model_calvetti_py/main.cpp +++ b/models/model_calvetti_py/main.cpp @@ -78,14 +78,14 @@ int main() { std::cout << "State sensitivities for timepoint " << rdata->ts[i_time] << std::endl; // nt x nplist x nx - for (int i_state = 0; i_state < rdata->nx; ++i_state) { + for (int i_state = 0; i_state < rdata->nx_rdata; ++i_state) { std::cout << "\td(" << state_ids[i_state] << ")/d(" << parameter_ids[model->plist(i_nplist)] << ") = "; // rdata->sx is a flat 3D array in row-major ordering std::cout << rdata->sx - [i_time * rdata->nplist * rdata->nx - + i_nplist * rdata->nx + i_state]; + [i_time * rdata->nplist * rdata->nx_rdata + + i_nplist * rdata->nx_rdata + i_state]; std::cout << std::endl; } diff --git a/models/model_dirac_py/main.cpp b/models/model_dirac_py/main.cpp index ecdff85a46..28b10015d2 100644 --- a/models/model_dirac_py/main.cpp +++ b/models/model_dirac_py/main.cpp @@ -78,14 +78,14 @@ int main() { std::cout << "State sensitivities for timepoint " << rdata->ts[i_time] << std::endl; // nt x nplist x nx - for (int i_state = 0; i_state < rdata->nx; ++i_state) { + for (int i_state = 0; i_state < rdata->nx_rdata; ++i_state) { std::cout << "\td(" << state_ids[i_state] << ")/d(" << parameter_ids[model->plist(i_nplist)] << ") = "; // rdata->sx is a flat 3D array in row-major ordering std::cout << rdata->sx - [i_time * rdata->nplist * rdata->nx - + i_nplist * rdata->nx + i_state]; + [i_time * rdata->nplist * rdata->nx_rdata + + i_nplist * rdata->nx_rdata + i_state]; std::cout << std::endl; } diff --git a/models/model_events_py/main.cpp b/models/model_events_py/main.cpp index ecdff85a46..28b10015d2 100644 --- a/models/model_events_py/main.cpp +++ b/models/model_events_py/main.cpp @@ -78,14 +78,14 @@ int main() { std::cout << "State sensitivities for timepoint " << rdata->ts[i_time] << std::endl; // nt x nplist x nx - for (int i_state = 0; i_state < rdata->nx; ++i_state) { + for (int i_state = 0; i_state < rdata->nx_rdata; ++i_state) { std::cout << "\td(" << state_ids[i_state] << ")/d(" << parameter_ids[model->plist(i_nplist)] << ") = "; // rdata->sx is a flat 3D array in row-major ordering std::cout << rdata->sx - [i_time * rdata->nplist * rdata->nx - + i_nplist * rdata->nx + i_state]; + [i_time * rdata->nplist * rdata->nx_rdata + + i_nplist * rdata->nx_rdata + i_state]; std::cout << std::endl; } diff --git a/models/model_jakstat_adjoint_py/main.cpp b/models/model_jakstat_adjoint_py/main.cpp index ecdff85a46..28b10015d2 100644 --- a/models/model_jakstat_adjoint_py/main.cpp +++ b/models/model_jakstat_adjoint_py/main.cpp @@ -78,14 +78,14 @@ int main() { std::cout << "State sensitivities for timepoint " << rdata->ts[i_time] << std::endl; // nt x nplist x nx - for (int i_state = 0; i_state < rdata->nx; ++i_state) { + for (int i_state = 0; i_state < rdata->nx_rdata; ++i_state) { std::cout << "\td(" << state_ids[i_state] << ")/d(" << parameter_ids[model->plist(i_nplist)] << ") = "; // rdata->sx is a flat 3D array in row-major ordering std::cout << rdata->sx - [i_time * rdata->nplist * rdata->nx - + i_nplist * rdata->nx + i_state]; + [i_time * rdata->nplist * rdata->nx_rdata + + i_nplist * rdata->nx_rdata + i_state]; std::cout << std::endl; } diff --git a/models/model_nested_events_py/main.cpp b/models/model_nested_events_py/main.cpp index ecdff85a46..28b10015d2 100644 --- a/models/model_nested_events_py/main.cpp +++ b/models/model_nested_events_py/main.cpp @@ -78,14 +78,14 @@ int main() { std::cout << "State sensitivities for timepoint " << rdata->ts[i_time] << std::endl; // nt x nplist x nx - for (int i_state = 0; i_state < rdata->nx; ++i_state) { + for (int i_state = 0; i_state < rdata->nx_rdata; ++i_state) { std::cout << "\td(" << state_ids[i_state] << ")/d(" << parameter_ids[model->plist(i_nplist)] << ") = "; // rdata->sx is a flat 3D array in row-major ordering std::cout << rdata->sx - [i_time * rdata->nplist * rdata->nx - + i_nplist * rdata->nx + i_state]; + [i_time * rdata->nplist * rdata->nx_rdata + + i_nplist * rdata->nx_rdata + i_state]; std::cout << std::endl; } diff --git a/models/model_neuron_py/main.cpp b/models/model_neuron_py/main.cpp index ecdff85a46..28b10015d2 100644 --- a/models/model_neuron_py/main.cpp +++ b/models/model_neuron_py/main.cpp @@ -78,14 +78,14 @@ int main() { std::cout << "State sensitivities for timepoint " << rdata->ts[i_time] << std::endl; // nt x nplist x nx - for (int i_state = 0; i_state < rdata->nx; ++i_state) { + for (int i_state = 0; i_state < rdata->nx_rdata; ++i_state) { std::cout << "\td(" << state_ids[i_state] << ")/d(" << parameter_ids[model->plist(i_nplist)] << ") = "; // rdata->sx is a flat 3D array in row-major ordering std::cout << rdata->sx - [i_time * rdata->nplist * rdata->nx - + i_nplist * rdata->nx + i_state]; + [i_time * rdata->nplist * rdata->nx_rdata + + i_nplist * rdata->nx_rdata + i_state]; std::cout << std::endl; } diff --git a/models/model_robertson_py/main.cpp b/models/model_robertson_py/main.cpp index ecdff85a46..28b10015d2 100644 --- a/models/model_robertson_py/main.cpp +++ b/models/model_robertson_py/main.cpp @@ -78,14 +78,14 @@ int main() { std::cout << "State sensitivities for timepoint " << rdata->ts[i_time] << std::endl; // nt x nplist x nx - for (int i_state = 0; i_state < rdata->nx; ++i_state) { + for (int i_state = 0; i_state < rdata->nx_rdata; ++i_state) { std::cout << "\td(" << state_ids[i_state] << ")/d(" << parameter_ids[model->plist(i_nplist)] << ") = "; // rdata->sx is a flat 3D array in row-major ordering std::cout << rdata->sx - [i_time * rdata->nplist * rdata->nx - + i_nplist * rdata->nx + i_state]; + [i_time * rdata->nplist * rdata->nx_rdata + + i_nplist * rdata->nx_rdata + i_state]; std::cout << std::endl; } diff --git a/models/model_steadystate_py/main.cpp b/models/model_steadystate_py/main.cpp index ecdff85a46..28b10015d2 100644 --- a/models/model_steadystate_py/main.cpp +++ b/models/model_steadystate_py/main.cpp @@ -78,14 +78,14 @@ int main() { std::cout << "State sensitivities for timepoint " << rdata->ts[i_time] << std::endl; // nt x nplist x nx - for (int i_state = 0; i_state < rdata->nx; ++i_state) { + for (int i_state = 0; i_state < rdata->nx_rdata; ++i_state) { std::cout << "\td(" << state_ids[i_state] << ")/d(" << parameter_ids[model->plist(i_nplist)] << ") = "; // rdata->sx is a flat 3D array in row-major ordering std::cout << rdata->sx - [i_time * rdata->nplist * rdata->nx - + i_nplist * rdata->nx + i_state]; + [i_time * rdata->nplist * rdata->nx_rdata + + i_nplist * rdata->nx_rdata + i_state]; std::cout << std::endl; } diff --git a/python/sdist/amici/numpy.py b/python/sdist/amici/numpy.py index f7b85d8b96..0f61a1305a 100644 --- a/python/sdist/amici/numpy.py +++ b/python/sdist/amici/numpy.py @@ -265,12 +265,12 @@ def __init__(self, rdata: ReturnDataPtr | ReturnData): ) self._field_dimensions = { "ts": [rdata.nt], - "x": [rdata.nt, rdata.nx], - "x0": [rdata.nx], - "x_ss": [rdata.nx], - "sx": [rdata.nt, rdata.nplist, rdata.nx], - "sx0": [rdata.nplist, rdata.nx], - "sx_ss": [rdata.nplist, rdata.nx], + "x": [rdata.nt, rdata.nx_rdata], + "x0": [rdata.nx_rdata], + "x_ss": [rdata.nx_rdata], + "sx": [rdata.nt, rdata.nplist, rdata.nx_rdata], + "sx0": [rdata.nplist, rdata.nx_rdata], + "sx_ss": [rdata.nplist, rdata.nx_rdata], # observables "y": [rdata.nt, rdata.ny], "sigmay": [rdata.nt, rdata.ny], diff --git a/src/hdf5.cpp b/src/hdf5.cpp index 32154f0672..50c5997793 100644 --- a/src/hdf5.cpp +++ b/src/hdf5.cpp @@ -421,7 +421,7 @@ void writeReturnData( if (!rdata.x.empty()) createAndWriteDouble2DDataset( - file, hdf5Location + "/x", rdata.x, rdata.nt, rdata.nx + file, hdf5Location + "/x", rdata.x, rdata.nt, rdata.nx_rdata ); if (!rdata.y.empty()) @@ -463,13 +463,13 @@ void writeReturnData( if (!rdata.sx0.empty()) createAndWriteDouble2DDataset( - file, hdf5Location + "/sx0", rdata.sx0, rdata.nplist, rdata.nx + file, hdf5Location + "/sx0", rdata.sx0, rdata.nplist, rdata.nx_rdata ); if (!rdata.sx.empty()) createAndWriteDouble3DDataset( file, hdf5Location + "/sx", rdata.sx, rdata.nt, rdata.nplist, - rdata.nx + rdata.nx_rdata ); if (!rdata.sy.empty()) @@ -699,7 +699,7 @@ void writeReturnDataDiagnosis( if (!rdata.J.empty()) createAndWriteDouble2DDataset( - file, hdf5Location + "/J", rdata.J, rdata.nx, rdata.nx + file, hdf5Location + "/J", rdata.J, rdata.nx_rdata, rdata.nx_rdata ); if (!rdata.x_ss.empty()) diff --git a/src/main.template.cpp b/src/main.template.cpp index ecdff85a46..28b10015d2 100644 --- a/src/main.template.cpp +++ b/src/main.template.cpp @@ -78,14 +78,14 @@ int main() { std::cout << "State sensitivities for timepoint " << rdata->ts[i_time] << std::endl; // nt x nplist x nx - for (int i_state = 0; i_state < rdata->nx; ++i_state) { + for (int i_state = 0; i_state < rdata->nx_rdata; ++i_state) { std::cout << "\td(" << state_ids[i_state] << ")/d(" << parameter_ids[model->plist(i_nplist)] << ") = "; // rdata->sx is a flat 3D array in row-major ordering std::cout << rdata->sx - [i_time * rdata->nplist * rdata->nx - + i_nplist * rdata->nx + i_state]; + [i_time * rdata->nplist * rdata->nx_rdata + + i_nplist * rdata->nx_rdata + i_state]; std::cout << std::endl; } diff --git a/src/rdata.cpp b/src/rdata.cpp index 5d88611e11..6fac090c09 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -37,8 +37,6 @@ ReturnData::ReturnData( ) : ModelDimensions(model_dimensions_) , ts(std::move(ts_)) - , nx(nx_rdata) - , nxtrue(nxtrue_rdata) , nplist(plist_.size()) , nmaxevent(nmaxevent_) , nt(ts.size()) @@ -135,7 +133,7 @@ void ReturnData::initializeFullReporting(bool quadratic_llh) { sigmaz.resize(nmaxevent * nz, 0.0); rz.resize(nmaxevent * nz, 0.0); - x.resize(nt * nx, 0.0); + x.resize(nt * nx_rdata, 0.0); w.resize(nt * nw, 0.0); preeq_numsteps.resize(3, 0); @@ -159,17 +157,17 @@ void ReturnData::initializeFullReporting(bool quadratic_llh) { } } - x0.resize(nx, getNaN()); - x_ss.resize(nx, getNaN()); + x0.resize(nx_rdata, getNaN()); + x_ss.resize(nx_rdata, getNaN()); if (sensi >= SensitivityOrder::first) { - sx0.resize(nx * nplist, getNaN()); - sx_ss.resize(nx * nplist, getNaN()); + sx0.resize(nx_rdata * nplist, getNaN()); + sx_ss.resize(nx_rdata * nplist, getNaN()); if (sensi_meth == SensitivityMethod::forward || sensi >= SensitivityOrder::second) { // for second order we can fill in from the augmented states - sx.resize(nt * nx * nplist, 0.0); + sx.resize(nt * nx_rdata * nplist, 0.0); sz.resize(nmaxevent * nz * nplist, 0.0); srz.resize(nmaxevent * nz * nplist, 0.0); } @@ -339,7 +337,7 @@ void ReturnData::getDataOutput( int it, Model& model, SolutionState const& sol, ExpData const* edata ) { if (!x.empty()) { - model.fx_rdata(slice(x, it, nx), sol.x); + model.fx_rdata(slice(x, it, nx_rdata), sol.x); } if (!w.empty()) model.getExpression(slice(w, it, nw), ts[it], sol.x); @@ -378,7 +376,7 @@ void ReturnData::getDataSensisFSA( ) { if (!sx.empty()) { model.fsx_rdata( - gsl::span(&sx.at(it * nplist * nx), nplist * nx), sol.sx, + gsl::span(&sx.at(it * nplist * nx_rdata), nplist * nx_rdata), sol.sx, sol.x ); } @@ -652,14 +650,14 @@ void ReturnData::invalidate(int const it_start) { invalidateSLLH(); if (!x.empty()) - std::fill(x.begin() + nx * it_start, x.end(), getNaN()); + std::fill(x.begin() + nx_rdata * it_start, x.end(), getNaN()); if (!y.empty()) std::fill(y.begin() + ny * it_start, y.end(), getNaN()); if (!w.empty()) std::fill(w.begin() + nw * it_start, w.end(), getNaN()); if (!sx.empty()) - std::fill(sx.begin() + nx * nplist * it_start, sx.end(), getNaN()); + std::fill(sx.begin() + nx_rdata * nplist * it_start, sx.end(), getNaN()); if (!sy.empty()) std::fill(sy.begin() + ny * nplist * it_start, sy.end(), getNaN()); } @@ -729,10 +727,10 @@ void ReturnData::applyChainRuleFactorToSimulationResults(Model const& model) { && sensi_meth == SensitivityMethod::adjoint) { if (!sx.empty() && !x.empty()) for (int ip = 0; ip < nplist; ++ip) - for (int ix = 0; ix < nxtrue; ++ix) + for (int ix = 0; ix < nxtrue_rdata; ++ix) for (int it = 0; it < nt; ++it) - sx.at(ix + nxtrue * (ip + it * nplist)) - = x.at(it * nx + nxtrue + ip * nxtrue + ix); + sx.at(ix + nxtrue_rdata * (ip + it * nplist)) + = x.at(it * nx_rdata + nxtrue_rdata + ip * nxtrue_rdata + ix); if (!sy.empty() && !y.empty()) for (int ip = 0; ip < nplist; ++ip) @@ -789,14 +787,14 @@ void ReturnData::applyChainRuleFactorToSimulationResults(Model const& model) { val *= pcoefficient[ip]; }; - chain_rule(sx, nxtrue, nx, nt); + chain_rule(sx, nxtrue_rdata, nx_rdata, nt); chain_rule(sy, nytrue, ny, nt); chain_rule(ssigmay, nytrue, ny, nt); chain_rule(sz, nztrue, nz, nmaxevent); chain_rule(ssigmaz, nztrue, nz, nmaxevent); chain_rule(srz, nztrue, nz, nmaxevent); - chain_rule(sx0, nxtrue, nx, 1); - chain_rule(sx_ss, nxtrue, nx, 1); + chain_rule(sx0, nxtrue_rdata, nx_rdata, 1); + chain_rule(sx_ss, nxtrue_rdata, nx_rdata, 1); } if (o2mode == SecondOrderMode::full) { @@ -843,7 +841,7 @@ void ReturnData::applyChainRuleFactorToSimulationResults(Model const& model) { } }; - chain_rule(sx, nxtrue, nx, nt); + chain_rule(sx, nxtrue_rdata, nx_rdata, nt); chain_rule(sy, nytrue, ny, nt); chain_rule(ssigmay, nytrue, ny, nt); chain_rule(sz, nztrue, nz, nmaxevent); @@ -883,7 +881,7 @@ void ReturnData::applyChainRuleFactorToSimulationResults(Model const& model) { } }; - chain_rule(sx, nxtrue, nx, nt); + chain_rule(sx, nxtrue_rdata, nx_rdata, nt); chain_rule(sy, nytrue, ny, nt); chain_rule(ssigmay, nytrue, ny, nt); chain_rule(sz, nztrue, nz, nmaxevent); diff --git a/tests/cpp/testfunctions.cpp b/tests/cpp/testfunctions.cpp index a75ec07113..220f89f56b 100644 --- a/tests/cpp/testfunctions.cpp +++ b/tests/cpp/testfunctions.cpp @@ -307,7 +307,7 @@ void verifyReturnDataSensitivities(H5::H5File const& file, std::string const& re expected = hdf5::getDoubleDataset3D(file, resultPath + "/sx", m, n, o); for(int ip = 0; ip < model->nplist(); ++ip) checkEqualArray(&expected[ip * model->nt() * model->nxtrue_rdata], - &rdata->sx[ip * model->nt() * rdata->nx], + &rdata->sx[ip * model->nt() * rdata->nx_rdata], model->nt() * model->nxtrue_rdata, atol, rtol, "sx"); } else { ASSERT_TRUE(rdata->sx.empty()); diff --git a/tests/cpp/unittests/testSerialization.cpp b/tests/cpp/unittests/testSerialization.cpp index d44e99611d..de4f0491dd 100644 --- a/tests/cpp/unittests/testSerialization.cpp +++ b/tests/cpp/unittests/testSerialization.cpp @@ -14,8 +14,8 @@ void checkReturnDataEqual( ASSERT_EQ(r.id, s.id); ASSERT_EQ(r.np, s.np); ASSERT_EQ(r.nk, s.nk); - ASSERT_EQ(r.nx, s.nx); - ASSERT_EQ(r.nxtrue, s.nxtrue); + ASSERT_EQ(r.nx_rdata, s.nx_rdata); + ASSERT_EQ(r.nxtrue_rdata, s.nxtrue_rdata); ASSERT_EQ(r.nx_solver, s.nx_solver); ASSERT_EQ(r.nx_solver_reinit, s.nx_solver_reinit); ASSERT_EQ(r.ny, s.ny);