diff --git a/.travis.yml b/.travis.yml index 4e66d1d315..e766f04ff6 100644 --- a/.travis.yml +++ b/.travis.yml @@ -23,10 +23,12 @@ matrix: - pyenv shell 2.7 3.6 compiler: gcc env: ENABLE_GCOV_COVERAGE=TRUE - after_success: + after_script: + # this can fail the build, after_success cannot - ./scripts/run-codecov.sh - bash <(curl -s https://codecov.io/bash) -f coverage.info -X fix -F cpp - bash <(curl -s https://codecov.io/bash) -f coverage_py.xml -F python + - os: osx osx_image: xcode9.3 @@ -58,6 +60,8 @@ install: - ./scripts/buildBNGL.sh - export BNGPATH=$(pwd)/ThirdParty/BioNetGen-2.3.2 - ./scripts/buildAmici.sh + - ./scripts/installAmiciArchive.sh + - ./scripts/installAmiciSource.sh script: # - cd $BASE_DIR/python/sdist diff --git a/include/amici/edata.h b/include/amici/edata.h index f0199337a0..a88903ac28 100644 --- a/include/amici/edata.h +++ b/include/amici/edata.h @@ -16,6 +16,9 @@ class ExpData { public: /** default constructor */ ExpData(); + + /** default copy constructor, needs to be declared to be generated in swig*/ + ExpData(const ExpData &) = default; /** * constructor that only initializes dimensions diff --git a/python/amici/__init__.py b/python/amici/__init__.py index 27a0d7f036..6fa1cf8ecf 100644 --- a/python/amici/__init__.py +++ b/python/amici/__init__.py @@ -35,7 +35,7 @@ from .amici import * hdf5_enabled = True has_clibs = True -except (ImportError, ModuleNotFoundError, AttributeError): +except (ImportError, ModuleNotFoundError, AttributeError): # pragma: no cover try: from . import amici_without_hdf5 as amici from .amici_without_hdf5 import * @@ -92,13 +92,11 @@ def runAmiciSimulation(model, solver, edata=None): return rdataToNumPyArrays(rdata) -def ExpData(rdata, sigma_y, sigma_z): - """ Convenience wrapper for ExpData constructor +def ExpData(*args): + """ Convenience wrapper for ExpData constructors Arguments: - rdata: rdataToNumPyArrays output - sigma_y: standard deviation for ObservableData - sigma_z: standard deviation for EventData + args: arguments Returns: ExpData Instance @@ -106,7 +104,16 @@ def ExpData(rdata, sigma_y, sigma_z): Raises: """ - return amici.ExpData(rdata['ptr'].get(), sigma_y, sigma_z) + if isinstance(args[0], dict): + # rdata dict created by rdataToNumPyArrays + return amici.ExpData(args[0]['ptr'].get(), *args[1:]) + elif isinstance(args[0], ExpDataPtr): + # the *args[:1] should be empty, but by the time you read this, + # the constructor signature may have changed and you are glad this + # wrapper did not break. + return amici.ExpData(args[0].get(), *args[:1]) + else: + return amici.ExpData(*args) def runAmiciSimulations(model, solver, edata_list, num_threads=1): diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index 4885eab55c..3461dfa8cb 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -1073,9 +1073,31 @@ def _computeEquation(self, name): [comp._dt for comp in self._states] ) - elif name in ['sx0', 'sx0_fixedParameters']: + elif name == 'sx0': self._derivative(name[1:], 'p', name=name) + elif name == 'sx0_fixedParameters': + # deltax = -x+x0_fixedParameters if x0_fixedParameters>0 else 0 + # deltasx = -sx+dx0_fixedParametersdx*sx+dx0_fixedParametersdp + # if x0_fixedParameters>0 else 0 + # sx0_fixedParameters = sx+deltasx = + # dx0_fixedParametersdx*sx+dx0_fixedParametersdp + self._eqs[name] = \ + self.eq('x0_fixedParameters').jacobian(self.sym('p')) + + for ip in range(self._eqs[name].shape[1]): + self._eqs[name][:,ip] += \ + self.eq('x0_fixedParameters').jacobian(self.sym('x')) \ + * self.sym('sx0') \ + + for index, formula in enumerate( + self.eq('x0_fixedParameters') + ): + if formula == 0 or formula == 0.0: + self._eqs[name][index, :] = \ + sp.zeros(1, self._eqs[name].shape[1]) + + elif name == 'JB': self._eqs[name] = self.eq('J').transpose() @@ -1155,7 +1177,11 @@ def _derivative(self, eq, var, name=None): # partial derivative if self.eq(eq).size and self.sym(var).size: - self._eqs[name] = self.eq(eq).jacobian(self.sym(var)) + if eq == 'Jy': + eq = self.eq(eq).transpose() + else: + eq = self.eq(eq) + self._eqs[name] = eq.jacobian(self.sym(var)) else: self._eqs[name] = sp.DenseMatrix([]) @@ -1643,7 +1669,8 @@ def _getFunctionBody(self, function, symbol): return lines if function == 'sx0_fixedParameters': - # here we specifically want to overwrite some of the values with 0 + # here we only want to overwrite values where x0_fixedParameters + # was applied lines.append(' ' * 4 + 'switch(ip) {') for ipar in range(self.model.np()): lines.append(' ' * 8 + f'case {ipar}:') @@ -1651,7 +1678,8 @@ def _getFunctionBody(self, function, symbol): self.model.eq('x0_fixedParameters') ): if formula != 0 and formula != 0.0: - lines.append(' ' * 12 + f'{function}[{index}] = 0.0;') + lines.append(' ' * 12 + f'{function}[{index}] = ' + f'{symbol[index, ipar]};') lines.append(' ' * 12 + 'break;') lines.append('}') elif function in sensi_functions: @@ -1724,11 +1752,8 @@ def _writeModelHeader(self): Raises: """ - if any([math != 0 and math != 0.0 for math in - self.model.eq('sx0_fixedParameters')]): - self.allow_reinit_fixpar_initcond = False - else: - self.allow_reinit_fixpar_initcond = True + + self.allow_reinit_fixpar_initcond = True templateData = { 'MODELNAME': str(self.modelName), diff --git a/scripts/buildAmici.sh b/scripts/buildAmici.sh index 7d2fb27b5c..19a5ca89c3 100755 --- a/scripts/buildAmici.sh +++ b/scripts/buildAmici.sh @@ -13,13 +13,5 @@ CPPUTEST_BUILD_DIR=${AMICI_PATH}/ThirdParty/cpputest-master/build/ CppUTest_DIR=${CPPUTEST_BUILD_DIR} cmake -DCMAKE_BUILD_TYPE=Debug .. make -# Disabled until cmake package is made compatible with updated setup.py -#make python-wheel -#pip3 install --user --prefix= `ls -t ${AMICI_PATH}/build/python/amici-*.whl | head -1` - make python-sdist set -x -python3 -m venv ${AMICI_PATH}/build/venv --clear -source ${AMICI_PATH}/build/venv/bin/activate -pip3 install --upgrade pip setuptools pkgconfig wheel numpy scipy matplotlib pysb -pip3 install $(ls -t ${AMICI_PATH}/build/python/amici-*.tar.gz | head -1) diff --git a/scripts/installAmiciArchive.sh b/scripts/installAmiciArchive.sh new file mode 100755 index 0000000000..4028843e2b --- /dev/null +++ b/scripts/installAmiciArchive.sh @@ -0,0 +1,24 @@ +#!/bin/bash +# +# Build libamici +# +set -e + +SCRIPT_PATH=$(dirname $BASH_SOURCE) +AMICI_PATH=$(cd $SCRIPT_PATH/.. && pwd) + +# Disabled until cmake package is made compatible with updated setup.py +#make python-wheel +#pip3 install --user --prefix= `ls -t ${AMICI_PATH}/build/python/amici-*.whl | head -1` + +rm -f ${AMICI_PATH}/python/sdist/amici/*.cxx +rm -f ${AMICI_PATH}/python/sdist/amici/*.so +rm -f ${AMICI_PATH}/python/sdist/amici/amici.py +rm -f ${AMICI_PATH}/python/sdist/amici/amici_without_hdf5.py + +# test install from archive +python3 -m venv ${AMICI_PATH}/build/venvArchive --clear +source ${AMICI_PATH}/build/venvArchive/bin/activate +pip3 install --upgrade pip setuptools pkgconfig wheel numpy +pip3 install $(ls -t ${AMICI_PATH}/build/python/amici-*.tar.gz | head -1) +deactivate diff --git a/scripts/buildAmiciDev.sh b/scripts/installAmiciSource.sh similarity index 56% rename from scripts/buildAmiciDev.sh rename to scripts/installAmiciSource.sh index 53e9373e79..930a56dd46 100755 --- a/scripts/buildAmiciDev.sh +++ b/scripts/installAmiciSource.sh @@ -7,14 +7,18 @@ set -e SCRIPT_PATH=$(dirname $BASH_SOURCE) AMICI_PATH=$(cd $SCRIPT_PATH/.. && pwd) -cd ${AMICI_PATH}/build -python3 -m venv ${AMICI_PATH}/build/venvDev --clear -rm -f ${AMICI_PATH}/python/sdist/amici/*.so +# Disabled until cmake package is made compatible with updated setup.py +#make python-wheel +#pip3 install --user --prefix= `ls -t ${AMICI_PATH}/build/python/amici-*.whl | head -1` + rm -f ${AMICI_PATH}/python/sdist/amici/*.cxx rm -f ${AMICI_PATH}/python/sdist/amici/*.so rm -f ${AMICI_PATH}/python/sdist/amici/amici.py rm -f ${AMICI_PATH}/python/sdist/amici/amici_without_hdf5.py -source ${AMICI_PATH}/build/venvDev/bin/activate -pip3 install --upgrade pip setuptools pkgconfig wheel numpy scipy matplotlib jupyter pysb -export ENABLE_AMICI_DEBUGGING=TRUE + +# test install from setup.py +python3 -m venv ${AMICI_PATH}/build/venv --clear +source ${AMICI_PATH}/build/venv/bin/activate +pip3 install --upgrade pip setuptools pkgconfig wheel numpy scipy matplotlib pysb pip3 install --verbose -e ${AMICI_PATH}/python/sdist +deactivate diff --git a/src/model_header.ODE_template.h b/src/model_header.ODE_template.h index 8a34577d57..0e168c9b24 100644 --- a/src/model_header.ODE_template.h +++ b/src/model_header.ODE_template.h @@ -38,6 +38,7 @@ extern void w_TPL_MODELNAME(realtype *w, const realtype t, const realtype *x, co extern void x0_TPL_MODELNAME(realtype *x0, const realtype t, const realtype *p, const realtype *k); extern void x0_fixedParameters_TPL_MODELNAME(realtype *x0, const realtype t, const realtype *p, const realtype *k); extern void sx0_TPL_MODELNAME(realtype *sx0, const realtype t,const realtype *x0, const realtype *p, const realtype *k, const int ip); +extern void sx0_fixedParameters_TPL_MODELNAME(realtype *sx0, const realtype t,const realtype *x0, const realtype *p, const realtype *k, const int ip); extern void xBdot_TPL_MODELNAME(realtype *xBdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *xB, const realtype *w, const realtype *dwdx); extern void xdot_TPL_MODELNAME(realtype *xdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w); extern void y_TPL_MODELNAME(double *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w); @@ -594,6 +595,18 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { sx0_TPL_MODELNAME(sx0, t, x0, p, k, ip); } + /** model specific implementation of fsx0_fixedParameters + * @param sx0 initial state sensitivities + * @param t initial time + * @param x0 initial state + * @param p parameter vector + * @param k constant vector + * @param ip sensitivity index + **/ + virtual void fsx0_fixedParameters(realtype *sx0, const realtype t,const realtype *x0, const realtype *p, const realtype *k, const int ip) override { + sx0_fixedParameters_TPL_MODELNAME(sx0, t, x0, p, k, ip); + } + /** model specific implementation of fsxdot * @param sxdot sensitivity rhs * @param t timepoint diff --git a/tests/testCoverage.py b/tests/testCoverage.py index 5d214d68b2..b1533ca2da 100755 --- a/tests/testCoverage.py +++ b/tests/testCoverage.py @@ -4,8 +4,22 @@ Generate coverage reports for the testModels and testSBML scripts exported format is cobertura xml """ - import coverage +# only consider amici module and ignore the swig generated amici.py +cov = coverage.Coverage( + source=['amici'], + omit=['*/amici.py','*/amici_without_hdf5.py'] +) + +# ignore code blocks containing import statements +cov.exclude('import') +cov.set_option('report:exclude_lines', ['pragma: no cover', ' raise ', + 'except:']) + +# coverage needs to be started before all other imports such that everything +# that happens during module import is also added to the report +cov.start() + import unittest import sys @@ -14,13 +28,6 @@ import testPandas import testPYSB -# only consider amici module and ignore the swig generated amici.py -cov = coverage.Coverage(source=['amici'],omit=['*/amici.py','*/amici_without_hdf5.py']) - -# ignore code blocks containing import statements -cov.exclude('import') -cov.start() - # build the testSuite from testModels and testSBML suite = unittest.TestSuite() suite.addTest(testSBML.TestAmiciSBMLModel()) diff --git a/tests/testModels.py b/tests/testModels.py index decc4651f4..0cd29baa88 100755 --- a/tests/testModels.py +++ b/tests/testModels.py @@ -8,6 +8,7 @@ import os import re import numpy as np +import copy class TestAmiciPregeneratedModel(unittest.TestCase): ''' @@ -55,11 +56,27 @@ def runTest(self): "/" + subTest + "/" + case + "/data", self.model.get()) rdata = amici.runAmiciSimulation(self.model, self.solver, edata) + + # todo: set higher tolerances in testcase options and + # regenerate results + if modelName == 'model_jakstat_adjoint_o2': + self.solver.setRelativeTolerance(1e-10) + self.solver.setAbsoluteTolerance(1e-10) + + if modelName in [ + 'model_jakstat_adjoint', 'model_nested_events', + 'model_steadystate' + ]: + self.solver.setRelativeTolerance(1e-12) + self.solver.setAbsoluteTolerance(1e-12) - if edata and self.solver.getSensitivityMethod() and self.solver.getSensitivityOrder(): - if not modelName.startswith('model_neuron'): - if(self.solver.getSensitivityOrder() and len(self.model.getParameterList())): - checkDerivatives(self.model, self.solver, edata) + if edata \ + and self.solver.getSensitivityMethod() \ + and self.solver.getSensitivityOrder() \ + and len(self.model.getParameterList()) \ + and not modelName.startswith('model_neuron') \ + and not case.endswith('byhandpreeq'): + checkDerivatives(self.model, self.solver, edata) if modelName == 'model_neuron_o2': self.solver.setRelativeTolerance(1e-12) @@ -81,10 +98,7 @@ def runTest(self): ) - - - -def checkDerivatives(model, solver, edata): +def checkDerivatives(model, solver, edata, atol=1e-4, rtol=1e-4, epsilon=1e-5): """Finite differences check for likelihood gradient Arguments: @@ -98,10 +112,10 @@ def func(x0, symbol='llh', x0full=None, plist=None, verbose=False): """Function of which the gradient is to be checked""" if plist is None: plist = [] - p = x0 + p = copy.deepcopy(x0) if len(plist): - p = x0full[:] - p[plist] = x0 + p = copy.deepcopy(x0full) + p[plist] = copy.deepcopy(x0) verbose and print('f: p=%s' % p) old_sensitivity_order = solver.getSensitivityOrder() @@ -124,18 +138,19 @@ def grad(x0, symbol='llh', x0full=None, plist=None, verbose=False): old_parameters = model.getParameters() old_plist = model.getParameterList() - p = x0 + p = copy.deepcopy(x0) if len(plist): model.setParameterList(plist) - p = x0full[:] - p[plist] = x0 + p = copy.deepcopy(x0full) + p[plist] = copy.deepcopy(x0) else: model.requireSensitivitiesForAllParameters() verbose and print('g: p=%s' % p) model.setParameters(p) + rdata = amici.runAmiciSimulation(model, solver, edata) - + model.setParameters(old_parameters) model.setParameterList(old_plist) @@ -148,13 +163,18 @@ def grad(x0, symbol='llh', x0full=None, plist=None, verbose=False): return res p = np.array(model.getParameters()) + + rdata = amici.runAmiciSimulation(model, solver, edata) for ip in range(model.np()): plist = [ip] - err_norm = check_grad(func, grad, p[plist], 'llh', p, [ip]) - print('sllh: p[%d]: |error|_2: %f' % (ip, err_norm)) + err_norm = check_grad(func, grad, copy.deepcopy(p[plist]), 'llh', + copy.deepcopy(p), [ip], epsilon=epsilon) + print(f'sllh: p[{ip}]: |error|_2: {err_norm}') + assert err_norm <= rtol * abs(rdata["sllh"][ip]) \ + or err_norm <= atol + - rdata = amici.runAmiciSimulation(model, solver, edata) leastsquares_applicable = solver.getSensitivityMethod() == amici.SensitivityMethod_forward @@ -166,26 +186,33 @@ def grad(x0, symbol='llh', x0full=None, plist=None, verbose=False): if leastsquares_applicable: for ip in range(model.np()): plist = [ip] - err_norm = check_grad(func, grad, p[plist], 'res', p, [ip]) + err_norm = check_grad(func, grad, copy.deepcopy(p[plist]), 'res', + copy.deepcopy(p), [ip], epsilon=epsilon) print('sres: p[%d]: |error|_2: %f' % (ip, err_norm)) + assert err_norm < atol \ + or err_norm < rtol * abs(rdata['sres'][:,ip].sum()) checkResults(rdata, 'FIM', np.dot(rdata['sres'].transpose(),rdata['sres']), 1e-8, 1e-4) checkResults(rdata, 'sllh', -np.dot(rdata['res'].transpose(),rdata['sres']), 1e-8, 1e-4) - ''' - print() - for ip in range(model.np()): - plist = [ip] - err_norm = check_grad(func, grad, p[plist], 'y', p, [ip]) - print('sy: p[%d]: |error|_2: %f' % (ip, err_norm)) - - print() - for ip in range(model.np()): - plist = [ip] - err_norm = check_grad(func, grad, p[plist], 'x', p, [ip]) - print('sx: p[%d]: |error|_2: %f' % (ip, err_norm)) - - ''' + + print() + for ip in range(model.np()): + plist = [ip] + err_norm = check_grad(func, grad, copy.deepcopy(p[plist]), 'y', + copy.deepcopy(p), [ip], epsilon=epsilon) + print('sy: p[%d]: |error|_2: %f' % (ip, err_norm)) + assert err_norm < atol \ + or err_norm < rtol * abs(rdata['sy'][:,ip].sum()) + + print() + for ip in range(model.np()): + plist = [ip] + err_norm = check_grad(func, grad, copy.deepcopy(p[plist]), 'x', + copy.deepcopy(p), [ip], epsilon=epsilon) + print('sx: p[%d]: |error|_2: %f' % (ip, err_norm)) + assert err_norm < atol \ + or err_norm < rtol * abs(rdata['sx'][:,ip].sum()) def verifySimulationResults(rdata, expectedResults, atol=1e-8, rtol=1e-4): @@ -206,7 +233,7 @@ def verifySimulationResults(rdata, expectedResults, atol=1e-8, rtol=1e-4): for field in expectedResults.keys(): if field == 'diagnosis': - for subfield in expectedResults[field].keys(): + for subfield in ['J', 'xdot']: checkResults(rdata, subfield, expectedResults[field][subfield][()], 0, 2) else: if field == 's2llh': diff --git a/tests/testPYSB.py b/tests/testPYSB.py index a615e4a6b3..58511d96d9 100755 --- a/tests/testPYSB.py +++ b/tests/testPYSB.py @@ -9,7 +9,7 @@ from pysb.simulator import ScipyOdeSimulator from pysb.examples import tyson_oscillator, robertson, \ - expression_observables, earm_1_3, bax_pore_sequential, bax_pore, \ + expression_observables, bax_pore_sequential, bax_pore, \ bngwiki_egfr_simple @@ -34,7 +34,6 @@ def runTest(self): self.compare_to_sbml_import() def compare_to_sbml_import(self): - constant_parameters = ['DRUG_0', 'KIN_0'] # -------------- PYSB ----------------- @@ -101,7 +100,7 @@ def compare_to_sbml_import(self): def compare_to_pysb_simulation(self): examples = [tyson_oscillator.model, robertson.model, - expression_observables.model, earm_1_3.model, + expression_observables.model, bax_pore_sequential.model, bax_pore.model, bngwiki_egfr_simple.model] for example in examples: diff --git a/tests/testPandas.py b/tests/testPandas.py index 9a774e5193..1c53f9257c 100755 --- a/tests/testPandas.py +++ b/tests/testPandas.py @@ -42,6 +42,8 @@ def setUp(self): self.solver = self.model.getSolver() rdata = amici.runAmiciSimulation(self.model, self.solver) self.edata = [amici.ExpData(rdata, 0.01, 0)] + # test copy constructor + self.edata_copy = amici.ExpData(self.edata[0]) def tearDown(self): os.chdir(self.resetdir) diff --git a/tests/testSBML.py b/tests/testSBML.py index 1337f62d75..c7f4c03ee8 100755 --- a/tests/testSBML.py +++ b/tests/testSBML.py @@ -5,6 +5,7 @@ import unittest import os import numpy as np +from testModels import checkDerivatives class TestAmiciSBMLModel(unittest.TestCase): ''' @@ -49,12 +50,16 @@ def test_presimulation(self): import test_model_presimulation as modelModule model = modelModule.getModel() solver = model.getSolver() + solver.setNewtonMaxSteps(0) model.setTimepoints(np.linspace(0, 60, 61)) + model.setSteadyStateSensitivityMode( + amici.SteadyStateSensitivityMode_simulationFSA + ) + solver.setSensitivityOrder(amici.SensitivityOrder_first) model.setReinitializeFixedParameterInitialStates(True) rdata = amici.runAmiciSimulation(model, solver) edata = amici.ExpData(rdata, 0.1, 0.0) - edata.fixedParametersPreequilibration = amici.DoubleVector([3, 0]) edata.fixedParameters = [10, 2] edata.fixedParametersPresimulation = [10, 2] edata.fixedParametersPreequilibration = [3, 0] @@ -62,6 +67,10 @@ def test_presimulation(self): amici.runAmiciSimulation(model, solver, edata), dict) + solver.setRelativeTolerance(1e-12) + solver.setAbsoluteTolerance(1e-12) + checkDerivatives(model, solver, edata, epsilon=1e-6) + def test_steadystate_scaled(self): ''' Test SBML import and simulation from AMICI python interface @@ -93,8 +102,9 @@ def test_steadystate_scaled(self): model = modelModule.getModel() model.setTimepoints(np.linspace(0, 60, 60)) solver = model.getSolver() + solver.setSensitivityOrder(amici.SensitivityOrder_first) rdata = amici.runAmiciSimulation(model, solver) - edata = [amici.ExpData(rdata, 0.01, 0)] + edata = [amici.ExpData(rdata, 1, 0)] rdata = amici.runAmiciSimulations(model, solver, edata) # check roundtripping of DataFrame conversion @@ -149,6 +159,11 @@ def test_steadystate_scaled(self): ) amici.getResidualsAsDataFrame(model, edata, rdata) + solver.setRelativeTolerance(1e-12) + solver.setAbsoluteTolerance(1e-12) + checkDerivatives(model, solver, edata[0], atol=1e-3, + rtol=1e-3, epsilon=1e-7) + if __name__ == '__main__': suite = unittest.TestSuite() diff --git a/version.txt b/version.txt index 88a7b22857..a48658c949 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.7.12 +0.7.13