bugfixes (#403)

* fix(python) fix edataToNumPyArrays function, fixes #402

* fix(core) fix stupid indexing error, lets sneak this fix in this PR

* feature(core) add solver unittests

* fix(python) fix logX replacement

* fix(core) fix segfault for steadystateproblem when timepoints are only specified through edata

* fix(core) fix restoring timepoints/fixedParameter logic and reduce abuse of rdata, fixes #405

* feature(core) more informative error codes from newton solver, fixes #404

* fix(core) fix rdata timepoints logic

* fix(core) rethrow instead of throw copy

* fix(serialization) fix rdata.ts serialization

* fix(serialization) fix r.ts const cast and add message when compiling python module with debug flags

* fix(serialization) r.ts is actually a vector

* fix(doc) adapt documentation for NewtonFailure

include/amici/defines.h

@@ -25,6 +25,7 @@ constexpr double pi = 3.14159265358979323846;
 #define AMICI_CONV_FAILURE              -4
 #define AMICI_ILL_INPUT                -22
 #define AMICI_ERROR                    -99
+#define AMICI_NOT_IMPLEMENTED         -999
 #define AMICI_SUCCESS                    0
 #define AMICI_DATA_RETURN                1
 #define AMICI_ROOT_RETURN                2

include/amici/exception.h

@@ -30,8 +30,8 @@ namespace amici {
     public:
         AmiException(char const* fmt, ...) : std::exception() {
             /** constructor with printf style interface
-             * @param[in] fmt error message with printf format
-             * @param[in] ... printf formating variables
+             * @param fmt error message with printf format
+             * @param ... printf formating variables
              */
             va_list ap;
             va_start(ap, fmt);
@@ -117,23 +117,23 @@ namespace amici {
     class CvodeException : public AmiException  {
     public:
         /** constructor
-         * @param[in] error_code error code returned by cvode function
-         * @param[in] function cvode function name
+         * @param error_code error code returned by cvode function
+         * @param function cvode function name
          */
         CvodeException(const int error_code, const char *function) :
-        AmiException("Cvode routine %s failed with error code (%i)",function,error_code){}
+        AmiException("Cvode routine %s failed with error code %i",function,error_code){}
     };
 
     /** @brief ida exception handler class
      */
     class IDAException : public AmiException  {
     public:
         /** constructor
-         * @param[in] error_code error code returned by ida function
-         * @param[in] function ida function name
+         * @param error_code error code returned by ida function
+         * @param function ida function name
          */
         IDAException(const int error_code, const char *function) :
-        AmiException("IDA routine %s failed with error code (%i)",function,error_code){}
+        AmiException("IDA routine %s failed with error code %i",function,error_code){}
     };
 
     /** @brief integration failure exception for the forward problem
@@ -148,8 +148,8 @@ namespace amici {
         /** time of integration failure */
         realtype time;
         /** constructor
-         * @param[in] code error code returned by cvode/ida
-         * @param[in] t time of integration failure
+         * @param code error code returned by cvode/ida
+         * @param t time of integration failure
          */
         IntegrationFailure(int code, realtype t) :
         AmiException("AMICI failed to integrate the forward problem") {
@@ -170,8 +170,8 @@ namespace amici {
         /** time of integration failure */
         realtype time;
         /** constructor
-         * @param[in] code error code returned by cvode/ida
-         * @param[in] t time of integration failure
+         * @param code error code returned by cvode/ida
+         * @param t time of integration failure
          */
         IntegrationFailureB(int code, realtype t) :
         AmiException("AMICI failed to integrate the backward problem") {
@@ -188,7 +188,7 @@ namespace amici {
     class SetupFailure : public AmiException {
     public:
         /** constructor, simply calls AmiException constructor
-         * @param[in] msg
+         * @param msg
          */
         SetupFailure(const char *msg) : AmiException(msg) {}
     };
@@ -200,10 +200,16 @@ namespace amici {
      */
     class NewtonFailure : public AmiException {
     public:
+        /** error code returned by solver */
+        int error_code;
         /** constructor, simply calls AmiException constructor
-         * @param[in] msg
+         * @param function name of the function in which the error occured
+         * @param code error code
          */
-        NewtonFailure(const char *msg) : AmiException(msg) {}
+        NewtonFailure(int code, const char *function) :
+        AmiException("NewtonSolver routine %s failed with error code %i",function,code) {
+            error_code = code;
+        }
     };
 
 

include/amici/rdata.h

@@ -45,7 +45,7 @@ class ReturnData {
     applyChainRuleFactorToSimulationResults(const Model *model);
 
     /** timepoints (dimension: nt) */
-    std::vector<realtype> ts;
+    const std::vector<realtype> ts;
 
     /** time derivative (dimension: nx) */
     std::vector<realtype> xdot;

include/amici/serialization.h

@@ -115,7 +115,7 @@ void serialize(Archive &ar, amici::ReturnData &r, const unsigned int version) {
     ar &const_cast<amici::SensitivityOrder &>(r.sensi);
     ar &const_cast<amici::SensitivityMethod &>(r.sensi_meth);
 
-    ar &r.ts;
+    ar &const_cast<std::vector<realtype> &>(r.ts);
     ar &r.xdot;
     ar &r.J;
     ar &r.z & r.sigmaz;

include/amici/solver.h

@@ -1000,31 +1000,31 @@ class Solver {
     /**
      * updates solver tolerances according to the currently specified member variables
      */
-    void setTolerances();
+    void applyTolerances();
 
     /**
      * updates FSA solver tolerances according to the currently specified member variables
      */
-    void setTolerancesFSA();
+    void applyTolerancesFSA();
 
     /**
      * updates ASA solver tolerances according to the currently specified member variables
      *
      * @param which identifier of the backwards problem
      */
-    void setTolerancesASA(int which);
+    void applyTolerancesASA(int which);
 
     /**
      * updates ASA quadrature solver tolerances according to the currently specified member variables
      *
      * @param which identifier of the backwards problem
      */
-    void setQuadTolerancesASA(int which);
+    void applyQuadTolerancesASA(int which);
 
     /**
      * updates all senstivivity solver tolerances according to the currently specified member variables
      */
-    void setSensitivityTolerances();
+    void applySensitivityTolerances();
 
     /** pointer to solver memory block */
     std::unique_ptr<void, std::function<void(void *)>> solverMemory;

python/amici/__init__.py

@@ -116,7 +116,13 @@ def edataToNumPyArrays(edata):
     """
     npExpData = {'ptr': edata}
 
-    fieldNames = ['my', 'sigmay', 'mz', 'sigmaz','fixedParameters','fixedParametersPreequilibration']
+    fieldNames = ['observedData', 'observedDataStdDev', 'observedEvents', 'observedEventsStdDev', 'fixedParameters',
+                  'fixedParametersPreequilibration']
+
+    edata.observedData = edata.getObservedData()
+    edata.observedDataStdDev = edata.getObservedDataStdDev()
+    edata.observedEvents = edata.getObservedEvents()
+    edata.observedEventsStdDev = edata.getObservedEventsStdDev()
 
     for field in fieldNames:
         npExpData[field] = getExpDataFieldAsNumPyArray(edata, field)
@@ -159,21 +165,21 @@ def getReturnDataFieldAsNumPyArray(rdata, field):
         Field Data as numpy array with dimensions according to model dimensions in rdata
 
     Raises:
-        
+
     """
 
     fieldDimensions = {'ts':  [rdata.nt],
                        'x':   [rdata.nt, rdata.nx],
                        'x0':  [rdata.nx],
                        'sx':  [rdata.nt, rdata.nplist, rdata.nx],
                        'sx0': [rdata.nplist, rdata.nx],
-                       
+
                        # observables
                        'y':       [rdata.nt, rdata.ny],
                        'sigmay':  [rdata.nt, rdata.ny],
                        'sy':      [rdata.nt, rdata.nplist, rdata.ny],
                        'ssigmay': [rdata.nt, rdata.nplist, rdata.ny],
-                       
+
                        # event observables
                        'z':       [rdata.nmaxevent, rdata.nz],
                        'rz':      [rdata.nmaxevent, rdata.nz],
@@ -189,7 +195,7 @@ def getReturnDataFieldAsNumPyArray(rdata, field):
                        'res':   [rdata.nt * rdata.nytrue],
                        'sres':  [rdata.nt * rdata.nytrue, rdata.nplist],
                        'FIM':   [rdata.nplist, rdata.nplist],
-                       
+
                        # diagnosis
                        'J':                       [rdata.nx, rdata.nx],
                        'xdot':                    [rdata.nx],
@@ -207,7 +213,7 @@ def getReturnDataFieldAsNumPyArray(rdata, field):
                        }
     if field == 't':
         field = 'ts'
-        
+
     return fieldAsNumpy(fieldDimensions, field, rdata)
 
 def getExpDataFieldAsNumPyArray(edata, field):
@@ -225,12 +231,12 @@ def getExpDataFieldAsNumPyArray(edata, field):
     """
 
     fieldDimensions = {# observables
-                       'my': [edata.nt, edata.nytrue],
-                       'sigmay': [edata.nt, edata.nytrue],
+                       'observedData': [edata.nt(), edata.nytrue],
+                       'observedDataStdDev': [edata.nt(), edata.nytrue],
 
                        # event observables
-                       'mz': [edata.nmaxevent, edata.nztrue],
-                       'sigmaz': [edata.nmaxevent, edata.nztrue],
+                       'observedEvents': [edata.nmaxevent, edata.nztrue],
+                       'observedEventsStdDev': [edata.nmaxevent, edata.nztrue],
 
                        # fixed parameters
                        'fixedParameters': [edata.fixedParameters.size()],

python/amici/sbml_import.py

@@ -789,7 +789,8 @@ def computeModelEquationsObjectiveFunction(self, observables={}, sigmas={}):
         # add user-provided observables or make all species observable
         if(observables):
             # Replace logX(.) by log(X,.) since symengine cannot parse the former
-            observables = { key: re.sub(r'(^|\W)log(\d+)\(', r'\1log(\2, ', formula) for key, formula in observables.items() }
+            for observable in observables:
+                observables[observable]['formula'] = re.sub(r'(^|\W)log(\d+)\(', r'\1log(\2, ', observables[observable]['formula'])
             self.observables = sp.DenseMatrix([observables[observable]['formula'] for observable in observables])
             observableNames = [observables[observable]['name'] if 'name' in observables[observable].keys()
                                else 'y' + str(index)

python/sdist/setup.py

@@ -101,6 +101,7 @@
 
 # Enable coverage?
 if 'ENABLE_GCOV_COVERAGE' in os.environ and os.environ['ENABLE_GCOV_COVERAGE'] == 'TRUE':
+    print("ENABLE_GCOV_COVERAGE was set to TRUE. Building AMICI with debug and coverage symbols.")
     cxx_flags.extend(['-g', '-O0',  '--coverage'])
     amici_module_linker_flags.append('--coverage')
 

src/amici.cpp

@@ -58,11 +58,7 @@ msgIdAndTxtFp warnMsgIdAndTxt = &printWarnMsgIdAndTxt;
  */
 std::unique_ptr<ReturnData> runAmiciSimulation(Solver &solver, const ExpData *edata, Model &model) {
 
-    auto rdata = std::unique_ptr<ReturnData>(new ReturnData(solver,&model));
-
-    if (model.nx <= 0) {
-        return rdata;
-    }
+    std::unique_ptr<ReturnData> rdata;
 
     auto originalFixedParameters = model.getFixedParameters(); // to restore after simulation
     auto originalTimepoints = model.getTimepoints();
@@ -79,7 +75,15 @@ std::unique_ptr<ReturnData> runAmiciSimulation(Solver &solver, const ExpData *ed
             model.setTimepoints(edata->getTimepoints());
         }
     }
+
     try{
+        rdata = std::unique_ptr<ReturnData>(new ReturnData(solver,&model));
+        if (model.nx <= 0) {
+            model.setFixedParameters(originalFixedParameters);
+            model.setTimepoints(originalTimepoints);
+            return rdata;
+        }
+
         auto fwd = std::unique_ptr<ForwardProblem>(new ForwardProblem(rdata.get(),edata,&model,&solver));
         fwd->workForwardProblem();
 

src/backwardproblem.cpp

@@ -46,8 +46,8 @@ void BackwardProblem::workBackwardProblem() {
      * workForwardProblem should be called before this is function is called
      */
 
-    if (model->nx <= 0 || rdata->sensi < SensitivityOrder::first ||
-        rdata->sensi_meth != SensitivityMethod::adjoint || model->nplist() == 0) {
+    if (model->nx <= 0 || solver->getSensitivityOrder() < SensitivityOrder::first ||
+        solver->getSensitivityMethod() != SensitivityMethod::adjoint || model->nplist() == 0) {
         return;
     }
 

src/edata.cpp

@@ -125,7 +125,7 @@ void ExpData::setObservedData(const std::vector<realtype> &observedData, int iy)
         throw AmiException("Input observedData did not match dimensions nt (%i), was %i", nt(), observedData.size());
 
     for (int it = 0; it < nt(); ++it)
-        this->observedData.at(iy + it*nytrue) = observedData.at(iy);
+        this->observedData.at(iy + it*nytrue) = observedData.at(it);
 }
 
 bool ExpData::isSetObservedData(int it, int iy) const {
@@ -164,7 +164,7 @@ void ExpData::setObservedDataStdDev(const std::vector<realtype> &observedDataStd
     checkSigmaPositivity(observedDataStdDev, "observedDataStdDev");
 
     for (int it = 0; it < nt(); ++it)
-        this->observedDataStdDev.at(iy + it*nytrue) = observedDataStdDev.at(iy);
+        this->observedDataStdDev.at(iy + it*nytrue) = observedDataStdDev.at(it);
 }
 
 void ExpData::setObservedDataStdDev(const realtype stdDev, int iy) {
@@ -203,7 +203,7 @@ void ExpData::setObservedEvents(const std::vector<realtype> &observedEvents, int
     }
 
     for (int ie = 0; ie < nmaxevent; ++ie)
-        this->observedEvents.at(iz + ie*nztrue) = observedEvents.at(iz);
+        this->observedEvents.at(iz + ie*nztrue) = observedEvents.at(ie);
 }
 
 bool ExpData::isSetObservedEvents(int ie, int iz) const {
@@ -242,7 +242,7 @@ void ExpData::setObservedEventsStdDev(const std::vector<realtype> &observedEvent
     checkSigmaPositivity(observedEventsStdDev, "observedEventsStdDev");
 
     for (int ie = 0; ie < nmaxevent; ++ie)
-        this->observedEventsStdDev.at(iz + ie*nztrue) = observedEventsStdDev.at(iz);
+        this->observedEventsStdDev.at(iz + ie*nztrue) = observedEventsStdDev.at(ie);
 }
 
 void ExpData::setObservedEventsStdDev(const realtype stdDev, int iz) {