diff --git a/@amidata/amidata.m b/@amidata/amidata.m index 03ee78f1aa..3ccc65a305 100644 --- 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 index bc49b23e7b..d9d930400c 100644 --- 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 index a87c70f70f..d1670197f0 100644 --- 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 index 969f6632fa..7e43a90f92 100644 --- 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 const& sigmay, std::vector const& mz, std::vector 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 fixedParameters; + /** condition-specific parameters for pre-equilibration of size Model::nk() or empty */ + std::vector fixedParametersPreequilibration; }; } // namespace amici -#endif /* _MY_EDATA */ +#endif /* AMICI_EDATA_H */ diff --git a/include/amici/model.h b/include/amici/model.h index 74f97aff58..000e64e0bf 100644 --- 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 index 9a62b87157..799e8850d5 100644 --- 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 index 2ff65781b4..452e900e0c 100644 --- 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 #include #include "amici/defines.h" diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index d81df60138..8777a452e1 100644 --- 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 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 index eed9119f0a..fc81b1b4c8 100644 --- 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(mxGetPr(dataPreeq), mxGetPr(dataPreeq) + m * n); + } + } + return edata; } diff --git a/src/rdata.cpp b/src/rdata.cpp index eee763a860..8c82365773 100644 --- 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