Use C++20

CMakeLists.txt

@@ -31,14 +31,14 @@ list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
 
 if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU" OR "${CMAKE_CXX_COMPILER_ID}"
                                                 STREQUAL "MINGW")
-  # require at least gcc 9.1 for proper C+17 support, e.g. std::reduce
-  if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 9.1)
-    message(FATAL_ERROR "GCC version must be at least 9.1!")
+  # require at least gcc 10.1 for std::ranges (C++20) support
+  if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 10.1)
+    message(FATAL_ERROR "GCC version must be at least 10.1!")
   endif()
 endif()
 
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 # Compiler flags

cmake/clang-tools.cmake

@@ -29,7 +29,7 @@ if(CLANG_TIDY)
     clang-tidy
     COMMAND
       sh -c
-      "${CLANG_TIDY} ${ALL_CXX_SOURCE_FILES} -- -std=c++17 -I${CMAKE_SOURCE_DIR}"
+      "${CLANG_TIDY} ${ALL_CXX_SOURCE_FILES} -- -std=c++20 -I${CMAKE_SOURCE_DIR}"
     WORKING_DIRECTORY ${CMAKE_SOURCE_DIR})
 else()
   message(STATUS "clang-tidy was not found")

doc/cpp_installation.rst

@@ -15,7 +15,7 @@ The following section describes building the AMICI C++ library:
 Prerequisites:
 
 * CBLAS compatible BLAS library
-* a C++17 compatible compiler
+* a C++20 compatible compiler
 * a C compiler
 * Optional:
 

include/amici/forwardproblem.h

@@ -456,8 +456,9 @@ class FinalStateStorer : public ContextManager {
                 auto final_time = fwd_->getFinalTime();
                 auto const timepoints = fwd_->model->getTimepoints();
                 if (!fwd_->timepoint_states_.count(final_time)
-                    && std::find(timepoints.cbegin(), timepoints.cend(), final_time)
-                           != timepoints.cend()) {
+                    && std::find(
+                           timepoints.cbegin(), timepoints.cend(), final_time
+                       ) != timepoints.cend()) {
                     fwd_->timepoint_states_[final_time] = fwd_->final_state_;
                 }
             } catch (std::exception const&) {
@@ -469,7 +470,7 @@ class FinalStateStorer : public ContextManager {
                 // `fwd_->{final_state_,timepoint_states_}` won't be set,
                 // and we assume that they are either not accessed anymore, or
                 // that there is appropriate error handling in place.
-                if(!std::uncaught_exceptions()) {
+                if (!std::uncaught_exceptions()) {
                     throw;
                 }
             }

matlab/@amimodel/compileAndLinkModel.m

@@ -34,9 +34,9 @@ function compileAndLinkModel(modelname, modelSourceFolder, coptim, debug, funs,
     end
 
     % compile flags
-    COPT = ['COPTIMFLAGS=''' coptim ' -DNDEBUG'' CXXFLAGS=''$CXXFLAGS -std=c++17'''];
+    COPT = ['COPTIMFLAGS=''' coptim ' -DNDEBUG'' CXXFLAGS=''$CXXFLAGS -std=c++20'''];
     if(debug)
-        DEBUG = ' -g CXXFLAGS=''$CXXFLAGS -Wall  -std=c++17 -Wno-unused-function -Wno-unused-variable'' ';
+        DEBUG = ' -g CXXFLAGS=''$CXXFLAGS -Wall  -std=c++20 -Wno-unused-function -Wno-unused-variable'' ';
         COPT = ''; % no optimization with debug flags!
     else
         DEBUG = '';

src/amici.cpp

@@ -237,12 +237,12 @@ std::unique_ptr<ReturnData> runAmiciSimulation(
     gsl_EnsuresDebug(rdata->posteq_cpu_timeB <= rdata->cpu_time_total);
     if (!posteq)
         gsl_EnsuresDebug(
-            std::is_sorted(rdata->numsteps.begin(), rdata->numsteps.end())
+            std::ranges::is_sorted(rdata->numsteps)
             || rdata->status != AMICI_SUCCESS
         );
     if (!preeq)
         gsl_EnsuresDebug(
-            std::is_sorted(rdata->numstepsB.begin(), rdata->numstepsB.end())
+            std::ranges::is_sorted(rdata->numstepsB)
             || rdata->status != AMICI_SUCCESS
         );
     rdata->messages = logger.items;

src/edata.cpp

@@ -129,7 +129,7 @@
 }
 
 void ExpData::setTimepoints(std::vector<realtype> const& ts) {
-    if (!std::is_sorted(ts.begin(), ts.end()))
+    if (!std::ranges::is_sorted(ts))
         throw AmiException(
             "Encountered non-monotonic timepoints, please order timepoints "
             "such that they are monotonically increasing!"
@@ -196,9 +196,7 @@
 
 void ExpData::setObservedDataStdDev(realtype const stdDev) {
     checkSigmaPositivity(stdDev, "stdDev");
-    std::fill(
-        observed_data_std_dev_.begin(), observed_data_std_dev_.end(), stdDev
-    );
+    std::ranges::fill(observed_data_std_dev_, stdDev);
 }
 
 void ExpData::setObservedDataStdDev(
@@ -292,9 +290,7 @@
 
 void ExpData::setObservedEventsStdDev(realtype const stdDev) {
     checkSigmaPositivity(stdDev, "stdDev");
-    std::fill(
-        observed_events_std_dev_.begin(), observed_events_std_dev_.end(), stdDev
-    );
+    std::ranges::fill(observed_events_std_dev_, stdDev);
 }
 
 void ExpData::setObservedEventsStdDev(
@@ -339,15 +335,10 @@
 }
 
 void ExpData::clear_observations() {
-    std::fill(observed_data_.begin(), observed_data_.end(), getNaN());
-    std::fill(
-        observed_data_std_dev_.begin(), observed_data_std_dev_.end(), getNaN()
-    );
-    std::fill(observed_events_.begin(), observed_events_.end(), getNaN());
-    std::fill(
-        observed_events_std_dev_.begin(), observed_events_std_dev_.end(),
-        getNaN()
-    );
+    std::ranges::fill(observed_data_, getNaN());
+    std::ranges::fill(observed_data_std_dev_, getNaN());
+    std::ranges::fill(observed_events_, getNaN());
+    std::ranges::fill(observed_events_std_dev_, getNaN());
 }
 
 void ExpData::applyDimensions() {

src/forwardproblem.cpp

@@ -87,9 +87,7 @@
     solver->setup(t0, model, x_, dx_, sx_, sdx_);
 
     if (model->ne
-        && std::any_of(roots_found_.begin(), roots_found_.end(), [](int rf) {
-               return rf == 1;
-           }))
+        && std::ranges::any_of(roots_found_, [](int rf) { return rf == 1; }))
         handleEvent(&t0, false, true);
 
     /* perform presimulation if necessary */
@@ -101,10 +99,9 @@
         t_ = model->t0();
         if (model->ne) {
             model->initEvents(x_, dx_, roots_found_);
-            if (std::any_of(
-                    roots_found_.begin(), roots_found_.end(),
-                    [](int rf) { return rf == 1; }
-                ))
+            if (std::ranges::any_of(roots_found_, [](int rf) {
+                    return rf == 1;
+                }))
                 handleEvent(&t0, false, true);
         }
     }
@@ -136,10 +133,10 @@
 
     // get list of trigger timepoints for fixed-time triggered events
     auto trigger_timepoints = model->get_trigger_timepoints();
-    auto it_trigger_timepoints = std::find_if(
-        trigger_timepoints.begin(), trigger_timepoints.end(),
-        [this](auto t) { return t > this->t_; }
-    );
+    auto it_trigger_timepoints
+        = std::ranges::find_if(trigger_timepoints, [this](auto t) {
+              return t > this->t_;
+          });
 
     /* loop over timepoints */
     for (it_ = 0; it_ < model->nt(); it_++) {
@@ -333,7 +330,7 @@
             }
         }
         if (initial_event) // t0 has no parameter dependency
-            std::fill(stau_.begin(), stau_.end(), 0.0);
+            std::ranges::fill(stau_, 0.0);
     } else if (solver->computingASA()) {
         /* store x to compute jump in discontinuity */
         x_disc_.push_back(x_);

src/hdf5.cpp

@@ -699,8 +699,8 @@ void writeReturnDataDiagnosis(
 // work-around for macos segfaults, use struct without std::string
 struct LogItemCStr {
     int severity;
-    const char* identifier;
-    const char* message;
+    char const* identifier;
+    char const* message;
 };
 
 void writeLogItemsToHDF5(
@@ -738,8 +738,7 @@ void writeLogItemsToHDF5(
         // ... therefore, as a workaround, we use a struct without std::string
         H5::CompType logItemType(sizeof(LogItemCStr));
         logItemType.insertMember(
-            "severity", HOFFSET(LogItemCStr, severity),
-            H5::PredType::NATIVE_INT
+            "severity", HOFFSET(LogItemCStr, severity), H5::PredType::NATIVE_INT
         );
         auto vlstr_type = H5::StrType(H5::PredType::C_S1, H5T_VARIABLE);
         logItemType.insertMember(

src/model.cpp

@@ -78,7 +78,7 @@
 };
 
 static void setNaNtoZero(std::vector<realtype>& vec) {
-    std::for_each(vec.begin(), vec.end(), [](double& val) {
+    std::ranges::for_each(vec, [](double& val) {
         if (std::isnan(val)) {
             val = 0.0;
         }
@@ -98,7 +98,7 @@
     std::vector<std::string> const& ids, std::vector<realtype> const& values,
     std::string const& id, char const* variable_name, char const* id_name
 ) {
-    auto it = std::find(ids.begin(), ids.end(), id);
+    auto it = std::ranges::find(ids, id);
     if (it != ids.end())
         return values.at(it - ids.begin());
 
@@ -121,7 +121,7 @@
     realtype value, std::string const& id, char const* variable_name,
     char const* id_name
 ) {
-    auto it = std::find(ids.begin(), ids.end(), id);
+    auto it = std::ranges::find(ids, id);
     if (it != ids.end())
         values.at(it - ids.begin()) = value;
     else
@@ -338,8 +338,8 @@
     std::vector<realtype> tmp_dvalues(allnodes, 0.0);
     std::vector<realtype> tmp_dslopes(allnodes, 0.0);
     for (int ip = 0; ip < nplist(); ip++) {
-        std::fill(tmp_dvalues.begin(), tmp_dvalues.end(), 0.0);
-        std::fill(tmp_dslopes.begin(), tmp_dslopes.end(), 0.0);
+        std::ranges::fill(tmp_dvalues, 0.0);
+        std::ranges::fill(tmp_dslopes, 0.0);
         fdspline_valuesdp(
             tmp_dvalues.data(), state_.unscaledParameters.data(),
             state_.fixedParameters.data(), plist(ip)
@@ -405,7 +405,7 @@
 ) {
     std::vector<realtype> rootvals(ne, 0.0);
     froot(simulation_parameters_.tstart_, x, dx, rootvals);
-    std::fill(roots_found.begin(), roots_found.end(), 0);
+    std::ranges::fill(roots_found, 0);
     for (int ie = 0; ie < ne; ie++) {
         if (rootvals.at(ie) < 0) {
             state_.h.at(ie) = 0.0;
@@ -796,7 +796,7 @@
 }
 
 void Model::setTimepoints(std::vector<realtype> const& ts) {
-    if (!std::is_sorted(ts.begin(), ts.end()))
+    if (!std::ranges::is_sorted(ts))
         throw AmiException("Encountered non-monotonic timepoints, please order"
                            " timepoints such that they are monotonically"
                            " increasing!");
@@ -813,9 +813,7 @@
 
 void Model::setStateIsNonNegative(std::vector<bool> const& nonNegative) {
     auto any_state_non_negative
-        = std::any_of(nonNegative.begin(), nonNegative.end(), [](bool x) {
-              return x;
-          });
+        = std::ranges::any_of(nonNegative, [](bool x) { return x; });
     if (nx_solver != nx_rdata) {
         if (any_state_non_negative)
             throw AmiException("Non-negative states are not supported with"
@@ -845,7 +843,7 @@
 
 void Model::setParameterList(std::vector<int> const& plist) {
     int np = this->np(); // cannot capture 'this' in lambda expression
-    if (std::any_of(plist.begin(), plist.end(), [&np](int idx) {
+    if (std::ranges::any_of(plist, [&np](int idx) {
             return idx < 0 || idx >= np;
         })) {
         throw AmiException("Indices in plist must be in [0..np]");
@@ -1112,7 +1110,7 @@
     std::vector<realtype> nllh(nJ, 0.0);
     for (int iyt = 0; iyt < nytrue; iyt++) {
         if (edata.isSetObservedData(it, iyt)) {
-            std::fill(nllh.begin(), nllh.end(), 0.0);
+            std::ranges::fill(nllh, 0.0);
             fJy(nllh.data(), iyt, state_.unscaledParameters.data(),
                 state_.fixedParameters.data(), derived_state_.y_.data(),
                 derived_state_.sigmay_.data(), edata.getObservedDataPtr(it));
@@ -1303,7 +1301,7 @@
     std::vector<realtype> nllh(nJ, 0.0);
     for (int iztrue = 0; iztrue < nztrue; iztrue++) {
         if (edata.isSetObservedEvents(nroots, iztrue)) {
-            std::fill(nllh.begin(), nllh.end(), 0.0);
+            std::ranges::fill(nllh, 0.0);
             fJz(nllh.data(), iztrue, state_.unscaledParameters.data(),
                 state_.fixedParameters.data(), derived_state_.z_.data(),
                 derived_state_.sigmaz_.data(),
@@ -1323,7 +1321,7 @@
     std::vector<realtype> nllh(nJ, 0.0);
     for (int iztrue = 0; iztrue < nztrue; iztrue++) {
         if (edata.isSetObservedEvents(nroots, iztrue)) {
-            std::fill(nllh.begin(), nllh.end(), 0.0);
+            std::ranges::fill(nllh, 0.0);
             fJrz(
                 nllh.data(), iztrue, state_.unscaledParameters.data(),
                 state_.fixedParameters.data(), derived_state_.rz_.data(),
@@ -1395,7 +1393,7 @@
     AmiVector const& x, AmiVectorArray const& sx
 ) {
 
-    std::fill(stau.begin(), stau.end(), 0.0);
+    std::ranges::fill(stau, 0.0);
 
     for (int ip = 0; ip < nplist(); ip++) {
         fstau(
@@ -1525,7 +1523,7 @@
 int Model::checkFinite(
     gsl::span<realtype const> array, ModelQuantity model_quantity, realtype t
 ) const {
-    auto it = std::find_if(array.begin(), array.end(), [](realtype x) {
+    auto it = std::ranges::find_if(array, [](realtype x) {
         return !std::isfinite(x);
     });
     if (it == array.end()) {
@@ -1624,7 +1622,7 @@
     gsl::span<realtype const> array, ModelQuantity model_quantity,
     size_t num_cols, realtype t
 ) const {
-    auto it = std::find_if(array.begin(), array.end(), [](realtype x) {
+    auto it = std::ranges::find_if(array, [](realtype x) {
         return !std::isfinite(x);
     });
     if (it == array.end()) {
@@ -1739,7 +1737,7 @@
     // check flat array, to see if there are any issues
     // (faster, in particular for sparse arrays)
     auto m_flat = gsl::make_span(m);
-    auto it = std::find_if(m_flat.begin(), m_flat.end(), [](realtype x) {
+    auto it = std::ranges::find_if(m_flat, [](realtype x) {
         return !std::isfinite(x);
     });
     if (it == m_flat.end()) {
@@ -1831,9 +1829,7 @@
 bool Model::getAlwaysCheckFinite() const { return always_check_finite_; }
 
 void Model::fx0(AmiVector& x) {
-    std::fill(
-        derived_state_.x_rdata_.begin(), derived_state_.x_rdata_.end(), 0.0
-    );
+    std::ranges::fill(derived_state_.x_rdata_, 0.0);
     /* this function  also computes initial total abundances */
     fx0(derived_state_.x_rdata_.data(), simulation_parameters_.tstart_,
         state_.unscaledParameters.data(), state_.fixedParameters.data());
@@ -1877,10 +1873,7 @@
     for (int ip = 0; ip < nplist(); ip++) {
         if (ncl() > 0)
             stcl = &state_.stotal_cl.at(plist(ip) * ncl());
-        std::fill(
-            derived_state_.sx_rdata_.begin(), derived_state_.sx_rdata_.end(),
-            0.0
-        );
+        std::ranges::fill(derived_state_.sx_rdata_, 0.0);
         fsx0(
             derived_state_.sx_rdata_.data(), simulation_parameters_.tstart_,
             computeX_pos(x), state_.unscaledParameters.data(),
@@ -2244,10 +2237,7 @@
             }
         }
     } else {
-        std::fill(
-            derived_state_.dJydy_matlab_.begin(),
-            derived_state_.dJydy_matlab_.end(), 0.0
-        );
+        std::ranges::fill(derived_state_.dJydy_matlab_, 0.0);
         for (int iyt = 0; iyt < nytrue; iyt++) {
             if (!edata.isSetObservedData(it, iyt))
                 continue;
@@ -2821,7 +2811,7 @@
 
 void Model::fw(realtype const t, realtype const* x, bool include_static) {
     if (include_static) {
-        std::fill(derived_state_.w_.begin(), derived_state_.w_.end(), 0.0);
+        std::ranges::fill(derived_state_.w_, 0.0);
     }
     fspl(t);
     fw(derived_state_.w_.data(), t, x, state_.unscaledParameters.data(),
@@ -3076,7 +3066,7 @@
     for (auto const& kv : state_independent_events_) {
         *(it++) = kv.first;
     }
-    std::sort(trigger_timepoints.begin(), trigger_timepoints.end());
+    std::ranges::sort(trigger_timepoints);
     return trigger_timepoints;
 }