feature(core,matlab) Allow specification of pre-equilibration paramet…

…ers in addition to existing fixed parameters. Closes #249
AMICI-dev · May 15, 2018 · 563edb7 · 563edb7
1 parent 4bc7cd2
commit 563edb7
Show file tree

Hide file tree

Showing 10 changed files with 95 additions and 16 deletions.
diff --git a/@amidata/amidata.m b/@amidata/amidata.m
@@ -32,6 +32,8 @@
         Sigma_Z = double.empty();
         % experimental condition
         condition = double.empty();
+        % experimental condition for preequilibration
+        conditionPreequilibration = double.empty();
     end
 
     methods
@@ -50,7 +52,7 @@
             %                 Z [ne,nz]
             %                 Sigma_Z [ne,nz]
             %                 condition [nk,1]
-            %                 
+            %                 conditionPreequilibration [nk,1]
             %                 if some fields are missing the function will try
             %                 to initialise them with NaNs with consistent
             %                 dimensions
@@ -115,6 +117,10 @@
                 else
                     D.nk = 0;
                 end
+                if(isfield(varargin{1},'conditionPreequilibration'))
+                    assert(D.nk == numel(varargin{1}.conditionPreequilibration));
+                    D.conditionPreequilibration = varargin{1}.conditionPreequilibration;
+                end
             elseif(nargin == 5)
                 D.nt = varargin{1};
                 D.ny = varargin{2};
@@ -169,6 +175,13 @@
             this.condition = double(value(:));
         end
 
+        function set.conditionPreequilibration(this,value)
+            assert(isnumeric(value),'AMICI:amimodel:condition:numeric','condition must have a numeric value!')
+            assert(ismatrix(value),'AMICI:amimodel:condition:ndims','condition must be a two dimensional matrix!')
+            assert(numel(value)==this.nk,'AMICI:amimodel:condition:ndims',['condition must have ' num2str(this.nk) ' (D.nk) elements!'])
+            this.conditionPreequilibration = double(value(:));
+        end
+
         function set.Y(this,value)
             assert(ismatrix(value),'AMICI:amimodel:Y:ndims','Y must be a two dimensional matrix!')
             assert(all(all(or(isnumeric(value),isnan(value)))),'AMICI:amimodel:Y:numeric','Y must have a numeric value!')

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -92,6 +92,10 @@ target_link_libraries(${PROJECT_NAME}
     -ldl -lz -lm
 )
 
+# Create targets to make the sources show up in IDEs for convenience
+add_custom_target(matlabInterface SOURCES src/interface_matlab.cpp src/returndata_matlab.cpp)
+add_custom_target(fileTemplates SOURCES src/CMakeLists.template.txt src/main.template.cpp)
+
 if($ENV{ENABLE_GCOV_COVERAGE})
 	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 --coverage")
 endif()

diff --git a/examples/example_adjoint/example_adjoint.m b/examples/example_adjoint/example_adjoint.m
@@ -29,6 +29,7 @@ function example_adjoint()
 options.maxsteps = 1e4;
 options.rtol = 1e-12;
 options.atol = 1e-12;
+
 % load mex into memory
 [~] = which('simulate_model_adjoint'); % fix for inaccessability problems
 sol = simulate_model_adjoint(t,log10(p),k,D,options);

diff --git a/include/amici/edata.h b/include/amici/edata.h
@@ -9,12 +9,13 @@ namespace amici {
 
 class Model;
 
-/** @brief struct that carries all information about experimental data */
+/** @brief ExpData carries all information about experimental or condition-specific data */
 class ExpData {
 
   public:
     /** default constructor */
     ExpData();
+
     /**
      * @brief ExpData
      * @param nytrue
@@ -23,6 +24,7 @@ class ExpData {
      * @param nmaxevent
      */
     ExpData(int nytrue, int nztrue, int nt, int nmaxevent);
+
     /**
      * @brief ExpData
      * @param nytrue
@@ -39,6 +41,7 @@ class ExpData {
             std::vector<realtype> const& sigmay,
             std::vector<realtype> const& mz,
             std::vector<realtype> const& sigmaz);
+
     /**
      * constructor that initializes with Model
      *
@@ -78,8 +81,13 @@ class ExpData {
     const int nt;
     /** maximal number of event occurences */
     const int nmaxevent;
+
+    /** condition-specific parameters of size Model::nk() or empty */
+    std::vector<realtype> fixedParameters;
+    /** condition-specific parameters for pre-equilibration of size Model::nk() or empty */
+    std::vector<realtype> fixedParametersPreequilibration;
 };
 
 } // namespace amici
 
-#endif /* _MY_EDATA */
+#endif /* AMICI_EDATA_H */
diff --git a/include/amici/model.h b/include/amici/model.h
@@ -932,7 +932,7 @@ namespace amici {
             throw AmiException("Requested functionality is not supported as (%s) is not implemented for this model!",__func__);
         }
 
-        /** model specific implementation of fdzdp
+        /** model specific implementation of fdrzdp
          * @param drzdp partial derivative of root output rz w.r.t. model parameters p
          * @param ie event index
          * @param t current time
@@ -1013,7 +1013,7 @@ namespace amici {
             throw AmiException("Requested functionality is not supported as (%s) is not implemented for this model!",__func__);
         }
 
-        /** model specific implementation of fdeltasx
+        /** model specific implementation of fdeltaqB
          * @param deltaqB sensitivity update
          * @param t current time
          * @param x current state
@@ -1204,7 +1204,7 @@ namespace amici {
          * @param h heavyside vector
          */
         virtual void fw(realtype *w, const realtype t, const realtype *x, const realtype *p,
-                        const realtype *k, const realtype *h) {};
+                        const realtype *k, const realtype *h) {}
 
         /** model specific implementation of dwdp
          * @param dwdp Recurring terms in xdot, parameter derivative
@@ -1216,7 +1216,7 @@ namespace amici {
          * @param w vector with helper variables
          */
         virtual void fdwdp(realtype *dwdp, const realtype t, const realtype *x, const realtype *p,
-                           const realtype *k, const realtype *h, const realtype *w) {};
+                           const realtype *k, const realtype *h, const realtype *w) {}
 
         /** model specific implementation of dwdx
          * @param dwdx Recurring terms in xdot, state derivative
@@ -1228,7 +1228,7 @@ namespace amici {
          * @param w vector with helper variables
          */
         virtual void fdwdx(realtype *dwdx, const realtype t, const realtype *x, const realtype *p,
-                           const realtype *k, const realtype *h, const realtype *w) {};
+                           const realtype *k, const realtype *h, const realtype *w) {}
 
         void getmy(const int it, const ExpData *edata);
 

diff --git a/include/amici/rdata.h b/include/amici/rdata.h
@@ -31,6 +31,11 @@ class ReturnData {
 
     ReturnData(Solver const& solver, const Model *model);
 
+    /**
+     * @brief initializeObjectiveFunction
+     */
+    void initializeObjectiveFunction();
+
     void invalidate(const realtype t);
     void invalidateLLH();
 

diff --git a/models/model_steadystate/wrapfunctions.h b/models/model_steadystate/wrapfunctions.h
@@ -1,6 +1,6 @@
 #ifndef _amici_wrapfunctions_h
 #define _amici_wrapfunctions_h
-/* Generated by AMICI 7e4e208f6b779f504c7f6213008012ff6006f710 */
+/* Generated by AMICI 4bc7cd2a972300a9064305c7713378fae209b017 */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"

diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp
@@ -86,22 +86,35 @@ void ForwardProblem::workForwardProblem() {
     } catch (...) {
         throw AmiException("AMICI setup failed due to an unknown error");
     }
-    /* initialise objective function values */
+
     if(edata){
-        rdata->llh = 0.0;
-        rdata->chi2 = 0.0;
-        std::fill(rdata->sllh.begin(),rdata->sllh.end(), 0.0);
-        std::fill(rdata->s2llh.begin(),rdata->s2llh.end(), 0.0);
+        rdata->initializeObjectiveFunction();
     }
 
     int ncheck = 0; /* the number of (internal) checkpoints stored so far */
     realtype tlastroot = 0; /* storage for last found root */
 
     /* if preequilibration is necessary, start Newton solver */
+    std::vector<realtype> originalFixedParameters; // to restore after pre-equilibration
     if (solver->getNewtonPreequilibration()) {
+        if(edata && edata->fixedParametersPreequilibration.size()) {
+            // Are there dedicated preequilibration parameters provided?
+            if(edata->fixedParametersPreequilibration.size() != (unsigned) model->nk())
+                throw AmiException("Number of fixed parameters (%d) in model does not match preequilibration parameters in ExpData (%zd).",
+                                   model->nk(), edata->fixedParametersPreequilibration.size());
+            originalFixedParameters = model->getFixedParameters();
+            model->setFixedParameters(edata->fixedParametersPreequilibration);
+        } else if(edata && edata->fixedParameters.size()) {
+            // ... or other condition parameters?
+            if(edata->fixedParameters.size() != (unsigned) model->nk())
+                throw AmiException("Number of fixed parameters (%d) in model does not match ExpData (%zd).",
+                                   model->nk(), edata->fixedParameters.size());
+            model->setFixedParameters(edata->fixedParameters);
+        }
+
+        // pre-equilibrate
         SteadystateProblem sstate = SteadystateProblem(&t,&x,&sx);
-        sstate.workSteadyStateProblem(rdata,
-                                       solver, model, -1);
+        sstate.workSteadyStateProblem(rdata, solver, model, -1);
     } else {
         for (int ix = 0; ix < model->nx; ix++) {
             rdata->x0[ix] = x[ix];
@@ -113,6 +126,18 @@ void ForwardProblem::workForwardProblem() {
         }
     }
 
+    if(edata && edata->fixedParameters.size()) {
+        // fixed parameter in model are superseded by those provided in edata
+        // Note: this changes those parameters of `model` permanently
+        if(edata->fixedParameters.size() != (unsigned) model->nk())
+            throw AmiException("Number of fixed parameters (%d) in model does not match ExpData (%zd).",
+                               model->nk(), edata->fixedParameters.size());
+        model->setFixedParameters(edata->fixedParameters);
+    } else if (originalFixedParameters.size()) {
+        // Restore original parameters if only preequilibration parameters but no regular condition parameters have been provided
+        model->setFixedParameters(originalFixedParameters);
+    }
+
     /* loop over timepoints */
     for (int it = 0; it < rdata->nt; it++) {
         if (rdata->sensi_meth == AMICI_SENSI_FSA &&

diff --git a/src/interface_matlab.cpp b/src/interface_matlab.cpp
@@ -268,6 +268,21 @@ ExpData *expDataFromMatlabCall(const mxArray *prhs[],
         throw AmiException("Field Sigma_Z not specified as field in data struct!");
 
     }
+
+    // preequilibration condition parameters
+    if (mxArray *dataPreeq = mxGetProperty(prhs[RHS_DATA], 0, "conditionPreequilibration")) {
+        int m = (int)mxGetM(dataPreeq);
+        int n = (int)mxGetN(dataPreeq);
+        if(m * n > 0) {
+            if (m * n != model.nk() || (m != 1 && n != 1)) {
+                throw AmiException("Number of preequilibration parameters (%dx%d) does "
+                                   "not match model (%d)", m, n, model.nk());
+            }
+            edata->fixedParametersPreequilibration =
+                    std::vector<realtype>(mxGetPr(dataPreeq), mxGetPr(dataPreeq) + m * n);
+        }
+    }
+
     return edata;
 }
 

diff --git a/src/rdata.cpp b/src/rdata.cpp
@@ -304,4 +304,12 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) {
     return;
 }
 
+void ReturnData::initializeObjectiveFunction()
+{
+    llh = 0.0;
+    chi2 = 0.0;
+    std::fill(sllh.begin(),sllh.end(), 0.0);
+    std::fill(s2llh.begin(),s2llh.end(), 0.0);
+}
+
 } // namespace amici