diff --git a/.travis.yml b/.travis.yml index d14f75c7da..ec0cf0acf8 100644 --- a/.travis.yml +++ b/.travis.yml @@ -58,13 +58,17 @@ install: - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then pyenv versions; fi - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then pyenv shell 2.7 3.6; fi - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export PATH=/Users/travis/Library/Python/3.6/bin:$PATH; fi - - pip3 install --user --upgrade pip setuptools wheel pkgconfig doxypypy coverage + - pip3 install --user --upgrade pip setuptools wheel pkgconfig doxypypy coverage scipy - ./scripts/buildSuiteSparse.sh - ./scripts/buildSundials.sh - ./scripts/buildCpputest.sh - ./scripts/buildAmici.sh script: + - cd $BASE_DIR/python/sdist + - python3 setup.py sdist --dist-dir=$BASE_DIR/build/python/ + - python3 -m pip install $BASE_DIR/build/python/amici-*.tar.gz --verbose + - cd $BASE_DIR - ./scripts/run-cpputest.sh - ./scripts/run-cppcheck.sh - cd $BASE_DIR/build && make python-tests && cd $BASE_DIR diff --git a/include/amici/rdata.h b/include/amici/rdata.h index 799e8850d5..40a8cda5e5 100644 --- a/include/amici/rdata.h +++ b/include/amici/rdata.h @@ -83,7 +83,7 @@ class ReturnData { /** state (dimension: nt x nx, row-major) */ std::vector x; - /** parameter derivative of state (dimension: nt x nx x nplist, + /** parameter derivative of state (dimension: nt x nplist x nx, * row-major) */ std::vector sx; @@ -93,12 +93,12 @@ class ReturnData { /** observable standard deviation (dimension: nt x ny, row-major) */ std::vector sigmay; - /** parameter derivative of observable (dimension: nt x ny x nplist, + /** parameter derivative of observable (dimension: nt x nplist x ny, * row-major) */ std::vector sy; /** parameter derivative of observable standard deviation (dimension: nt x - * ny x nplist, row-major) */ + * nplist x ny, row-major) */ std::vector ssigmay; /** observable (dimension: nt x ny, row-major) */ diff --git a/matlab/mtoc/config/Doxyfile.template b/matlab/mtoc/config/Doxyfile.template index 9dd3eed592..ddf6cc10c6 100644 --- a/matlab/mtoc/config/Doxyfile.template +++ b/matlab/mtoc/config/Doxyfile.template @@ -697,6 +697,7 @@ EXCLUDE_PATTERNS = "_SourceDir_/models/*" \ "_SourceDir_/matlab/auxiliary/*" \ "_SourceDir_/matlab/SBMLImporter/*" \ "_SourceDir_/python/test/*" \ + "_SourceDir_/python/sdist/*" \ "_SourceDir_/build/*" \ "_SourceDir_/build_xcode/*" \ "_SourceDir_/swig/*" diff --git a/models/model_dirac/model_dirac.h b/models/model_dirac/model_dirac.h index 4ee8326c2b..e8d2601101 100644 --- a/models/model_dirac/model_dirac.h +++ b/models/model_dirac/model_dirac.h @@ -1,6 +1,6 @@ #ifndef _amici_model_dirac_h #define _amici_model_dirac_h -/* Generated by amiwrap (R2017b) 2632f34b3b684c381e19069e73b4b1277492fcf1 */ +/* Generated by amiwrap (R2017b) 53e279f979516a66d9b07a42c171c4fa35b6f26d */ #include #include #include "amici/defines.h" @@ -62,7 +62,7 @@ class Model_model_dirac : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_dirac(*this); }; - const char* getAmiciVersion() const { return "2632f34b3b684c381e19069e73b4b1277492fcf1"; }; + const char* getAmiciVersion() const { return "53e279f979516a66d9b07a42c171c4fa35b6f26d"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_dirac(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_events/model_events.h b/models/model_events/model_events.h index 19a38d5160..9292d758bd 100644 --- a/models/model_events/model_events.h +++ b/models/model_events/model_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_events_h #define _amici_model_events_h -/* Generated by amiwrap (R2017b) 2632f34b3b684c381e19069e73b4b1277492fcf1 */ +/* Generated by amiwrap (R2017b) 53e279f979516a66d9b07a42c171c4fa35b6f26d */ #include #include #include "amici/defines.h" @@ -76,7 +76,7 @@ class Model_model_events : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_events(*this); }; - const char* getAmiciVersion() const { return "2632f34b3b684c381e19069e73b4b1277492fcf1"; }; + const char* getAmiciVersion() const { return "53e279f979516a66d9b07a42c171c4fa35b6f26d"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_events(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_jakstat_adjoint/model_jakstat_adjoint.h b/models/model_jakstat_adjoint/model_jakstat_adjoint.h index 7c733dcd6a..a060ed3171 100644 --- a/models/model_jakstat_adjoint/model_jakstat_adjoint.h +++ b/models/model_jakstat_adjoint/model_jakstat_adjoint.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_h #define _amici_model_jakstat_adjoint_h -/* Generated by amiwrap (R2017b) 2632f34b3b684c381e19069e73b4b1277492fcf1 */ +/* Generated by amiwrap (R2017b) 53e279f979516a66d9b07a42c171c4fa35b6f26d */ #include #include #include "amici/defines.h" @@ -65,7 +65,7 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_jakstat_adjoint(*this); }; - const char* getAmiciVersion() const { return "2632f34b3b684c381e19069e73b4b1277492fcf1"; }; + const char* getAmiciVersion() const { return "53e279f979516a66d9b07a42c171c4fa35b6f26d"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_jakstat_adjoint(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h index 85debf6dc8..c337e493a4 100644 --- a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h +++ b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_o2_h #define _amici_model_jakstat_adjoint_o2_h -/* Generated by amiwrap (R2017b) 2632f34b3b684c381e19069e73b4b1277492fcf1 */ +/* Generated by amiwrap (R2017b) 53e279f979516a66d9b07a42c171c4fa35b6f26d */ #include #include #include "amici/defines.h" @@ -65,7 +65,7 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_jakstat_adjoint_o2(*this); }; - const char* getAmiciVersion() const { return "2632f34b3b684c381e19069e73b4b1277492fcf1"; }; + const char* getAmiciVersion() const { return "53e279f979516a66d9b07a42c171c4fa35b6f26d"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_jakstat_adjoint_o2(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_nested_events/model_nested_events.h b/models/model_nested_events/model_nested_events.h index 92a18b5861..fb81aa82a6 100644 --- a/models/model_nested_events/model_nested_events.h +++ b/models/model_nested_events/model_nested_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_nested_events_h #define _amici_model_nested_events_h -/* Generated by amiwrap (R2017b) 2632f34b3b684c381e19069e73b4b1277492fcf1 */ +/* Generated by amiwrap (R2017b) 53e279f979516a66d9b07a42c171c4fa35b6f26d */ #include #include #include "amici/defines.h" @@ -65,7 +65,7 @@ class Model_model_nested_events : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_nested_events(*this); }; - const char* getAmiciVersion() const { return "2632f34b3b684c381e19069e73b4b1277492fcf1"; }; + const char* getAmiciVersion() const { return "53e279f979516a66d9b07a42c171c4fa35b6f26d"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_nested_events(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_neuron/model_neuron.h b/models/model_neuron/model_neuron.h index 49089a5901..6a0156e238 100644 --- a/models/model_neuron/model_neuron.h +++ b/models/model_neuron/model_neuron.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_h #define _amici_model_neuron_h -/* Generated by amiwrap (R2017b) 2632f34b3b684c381e19069e73b4b1277492fcf1 */ +/* Generated by amiwrap (R2017b) 53e279f979516a66d9b07a42c171c4fa35b6f26d */ #include #include #include "amici/defines.h" @@ -79,7 +79,7 @@ class Model_model_neuron : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_neuron(*this); }; - const char* getAmiciVersion() const { return "2632f34b3b684c381e19069e73b4b1277492fcf1"; }; + const char* getAmiciVersion() const { return "53e279f979516a66d9b07a42c171c4fa35b6f26d"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_neuron(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_neuron_o2/model_neuron_o2.h b/models/model_neuron_o2/model_neuron_o2.h index 12ce6d9d7c..11d3d00167 100644 --- a/models/model_neuron_o2/model_neuron_o2.h +++ b/models/model_neuron_o2/model_neuron_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_o2_h #define _amici_model_neuron_o2_h -/* Generated by amiwrap (R2017b) 2632f34b3b684c381e19069e73b4b1277492fcf1 */ +/* Generated by amiwrap (R2017b) 53e279f979516a66d9b07a42c171c4fa35b6f26d */ #include #include #include "amici/defines.h" @@ -81,7 +81,7 @@ class Model_model_neuron_o2 : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_neuron_o2(*this); }; - const char* getAmiciVersion() const { return "2632f34b3b684c381e19069e73b4b1277492fcf1"; }; + const char* getAmiciVersion() const { return "53e279f979516a66d9b07a42c171c4fa35b6f26d"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_neuron_o2(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_robertson/model_robertson.h b/models/model_robertson/model_robertson.h index 408d40a868..d6f9ebcf17 100644 --- a/models/model_robertson/model_robertson.h +++ b/models/model_robertson/model_robertson.h @@ -1,6 +1,6 @@ #ifndef _amici_model_robertson_h #define _amici_model_robertson_h -/* Generated by amiwrap (R2017b) 2632f34b3b684c381e19069e73b4b1277492fcf1 */ +/* Generated by amiwrap (R2017b) 53e279f979516a66d9b07a42c171c4fa35b6f26d */ #include #include #include "amici/defines.h" @@ -63,7 +63,7 @@ class Model_model_robertson : public amici::Model_DAE { virtual amici::Model* clone() const override { return new Model_model_robertson(*this); }; - const char* getAmiciVersion() const { return "2632f34b3b684c381e19069e73b4b1277492fcf1"; }; + const char* getAmiciVersion() const { return "53e279f979516a66d9b07a42c171c4fa35b6f26d"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override { J_model_robertson(J, t, x, p, k, h, cj, dx, w, dwdx); diff --git a/models/model_steadystate/model_steadystate.h b/models/model_steadystate/model_steadystate.h index c6acd85107..8556625eb1 100644 --- a/models/model_steadystate/model_steadystate.h +++ b/models/model_steadystate/model_steadystate.h @@ -1,6 +1,6 @@ #ifndef _amici_model_steadystate_h #define _amici_model_steadystate_h -/* Generated by amiwrap (R2017b) 2632f34b3b684c381e19069e73b4b1277492fcf1 */ +/* Generated by amiwrap (R2017b) 53e279f979516a66d9b07a42c171c4fa35b6f26d */ #include #include #include "amici/defines.h" @@ -62,7 +62,7 @@ class Model_model_steadystate : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_steadystate(*this); }; - const char* getAmiciVersion() const { return "2632f34b3b684c381e19069e73b4b1277492fcf1"; }; + const char* getAmiciVersion() const { return "53e279f979516a66d9b07a42c171c4fa35b6f26d"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_steadystate(J, t, x, p, k, h, w, dwdx); diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 3f218da684..b802b87c18 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -2,11 +2,11 @@ find_package(PythonInterp 3.6 REQUIRED) # Configure setup.py and copy to output directory set(AMICI_VERSION "@AMICI_VERSION@") # to be replaced later -SET(SETUP_PY_IN ${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in) -SET(SETUP_PY_OUT ${CMAKE_CURRENT_BINARY_DIR}/setup.py) -CONFIGURE_FILE(${SETUP_PY_IN} ${SETUP_PY_OUT}) +set(SETUP_PY_IN ${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in) +set(SETUP_PY_OUT ${CMAKE_CURRENT_BINARY_DIR}/setup.py) +configure_file(${SETUP_PY_IN} ${SETUP_PY_OUT}) -CONFIGURE_FILE(${PROJECT_SOURCE_DIR}/python/amici/setup.template.py ${CMAKE_CURRENT_BINARY_DIR}/setup.template.py) +configure_file(${PROJECT_SOURCE_DIR}/python/amici/setup.template.py ${CMAKE_CURRENT_BINARY_DIR}/setup.template.py) # Copy further amici-python-module files to binary_dir add_custom_target(python-module-files @@ -53,7 +53,9 @@ add_custom_target(python-wheel COMMAND ${PYTHON_EXECUTABLE} setup.py bdist_wheel --dist-dir=${CMAKE_CURRENT_BINARY_DIR} ) + add_custom_command( OUTPUT always_rebuild COMMAND cmake -E echo ) + diff --git a/python/amici/__init__.py b/python/amici/__init__.py index 9344ba4020..cce06bf6ab 100644 --- a/python/amici/__init__.py +++ b/python/amici/__init__.py @@ -69,8 +69,9 @@ def runAmiciSimulation(model, solver, edata=None): Raises: """ - if edata: + if edata and edata.__class__.__name__ == 'ExpDataPtr': edata = edata.get() + rdata = amici.runAmiciSimulation(solver.get(), edata, model.get()) return rdataToNumPyArrays(rdata) diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py index 735361ad97..4cacd9fcac 100644 --- a/python/amici/sbml_import.py +++ b/python/amici/sbml_import.py @@ -133,7 +133,9 @@ def __init__(self, SBMLFile): ' const realtype *k, const realtype *h)'}, 'dydp': { 'signature': '(double *dydp, const realtype t, const realtype *x, const realtype *p,' - ' const realtype *k, const realtype *h)'}, + ' const realtype *k, const realtype *h, const int ip)', + 'sensitivity': True}, + 'qBdot': { 'signature': '(realtype *qBdot, const int ip, const realtype t, const realtype *x, const realtype *p,' ' const realtype *k, const realtype *h, const realtype *xB, const realtype *w,' @@ -222,7 +224,7 @@ def loadSBMLFile(self, SBMLFile): self.sbml = self.sbml_doc.getModel() - def sbml2amici(self, modelName, output_dir=None, observables={}, constantParameters=[], sigmas={}): + def sbml2amici(self, modelName, output_dir=None, observables={}, constantParameters=[], sigmas={}, verbose=False): """Generate AMICI C++ files for the model provided to the constructor. Arguments: @@ -231,6 +233,7 @@ def sbml2amici(self, modelName, output_dir=None, observables={}, constantParamet observables: dictionary(observableName:formulaString) to be added to the model sigmas: dictionary(observableName: sigma value or (existing) parameter name) constantParameters: list of SBML Ids identifying constant parameters + verbose: more verbose output if True Returns: Raises: @@ -243,7 +246,7 @@ def sbml2amici(self, modelName, output_dir=None, observables={}, constantParamet self.computeModelEquations(observables, sigmas) self.prepareModelFolder() self.generateCCode() - self.compileCCode() + self.compileCCode(verbose) def setName(self, modelName): @@ -757,7 +760,7 @@ def computeModelEquationsObjectiveFunction(self, observables={}, sigmas={}): self.symbols['my']['sym'] = sp.DenseMatrix([sp.sympify('m' + str(symbol)) for symbol in observableSyms]) - loglikelihoodString = lambda strSymbol: '0.5*sqrt(2*pi*sigma{symbol}**2) + 0.5*(({symbol}-m{symbol})/sigma{symbol})**2'.format(symbol=strSymbol) + loglikelihoodString = lambda strSymbol: '0.5*log(2*pi*sigma{symbol}**2) + 0.5*(({symbol}-m{symbol})/sigma{symbol})**2'.format(symbol=strSymbol) self.functions['Jy']['sym'] = sp.DenseMatrix([sp.sympify(loglikelihoodString(str(symbol))) for symbol in observableSyms]) self.functions['dJydy']['sym'] = self.functions['Jy']['sym'].jacobian(observableSyms) self.functions['dJydsigma']['sym'] = self.functions['Jy']['sym'].jacobian(sigmaYSyms) @@ -782,7 +785,7 @@ def computeModelEquationsSensitivitesCore(self): self.functions['dydp']['sym'] = self.functions['y']['sym']\ .jacobian(self.symbols['parameter']['sym']) self.functions['dydx']['sym'] = self.functions['y']['sym']\ - .jacobian(self.symbols['species']['sym']) + .jacobian(self.symbols['species']['sym']).transpose() self.functions['dwdp']['sym'] = self.fluxVector.jacobian(self.symbols['parameter']['sym']) @@ -882,26 +885,43 @@ def generateCCode(self): os.path.join(self.modelPath, 'main.cpp')) - def compileCCode(self): + def compileCCode(self, verbose=False): """Compile the generated model code Arguments: - + verbose: Make model compilation verbose Returns: Raises: """ - + + # setup.py assumes it is run from within the model directory moduleDir = self.modelPath - oldCwd = os.getcwd() - os.chdir(moduleDir) # setup.py assumes it is run from within the model dir - from distutils.core import run_setup - run_setup('setup.py', - script_args=["build_ext", - '--build-lib=%s' % moduleDir, - ]) - os.chdir(oldCwd) + + script_args = [sys.executable, '%s%ssetup.py' % (moduleDir, os.sep)] + + if verbose: + script_args.append('--verbose') + else: + script_args.append('--quiet') + + script_args.extend(['build_ext', '--build-lib=%s' % moduleDir]) + + # distutils.core.run_setup looks nicer, but does not let us check the + # result easily + try: + result = subprocess.run(script_args, + cwd=moduleDir, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + check=True) + except subprocess.CalledProcessError as e: + print(e.output.decode('utf-8')) + raise + + if verbose: + print(result.stdout.decode('utf-8')) def writeIndexFiles(self,name): """Write index file for a symbolic array. diff --git a/python/amici/setup.template.py b/python/amici/setup.template.py index 1b86d32152..66bb3a26e8 100644 --- a/python/amici/setup.template.py +++ b/python/amici/setup.template.py @@ -26,17 +26,18 @@ def getAmiciLibs(): import pkgconfig h5pkgcfg = pkgconfig.parse("hdf5") -cxx_flags = ['-std=c++0x'] +cxx_flags = ['-std=c++11'] #linker_flags = ['${BLAS_LIBRARIES}'] linker_flags = [] if 'ENABLE_GCOV_COVERAGE' in os.environ and os.environ['ENABLE_GCOV_COVERAGE'] == 'TRUE': cxx_flags.extend(['-g', '-O0', '--coverage']) linker_flags.append('--coverage') -libraries = ['cblas',# TODO generic BLAS +libraries = [*getAmiciLibs(), + 'cblas',# TODO generic BLAS 'hdf5_hl_cpp', 'hdf5_hl', 'hdf5_cpp', 'hdf5'] sources = ['swig/TPL_MODELNAME.i', *getModelSources()] -libraries.extend(getAmiciLibs()) + # Remove the "-Wstrict-prototypes" compiler option, which isn't valid for diff --git a/python/examples/example_steadystate/ExampleSteadystate.ipynb b/python/examples/example_steadystate/ExampleSteadystate.ipynb new file mode 100644 index 0000000000..b31dfb584f --- /dev/null +++ b/python/examples/example_steadystate/ExampleSteadystate.ipynb @@ -0,0 +1,1157 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# AMICI Python example \"steadystate\"\n", + "Example using [model_steadystate_scaled.sbml] model to demonstrate and test SBML import and Amici Python interface." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "sbml_file = 'model_steadystate_scaled.sbml'\n", + "model_name = 'model_steadystate_scaled'\n", + "model_output_dir= 'model_steadystate_scaled'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The example model" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Species: ['x1', 'x2', 'x3']\n", + "\n", + "Reactions:\n", + " r1: 2 x1 -> x2\t\t[p1 * x1^2]\n", + " r2: x1 + x2 -> x3\t\t[p2 * x1 * x2]\n", + " r3: x2 -> 2 x1\t\t[p3 * x2]\n", + " r4: x3 -> x1 + x2\t\t[p4 * x3]\n", + " r5: x3 -> \t\t[k4 * x3]\n", + " r6: -> x1\t\t[p5]\n" + ] + } + ], + "source": [ + "import libsbml\n", + "SBMLreader = libsbml.SBMLReader()\n", + "sbml_doc = SBMLreader.readSBML(sbml_file)\n", + "sbml_model = sbml_doc.getModel()\n", + "\n", + "print('Species: ', [s.getId() for s in sbml_model.getListOfSpecies()])\n", + "\n", + "print('\\nReactions:')\n", + "for reaction in sbml_model.getListOfReactions():\n", + " reactants = ' + '.join(['%s %s'%(int(r.getStoichiometry()) if r.getStoichiometry() > 1 else '', r.getSpecies()) for r in reaction.getListOfReactants()])\n", + " products = ' + '.join(['%s %s'%(int(r.getStoichiometry()) if r.getStoichiometry() > 1 else '', r.getSpecies()) for r in reaction.getListOfProducts()])\n", + " reversible = '<' if reaction.getReversible() else ''\n", + " print('%3s: %10s %1s->%10s\\t\\t[%s]' % (reaction.getId(), \n", + " reactants,\n", + " reversible,\n", + " products,\n", + " libsbml.formulaToL3String(reaction.getKineticLaw().getMath())))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Importing an SBML model, compiling and generating an AMICI module" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'observable_x1': 'x1', 'observable_x2': 'x2', 'observable_x3': 'x3', 'observable_x1_scaled': 'scaling_x1 * x1', 'observable_x2_offsetted': 'offset_x2 + x2', 'observable_x1withsigma': 'x1'}\n", + "running build_ext\n", + "building 'model_steadystate_scaled/_model_steadystate_scaled' extension\n", + "swigging swig/model_steadystate_scaled.i to swig/model_steadystate_scaled_wrap.cpp\n", + "swig -python -c++ -modern -outdir model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/swig -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -o swig/model_steadystate_scaled_wrap.cpp swig/model_steadystate_scaled.i\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c swig/model_steadystate_scaled_wrap.cpp -o build/temp.linux-x86_64-3.6/swig/model_steadystate_scaled_wrap.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_JvB.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JvB.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_Jy.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_Jy.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_sx0.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_sx0.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_sigma_y.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_sigma_y.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_Jv.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_Jv.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_J.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_J.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_xdot.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_xdot.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_dxdotdp.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dxdotdp.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_dwdx.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dwdx.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c wrapfunctions.cpp -o build/temp.linux-x86_64-3.6/wrapfunctions.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_dydx.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dydx.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_dwdp.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dwdp.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_JSparseB.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JSparseB.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_JSparse.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JSparse.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_dydp.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dydp.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_sxdot.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_sxdot.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_w.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_w.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_JB.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JB.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_JDiag.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JDiag.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_qBdot.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_qBdot.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_x0.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_x0.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_xBdot.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_xBdot.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_y.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_y.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_dJydsigma.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dJydsigma.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c model_steadystate_scaled_dJydy.cpp -o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dJydy.o -std=c++11\n", + "x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled -I/home/dweindl/.local/lib/python3.6/site-packages/amici/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/sundials/include -I/home/dweindl/.local/lib/python3.6/site-packages/amici/ThirdParty/SuiteSparse/include -I/usr/include/hdf5/serial -I/usr/include/python3.6m -c main.cpp -o build/temp.linux-x86_64-3.6/main.o -std=c++11\n", + "x86_64-linux-gnu-g++ -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -specs=/usr/share/dpkg/no-pie-link.specs -Wl,-z,relro -Wl,-Bsymbolic-functions -specs=/usr/share/dpkg/no-pie-link.specs -Wl,-z,relro -g -fdebug-prefix-map=/build/python3.6-sXpGnM/python3.6-3.6.3=. -specs=/usr/share/dpkg/no-pie-compile.specs -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 build/temp.linux-x86_64-3.6/swig/model_steadystate_scaled_wrap.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JvB.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_Jy.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_sx0.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_sigma_y.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_Jv.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_J.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_xdot.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dxdotdp.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dwdx.o build/temp.linux-x86_64-3.6/wrapfunctions.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dydx.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dwdp.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JSparseB.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JSparse.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dydp.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_sxdot.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_w.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JB.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_JDiag.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_qBdot.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_x0.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_xBdot.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_y.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dJydsigma.o build/temp.linux-x86_64-3.6/model_steadystate_scaled_dJydy.o build/temp.linux-x86_64-3.6/main.o -L/usr/lib/x86_64-linux-gnu/hdf5/serial -L/home/dweindl/.local/lib/python3.6/site-packages/amici/libs -lamici -lsundials -lsuitesparse -lcblas -lhdf5_hl_cpp -lhdf5_hl -lhdf5_cpp -lhdf5 -o /home/dweindl/src/AMICI-devel/python/examples/example_steadystate/model_steadystate_scaled/model_steadystate_scaled/_model_steadystate_scaled.cpython-36m-x86_64-linux-gnu.so\n", + "\n" + ] + } + ], + "source": [ + "import amici\n", + "\n", + "import os\n", + "import sys\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "def createModule(sbml_file, model_name, model_output_dir):\n", + " \"\"\"Create Python module from SBML model\"\"\"\n", + " sbmlImporter = amici.SbmlImporter(sbml_file)\n", + " sbml = sbmlImporter.sbml\n", + " \n", + " observables = amici.assignmentRules2observables(sbml, filter=lambda variableId: \n", + " variableId.startswith('observable_') and not variableId.endswith('_sigma'))\n", + " \n", + " print(observables)\n", + " \n", + " sbmlImporter.sbml2amici(model_name, model_output_dir, verbose=True,\n", + " observables=observables,\n", + " constantParameters=['k4'],\n", + " sigmas={'observable_x1withsigma': 'observable_x1withsigma_sigma'})\n", + "\n", + "createModule(sbml_file, model_name, model_output_dir)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To use python-generated model in matlab, run:\n", + "```\n", + "modelName = '';\n", + "modelDir = '';\n", + "amimodel.compileAndLinkModel(modelName, modelDir, [], [], [], []);\n", + "amimodel.generateMatlabWrapper(3, 6, 8, 1, 0, 0, [], [ modelDir '/simulate_test.m'], modelName, 'lin', 1, 1);\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Running simulations and analyzing results" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Parameters: [1. 0.5 0.4 2. 0.1 1. 1. 1. ]\n", + "\n", + " t: [ 0. 1.01694915 2.03389831 3.05084746 4.06779661 5.08474576\n", + " 6.10169492 7.11864407 8.13559322 9.15254237 10.16949153 11.18644068\n", + " 12.20338983 13.22033898 14.23728814 15.25423729 16.27118644 17.28813559\n", + " 18.30508475 19.3220339 20.33898305 21.3559322 22.37288136 23.38983051\n", + " 24.40677966 25.42372881 26.44067797 27.45762712 28.47457627 29.49152542\n", + " 30.50847458 31.52542373 32.54237288 33.55932203 34.57627119 35.59322034\n", + " 36.61016949 37.62711864 38.6440678 39.66101695 40.6779661 41.69491525\n", + " 42.71186441 43.72881356 44.74576271 45.76271186 46.77966102 47.79661017\n", + " 48.81355932 49.83050847 50.84745763 51.86440678 52.88135593 53.89830508\n", + " 54.91525424 55.93220339 56.94915254 57.96610169 58.98305085 60. ]\n", + " x: [[0.1 0.4 0.7 ]\n", + " [0.57995051 0.73365809 0.0951589 ]\n", + " [0.55996496 0.71470091 0.0694127 ]\n", + " [0.5462855 0.68030366 0.06349394]\n", + " [0.53561883 0.64937432 0.05923555]\n", + " [0.52636487 0.62259567 0.05568686]\n", + " [0.51822014 0.59943346 0.05268079]\n", + " [0.51103767 0.57935661 0.05012037]\n", + " [0.5047003 0.56191593 0.04793052]\n", + " [0.49910666 0.54673518 0.0460508 ]\n", + " [0.49416809 0.53349812 0.04443206]\n", + " [0.48980688 0.52193768 0.043034 ]\n", + " [0.48595477 0.51182733 0.04182339]\n", + " [0.48255177 0.50297415 0.04077267]\n", + " [0.47954512 0.49521321 0.03985882]\n", + " [0.47688834 0.48840307 0.03906254]\n", + " [0.4745405 0.48242201 0.03836756]\n", + " [0.47246548 0.47716503 0.0377601 ]\n", + " [0.47063147 0.4725413 0.03722844]\n", + " [0.46901037 0.46847203 0.03676259]\n", + " [0.4675774 0.46488882 0.03635397]\n", + " [0.46631066 0.46173208 0.03599523]\n", + " [0.46519082 0.45894988 0.03568002]\n", + " [0.46420083 0.45649685 0.03540285]\n", + " [0.4633256 0.45433333 0.03515899]\n", + " [0.46255181 0.45242458 0.0349443 ]\n", + " [0.46186769 0.45074017 0.0347552 ]\n", + " [0.46126283 0.44925338 0.03458856]\n", + " [0.46072804 0.44794076 0.03444166]\n", + " [0.46025521 0.44678169 0.03431212]\n", + " [0.45983714 0.44575805 0.03419785]\n", + " [0.4594675 0.44485389 0.03409701]\n", + " [0.45914066 0.44405515 0.03400802]\n", + " [0.45885167 0.44334948 0.03392946]\n", + " [0.45859615 0.44272595 0.03386009]\n", + " [0.45837021 0.44217497 0.03379883]\n", + " [0.45817044 0.44168806 0.03374473]\n", + " [0.45799379 0.44125773 0.03369693]\n", + " [0.4578376 0.44087739 0.03365471]\n", + " [0.45769949 0.44054121 0.0336174 ]\n", + " [0.45757737 0.44024405 0.03358444]\n", + " [0.45746939 0.43998138 0.03355531]\n", + " [0.45737391 0.43974917 0.03352956]\n", + " [0.45728949 0.4395439 0.03350681]\n", + " [0.45721483 0.43936243 0.0334867 ]\n", + " [0.45714882 0.43920199 0.03346892]\n", + " [0.45709046 0.43906015 0.03345321]\n", + " [0.45703884 0.43893475 0.03343932]\n", + " [0.45699321 0.43882387 0.03342704]\n", + " [0.45695285 0.43872585 0.03341618]\n", + " [0.45691717 0.43863918 0.03340658]\n", + " [0.45688562 0.43856255 0.0333981 ]\n", + " [0.45685772 0.43849479 0.0333906 ]\n", + " [0.45683305 0.43843488 0.03338397]\n", + " [0.45681123 0.43838191 0.0333781 ]\n", + " [0.45679194 0.43833508 0.03337292]\n", + " [0.45677489 0.43829366 0.03336833]\n", + " [0.4567598 0.43825705 0.03336428]\n", + " [0.45674647 0.43822467 0.0333607 ]\n", + " [0.45673468 0.43819604 0.03335753]]\n", + " x0: [0.1 0.4 0.7]\n", + " sx: None\n", + " sx0: None\n", + " y: [[0.1 0.4 0.7 0.1 1.4 0.1 ]\n", + " [0.57995051 0.73365809 0.0951589 0.57995051 1.73365809 0.57995051]\n", + " [0.55996496 0.71470091 0.0694127 0.55996496 1.71470091 0.55996496]\n", + " [0.5462855 0.68030366 0.06349394 0.5462855 1.68030366 0.5462855 ]\n", + " [0.53561883 0.64937432 0.05923555 0.53561883 1.64937432 0.53561883]\n", + " [0.52636487 0.62259567 0.05568686 0.52636487 1.62259567 0.52636487]\n", + " [0.51822014 0.59943346 0.05268079 0.51822014 1.59943346 0.51822014]\n", + " [0.51103767 0.57935661 0.05012037 0.51103767 1.57935661 0.51103767]\n", + " [0.5047003 0.56191593 0.04793052 0.5047003 1.56191593 0.5047003 ]\n", + " [0.49910666 0.54673518 0.0460508 0.49910666 1.54673518 0.49910666]\n", + " [0.49416809 0.53349812 0.04443206 0.49416809 1.53349812 0.49416809]\n", + " [0.48980688 0.52193768 0.043034 0.48980688 1.52193768 0.48980688]\n", + " [0.48595477 0.51182733 0.04182339 0.48595477 1.51182733 0.48595477]\n", + " [0.48255177 0.50297415 0.04077267 0.48255177 1.50297415 0.48255177]\n", + " [0.47954512 0.49521321 0.03985882 0.47954512 1.49521321 0.47954512]\n", + " [0.47688834 0.48840307 0.03906254 0.47688834 1.48840307 0.47688834]\n", + " [0.4745405 0.48242201 0.03836756 0.4745405 1.48242201 0.4745405 ]\n", + " [0.47246548 0.47716503 0.0377601 0.47246548 1.47716503 0.47246548]\n", + " [0.47063147 0.4725413 0.03722844 0.47063147 1.4725413 0.47063147]\n", + " [0.46901037 0.46847203 0.03676259 0.46901037 1.46847203 0.46901037]\n", + " [0.4675774 0.46488882 0.03635397 0.4675774 1.46488882 0.4675774 ]\n", + " [0.46631066 0.46173208 0.03599523 0.46631066 1.46173208 0.46631066]\n", + " [0.46519082 0.45894988 0.03568002 0.46519082 1.45894988 0.46519082]\n", + " [0.46420083 0.45649685 0.03540285 0.46420083 1.45649685 0.46420083]\n", + " [0.4633256 0.45433333 0.03515899 0.4633256 1.45433333 0.4633256 ]\n", + " [0.46255181 0.45242458 0.0349443 0.46255181 1.45242458 0.46255181]\n", + " [0.46186769 0.45074017 0.0347552 0.46186769 1.45074017 0.46186769]\n", + " [0.46126283 0.44925338 0.03458856 0.46126283 1.44925338 0.46126283]\n", + " [0.46072804 0.44794076 0.03444166 0.46072804 1.44794076 0.46072804]\n", + " [0.46025521 0.44678169 0.03431212 0.46025521 1.44678169 0.46025521]\n", + " [0.45983714 0.44575805 0.03419785 0.45983714 1.44575805 0.45983714]\n", + " [0.4594675 0.44485389 0.03409701 0.4594675 1.44485389 0.4594675 ]\n", + " [0.45914066 0.44405515 0.03400802 0.45914066 1.44405515 0.45914066]\n", + " [0.45885167 0.44334948 0.03392946 0.45885167 1.44334948 0.45885167]\n", + " [0.45859615 0.44272595 0.03386009 0.45859615 1.44272595 0.45859615]\n", + " [0.45837021 0.44217497 0.03379883 0.45837021 1.44217497 0.45837021]\n", + " [0.45817044 0.44168806 0.03374473 0.45817044 1.44168806 0.45817044]\n", + " [0.45799379 0.44125773 0.03369693 0.45799379 1.44125773 0.45799379]\n", + " [0.4578376 0.44087739 0.03365471 0.4578376 1.44087739 0.4578376 ]\n", + " [0.45769949 0.44054121 0.0336174 0.45769949 1.44054121 0.45769949]\n", + " [0.45757737 0.44024405 0.03358444 0.45757737 1.44024405 0.45757737]\n", + " [0.45746939 0.43998138 0.03355531 0.45746939 1.43998138 0.45746939]\n", + " [0.45737391 0.43974917 0.03352956 0.45737391 1.43974917 0.45737391]\n", + " [0.45728949 0.4395439 0.03350681 0.45728949 1.4395439 0.45728949]\n", + " [0.45721483 0.43936243 0.0334867 0.45721483 1.43936243 0.45721483]\n", + " [0.45714882 0.43920199 0.03346892 0.45714882 1.43920199 0.45714882]\n", + " [0.45709046 0.43906015 0.03345321 0.45709046 1.43906015 0.45709046]\n", + " [0.45703884 0.43893475 0.03343932 0.45703884 1.43893475 0.45703884]\n", + " [0.45699321 0.43882387 0.03342704 0.45699321 1.43882387 0.45699321]\n", + " [0.45695285 0.43872585 0.03341618 0.45695285 1.43872585 0.45695285]\n", + " [0.45691717 0.43863918 0.03340658 0.45691717 1.43863918 0.45691717]\n", + " [0.45688562 0.43856255 0.0333981 0.45688562 1.43856255 0.45688562]\n", + " [0.45685772 0.43849479 0.0333906 0.45685772 1.43849479 0.45685772]\n", + " [0.45683305 0.43843488 0.03338397 0.45683305 1.43843488 0.45683305]\n", + " [0.45681123 0.43838191 0.0333781 0.45681123 1.43838191 0.45681123]\n", + " [0.45679194 0.43833508 0.03337292 0.45679194 1.43833508 0.45679194]\n", + " [0.45677489 0.43829366 0.03336833 0.45677489 1.43829366 0.45677489]\n", + " [0.4567598 0.43825705 0.03336428 0.4567598 1.43825705 0.4567598 ]\n", + " [0.45674647 0.43822467 0.0333607 0.45674647 1.43822467 0.45674647]\n", + " [0.45673468 0.43819604 0.03335753 0.45673468 1.43819604 0.45673468]]\n", + " sigmay: [[1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]]\n", + " sy: None\n", + " ssigmay: None\n", + " z: None\n", + " rz: None\n", + " sigmaz: None\n", + " sz: None\n", + " srz: None\n", + " ssigmaz: None\n", + " sllh: []\n", + " s2llh: None\n", + " J: [[-2.04603672 0.69437133 0.21909802]\n", + " [ 0.57163266 -0.62836734 0.22836734]\n", + " [ 2. 2. -3. ]]\n", + " xdot: [-1.08973588e-05 -2.64550670e-05 -2.92781904e-06]\n", + " status: 0.0\n", + " llh: nan\n", + " chi2: nan\n", + "newton_numlinsteps: None\n", + "newton_numsteps: [[0 0]]\n", + " numsteps: [ 0 100 145 172 188 199 206 210 214 217 221 224 226 229 231 233 236 238\n", + " 240 243 245 248 250 252 255 256 258 259 261 262 264 265 267 269 270 272\n", + " 273 275 276 278 279 281 282 283 284 285 286 287 288 289 290 291 292 293\n", + " 294 295 296 297 298 299]\n", + " numrhsevals: [ 0 115 161 189 206 218 226 232 236 239 243 247 250 253 255 257 260 262\n", + " 264 267 269 272 274 276 279 281 283 284 286 287 289 290 292 294 295 297\n", + " 298 300 301 303 304 307 308 309 310 311 312 313 314 315 316 317 318 319\n", + " 320 321 322 323 324 325]\n", + "numerrtestfails: [0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n", + " 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]\n", + "numnonlinsolvconvfails: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", + " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", + " order: [0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5\n", + " 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5]\n", + " numstepsB: []\n", + "numrhsevalsB: []\n", + "numerrtestfailsB: []\n", + "numnonlinsolvconvfailsB: []\n" + ] + } + ], + "source": [ + "sys.path.insert(0, os.path.abspath(model_output_dir))\n", + "import model_steadystate_scaled as modelModule\n", + "\n", + "model = modelModule.getModel()\n", + "\n", + "# show default parameters\n", + "print('Parameters:', np.array(model.getParameters()))\n", + "\n", + "# simulation using default parameters\n", + "model.setTimepoints(amici.DoubleVector(np.linspace(0, 60, 60))) \n", + "solver = model.getSolver()\n", + "rdata = amici.runAmiciSimulation(model, solver)\n", + "\n", + "print()\n", + "#np.set_printoptions(threshold=8, edgeitems=2)\n", + "for key, value in rdata.items():\n", + " print('%12s: ' % key, value)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plotting tractories" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEZCAYAAABiu9n+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8HWW9+PHP9+xJmnRNS9ukq4VS2S0FBKplu4CyXARlU1AuyBXEjSugyHa5PwSvoCheQGW5ChQuawUEWQpKFWiRCrSltNRCQ1uapkvaNGf//v6YOcnJyTnJOWnOSU7P9/1yXjPzzDMzz5SYb57nmXkeUVWMMcaYbDwDXQBjjDGDlwUJY4wxOVmQMMYYk5MFCWOMMTlZkDDGGJOTBQljjDE5WZAwpshEZImIfHYQlON2EfnRQJfDlBcLEmZQEZHDROSvIrJVRDaJyAIROdA9dq6IvFLAtSaJiIqIr49l2anzU1T1k6r60s5cQ0TuEZHrd7IcF6rqf+7MNUzl2akffmP6k4jUAU8C/w48BASAw4HIQJarJyLiU9X4QJejNyLiVdXEQJfDlB+rSZjBZHcAVX1AVROq2q6qf1LVt0RkT+B24BAR2S4iWwBE5HMi8qaItIrIGhG5Ju16f3bXW9xzDnHP+ZqILBORzSLyrIhMzFGebue7tZkFInKLiGwCrhGRqSLyooi0iMhGEblPRIalLiIiq0XkKHfbIyKXi8j7bv6HRGREWt5UTWqL+zznisgFwFnA991y/MHNu6eIvOTmXSIiJ6Zd5x4R+R8ReVpE2oA5mbUREfm8iCx2z/+riOyTduwyEflIRLaJyHIRObKA/45mV6KqttgyKBagDmgB7gWOA4ZnHD8XeCUj7bPA3jh/8OwDfAyc7B6bBCjgS8t/MrAS2BOnJn0l8Ncc5cl2/rlAHPime34V8AngaCAI1OMEl5+lnbMaOMrd/jbwKtDg5r8DeMA9NgHYBpwB+IGRwH7usXuA69Ou6Xef4wc4Na4j3HP3SMu/FTjU/bcJpV8DOADYABwEeIFz3HIGgT2ANcC4tH+HqQP982HLwCxWkzCDhqq2Aofh/GL+NdAsIvNEZEwP57ykqm+ralJV3wIeAD7Tw22+DtygqsvUaSb6f8B+PdQmslmrqr9Q1bg6tZ2VqvqcqkZUtRm4uYcyfB34oao2qWoEuAY41e33OAt4Xp2aVExVW1R1cY7rHAwMAX6sqlFVfRGnqe6MtDxPqOoC998mnHH++cAdqvqaOrW2e3Ga9Q4GEjjBYoaI+FV1taq+X8C/j9mFWJAwg4r7y/tcVW0A9gLGAT/LlV9EDhKR+SLSLCJbgQuBUT3cYiLwc7eJZQuwCRBgfAHFXJNRhtEiMtdtnmkFft9DGSYCj6XdfxnOL+UxQCOQ7y/jccAaVU2mpX2Q8RxryG0i8L1UOdyyNOLUHlbi1HiuATa4zzYuz3KZXYwFCTNoqeq7OE0ke6WSsmS7H5gHNKrqUJx+C+kh/xrg66o6LG2pUtW/ZitCrqJl7N/gpu2jqnXA2WllyHb/4zLuH1LVj9xjU/O851qgUUTS/z88Afgoj/KnyvFfGeWoVtUHAFT1flU9DCeYKHBjD9cyuzALEmbQEJHpIvI9EWlw9xtxmk9edbN8DDSISCDttFpgk6qGRWQWcGbasWYgCUxJS7sduEJEPuneY6iInJajSNnOz6YW2I7TwT0e+I8e8t4O/FeqeUtE6kXkJPfYfcBRIvJFEfGJyEgR2c899nFGOV4D2nA6s/3udxgnAHN7KWvKr4EL3ZqYiEiN+xJArYjsISJHiEgQCAPtOLUdU4EsSJjBZBtOR+pr7hs5rwLvAN9zj78ILAHWi8hGN+0bwHUisg24CufVWQBUdQfwX8ACt0nlYFV9DOev4rlu09A7OJ3k3WQ7P0e5r8XpCN4KPAU82sMz/hyn5vMnt8yvus+Mqn4IHO8+7yZgMbCve95vcfoItojI46oaBU50y74R+BXwFbf21StVXYTTL/FLYDNOJ/i57uEg8GP3uuuB0Tgd5KYCiapNOmRMMYnIh8DZqvrnXjMbM8hYTcKYIhKRepzXYlcPcFGM6RMLEsYUiTjDiawAfuE2JRlTdqy5yRhjTE5WkzDGGJNT2Q/wN2rUKJ00adJAF8MYY8rGG2+8sVFV6/PJW/ZBYtKkSSxatGigi2GMMWVDRD7IN681NxljjMnJgoQxxpicLEgYY4zJqez7JIwxpphisRhNTU2Ew5mjrQ9+oVCIhoYG/H5/n69hQcIYY3rQ1NREbW0tkyZNQiTX4L6Dj6rS0tJCU1MTkydP7vN1rLnJGGN6EA6HGTlyZFkFCAARYeTIkTtdA7IgYYwxvSi3AJHSH+W2IJHuo7/DiucHuhTGGDNoWJBI99INcN+p8NZDvec1xpgKYEEiXXgroPDYhbD0iYEujTHGDDgLEuki22DKHBj/KXj4PHjv2YEukTHGADBnzhyee+45AK688kouueSSktzXXoFNF9kGY/eD434M954ID34ZznwQps4Z6JIZYwaBa/+whKVrW/v1mjPG1XH1CZ/s/d7XXstVV13Fhg0bePPNN5k3b16/liMXq0mki7RCsBZCQ+HLj8HIT8ADZ8DqBQNdMmNMhZs9ezaqys0338zcuXPxer20tbVxzjnncP7553PfffcV5b5Wk0hRdWoSoTpnv3oEfOUJuOd4mHsGXLLYSTPGVKx8/uIvlrfffpt169YxatQoamtrAXj00Uc59dRTOeGEE/jSl77EWWed1e/3tZpESmwHaNKpSaQMqYfT7oVwK7xyy8CVzRhT0datW8dZZ53FE088QU1NDc8+6/SXNjU10djYCIDX6y3KvUsWJETkWBFZLiIrReTyLMdvEZHF7vKeiGwpVdkApxYBXYMEwJgZsO/p8PqdsPWjkhbJGGN27NjBKaecwk9/+lP23HNPfvSjH3HNNdcA0NDQQFNTEwDJZLIo9y9Jc5OIeIHbgKOBJmChiMxT1aWpPKr6nbT83wT2L2aZVJWEJvB53H+CjiBR1z3zZ6+Atx+Gl2+EE28tZrGMMaaL6upq/va3v3Xsz549u2P/lFNO4eKLL+app57ihBNOKMr9S9UnMQtYqaqrAERkLnASsDRH/jOAq4tZoOMfPZ79Ru/HDYff4CRE3DcWMmsSAMMnwoHnweu/hk9/E0ZNK2bRjDEmLzU1Ndx9991FvUepmpvGA2vS9pvctG5EZCIwGXgx18VE5AIRWSQii5qbm/tUoIA3QDQR7UzI1dyUcvil4AvBi9f36X7GGFOOShUkso0ypTnyng48rKqJXBdT1TtVdaaqzqyvz2su726C3mBhQWJIPXz6Ylj6uDPGkzHGVIBSBYkmoDFtvwFYmyPv6cADxS6Q3+snkoh0JvQWJAAOuRiqRsAL1xW3cMYYM0iUKkgsBKaJyGQRCeAEgm6fC4rIHsBw4G+Zx/pb0BvMESSydFynhOpg9qWwaj6seqmo5TPGmMGgJEFCVePAxcCzwDLgIVVdIiLXiciJaVnPAOaqaq6mqH4T8AaIJWOdCamO68CQnk+ceR7UNcDz1zof4BljzC6sZF9cq+rTwNMZaVdl7F9TqvIEPVlqEr4Q+AI9n+gPwZwfwBPfgLf/D/b5YnELaowxA6hiv7jO+nZTT/0R6fY9A8btD3/6UWczlTHG7IIqOkh0q0nkGyQ8Hjj+v2H7evjzT4pTQGOMGQQqNkhkfQU23yAB0DAT9jsb/vYr2Lii/wtojDGDQMUGiezNTT282ZTNUVeDvwr+eJl1YhtjisomHSqx7s1NrTC0MfcJ2QwZ7XRiP3M5LH8apn+ufwtpjBlc/ng5rH+7f6+5297ORGe9sEmHSizoDRJNRul427bQ5qaUA/8N6veEZ66AWHv/FtIYY1zZJh1atWoV5513HqeeemrR7luxNYmgNwhALBkj4A30PUh4/XD8TXDvCbDgVvjsZf1cUmPMoJHHX/zFkm3SoSlTpvDb3/62qEGiYmsSfo8foLPJqa9BAmDybJhxMrxyM2xe3T8FNMYYV65Jh0qhYoNEqiYRSUQgHoFEtO9BAuBf/gs8fnj8IijS5B/GmMrT06RDpVDxQSKWiOU3blNvhjbAsTfAB6/A63f0QwmNMaZz0qGjjz4a6DrpUEtLCxdeeCFvvvkmN9xwQ1HuX7F9EgGvM/xGJBGBhNvhvDM1CYD9z4Zlf4Dnr4GpR0L97jt3PWOM6cHIkSO5/fbbi3qPiq1JdAkS4R5mpSuEiDO9qb8KHr8QEvGdLKUxxgysig0SqeamaCKa31wS+ardDT53M3z0Brxyy85fzxhjBlDFBolUTSKa7OcgAbDXKbDXF+DlH8O6f/TPNY0xZgBUbpDwpDU39UfHdabj/xuqR8JjFzpvTxljTBmq2CDRtbmpn/ok0lWPgBN/CRuWOkOKG2NMGarYINGl47q/m5tSdj8GDr7IeSX2Hw/277WNMaYEKjZIdOu49vicmen629HXwsTD4A/fgnVv9f/1jTGmiCo2SHR0XKeCRLDWeYW1v3n9cNrdUDUcHjwbdmzq/3sYY0yRlCxIiMixIrJcRFaKyOU58nxRRJaKyBIRub+Y5enW3NTfTU3phoyGL/0Otq2DR/4Nkoni3csYY/pRSb64FhEvcBtwNNAELBSReaq6NC3PNOAK4FBV3Swio4tZpvRRYPs04VChGmbCcTfBk9+Gl26AI64s7v2MMf3uxtdv5N1N7/brNaePmM5ls3ofPXrOnDn84Ac/4Oijj+bKK6+ktbWVW2+9tV/Lkk2phuWYBaxU1VUAIjIXOAlYmpbnfOA2Vd0MoKobilmgrjWJ1uLWJFI+da7zkd2ffwJj94U9Tyj+PY0xu4SBmnSoVEFiPLAmbb8JOCgjz+4AIrIA8ALXqOoz2S4mIhcAFwBMmDChTwXyiQ9BOpubhhS14uIQcb6f2LDMaXb68uMw8ZDi39cY0y/y+Yu/WNInHXrppZfwer08/vjjPPXUU2zYsIGLLrqIY445pt/vW6o+iWw9wpmTQvuAacBngTOA34jIsGwXU9U7VXWmqs6sr6/vW4FEnNnp0juuS8EfgjMfdEaNfeBL8PGS0tzXGFPWUpMOBYPBjkmHTj75ZH79619zzz338OCDxXnNvlRBoglIn0C6AVibJc8TqhpT1X8Cy3GCRtEEvIHSBwmAmlHw5cfAXw2//wJs+bB09zbGlJ3eJh26/vrrueiii4py71IFiYXANBGZLCIB4HQgs0HtcWAOgIiMwml+WlXMQgW9wdK83ZTNsAlw9qMQ2wG/+1do21ja+xtjykJPkw6pKpdddhnHHXccBxxwQFHuX5Igoapx4GLgWWAZ8JCqLhGR60TkRDfbs0CLiCwF5gP/oaotxSxXwBsgGg9DvB2CdSSSSjRewlnlxsyAMx6ErU1w32kQ2V66extjykJPkw794he/4Pnnn+fhhx8u2rwSJZt0SFWfBp7OSLsqbVuB77pLSQS8ASKxNmcnWMsVj77FwtWbeejrh1BfGyxNISYeAqfdA3PPgrlnwBlzIVBTmnsbY8raJZdcwiWXXFLUe1TsF9fgNDdF452z0n20pZ1/bmzjnLtepzUcK11B9jgO/vV2WP0K/O4UCG8t3b2NMaYHFR0kAp4A0VhnkGiPJhhTF+S9j7dxwf8uIhwr4ZfR+3zRqVF89AbcewK0FbWlzRhj8lLZQcIbIJJWk2iPJdl7/DB++sV9eXXVJr41900Sycw3dYtoxklw+v2w4V2453Ow7ePS3dsYY7Ko6CDhfCfhTggUrCMcS1AV8HLSfuO5+oQZPLvkY658/G2c7pIS2f0YOPth57XYu4+FLWt6P8cYU1Ql/R3Qj/qj3BUdJJzvJFJBotYJEn7nn+Srh07m4jmf4IHX1/DjZ94t7Q/J5NnwlcedJqe7/sWGGDdmAIVCIVpaWsouUKgqLS0thEI7NwVCyd5uGowC3gCRRNTZCdbSHktQ5fd2HP/eMbuzeUeUO15exfsbtvPfp+3LsOpAaQrXOAu++hTc/yUnUPzrHTDjxN7PM8b0q4aGBpqammhubh7oohQsFArR0NCwU9eo6CAR9AaJJt23mNyO61BakBARrj95L3YfU8v1Ty3l+J//hV+ceQCfmji8NAXcbW84fz7MPRMe+jLMuRJmX1qceS+MMVn5/X4mT5480MUYMNbclIwDQtJXQySe7BIkwAkU53x6Eo/8+6fxeoUv3fE37nj5fZKl6tCuHQPnPgX7fAnmX+8MDJh6I8sYY4qsooNE0BskonEI1hFJOL/0qwLerHn3aRjGk988nKNnjOGGP77Lefcu5KMtJfpl7Q85zU1HXg3vPAJ3Hw+bV5fm3saYilbRQSLgCRDVZEd/BNClTyLT0Co/vzrrAK498ZP89f0WjvzpS/z8+RWl+Z5CBA7/rvOKbMv7cPvh8PbDxb+vMaaiVXaQ8AaIkESDQzqCRMjf8z9Jqvnphe99hiOnj+GW59/jyJ++zNNvryvN2w/Tj4cL/wL10+GR8+Dxb9iYT8aYoqnoIJGawjTudloD3fokcmkYXs1tZx3AA+cfTG3Ixzfu+ztn/vo1Fq/ZUrTydhg+Eb76R5j9fVh8P9wxG9YuLv59jTEVp6KDRMcUpoGajiajnpqbsjlk6kie/OZh/OdJn+Td9a2cfNsCzrtnIW83FXn8Ja8PjvghnPMHpyP7N0fB/P8H8Uhx72uMqSgWJIBIsDNI5FuTSOfzevjyIZP48/fncOkxu7Pog82c8MtXOP9/F7FkbZGDxeTD4d8XwCdPhpdvhP85FFYvKO49jTEVo6KDRKq5KRqo7uy4zvF2Uz5qQ34uPmIaf7lsDt85andeXdXC5259ha/ds5C/rtxYvD6L6hHwhd/A2Y9AIgL3HA/zvgntm4tzP2NMxajoIJGqSUT91R19EoU2N2VTF/LzraOm8cplR/Cdo3bnraYtnPmb1zj+1ld45I2m4k1s9Imj4BuvwqcvgTfvg1/OctbJEo5ma4zZpVR2kBDng/OIv4qw+4u7L81NuQyt6gwWN31hHxLJJN/7v39w2I0v8rPn32NtMb6zCNTAMf8JF8x3pkh94htOx/bKF/r/XsaYXV5FB4mgOoEh6g8Sjub3CmxfhPxevnhgI89+ezb/+7VZTB9bx89fWMGhN77IuXe/zjPvrCeW6Ofaxdh94d+eh1Pvgkgr/P4UZ0Kj9e/0732MMbu0ih67KZBwAkPUF8zrY7qdJSLM3r2e2bvXs2bTDh5atIaHFq3hwt+/waghQU45YDwn7juOT46rQ/pjfCYR2OsLMP3z8Pqv4c83we2Hwd6nOWNA1e+x8/cwxuzSCv6zWURqRKTg36QicqyILBeRlSJyeZbj54pIs4gsdpd/K/QehQom4gBEfIF+6bguROOIar53zB4suOwIfnvOTPZrHMZdr/yTz//iFY786cvc/Nx7rNzQTx/J+YLw6YvhksXO+t0n4baD4KGv2DDkxpge9VqTEBEPcDpwFnAgEAGCItIMPA3cqaorermGF7gNOBpoAhaKyDxVXZqR9UFVvbjwx+ibgDtMeNQbILzDbW7ylSZIpPi8Ho7ccwxH7jmGLTui/PGd9cxbvJZfvLiCW19YwfTdajlmxhiOnrEbe43fyRpG9Qg45no49Nvw6q+c2sXSJ2D3Y+Gw70DjQTbCrDGmi3yam+YDzwNXAO+oOg35IjICmAP8WEQeU9Xf93CNWcBKVV3lnjsXOAnIDBIlFXA/PIt4fbTHEgR8HjyegfslOaw6wBmzJnDGrAl83BrmybfW8ew76/nl/JXc+uJKdqsLcdSM0Ry55xgOnjyy77WemlFw5FXOW1Cv3+kEjLv+xenHmHWB00Tlr+rfhzPGlKV8gsRRqhrLTFTVTcAjwCMi4u/lGuOB9Hk4m4CDsuT7gojMBt4DvqOqWefuFJELgAsAJkyY0PsT5BCMp2oSPsLRRFH7Iwo1pi7EeYdN5rzDJtOyPcL85c08v/RjHv37R/z+1Q8JeD3MnDScw6aN4vBP1PPJcXWFB7iqYfCZ78PB34C3HnRqFk9cBH/6ERzwFTjwPOcNKWNMxeo1SGQLEH3Ik+23V+aXZX8AHlDViIhcCNwLHJHjfncCdwLMnDmzz1+oBWNhAKIeb7dZ6QaTkUOCnPqpBk79VAPhWILX/rmJV1Y085cVG7npmeXcxHKGV/uZOWkEsyaN4MDJI/jkuDr83jy7nIJDnIAw82uw+i/w2h3w11thwc+dL7r3PdOZFS9QU9wHNcYMOvn0SWyj+y90cH7xq6rW5XGfJqAxbb8BWJueQVVb0nZ/DdyYx3V3it8NEhHxEI4lS9ZpvTNCfi+f2b2ez+xeD0DztggLVm5kwcqNvL56E88t/RiA6oCX/ScM44AJw9m3YRj7NA5ldG0vc92KOPNrT54NW9Y4gwf+4354/EJ4+lKYcRLsezpMPBQ8g//fyhiz8/KpSdT2w30WAtNEZDLwEU5H+JnpGURkrKquc3dPBJb1w317FHRneIt6hPZYgqCv/D4bqa8NcvL+4zl5//EAfNwaZuHqTSz85yZeX72Z2+avJDWJ3vhhVezTMJS9xg9lxrg6ZoytY3RtMHtn+LBG+OxlTnPUh39zAsaSx2HxfVBT77xWO+MkmHS4M9igMWaXVJL/d6tqXEQuBp4FvMBdqrpERK4DFqnqPOASETkRiAObgHOLXa5gdAcA0WSccCxRFjWJ3oypC/H5fcbx+X3GAbAjGmfJ2lb+sWYLi9ds4R9NW/jjO+s78o+sCbDn2Dqm71bLtDFD+MRoZ10XcruZRGDip53luJvgvWdg2Tx46yF4426oGuHMcbH7sTD5MxDKp2JpjCkXeQcJEZkJ/BCY6J6Xam7aJ5/zVfVpnFdm09OuStu+AucNqpLxRbYhqkQSEdoHWcd1f6kO+Dhw0ggOnDSiI21re4x317WydF0ry9a1smRtK7979QMiaWNK7VYXYuroGiaPqmHSSGc9eVQNjXuejH+vU5zhyVe+4LxCu3QevPl78PhgwiEw7RiYdrQzMZK9UmtMWSukJnEf8B/A20CRRqgrLYluJ4AQTUQJxxMMrertJa1dw9AqPwdNGclBU0Z2pCWSStPmHaz4eDsrNmxnxcfbeL95O/MWr6U1HO/I5xEYO7SKhuFVNI5opHH492k86gfsHltK48YF1DbNx/Pcj+C5H8GQMTDpMKdJavJsGDHFgoYxZaaQINHsNgvtOiLbCCAdNYn+HNyv3Hg9wsSRNUwcWcNRM8Z0pKsqm3fE+OfGNlZvbOODljbWbG5nzaYdvLJiIx9vC9M5AvqhwKHMqG7l2KolHKTvsOe7L1H3ziMAhKvGEB47C2mcRdWUgwmM3w98gZI/qzEmf4UEiatF5DfACzhfXQOgqo/2e6lKJbKNoM9DNBklHEtWdJDIRUQYURNgRE2AT00c3u14JJ7go83trNsaZu0WZ71uaztvbJnKH7f9Cxva2xka+YBDPEs5JLGUA95fwLhVf4CXIYyfFZ5P8EFoOutrprN12AwSI6YyvKaKYdV+hlUHGFrld7ar/NRV+e2/kTElVkiQ+CowHfDT2dykQFkHiYDPQzQRdb6TCJTf200DLejzMqV+CFPqh+TME40n2bg9wvrWMEu2R3m9+QOC695g2KY3GdP6NsfseIrAjsegGXZokKU6kaXJiSzWRpYnG3hPG2jFuX7A56Eu5KMu5Ke2yk9dyMeQoLuEfNS665qgj5qAj+qAl5qgs65296sCXqoDXkI+74B+YW9MOSgkSOyrqnsXrSQDIbKNQG01kUTEebvJ/kotioDPw7hhVYwblhrqYwzOSC2uRBw2vgfrFlO1djH7rV3MAR//DU/suY4sbcHRtFRPYUOgkY+8DXwg41mlY1ndPox1W8O0ReJsD8fZHo1TyASAIb+HKr+XkN9Lld9L0O8l5PcQ8nkJZqwDPg9Bn8ddO/sBn4eAVzq2/V5nCXhT24LPXXfsezz43H2fp3Pf63HSPEL/jAJsTD8oJEi8KiIzsgzKV74irQSlrqMmYU0ZA8TrgzEzYMwMZL8znR9KVdjaBBuWwYal1GxYRk3zu0xo+SMzo2mj4/qrYfgkGDsJhk8mOWwikdoJ7Kgex7bgWLZrkLZInLZonB3RBDuiCdqjCdpjznY45uyHY05aOJYgEk8SjiVoa4t32Y/Gk86SSBJLFGkqWpfP4wSN1NpZPBn7gkdw1+lpXdNTxzyp/CKIm8cjgsdDR75UgJLUMXfdmebu0zWfAB6Ps6ZLHhCcfJItLbXvBsXUsc7trvlTUtdKz0damdxidGRIz+sck27nZjve9UJd79M9rft10nW5Zg/5csn2h8MJ+4wt+h8UhQSJw4BzROSfOH0SBb0CO+ioOn0SXj/t8Qiq/TsrndlJIs4HfcMaYfdjOtNVYds62LgCWlZAy/uwebWzrHoJT2wHVUAVMBIgNAyGutepHQt1Y5316LFQN855A6tqeMFvXSWTSjSRJOIGjliiM4BE40niSSWWSBJz0+IJdz+pxBNO/nhSSSSVWEJJJJPuWt1093hCO/IlNH0/SUKdcqSOJVPrjPR4Mkkk7qark67utiod5ypOWipd1XnrTXHO1bT8mrafOq/jmgBp29qx7oefC9PF5/ceW/QXBvMZluMQ4FXg2OIWpcRiO0CT+D0Bwu5osNbcVAZEnF/udeNgyme6HlOF7RucgLF1jVMTSa03f+B8Od6+ufs1PX4YMtpZakbDkHqoHgnVo5wRc6tHufvDnYASHIrH4yHk8dofFgVKDyzpwUPRjiDSsaZrXifNPZ6Rp/NYZz7Srp2emPVYl/t2LW/m8a7P01mOHvN1T+py7Z7y5bomlOaN8nxqEufgzAXxHvAM8Iyqru/5lDIQ2QZA0BtgUypI7AJfXFc0Eagd4yxZBxnG+Qhw23qnNtK61gkq2z921m0bnLT1b0HbRkjmGLdSPE4NpWo4hIY6o+mGhnYuwTp3qU1bhkAgtXYXT+W9KJHeHJR93E8z2OQzdtOFACIyHTgOuEdEhuLMM/EMsEBVE0UtZTG4QSLgCxFxJx8qxvzWZpDxV8GIyc7SE1VnbvAdLdDW4qzbN7vLJme9Y5OTp32LMyBieCuEt4D789RRW8djAAAZjklEQVQrXxUEqsFf466rnZF2/VXuUu2sfVXgDzlrX9BNC7lLwF0HwRt09r1Bd9+fsR1wak0VGJxM3+XdJ6Gq7wLvAreISBXOhEOnATcDM4tTvCKKtAIQ9FYRiWwBrLnJpBHprBmMmFLYufGI80dIpBXCrc462gaR7RB1l8h2iLVBdIfT9Bnd7m63O7WaWLu7uGnxMPTX32IenxMsvAHnpYEu26n9tG2PL20/tXidtbhrjydt2+tue51aV2o/fdvjcfa7pHvyXMRd3H0k41h6Gmnb0st2xhrSzs+WJ9sxslwj8zhd07vsk3Fetv20tLpxRW9z6tMAf6rajjMO09O95R20UjUJf4joDqe5ydqXTb/wuX+914zq3+smYk6wiIUh3g7xqLMfj0Ai4gSTRMzZjkfddQSScad2k4i614g4TWkJNz0Zc8+LOXmT8bTtGCQTzrU70hNOwErlTSadtSacY8k4aDItn7vWXWI0n8Hl6i1Fv0Wh80lkhqx855MYfDqam6qJus0DVpMwg5rX7yzB/hi9f4Akk50BI+muU/uqaWmpdE3bT3bm0ySgnddBM9Kz7WfbTnb2hKdey+roJda0YxlryCONrmkd25n3IfvxbPtd0kqjVPNJDD6pjutADXG3g9JqEsYUmccDWJ9IOSmouUlE9gUOd3f/rKpv9X+RSiQVJPxDiCXdmoS93WSMMV3kHdJF5Fs4w4WPdpf7ROSbxSpY0bkd1/5ADTGNAmrNTcYYk6GQmsR5wEGq2gYgIjcCfwN+UYyCFV1kG3iDBP3VuN+MWnOTMcZkKKRxUID0d/ASlPPXMJFtEKwl6A06+564fSdhjDEZCqlJ3A28JiKPufsnA7/t/yKVSLgVgrX4Pc5sdCIxq0kYY0yGvP90VtWbga8Bm4DNwFdV9WeF3ExEjhWR5SKyUkQu7yHfqSKi7rzaxZFRk/D5Evi9VpMwxph0Bb3dpKpvAG/05UYi4sUZA+pooAlYKCLzMoceF5Fa4BLgtb7cJ2+RbRCsI+B1ps8M+W2ISmOMyVTI200zReQxEfm7iLwlIm+LSCGvwM4CVqrqKlWNAnOBk7Lk+0/gJiBcwLULl1GTCPjta1BjjMlUSE3iPuA/gLfpnL60EOOBNWn7TWQM1Ski+wONqvqkiFya60IicgFwAcCECRP6UBScV2CDtR01iaAFCWOM6aaQINGsqvN24l7Z3oTqaOMREQ9wC3BubxdS1TuBOwFmzpzZt3aiLz8GXj+BHU0ABPzlN5CtMcYUWyFB4moR+Q3wAs7MdACo6qN5nt8ENKbtNwBr0/Zrgb2Al9zp+HYD5onIiaq6qIBy5mfkVACCkWYAAtYnYYwx3RQSJL4KTAf8dDY3KZBvkFgITBORycBHwOnAmamDqroV6Bg2U0ReAi4tSoBI0/F2kzdezNsYY0xZKiRI7Kuqe/f1RqoaF5GLgWcBL3CXqi4RkeuARTvZlNVnqe8k/D5rbjLGmEyFBIlXRWRG5iurhVDVbnNQqOpVOfJ+tq/3KURnTcKChDHGZCokSBwGnCsiq3D6JARnPol9ilKyEukMEvZ2kzHGZCokSBybJa3se3tTr8B6rCZhjDHdFBIk6oEfAhMzzivrmkQqSHg91nFtjDGZSvkx3aCUam6ymoQxxnRXyo/pBid1Rn4Vq0kYY0w3pfyYblAKx5No0oeIBQljjMlUyo/pBqVwNAHqA4kNdFGMMWbQKdnHdINVeyyBqg/E+iSMMSZTIbPsvCoiM4pWkgESjiUhaTUJY4zJptCP6c4RkX+yC31M1x5zmpuSWJAwxphMO/sxXdlrjzrNTRYkjDGmu16DhIiIOj7oLU//Fq00wrEEqB/F3m4yxphM+fRJzBeRb4pIlyngRCQgIkeIyL3AOcUpXvGFYwk06SOh0YEuijHGDDr5NDcdC3wNeMCdC2ILEMIZ7vtPwC2qurh4RSyuVJ9EXK25yRhjMvUaJFQ1DPwK+JWI+HEmBmpX1S3FLlwpdAYJq0kYY0ymQjquUdUYsK5IZRkQqY7rhIYHuijGGDPoFPKdxC4pEne+k4glrSZhjDGZKj5ItEcTgI9owoKEMcZkKjhIiEiNiHiLUZiB0B5L4JMAUatJGGNMN70GCRHxiMiZIvKUiGwA3gXWicgSEfmJiEzL92YicqyILBeRlSJyeZbjF4rI2yKyWEReKcUwIE6Q8FtNwhhjssjrOwlgKnAFsJuqNqrqaOBw4FXgxyJydm8XcWsftwHHATOAM7IEgftVdW9V3Q+4Cbg5/0fpm3Asgc8TIJKIUKbfAxpjTNHk83bTUaoaE5GJqtoxI52qbgIeAR5xX43tzSxgpaquAhCRucBJwNK0a7am5a+hBHNodwQJTRLXOP68HsUYYypDrzUJ97VXgMcyj4nIwRl5ejIeWJO23+SmZV7zIhF5H6cmcUke190p7dEEAY8zz3UsYR/UGWNMunz6JL4oIj8GakVkz4xO6zsLuJdkSetWU1DV21R1KnAZcGWOMl0gIotEZFFzc3MBReiuPdYZJCKJSC+5jTGmsuTTJ7EAp0loOE4fwQoR+buIPAm0F3CvJqAxbb8BWNtD/rnAydkOqOqdqjpTVWfW19cXUITuwrEkAZ8FCWOMySafYTk+Av5XRN5X1QUAIjICmIzzplO+FgLT3PGfPgJOB85MzyAi01R1hbv7OWAFRRaOJajyBAHsDSdjjMlQyFDhC1Jpbqf1psw8PV1HVeMicjHwLM7ggHep6hIRuQ5YpKrzgItF5CggBmymBKPLtscSDPcFIW5BwhhjMuXzdtN8EXkEeEJVP0wlikgAd7Y6nNdk7+ntQqr6NPB0RtpVadvfyq/Y/ac9miDkBolI0pqbjDEmXV+HCq/C6c8o+6HCw7EEIbdPwmoSxhjTVcUPFR6OJQn5QoAFCWOMyVTwUOEi8u+AV0QWA4vTOprLTjyRJJpIUu13goS93WSMMV0VFCTA6UMQkTHA/sAXRGSqqp7f/0UrvnDc+YC82m/NTcYYk03eQUJEfgZ8x33T6WPgGXcpW+FYAoCagNUkjDEmm0KGCt8OzBORGgAROUZEFvRyzqDmzCVBR3OT1SSMMaarvGsSqnqliJwJvCQiEaAN6DbcdzlJ1SSGBKoACxLGGJOpkOamI4HzcYLDWOA8VV1erIKVQjjm9EnUWnOTMcZkVUhz0w+BH6nqZ4FTgQdF5IiilKpE2t2aRG3IrUnY7HTGGNNFIc1NR6Rtvy0ix+HMJ/HpYhSsFNo7mpts7CZjjMmm4DmuU1R1HXBkP5al5FId11UBHwF3djpjjDGd+hwkAFS1kKHCB51I3A0Sfi8Bb8BqEsYYk2GngkS566xJOEHCahLGGNNVZQcJt08i5PMS9AatJmGMMRksSODUJCxIGGNMdxUdJFLfSQR9HvxevzU3GWNMhgoPEglCfg8iQtATtO8kjDEmQ0UHifZogiq/F8DebjLGmCwqO0jEugYJa24yxpiuKjpIhGMJQgEnSFjHtTHGdFeyICEix4rIchFZKSLdRo8Vke+KyFIReUtEXhCRicUukzO/tTU3GWNMLiUJEiLiBW4DjgNmAGeIyIyMbG8CM1V1H+Bh4KZil6s9lqAqrSZhzU3GGNNVqWoSs4CVqrpKVaPAXOCk9AyqOl9Vd7i7rwINxS6UdVwbY0zPShUkxgNr0vab3LRczgP+mOugiFwgIotEZFFzc3OfCxWOJQmlgoQnQCRpNQljjElXqiAhWdI0a0aRs4GZwE9yXUxV71TVmao6s76+vs+FSn0nAdZxbYwx2eQ9n8ROagIa0/YbgLWZmUTkKJzJjT6jqkX/sz7zFVgLEsYY01WpahILgWkiMllEAsDpwLz0DCKyP3AHcKKqbihFodI7rgPeAAlNEE/GS3FrY4wpCyUJEqoaBy4GngWWAQ+p6hIRuU5ETnSz/QQYAvyfiCwWkXk5Ltdvwmk1iaDXZqczxphMpWpuQlWfBp7OSLsqbfuoUpUFIJlUwrEkwbTmJnCCRLW/upRFMcaYQativ7iOxJ0RYDNrEvathDHGdKrYINExl4T7dlN6TcIYY4yjYoNEOG3CIUgLEjZcuDHGdKjYINExdWmqucljzU3GGJOpcoNEtGuQsOYmY4zprmKDREdzU0aQsJqEMcZ0quAg4b7dFLDvJIwxJpeKDRIdfRI+CxLGGJNLxQeJqoDzT+D3+gFrbjLGmHQVGyTCGR3X9jGdMcZ0V7lBIt614zoVJGLJ2ICVyRhjBpuKDRK5XoG1moQxxnSq3CCR8TFdwGNBwhhjMlV0kAj4PHg9zqR5qZpELGHNTcYYk1KxQSISS3b0RwB4xIPf47eahDHGpKnYINEe7ZzfOiXoDVqQMMaYNJUbJNJmpUuxea6NMaarig4SoWxBwoYKN8aYDhUbJMKxRMe4TSnW3GSMMV2VLEiIyLEislxEVorI5VmOzxaRv4tIXEROLXZ5wrFEx7hNKX6P35qbjDEmTUmChIh4gduA44AZwBkiMiMj24fAucD9pShTu9UkjDGmV74S3WcWsFJVVwGIyFzgJGBpKoOqrnaPJUtRoPZo947roDdo30kYY0yaUjU3jQfWpO03uWl9IiIXiMgiEVnU3Nzcp2uMqQsxdmioS1rAG7CahDHGpClVTUKypGlfL6aqdwJ3AsycObNP17n//IO7pQ0NDmXZpmW0tLcwsmpkX4tnjDG7jFLVJJqAxrT9BmBtie6dt6/v83XaY+1c/derUe1zDDPGmF1GqYLEQmCaiEwWkQBwOjCvRPfO27Th0/juzO/yctPLPLT8oYEujjHGDLiSBAlVjQMXA88Cy4CHVHWJiFwnIicCiMiBItIEnAbcISJLSlG2TGdOP5NDxx/KTxb9hFVbVg1EEYwxZtCQcm9WmTlzpi5atKhfr7mxfSOnPHEKY2rGcN/x93WMEGuMMbsCEXlDVWfmk7div7juyaiqUVz76Wt5d9O7/PLNXw50cYwxZsBYkMhhzoQ5fHH3L3L3krt5dd2rA10cY4wZEBYkenDpgZcyqW4S33rxW1zxlyt44cMXCMfDA10sY4wpmVJ9J1GWqnxV3Hbkbdz51p3MXzOfJ1c9SZWvisPHH86RE45k5m4zGV09eqCLaYwxRWMd13mKJWMsXL+Q5z94nhc+fIFN4U0AjB8yngNGH8D+Y/Zn//r9mTx0Ml6Pt5erGWPMwCmk49qCRB8kkgmWtizlzQ1vdiwt4RbAqX3sMXwPpo+YzoyRM5g+YjpThk0h6A2WtIzGGJOLBYkSU1XWbFvD4ubFLGtZxtKWpSzfvJy2WBvgzJ/dWNvIlKFTmDpsKlOGTmHy0Mk01jYyNDh0QMtujKk8hQQJ65PoByLChLoJTKibwIlTTwQgqUnWbFvDspZlvL/1fd7f4ix/afoLcY13nDs0OJQJtRNorG1k/JDxjBsyjnE14xg7ZCy71exGla9qoB7LGGMsSBSLRzxMrJvIxLqJXdJjiRgfbvuQD1o/YM22NXzY+iEfbvuQxRsW88zqZ0hq15HShweHU19dT311PWOqx1BfVU99VT0jq0Y6S2gkI0IjqPHXIJJtHEVjjOk7CxIl5vf6mTpsKlOHTe12LJ6Ms2HHBta1rWPt9rWsb1vP+rb1bGjfQPOOZt7b9B4t4ZZugQScuTCGBYcxPDScocGhDA8666HBodQF6jqXYB21gVqG+Ic4S2AIPo/9GBhjsrPfDoOIz+NzmpuGjONTYz6VNU88GWdTeBObwptoaW+hJdzirNtb2BLZwtbIVjZHNrOsbRlbIlvYFt2WNaikq/JVUeOvodpX7az91R37Vb6qbkvIF+pYh7whQr4QQW+wYwn5QgS8AYLeIAFPAJ/HZ7UcY8qUBYky4/P4GF09Ou/vM5KapC3WRmu0ldZIK63RVrZHt7Mtto22WBvbotvYHt3O9th2dsR20BZvoy3WRvOOZnbEd9Aea6c97izpfSmFEISAN0DAE8Dv9RPwBvB7/J37biDxe/xd1t0W6brvFW+XbY94um17xNOx7/V4O7cz1qklc19E8IoXQbqkeUjbFg8enO1UvtR2Kq+IdKa5eYCO81N5xZ16JX3fAqwZSBYkdnEe8VAbqKU2UMv4IX2eDBBw+lN2xHcQSUQIx8OEE2FnHQ8TSUS6LOF4mFgyRjQRJZKIEE1GiSacJZaMdRyLJWLENEY8ESeWdK4fTzrb6et4Mk5c4ySSiY7teLJvQatcdQkcCM7/JGdAyRZ43I3ONEkLShlp6dfJTO92Tfe6mWm9npMlX7Z795a3kGvmc24hdiaI7+y9HzvpsY4/OIrFgoTJm9/rZ6h38Lyyq6okNEFSk8STcRKacIKIxklqsiM9qcku+VLHUmmZ26lFURLJBEncfVWSuOu0PB3HNNlxXFU7j9F1H+iWnjofwE3B+Z92y9eRR7V73vTjafsd/2YZeVL/jqnt9H/bLnnSjmd7bT79Wplp2a7X03/TbOf3mNZDeXqVJVve5xaYt9u5ZfL5gQUJU7ZEBJ84P8I2nLsxxWED/BljjMnJgoQxxpicLEgYY4zJyYKEMcaYnEoaJETkWBFZLiIrReTyLMeDIvKge/w1EZlUyvIZY4zpqmRBQkS8wG3AccAM4AwRmZGR7Txgs6p+ArgFuLFU5TPGGNNdKWsSs4CVqrpKVaPAXOCkjDwnAfe62w8DR4p9bmqMMQOmlEFiPLAmbb/JTcuaR1XjwFZgZOaFROQCEVkkIouam5uLVFxjjDGl/JguW40g85PDfPKgqncCdwKISLOIfNDHMo0CNvbx3MFmV3mWXeU5wJ5lsNpVnmVnnmNi71kcpQwSTUBj2n4DsDZHniYR8QFDgU09XVRV6/taIBFZlO/sTIPdrvIsu8pzgD3LYLWrPEupnqOUzU0LgWkiMllEAsDpwLyMPPOAc9ztU4EXtVwGODHGmF1QyWoSqhoXkYuBZwEvcJeqLhGR64BFqjoP+C3wOxFZiVODOL1U5TPGGNNdSQf4U9Wngacz0q5K2w4Dp5WwSHeW8F7Ftqs8y67yHGDPMljtKs9SkucQa80xxhiTiw3LYYwxJicLEsYYY3KqyCDR2xhSg5mI3CUiG0TknbS0ESLynIiscNfDB7KM+RKRRhGZLyLLRGSJiHzLTS+75xGRkIi8LiL/cJ/lWjd9sjsO2Qp3XLKymB1JRLwi8qaIPOnul+tzrBaRt0VksYgsctPK7ucLQESGicjDIvKu+/+ZQ0rxLBUXJPIcQ2owuwc4NiPtcuAFVZ0GvODul4M48D1V3RM4GLjI/W9Rjs8TAY5Q1X2B/YBjReRgnPHHbnGfZTPO+GTl4FvAsrT9cn0OgDmqul/aNwXl+PMF8HPgGVWdDuyL89+n+M/SMR9vhSzAIcCzaftXAFcMdLkKfIZJwDtp+8uBse72WGD5QJexj8/1BHB0uT8PUA38HTgI54tYn5ve5WdvsC44H7q+ABwBPIkzEkLZPYdb1tXAqIy0svv5AuqAf+K+bFTKZ6m4mgT5jSFVbsao6joAdz16gMtTMHdY+P2B1yjT53GbaBYDG4DngPeBLeqMQwbl87P2M+D7QNLdH0l5Pgc4w/r8SUTeEJEL3LRy/PmaAjQDd7vNgL8RkRpK8CyVGCTyGh/KlI6IDAEeAb6tqq0DXZ6+UtWEqu6H85f4LGDPbNlKW6rCiMjngQ2q+kZ6cpasg/o50hyqqgfgNC9fJCKzB7pAfeQDDgD+R1X3B9ooUTNZJQaJfMaQKjcfi8hYAHe9YYDLkzcR8eMEiPtU9VE3uWyfB0BVtwAv4fSzDHPHIYPy+Fk7FDhRRFbjDOd/BE7NotyeAwBVXeuuNwCP4QTvcvz5agKaVPU1d/9hnKBR9GepxCCRzxhS5SZ9zKtzcNr2Bz13rpDfAstU9ea0Q2X3PCJSLyLD3O0q4CicjsX5OOOQQRk8i6peoaoNqjoJ5/8bL6rqWZTZcwCISI2I1Ka2gWOAdyjDny9VXQ+sEZE93KQjgaWU4lkGukNmgDqBjgfew2kz/uFAl6fAsj8ArANiOH9dnIfTZvwCsMJdjxjocub5LIfhNFu8BSx2l+PL8XmAfYA33Wd5B7jKTZ8CvA6sBP4PCA50WQt4ps8CT5brc7hl/oe7LEn9f70cf77ccu8HLHJ/xh4HhpfiWWxYDmOMMTlVYnOTMcaYPFmQMMYYk5MFCWOMMTlZkDDGGJOTBQljjDE5WZAwxhiTkwUJY4wxOVmQMKYPRKRBRL6U41iViLzsDkuf7XhARP6cNsyFMYOWBQlj+uZInLFzsvka8KiqJrIdVNUoztexWYOMMYOJBQljCiQihwE3A6e6M55NzshyFu4YOu74QU+5M9a9k1b7eNzNZ8ygZtVdYwqkqq+IyELgUlV9J/2YO2jkFFVd7SYdC6xV1c+5x4e66e8AB5aoyMb0mdUkjOmbPXBmBcs0CtiStv82cJSI3Cgih6vqVnDmngCiqVFKjRmsLEgYUyARGQlsVdVYlsPtQCi1o6rvAZ/CCRY3iMhVaXmDQLiYZTVmZ1lzkzGFm0yOSXdUdbM7jWlIVcMiMg7YpKq/F5HtwLnQEWiacwQaYwYNq0kYU7h3gVFuR/Snsxz/E85cGQB7A6+7c1//ELjeTZ8DPF30khqzk2w+CWP6mYjsD3xXVb/cQ55HgStUNVu/hjGDhtUkjOlnqvomML+nj+mAxy1AmHJgNQljjDE5WU3CGGNMThYkjDHG5GRBwhhjTE4WJIwxxuRkQcIYY0xOFiSMMcbk9P8BP8AfRfh7qiwAAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEZCAYAAAC5AHPcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8HMWZ8PHfM6fuw5J8yraMMdiYcBgBTsxhDCHcV9hgcsFuEhYWwm5Idgk5d9kckGyS5U1CCMmyJBsCJNiAkxDM6YQjgA+cGNuAD4wtfMiSrfsYzczz/tE90mg8kkZjzWhkP9/PZ+juqurualnMo+rqrhJVxRhjjBkuz2hXwBhjzNhkAcQYY0xaLIAYY4xJiwUQY4wxabEAYowxJi0WQIwxxqTFAogxgIj8u4j8arTrMZih6igi20TknGzWyRzeLICYw4aIXCsi60SkQ0R2i8hPRKRstOtlzFhlAcQcFkTk88CdwL8CpcB8YDrwtIgEslQHXzbOY0y2WAAxhzwRKQH+A/isqj6pqj2qug34CE4Q+bhbNE9EHhaRVhFZIyLHxx3jVhF5z817S0TOdtM9IvJFEdkiIo0i8hsRGefm1YiIisinRGQ78JyIPCkiNyXU768icoW7fpeI7BCRFhFZLSKnJ1zOgHVMOOZg9coTkV+56U0islJEJhzkj9kchiyAmMPBB4A8YGl8oqq2AX8EPugmXQr8FhgH/Bp4TET8InI0cBNwsqoWAx8Ctrn73AxcBpwJTAb2Az9OOP+ZwBx3v18DV8cyROQYnCD2BzdpJXBCXB1+KyJ5ccdKWsck1zxYva7BaYVNBSqA64HOJMcwZlAWQMzhoBJoUNVwkrxdbj7AalV9RFV7gO/jBJ35QAQIAseIiF9Vt6nqFneffwS+rKp1qtoN/DtwZcLtqn9X1XZV7QQeBU4Qkelu3seApe6+qOqvVLVRVcOq+j33vEfHHWugOiYarF49OIHjSFWNqOpqVW0Z+sdoTH8WQMzhoAGoHKAPYpKbD7AjlqiqUaAOmKyqm4F/wfkSrheRh0Rkslt0OvCoeyuoCdiIE3DibwnFH7cVp7Wx2E1aDDwQyxeRz4vIRhFpdo9XSl+AG7COSa5rsHr9H7AceEhEdorIdwZoxRgzKAsg5nDwF6AbuCI+UUQKgfOBZ92kqXF5HqAa2Amgqr9W1dNwvpgVp0MenC/081W1LO6Tp6rvxZ0qccjrB4GrReT9QD7wvHvO04FbcfpmylW1DGgGJG7fAeuYYMB6uX1A/6Gqx+Dc3rsI+GTyH50xA7MAYg55qtqM04n+QxE5z+3XqMHpS6jD+Ysc4CQRucJtqfwLTtB5RUSOFpFFIhIEunD6CyLuPvcA34zdkhKRKhG5dIgqPYETiG4HHnZbEgDFQBjYC/hE5GtAScK+SeuY5BwD1ktEzhKR94mIF2jBuaUVSXIMYwZlAcQcFlT1O8CXgP/C+dJ8Feev9LNj/Q/A48BVOB3OnwCucPsagsAdOLe6dgPj3WMB3AUsA54SkVacL/NTh6hLN06H/jk4HeExy3E69d8G3sUJVjsSdh+ojokGq9dE4BH357AR+BOQ0y9RmtwkNqGUMcaYdFgLxBhjTFosgBhjjEmLBRBjjDFpsQBijDEmLYf04G6VlZVaU1Mz2tUwxpgxY/Xq1Q2qWpVK2UM6gNTU1LBq1arRroYxxowZIvJuqmXtFpYxxpi0WAAxxhiTFgsgxhhj0pK1PhARuQ9n0LZ6VT02Sf6/4gxtHavXHKBKVfeJyDagFWe8nrCq1man1sYYM7Cenh7q6uro6uoa7aoMW15eHtXV1fj96Q/EnM1O9PuBHwG/TJapqt8FvgsgIhcDn1PVfXFFzlLVhmT7GmPMaKirq6O4uJiamhpEZOgdcoSq0tjYSF1dHTNmzEj7OFm7haWqfwb2DVnQcTXOkNfGGJOzurq6qKioGFPBA0BEqKioOOiWU871gYhIAXAesCQuWXFGFV0tItcNsf91IrJKRFbt3bs3k1U1xpgxFzxiRqLeORdAgIuBlxJuXy1Q1Xk4k//cKCJnDLSzqt6rqrWqWltVldK7MCnZubmJLWvqiUSiQxc2xpjDQC6+SLiYhNtXqhqbFa5eRB4FTgH+nM1KrfjVm+zf3UFhaYBjz5zCMadNoaAkkM0qGGNMTsmpACIipcCZwMfj0goBj6q2uuvn4szkljUaVZobOpk2dxwovLrsHVY+sY1ZJ03ghA9Oo7K6KJvVMcaYnJC1W1gi8iDO3NRHi0idiHxKRK4Xkevjil0OPKWq7XFpE4AXReSvwGvAH1T1yWzVG6C9uZtoWJlxfBUX33wCH/33U5l72hS2rt3Lb7+1kr89vwObmMsYMxrWrVvHggULerfXrFnDokWLsnLurLVAVPXqFMrcj/O4b3zaVuD4zNQqNS0NnQCUVOYBUD6xkDMWH8UpF8/g2V9s5IWHN7F7awsLP3Y0gbycatQZY7LkP363ng07W0b0mMdMLuHrF88dtMzcuXPZsmULkUgEr9fL5z//eb73ve+NaD0GYt92KWje6zzqVlKZ3y89r9DPBde/j9XL3+W1ZVtpqGvj/H88lvKJhaNRTWPMYcjj8TB37lzWr1/Ppk2bmDZtGvPmzaO9vZ1/+qd/IhAIsHDhQj72sY8NfbBhsgCSgpaGTkSgeFzeAXniEWrPr2FCTQlP/c96fvvtVSz65ByOPGn8KNTUGDNahmopZNL8+fN56aWXuPvuu3nySecO/9KlS7nyyiu5+OKLueqqqzISQHLxMd6c09LQSVF5Hl7fwD+uqXPGcdWXT2bc5EKW/+wNnvu/jYS6wlmspTHmcDV//ny+8pWvcPnllzNlyhTAeUt+6tSpAHi93oyc1wJICloaOimpOrD1kaioPI/LPz+PeedN582Xd/HwN15j15bmLNTQGHM4mz17NsFgkFtvvbU3rbq6mrq6OgCi0cy8v2YBJAXNDV0H9H8MxOvz8P7LZnLZ5+cB8Oh/reaVx7YQCdsLiMaYzLjrrrv49re/TWFhX//rFVdcwZIlS7jhhhu4+OKLM3Je6wMZQk93hM6WUMoBJGbykWVc9eVTePG3m1j95Lts37CPhR87mvHTSzJUU2PM4WbLli1ceOGFLFiwgGuuuaZfXmFhIf/7v/+b0fNbABlC7BHe0mEGEIBAvo9Fn5xDzXGVrPj1W/z2jlUce8YUTr3kCPIK0x9C2RhjAGbOnMmbb745aue3ADKEvndAhh9AYo44oYopR5fz2rKtrFtRx+bV9XzgipnMnj8J8YzNgdiMMcb6QIbQ0uC+A5JCJ/pggvk+Tr/qKD7y5ZMpG1/Ac798k6X/tYbdW62T3RgzNlkAGUJzQyf+PO+I3XKqrC7mii/MY9En59C8t4Ml31nNH+7+Gw11bSNyfGOMyRa7hTWEloZOSirzR3TMf/EIcz4wiZnzqvjb83W8/tR2Hv7Ga8yqHc8pFx9B2YSCETuXMcZkigWQIbTs7czY0CSBPB+159dw7BlTWPv0dv763A42r9nLrJPHc8I506iaWpyR8xpjzEiwADIIjSotjV1MP7Yio+fJK/Qz/7KZHLdoKmuWv8v6F3fy9qt7qJ5dzokfnMbUY8aN2VnPjDGHLgsgg+hoCRHpiR7UE1jDUVAS4LS/m0XtBTVseHEnf31uB7/74V+pmFLI+xZWM+vkCTbarzEmZ9i30SCaY++AVGUngMTkFfqZ96HpHH/2VDat3MPaZ7az4oG3eOmRzcyqHc+c0yYzoabEWiXGmFFlT2ENYiTeATkYXp+H2e+fxFVfOYUP/9tJzDxpPG+v3MOSO1fz8DdeY+0z22nb3z0qdTPG5IbDYkKpsah5bycIFFcc3DsgB0tEmHhEKROPKOX0v5vF2yv3sOHFnbz0yGZeWrKZyUeWMevkCcw8sYr8Ypun3ZhR8ccvwu51I3vMie+D8+8YtIhNKJWjnGHcg4MO455tgXwfx54xhWPPmML+3e1sXl3PppV7+NOv3+LPD73NlKPKmHF8JTXvqxy1lpMxJnsGmlBq69atfPOb36S5uZlHHnkkI+e2ADKIlr1daY2BlS3lEws5+cIZ1F5Qw76d7WxauYeta/fywsObeOHhTZRPKmTGcRVMP7aSCUeU4PXmTiA05pAzREshk5JNKHXEEUfwP//zP1x55ZUZO68FkEG0NHRm/BHekSAiVEwpomJKEfMvm0lTfQfvrmtk27oG1j69gzXLt+MLepl8ZBnVs8upPrqcyuoiG4fLmEPE/Pnzufbaa7nxxht7J5TKhqwFEBG5D7gIqFfVY5PkLwQeB95xk5aq6u1u3nnAXYAX+LmqZjzU94QidKQxjHsuKBtfQNnZBRx/9lS6O8O89+Z+6t7cR91b+3l5SSMAwUIfk2aWMfGIEibNLGP89GJ8gczMWmaMyaxkE0plQzZbIPcDPwJ+OUiZF1T1ovgEEfECPwY+CNQBK0VkmapuyFRFIe4JrIMcRHG0BfN9HHFiFUecWAVA2/5u3ntrH3VvN7F7SzPb/tYAgMcrVE0rZvz0EndZTPnEAjx228uYnJdsQqnGxka+/OUv8/rrr/Ptb3+b2267bcTPm7UAoqp/FpGaNHY9BdisqlsBROQh4FIgwwHEHYV3DLZABlNUHuTo+ZM4ev4kADrbQuze2sLuLc3s3trMm3/ZxboVzjSYPr+HiuoiKqcWUzG5kIophYybXGRzmRiTIwabUKqiooJ77rkno+fPtT6Q94vIX4GdwBdUdT0wBdgRV6YOOHWgA4jIdcB1ANOmTUu7Ii17059IaizJLwow47hKZhxXCUA0qjTXd1D/bit7322lfnsLm1buYX1nuHefwrIg4yYXUjahgPIJBZRNdJaFZUF7udGYLLIJpfqsAaarapuIXAA8BswCkn0j6UAHUdV7gXsBamtrByw3lJaGTvxBL3lFh9df2x6PUD6xkPKJhRx96kQAVJX2pm4a32un8b02Gne2sX9XB29u2UVPd6R3X1/AQ0llPqVV+ZRUxj55FFfkUTwuz4ZhMeYQkzP/R6tqS9z6EyJyt4hU4rQ4psYVrcZpoWRUJoZxH6tEhKLyPIrK8/o9leYElhBNe9pp2tNB055Omhs6ad7byY4N+wj3RPsdJ1jgo2icE0wKy4IUlQUoLAtSWBrsXQYLfPZ0mDFjRM4EEBGZCOxRVRWRU3CGWWkEmoBZIjIDeA9YDHw00/VpbuiibPyhffvqYDmBJUhReZDq2eP65akqHS0hWhq6aNvXRav7advXRWtjF7u3NNPV3nPAMT0eIb8kQEFJgPziAAXFfvKKA+QX+ckr8rvLAHmFPvIK/QQLfNbRb8woyeZjvA8CC4FKEakDvg74AVT1HuBK4AYRCQOdwGJVVSAsIjcBy3Ee473P7RvJGFWlpaGTaXPHDV3YJCUiTsuiNAgzS5OWCfdE6GgO0dbUTXtTNx0tITpaQnS6y46WEPt2tdHV2nNAayZeIM9L0A0mwXwfgXyfu+7Hn+8lkOcjkOclkO/DH3S2/XnevvWgF1/AY61NY4Ypm09hXT1E/o9wHvNNlvcE8EQm6pVMbBj3Q70DfbT5/N7evpKh9HRH6GwN0dnWQ3d7D13tPXS1h+nucNa728N0d4YJdYZpaeiku8PZ7umODNJjFkfAF/DiD3jcgOLt3fYFvPj8ztIb8Djrfg9ev5Pu9Xnwuksn3Vn3+MTJ83nwuuser7Peu/R58HgFj0csgJkxJ2duYeWS2BNYh9ojvGOZP+jFH0wt2MTTqNLTHSHUFSbU6Sx7uiP0dEXo6Q67eRF6QhHC3RF6QlF3GSEcihAORelsDRHuifZuh3uihHsiRMNpP6ORlMcjeNzgEgsqHm/s46SJJy7d42yLuy0SW9IvTzzgkfhtwSP0bcfWxSkrkjwNictL3Hb7rfrn9S8HOMdCeh+Nie3vJMWnO/8R6Fe2t0gsLZbQu51sPxIKxD2ZEzt/vMRj9xU9QCQc7fcgSf+DDy6jfy6I8wdaplkASaJvGPex/RKhcb7YAu5tLcpH9tgaVcLhKJFQlEjY+YR73PWeKNFIlEiPEon05UcjSjSsvevOMkokom6eWyaqcUsnTd00jcblR5z8cKgv3Vk69dOoohpL58BtjZWJKw+ptdoMJ3+8nP2720e7GgfweIXK6sxPiW0BJInmhi4QKKmwFogZmHgEf8CL/xAcAkbVDSpugEHp23aDDQpKXyAiLh/cMnHH6pfv7hsLVPHpsW1nmaxsbyUTgl38sRPKxJfT3tJx6wl5vSdJ3KffD4k29gxrwrkD6pIh2bobagEkiZaGTorKgnj99nSPOTz13lpCnEdXTFIbNzYQLBjdd8XWrVvH9ddfz0svvQQ4E0p94Qtf4Lnnnsv4uS2AJBF7B8QYY1J152t38ua+kX0rfPa42dx6yuADJNqEUjmmZW8nU+fm/jDuxhgz0IRSjz32GH/4wx+or6/nxhtv5Nxzzx3xc1sASRCNKmUTCqiamvkOKGPMoWOolkImJZtQ6rLLLuOyyy5j//79fOELX7AAkg0ej3DZLfNGuxrGGJOywSaU+sY3vsGNN96YkfNaL7ExxoxxySaUUlVuvfVWzj//fObNy8wfxdYCMcaYMS7ZhFI//OEPeeaZZ2hubmbz5s1cf/31I35eCyDGGDNGDTah1M0338zNN9+c0fNbADHGmDFqtCeUsj4QY4wxabEAYowxJi0WQIwxxqTFAogxxpi0WAAxxhiTFgsgxhhj0mIBxBhjTFosgBhjjElL1gKIiNwnIvUi8sYA+R8Tkb+5n5dF5Pi4vG0isk5E1orIqmzV2Rhjct26detYsGBB7/aaNWtYtGhRVs6dzTfR7wd+BPxygPx3gDNVdb+InA/cC5wal3+WqjZktorGGJOe3d/6Ft0bR/at8OCc2Uz80pcGLXNYTCilqn8WkZpB8l+O23wFqM50nYwxZqwbaEKpjRs3ctddd9HQ0MDZZ5/NDTfcMOLnztWxsD4F/DFuW4GnRESBn6rqvQPtKCLXAdcBTJs2LaOVNMaYmKFaCpmUbEKpOXPmcM899xCNRvnMZz6TkfPmXCe6iJyFE0Dip/daoKrzgPOBG0XkjIH2V9V7VbVWVWurqqoyXFtjjBl98+fP5ytf+QqXX355vwmlli1bxmmnncbZZ5+dkfPmVAARkeOAnwOXqmpjLF1Vd7rLeuBR4JTRqaExxuSeZBNKAVxyySW8/PLLPPDAAxk5b87cwhKRacBS4BOq+nZceiHgUdVWd/1c4PZRqqYxxuScZBNKrVixgqVLl9Ld3c0FF1yQkfNmLYCIyIPAQqBSROqArwN+AFW9B/gaUAHcLSIAYVWtBSYAj7ppPuDXqvpktuptjDG5arAJpRYuXMjChQszev5sPoV19RD5nwY+nSR9K3D8gXsYY8zhzSaUMsYYMyZZADHGGJMWCyDGGGPSYgHEGGNMWiyAGGOMSYsFEGOMMWmxAGKMMSYtFkCMMcakxQKIMcaMYYfLhFLGGHPIeuE3b9Owo21Ej1k5tYjTP3LUoGVGc0Ipa4EYY8wYFj+h1JIlS3onlAJob2/npJNO4ve//31Gzm0tEGOMGQFDtRQyKdmEUgB33nknH/nIRzJ2Xgsgxhgzxs2fP59rr72WG2+8sXdCqWeeeYZjjjmGrq6ujJ3XAogxxoxxySaUev7552lvb2fDhg3k5+dzwQUX4PGMbK+FBRBjjBnjkk0o9c1vfhOA+++/n8rKyhEPHmCd6MYYM2Zt2bKF2bNn09nZecCEUjHXXnstF110UUbOP+wWiDutbJeqRjJQH2OMMSnK+QmlRMQjIh8VkT+ISD3wJrBLRNaLyHdFZFbmq2mMMSbXpHIL63lgJnAbMFFVp6rqeOB04BXgDhH5eAbraIwxJgelcgvrHFXtSUxU1X3AEmCJiPhHvGbGGGNyWioB5LMiEr+tQAPwoqq+A5AswCQjIvcBFwH1qnpsknwB7gIuADqAa1V1jZt3DfAVt+g3VPUXqZzTGGNMZqRyC6s44VMC1AJ/FJHFwzzf/cB5g+SfD8xyP9cBPwEQkXHA14FTgVOAr4tI+TDPbYwxZgQN2QJR1f9Ilu5+qT8DPJTqyVT1zyJSM0iRS4FfqqoCr4hImYhMAhYCT7u3zRCRp3EC0YOpntsYY8zISvs9EPfLXIYsODxTgB1x23Vu2kDpxhhjRknaAUREFgH7R7AukDwg6SDpBx5A5DoRWSUiq/bu3TuilTPGGNMnlfdA1onI3xI+dcAdwI0jXJ86YGrcdjWwc5D0A6jqvapaq6q1VVVVI1w9Y4zJLbk+oVTiO/AKNKpqewbqswy4SUQewukwb1bVXSKyHPhWXMf5uTjvpRhjTE54/v57qX9364gec/z0Izjr2usGLTOaE0ql0on+brJ0EVkAfFRVU26FiMiDOB3ilW4r5uuA3z3PPcATOI/wbsZ5jPfv3bx9IvKfwEr3ULfHOtSNMeZwFj+h1KZNm3onlFqxYgVf/epXmTt3LosXL2bhwoUjfu5hjYUlIicAHwU+ArwDLB3O/qp69RD5ygC3xVT1PuC+4ZzPGGOyZaiWQiYlm1BKRCgqKqKrq4vq6uqMnHfIACIiRwGLgauBRuBhQFT1rIzUyBhjzLAkm1Dq9NNP58wzz2TPnj3ccsstPPDAAyN+3lRaIG8CLwAXq+pmABH53IjXxBhjTFqSTSgVm/+jvLyc7u7ujJw3lQDyYZwWyPMi8iTOi4Mj/f6HMcaYNCWbUGrp0qUsX76cpqYmbrrppoycN5VO9EeBR915QC4DPgdMEJGfAI+q6lMZqZkxxphBbdmyhQsvvJAFCxYcMKHUFVdcwRVXXJHR86fcie4+tvsA8IA7jMnfAV8ELIAYY8woGAsTSh1wu0pV96nqT1V10UBljDHGHNpSmlBKRD4rItPiE0UkICKLROQXQPLJeI0xxhyyUrmFdR7wD8CDIjIDaALyAC/O7asfqOrazFXRGGNyl6oyFm/COK/dHZxUOtG7gLuBu92ZByuBTlVtOuizG2PMGJaXl0djYyMVFRVjKoioKo2NjeTl5R3UcYb1Jro78+CugzqjMcYcIqqrq6mrq2Msjvydl5d30G+oDyuAGGOM6eP3+5kxY8ZoV2PUpD0fiDHGmMPbsAOIiBSKiDcTlTHGGDN2pPIeiEdEPioifxCRepyxsXaJyHoR+a6IzMp8NY0xxuSalN4DAWbiTOA0UVWnqup44HTgFeAOEfl4ButojDEmB6XSiX6OqvaIyHRVjcYS3QmdlgBL3Md7jTHGHEaGbIG4j+4CPJqYJyLzE8oYY4w5TKTSB/IREbkDKBaROQkd6PdmrmrGGGNyWSq3sF7CGbrk08D3gaNFpAnYCXRmsG7GGGNyWCpDmbwH/FJEtqjqSwDucO4zcJ7IMsYYcxhKeTj3WPBw1/ep6mp3jpCUh3MXkfNE5C0R2SwiX0yS/wMRWet+3nZbOrG8SFzeslTOZ4wxJnNSuYX1vIgsAR5X1e2xRBEJAKfhDOX+PHD/YAdx+05+DHwQqANWisgyVd0QK6Oqn4sr/1ngxLhDdKrqCSnU1xhjTBak8h7IeUAEZzj3nSKyQUTeATYBV+MM535/Csc5BdisqltVNYQzt/qlg5S/GngwheMaY4wZBdkczn0KsCNuuw44NVlBEZmO08fyXFxynoisAsLAHar62DDPb4wxZgQNezh3EbkB8IrIWmCtqm5Kcfdk/SQDzWiyGHhEVSNxadNUdaeIHAE8JyLrVHXLAScRuQ64DmDatGmJ2elrbwSNQlHVyB3TGGPGsGEPpqiqXwP+H9AKfFhEfpbirnXA1LjtapxHgZNZTMLtK1Xd6S63Aivo3z8SX+5eVa1V1dqqqhH8sn/oo/D92fDbv4ftr8IIzOZljDFjWcoBRESeEZHjAVR1j6o+qap3qOpnUjzESmCWiMxwO+AXAwc8TSUiRwPlwF/i0spFJOiuVwILgA2J+2aMKtRvgPIa2Pws3Hcu3HsmvP4A9HRlrRrGGJNLhtMC+TfgByLyvyIyabgnUtUwcBOwHNgI/EZV14vI7SJySVzRq4GHtP+EvXOAVSLyV5wnvu6If3or4zr3Q3cL1P4D3LIBLvw+hLvh8X+CH58MO1/PWlWMMSZXyHAnVheRDwNfA5YC31HVnH0bvba2VletWjWsfVSVZ7c/S3VxNbPHzXYS31sNP1sEi38Nsy+MFYQtz8Gym6G9Hi74Lsy7BsbQvMjGGJNIRFaram0qZYfVB+K+MPgW8BPgs8AmEfnE8KuYu0SEL734JX635Xd9ifu3OcvymviCcOTZ8I9/hukL4Hf/DI/fBD05G0+NMWZEDacP5EXgPeAHOI/kXgssBE4RkUNqUMXSYCnN3c19CbEAUjb9wMKFFfDxJXDmrbD2V/DzD8K+rVmppzHGjKbhPMZ7PbBeD7zn9VkR2TiCdRp1ZcGyAwNIYRUEi5Lv4PHCWV+CKbWw9DPw04Vw4X/B+/7ObmkZYw5ZKbdAVPWNJMEj5sIRqk9OKA2W0tQd957k/m39b18N5Khz4R//BONnO4HkkX+Ajn2ZqqYxxoyqYb8Hkoz7bsYhoyxYll4AAafctU/Aoq/CxmXwkwWw5fkM1NIYY0bXiASQQ02/W1iRHmiuSz2AAHh9cMYX4NPPQrAY/u8y+OOt0N2WkfoaY8xosACSREmghOZQM1GNQvMOZwiT4QSQmMknOLe0Tr0eXr0HfnwKbHjc3mI3xhwSLIAkURYsI6pR2nrakj/COxz+fDj/TviHpyB/HPzmk/CrK6Bh80hV1xhjRoUFkCTK8soAaO5qPvgAEjPtVLhuBZz/HahbBT95Pzx7u93WMsaMWRZAkigLOgGkqbvJCSDeABQPe/SWA3l9cOo/wk2rYO4V8ML34K7j4S9325haxpgxxwJIEiWBEiAugJRNd971GCnFE+CKnzqd7BPmwvLb4IfzYPX9Tqe9McaMARZAkoi1QJpDzbDvnYO/fTWQ6lq4Zhl88nGnhfO7f3Y62l//lTNYozHG5DALIEn0BpCupuG9A5KuIxbCp5+Bqx8CfyE8fiP89/ucW1yd+zN7bmOMSZMFkCSKA8UIQlP7bmcY9/Ia1j69gpd/+3jmTioCR58P178An3hCSSmrAAAaMUlEQVTUubX17O3w/bnOOySNB0y+aIwxo2pYU9oeLrweLyXBEppa6pyE8hpe/u536QwIf33wZ0w59Rwu+NyN+AL+kT+5CMxc5Hx2vwF/+RGs/LnzHknN6XDStTD7IvDnjfy5jTFmGKwFMoDSQCnN7fXORnkNYY+SF4oSJsKmNc9wz+JL+fUtt9JU35i5Skw8Fi6/B/7lDVj0FWjaDks+5Uyt+8cvwu519lKiMWbUWAAZQFmwjOYuJzh0B8fT4xUKyqdw3a+WUnPMB/CosOu99dz/T5/gp1dewW+/+p/s3ro9M5UpmQRn/CvcvBY+8ZjTZ7Ly53DPafDjU2HFndCwKTPnNsaYAdgtrAGUBktp3P8OFFSyc8deECFYWk6wIJ8Pf/1LRCIRnr/vV7zz3JN0hprY/var/PqLr5Af9lA0eSYzzzyL2kvOIxAMjlylPB6YeZbzaW+EDY/CG4/Cim/Dim/BxPfB3MvhqPNh/BwbSt4Yk1EWQAZQFixja6QTymvYvckZbLigsqo33+v1cs5nroHPXEM41MNLDz/KlmefpLN5N/V7N1P/yGZee/in5Kmf/Mpqqk+dz0kXXUDZ+IqRqWBhBZz8aefTsgs2PAZvLHE63p+9HcqmwVHnwVEfcvpOfCMYyIwxBgsgAyoNltKkYSivYd8O59ZUyeTJScv6An7O/MRHOPMTHwHg7dfW8Lff/4H9mzcQ6mymoeldGpa/y1+ffIhgjxIIllA0dQY18z/A8ecupKC4+OAqWzIJ5t/gfFp2wqan4O3lsOb/4LV7wV8A094PR5wJM86Eicc5rRljjDkIFkAGUOovpl2gp2wabWuczvSqGTUp7XvUKfM46pR5vdt1b25m9WPLaHx7A6GmejrCLbRsX8fO7ev4y0P3EAwr/kAh+ZWTGT/nWI464wNMmzsbrzeNt99LJjtPap10rTM/+zsvwOanYeuf4OmvOWXyy51WybT3w9RTYdJx4M3AE2XGmENaVgOIiJwH3AV4gZ+r6h0J+dcC38WZex3gR6r6czfvGuArbvo3VPUXmaxrmft0U3PxBLr2vwHAlKNnpXWs6tlHUv3FW3q3I5EI6//0Fzb9+U80v7OZUFsjXeE2Whu2UP/CFt544XF8kSj+iAd/sJC8yomUz5hJ9XHHcdT8k1JvsfjznVkSjzrX2W7ZBe/8Gd75kxNYNi5z0n35MGWeE0ymzINJJ0BptfWhGGMGlbUAIiJe4MfAB4E6YKWILFPVDQlFH1bVmxL2HQd8HagFFFjt7pux17TLQs7ghs2F4+hpa8EbiTJu0vgRObbX6+W4Radx3KLT+qXv2Pg2G577Ew1vv0nH3p2Eu1rp7GmlZW879Xu38NZrT/Hsz5RgOIoPP968IoKl4yieXE3lkUdSc8JxTD5q5sAtl5JJcPxVzgec2107XoUdrznLl/8fRMNOXkGFE0gmn+B0zo8/BsbNdAaENMYYstsCOQXYHJv+VkQeAi4FEgNIMh8CnlbVfe6+TwPnAQ9mqK6UdLUC0JRfTKSrA180U2fqM3XOUUydc9QB6XVvbmbzK69Rv+lt2nfVEWrdTyTcRXeoiZZ9Lezdt42tb7zIa4+BJxrFHwEvXrz+fPyFJeSNq6SwaiLjpk9l4qxZTJ0zi2BBvnO7a+7lzgecW1571sPO12HnWti1Fl78b9CIk+8NQtVRMH6uM+97xSyoPMoZ6sUXyPwPyBiTU7IZQKYAO+K264BTk5T7sIicAbwNfE5Vdwyw75RkJxGR64DrAKZNm5Z2Zcva9wHQ5PURjXTj1dHrdK6efSTVs49Mmle/vY6tq9ey5623ad25g679jUS62oiEu+nuaaW1rR3ad8OON2CNu5Mq/ojijYJHvHi8Qbx5BfgLigmUlJBfNo6i8adQPv0yKt5fxYTiLvJb34H69bBng3ML7G8P9VVCvE4QqTgSxs2A8hnOdnkNlE93bqUZYw452QwgyW6oJ75G/TvgQVXtFpHrgV8Ai1Lc10lUvRe4F6C2tjbt17TLWp2O85aeNqIawefNzcdgx0+rZvy06gHzuzs62b7+TXZv3sL+7XW0791N1/5Gwh2tRHq6iEbC9EQ66OzsIBLaB01AkvchvZEo3qjiUcEjs/B4ZuP3QNCr5Hsi5Hu7KfLUUchG8j1dBP09BH1h8gIR/MXFBKom4a+sxls+3elfKZkExZOdZdFEa8EYMwZlM4DUAVPjtquBnfEFVDV+XJCfAXfG7bswYd8VI17DOGXNOyHPmRMk7IFgoCCTp8uYYEE+s04+kVknnzhk2ZbG/eze8g6NO96jZddu2hr30tW0n562VsJdHURDXUTDIVTDhCMhQlGlJSpEvB4IB4EgMC75wd8C0SZ80X14o6vxaRQvik+j+IjgBTweQcSDeL2I1we+AOIPIoE8CBZCsAgpKMGbX4QvLw9/fj7+/DwCefkECwsIFBQQLMgnWFRIsCCfguJigvl5+IKB9J5oM8YMKpsBZCUwS0Rm4DxltRj4aHwBEZmkqrvczUuAje76cuBbIlLubp8L3JbJyubvfxf/pFL27asHnwdf4UG+qzEGlFSUU1JRDnGPIKci1N3Nvl172L9zD20N+2jf10hHczPdLc10t7YS7uok0t1FpLuLaE8Iwl3OxFnRMD1ECKFEwf0IkagQVQ/RaARCHdDeAexL/8JU8SiIKqIgxJYAgrjtVEGc/0psS9wn0QSR2LoH8cTSPCDu0iMgHicAioDH05fv8SCevrLidct5PG45QTxexBM7lpsHfemxcu5+seN5PF7E69bLK73n93i8ONWXfnX0ePvq4Yldk8c9jkfcc3j7HsDrrb/703KvzYMHce/qxuoWS4gdJ5bn7Bef56x73GMSlxY7j5Mf27cvTdyJ3WSAO8qeWEbccZznd/pzfmY4P4PecgceVJK8L+XxJLshkpwkm4guxbvhnhTf1Upab6+HyskTUzvRQchaAFHVsIjchBMMvMB9qrpeRG4HVqnqMuBmEbkECON8Y1zr7rtPRP4TJwgB3B7rUM+Izv1IVxOl3gm0ba+jCAiWlg+52+EqEAwysWYaE2vS73M6gCqh1kY692ynY88OQo276GnaS09TA9GOFiLtrUQ624mGOtHubqKhEBrpQSNRIupxPlEPUYWIepzApB6iKkQQInjddA8RnPzej4I6YQbVvsCmgAoQcfPFuY+qANKXrwhIX7ox2eYPR7h5yR8zfp6sPpOpqk8ATySkfS1u/TYGaFmo6n3AfRmtYMz+dwEoCxQTfrcBgMKqCVk5tXGJECipJFBSSemsYbSIIj3Q1QJdTdDVDKE26G510rpbnfldQu1OemzZ3QY9Hc52T4fT6ulpd55Ki4QO6jIi+FFvEJUAUU+QiPiIip8oAaL4iHj8KF5UvETVR1S8RJybe0TxEBUvqh4iUUHFi4oHjQpRcVtp6kERoiqoCBp1Al/UXY+6QVDd1pUiRKPqHAdQ7ctXVWfb+QdAcbZRJ1gCbtnYINBOGVSIuntpX0R18pC+zkrtGzy6X1rvurr1jCvj7hCf3m8EaumrL+o0KzUa3/WZpBu0tw5xeclGtU6SpsMa/Xrgcw+5Z8rnSV7OG8jOdA/2UH8y+7cBUJo3Dml0Huctq0760JfJNV6/M05Y4QiNORaNOIGkp9MJLuEuZz3c5a53QbgTwiGIdDtTEYe7e2/TeSPdTlALdzv5kbATlCIhJz2WHw276yFnPRp20yMQ7QHtcR6njkT68qNh0BF+vjyxwTRmG1AS1/pLWIc08pLlx6X1W41PS/YDHOI4qew7VLnCqgPTMsACSDJuACkrGI+v+U0gP+VhTMwhxuOFYJHzyUXRqBNYomE32LhLjSRZRpOkR/s+sXRVN83Nj0ZxmgoJZWNNCo3GLaP900nI783TvvX4MiTLS1zGLj4uDQ4sF5/Wu84w8uIk5h9QJtl+g7Vqhmj9DLrvEOUAAtn5fbUAksz+bVBQQWlBJYEO583sKXPSG8bEmIzyeACPjWVmRoUFkGT2vwPlNZQGS9FuxReJUlJeNtq1MsaYnGJjeiezfxuU11AWLMMf9uCNDtBMNMaYw5gFkETRiDNqrRtAPFEvHvsxGWPMAewWViKPF26rg0iI0j2vshP3rWhjjDH92DdjMr4A+AKUBkoJe3z4AtZBaYwxiezezCD87RDxetB8G+jPGGMSWQAZROc2Z2zHaJEFEGOMSWQBZBDtO/YCEC6zW1jGGJPIAsgg2nbtAaCr3LqKjDEmkQWQQbTvdSaVaq0YswMCGWNMxlgAGUR3834A9o+LjHJNjDEm91gAGUS4sxVfOEITbaNdFWOMyTkWQAYRCXXh1SjN3c2jXRVjjMk5FkAGEY2E8ODMi26MMaY/CyCDiEoU9Xpo6W4hErV+EGOMiWcBZACRSIQer0AgiKK0hlpHu0rGGJNTLIAMYN/OPUQ9HjyFBQA0h6wfxBhj4mU1gIjIeSLylohsFpEvJsm/RUQ2iMjfRORZEZkelxcRkbXuZ1mm67rz7c0A+MrKAesHMcaYRFl7xVpEvMCPgQ8CdcBKEVmmqhviir0O1Kpqh4jcAHwHuMrN61TVE7JV34Zt2wAoqJoIYE9iGWNMgmy2QE4BNqvqVlUNAQ8Bl8YXUNXnVbXD3XwFqM5i/fpp2bkTgMppTiPIWiDGGNNfNgPIFGBH3HadmzaQTwF/jNvOE5FVIvKKiFw20E4icp1bbtXevXvTrmxHYwMA046eC0BTlwUQY4yJl80AkmxAqaSTjYvIx4Fa4LtxydNUtRb4KPDfIjIz2b6qeq+q1qpqbVVVVdqV7W5pAlVmzD4Gj3isE90YYxJkM4DUAVPjtquBnYmFROQc4MvAJaraHUtX1Z3uciuwAjgxk5WNdLTij0TJy8unJFBifSDGGJMgmwFkJTBLRGaISABYDPR7mkpETgR+ihM86uPSy0Uk6K5XAguA+M73ERfp6cIbdRpNZcEy6wMxxpgEWXsKS1XDInITsBzwAvep6noRuR1YparLcG5ZFQG/FRGA7ap6CTAH+KmIRHGC3h0JT2+NuGg0jNfj/HhKg6UWQIwxJkFWZ0pS1SeAJxLSvha3fs4A+70MvC+ztesvIlF8vkLACSD1HfVD7GGMMYcXexM9CWcYEw++PCeAlAXLrA/EGGMSWABJov6d7ahH8BeXAnYLyxhjkrEAksSuTVsAyK+oBJwWSGe4k1AkNJrVMsaYnGIBJImGbe8CUDJxEuAEELC30Y0xJp4FkCRadu8CoGL6NABKgiWAjYdljDHxLIAk0bWvEYDJRx0JWAvEGGOSsQCSRKilCYkq42c4LZBYALEWiDHG9LEAkkS4qx1/JIrX6wWsBWKMMclYAEkiGu7Cq30/mpKA0wdiAcQYY/pYAEkiEg3j8fS9pJ/vyyfgCdDS3TKKtTLGmNxiASSJiEfx+vN6t0XEBlQ0xpgEFkASRCIRRMFXWNovfVz+OF7d9Sqrdq8apZoZY0xuEdWkczodEmpra3XVqpH5wn9l1yt89aWvsrt9Nx+q+RC3nHQLk4smj8ixjTEmV4jIanfyviFZCyRF8yfNZ9lly7jh+BtYsWMFlzx2CXevvZvOcOdoV80YY0aFtUDSsKttF99b/T2Wb1uOz+NjStEUqourqS6qZmrxVCYXTaYqv4rK/Eoq8yvJ8+UNfVBjjMkBw2mBWAA5CGv2rGFF3Qrea32PurY66lrraAkd+KRWsb+YivwKxuWNozyvnLJgGePyxlEWLKM0WEpJoKR3WRIsoThQTJ43D3dSLWOMyZrhBJCsTih1qJk3YR7zJszrl9bc3cyu9l00dDawt2Ovs+x0lk3dTbzb8i5ru9bS1N1ERCMDHtsnPgoDhRT5iygOFFPgK6DQX0iB3136Csj35VPgd5axT543jzyf8wl6g866N4+AN+Bse/PweXwWnIwxB80CyAgrDZZSGiwdslxUo7SGWmnpbqEl1EJzqJmWUAst3S20hlpp72mnNdRKW08bbaE22sPtNHY1sqN1Bx09HbSH2+no6UAZfgtSEILeIH6vn4AnQMDrfPwev/Px+vvWPX58Hl/vMnHdJz68Hu8B215x0rzixSMefB4fHvHgFSfP4/H05sWWsU/8tle8iEjv0oOnX9n4NBFxlvQtB0oTkX5pwIFp4vysevdxl8YYhwWQUeIRT8rBZiCqSnekm85wJ53hTrrCXc4y0kV3uJvOSCfd4W66I32fUCTUbz0UDRGKhOiJ9BCKhghHw/REe+iJ9tAV7qIl2tKbFo6GD/xo33o6wWwsig8ksQCTGGwSy8Vv9+7jJB4QwGJle/dLCFrx5fttx+1DkjiXuN8B+yRJGyxgJtt3oH0GKjuS+w5HJv4QyEQ901UWLOMX5/8i4+exADKGiUjv7apyyke7OkSiEaIaJaxhItEIEY0QjoaJapSIOtuRaISwholGnbSoRnvz49dVtf8S7ZcX1ShRon3rGu0to6ooffvHtmNlYvvEAl58er80d1tRUHrPFysT6z+M37c3iGpfem+ZuO34/ZKtJ0rMS6xDv7Ip7J/OvkOVG87+qe47nL9JUv0DJhN/6ORaX3JxoDgr58lqABGR84C7AC/wc1W9IyE/CPwSOAloBK5S1W1u3m3Ap4AIcLOqLs9i1U0KvB4vXrz48Y92VYwxWZC190BExAv8GDgfOAa4WkSOSSj2KWC/qh4J/AC40933GGAxMBc4D7jbPZ4xxphRks0XCU8BNqvqVlUNAQ8BlyaUuRSI3bh7BDhbnJuVlwIPqWq3qr4DbHaPZ4wxZpRkM4BMAXbEbde5aUnLqGoYaAYqUtwXABG5TkRWiciqvXv3jlDVjTHGJMpmAEn2iEJiz9NAZVLZ10lUvVdVa1W1tqqqaphVNMYYk6psBpA6YGrcdjWwc6AyIuIDSoF9Ke5rjDEmi7IZQFYCs0RkhogEcDrFlyWUWQZc465fCTynzvNxy4DFIhIUkRnALOC1LNXbGGNMEll7jFdVwyJyE7Ac5zHe+1R1vYjcDqxS1WXA/wD/JyKbcVoei91914vIb4ANQBi4UXWQcUCMMcZknA2maIwxppeNxusSkb3Au2nuXgk0jGB1RtOhci2HynWAXUsuOlSuAw7uWqarakpPIB3SAeRgiMiqVKNwrjtUruVQuQ6wa8lFh8p1QPauxWYkNMYYkxYLIMYYY9JiAWRg9452BUbQoXIth8p1gF1LLjpUrgOydC3WB2KMMSYt1gIxxhiTFgsgxhhj0mIBJIGInCcib4nIZhH54mjXZzhE5D4RqReRN+LSxonI0yKyyV2O/tSFKRCRqSLyvIhsFJH1IvLPbvqYuh4RyROR10Tkr+51/IebPkNEXnWv42F3eJ8xQUS8IvK6iPze3R6T1yIi20RknYisFZFVbtqY+v2KEZEyEXlERN50/595fzauxQJInBQnvcpl9+NMuBXvi8CzqjoLeNbdHgvCwOdVdQ4wH7jR/bcYa9fTDSxS1eOBE4DzRGQ+zmRpP3CvYz/OZGpjxT8DG+O2x/K1nKWqJ8S9MzHWfr9i7gKeVNXZwPE4/z6Zv5beOaPtA/B+YHnc9m3AbaNdr2FeQw3wRtz2W8Akd30S8NZo1zHN63oc+OBYvh6gAFgDnIrzlrDPTe/3e5fLH5yRsJ8FFgG/x5lqYaxeyzagMiFtzP1+ASXAO7gPRWXzWqwF0l/KE1eNIRNUdReAuxw/yvUZNhGpAU4EXmUMXo97y2ctUA88DWwBmtSZNA3G1u/ZfwP/BkTd7QrG7rUo8JSIrBaR69y0Mff7BRwB7AX+1721+HMRKSQL12IBpL+UJ64y2SEiRcAS4F9UtWW065MOVY2o6gk4f72fAsxJViy7tRo+EbkIqFfV1fHJSYrm/LW4FqjqPJxb1jeKyBmjXaE0+YB5wE9U9USgnSzderMA0t+hOHHVHhGZBOAu60e5PikTET9O8HhAVZe6yWP2elS1CViB06dT5k6aBmPn92wBcImIbAMewrmN9d+MzWtBVXe6y3rgUZzgPhZ/v+qAOlV91d1+BCegZPxaLID0l8qkV2NN/CRd1+D0JeQ8ERGc+WE2qur347LG1PWISJWIlLnr+cA5OB2cz+NMmgZj4DoAVPU2Va1W1Rqc/zeeU9WPMQavRUQKRaQ4tg6cC7zBGPv9AlDV3cAOETnaTTobZ+6kjF+LvYmeQEQuwPmrKjbp1TdHuUopE5EHgYU4QznvAb4OPAb8BpgGbAf+TlX3jVYdUyUipwEvAOvou9/+JZx+kDFzPSJyHPALnN8nD/AbVb1dRI7A+St+HPA68HFV7R69mg6PiCwEvqCqF43Fa3Hr/Ki76QN+rarfFJEKxtDvV4yInAD8HAgAW4G/x/19I4PXYgHEGGNMWuwWljHGmLRYADHGGJMWCyDGGGPSYgHEGGNMWiyAGGOMSYsFEGOMMWmxAGKMMSYtFkCMGWEiUi0iVw2Qly8if3KnDkiWHxCRP8cNDWJMzrIAYszIOxtnLKJk/gFYqqqRZJmqGsIZLj1pADIml1gAMWYEuUOwfB+40p3pbkZCkY/hjknkjsf0B3e2wjfiWi2PueWMyWnWTDZmBKnqiyKyEmecqDfi89wBOo9Q1W1u0nnATlW90M0vddPfAE7OUpWNSZu1QIwZeUfjzAaXqBJoitteB5wjIneKyOmq2gzO/CFAKDZarDG5ygKIMSPIHc21WVV7kmR3AnmxDVV9GzgJJ5B8W0S+Flc2CHRlsq7GHCy7hWXMyJrBABMqqep+d3rbPFXtEpHJwD5V/ZWItAHXQm8Q2jtAEDImZ1gLxJiR9SZQ6XaKfyBJ/lPAae76+4DX3PnSvwx8w00/C3gi4zU15iDZfCDGZJGInAjcoqqfGKTMUuA2VU3Wj2JMzrAWiDFZpKqvA88P9iIh8JgFDzMWWAvEGGNMWqwFYowxJi0WQIwxxqTFAogxxpi0WAAxxhiTFgsgxhhj0mIBxBhjTFr+P6m17jZSjeKJAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def plotStateTrajectories(rdata):\n", + " for ix in range(rdata['x'].shape[1]):\n", + " plt.plot(rdata['t'], rdata['x'][:, ix], label='$x_%d$' % ix)\n", + " plt.xlabel('$t$ (s)')\n", + " plt.ylabel('$x_i(t)$ (mmol/ml)')\n", + " plt.legend()\n", + " plt.title('State trajectories')\n", + " plt.show()\n", + " \n", + "def plotObservableTrajectories(rdata):\n", + " for iy in range(rdata['y'].shape[1]):\n", + " plt.plot(rdata['t'], rdata['y'][:, iy], label='$y_%d$' % iy)\n", + " plt.xlabel('$t$ (s)')\n", + " plt.ylabel('$y_i(t)$ (AU)')\n", + " plt.legend()\n", + " plt.title('Observables')\n", + " \n", + " plt.show()\n", + " sys.path.insert(0, 'test')\n", + "\n", + "plotStateTrajectories(rdata)\n", + "plotObservableTrajectories(rdata)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing likelihood" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Loglikelihood -60.649943\n" + ] + } + ], + "source": [ + "model = modelModule.getModel()\n", + "model.setTimepoints(amici.DoubleVector(np.linspace(0, 10, 11))) \n", + "solver = model.getSolver()\n", + "rdata = amici.runAmiciSimulation(model, solver)\n", + "\n", + "edata = amici.ExpData(model.get())\n", + "edata.my = amici.DoubleVector(rdata['y'].flatten())\n", + "edata.sigmay = amici.DoubleVector(np.ones(shape=rdata['y'].shape).flatten())\n", + "rdata = amici.runAmiciSimulation(model, solver, edata)\n", + "\n", + "print('Loglikelihood %f' % rdata['llh'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Forward sensitivity analysis" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " sx: [[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-2.00747250e-01 1.19873139e-01 -9.44167985e-03]\n", + " [-1.02561396e-01 -1.88820454e-01 1.01855972e-01]\n", + " [ 4.66193077e-01 -2.86365372e-01 2.39662449e-02]\n", + " [ 4.52560294e-02 1.14631370e-01 -3.34067919e-02]\n", + " [ 4.00672911e-01 1.92564093e-01 4.98877759e-02]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-2.23007240e-01 1.53979022e-01 -1.26885280e-02]\n", + " [-1.33426939e-01 -3.15955239e-01 9.49575030e-02]\n", + " [ 5.03470377e-01 -3.52731535e-01 2.81567412e-02]\n", + " [ 3.93630714e-02 1.10770683e-01 -1.05673869e-02]\n", + " [ 5.09580304e-01 4.65255489e-01 9.24843702e-02]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-2.14278104e-01 1.63465064e-01 -1.03268418e-02]\n", + " [-1.60981967e-01 -4.00490452e-01 7.54810648e-02]\n", + " [ 4.87746419e-01 -3.76014315e-01 2.30919334e-02]\n", + " [ 4.28733680e-02 1.15473583e-01 -6.63571687e-03]\n", + " [ 6.05168647e-01 7.07226039e-01 1.23870914e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-2.05888038e-01 1.69308689e-01 -7.93085660e-03]\n", + " [-1.84663809e-01 -4.65451966e-01 5.95026117e-02]\n", + " [ 4.66407064e-01 -3.87612079e-01 1.76410128e-02]\n", + " [ 4.52451104e-02 1.19865712e-01 -4.73313094e-03]\n", + " [ 6.90798449e-01 9.20396633e-01 1.49475827e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.98803165e-01 1.73327268e-01 -6.03008179e-03]\n", + " [-2.04303740e-01 -5.16111388e-01 4.68785776e-02]\n", + " [ 4.47070326e-01 -3.94304029e-01 1.32107437e-02]\n", + " [ 4.69732048e-02 1.22961727e-01 -3.35899442e-03]\n", + " [ 7.68998995e-01 1.10844286e+00 1.70889328e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.92789113e-01 1.75978657e-01 -4.54517629e-03]\n", + " [-2.20500138e-01 -5.55540705e-01 3.68776526e-02]\n", + " [ 4.30424855e-01 -3.97907706e-01 9.75257113e-03]\n", + " [ 4.82793652e-02 1.24952071e-01 -2.30991637e-03]\n", + " [ 8.40805131e-01 1.27504628e+00 1.89020151e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.87672774e-01 1.77588334e-01 -3.38318222e-03]\n", + " [-2.33807210e-01 -5.86081383e-01 2.89236334e-02]\n", + " [ 4.16201399e-01 -3.99295277e-01 7.06598588e-03]\n", + " [ 4.92546648e-02 1.26089711e-01 -1.50412006e-03]\n", + " [ 9.06806543e-01 1.42334018e+00 2.04522708e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.83320440e-01 1.78410042e-01 -2.47240692e-03]\n", + " [-2.44690164e-01 -6.09568485e-01 2.25774266e-02]\n", + " [ 4.04061655e-01 -3.99063012e-01 4.97908386e-03]\n", + " [ 4.99612484e-02 1.26581014e-01 -8.85891342e-04]\n", + " [ 9.67473970e-01 1.55589415e+00 2.17895305e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.79620591e-01 1.78640114e-01 -1.75822439e-03]\n", + " [-2.53540123e-01 -6.27448857e-01 1.75019839e-02]\n", + " [ 3.93704970e-01 -3.97656641e-01 3.35895484e-03]\n", + " [ 5.04492282e-02 1.26586733e-01 -4.13401240e-04]\n", + " [ 1.02322336e+00 1.67481439e+00 2.29524046e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.76478441e-01 1.78430281e-01 -1.19867662e-03]\n", + " [-2.60686971e-01 -6.40868686e-01 1.34365068e-02]\n", + " [ 3.84873835e-01 -3.95414931e-01 2.10369522e-03]\n", + " [ 5.07601805e-02 1.26231631e-01 -5.46465317e-05]\n", + " [ 1.07443160e+00 1.78183962e+00 2.39710937e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]]\n", + " sx0: [[0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0.]]\n", + " sigmay: [[1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1.]]\n", + " sy: [[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-2.00747250e-01 1.19873139e-01 -9.44167985e-03 -2.00747250e-01\n", + " 1.19873139e-01 -2.00747250e-01]\n", + " [-1.02561396e-01 -1.88820454e-01 1.01855972e-01 -1.02561396e-01\n", + " -1.88820454e-01 -1.02561396e-01]\n", + " [ 4.66193077e-01 -2.86365372e-01 2.39662449e-02 4.66193077e-01\n", + " -2.86365372e-01 4.66193077e-01]\n", + " [ 4.52560294e-02 1.14631370e-01 -3.34067919e-02 4.52560294e-02\n", + " 1.14631370e-01 4.52560294e-02]\n", + " [ 4.00672911e-01 1.92564093e-01 4.98877759e-02 4.00672911e-01\n", + " 1.92564093e-01 4.00672911e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.80072436e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-2.23007240e-01 1.53979022e-01 -1.26885280e-02 -2.23007240e-01\n", + " 1.53979022e-01 -2.23007240e-01]\n", + " [-1.33426939e-01 -3.15955239e-01 9.49575030e-02 -1.33426939e-01\n", + " -3.15955239e-01 -1.33426939e-01]\n", + " [ 5.03470377e-01 -3.52731535e-01 2.81567412e-02 5.03470377e-01\n", + " -3.52731535e-01 5.03470377e-01]\n", + " [ 3.93630714e-02 1.10770683e-01 -1.05673869e-02 3.93630714e-02\n", + " 1.10770683e-01 3.93630714e-02]\n", + " [ 5.09580304e-01 4.65255489e-01 9.24843702e-02 5.09580304e-01\n", + " 4.65255489e-01 5.09580304e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.60534516e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-2.14278104e-01 1.63465064e-01 -1.03268418e-02 -2.14278104e-01\n", + " 1.63465064e-01 -2.14278104e-01]\n", + " [-1.60981967e-01 -4.00490452e-01 7.54810648e-02 -1.60981967e-01\n", + " -4.00490452e-01 -1.60981967e-01]\n", + " [ 4.87746419e-01 -3.76014315e-01 2.30919334e-02 4.87746419e-01\n", + " -3.76014315e-01 4.87746419e-01]\n", + " [ 4.28733680e-02 1.15473583e-01 -6.63571687e-03 4.28733680e-02\n", + " 1.15473583e-01 4.28733680e-02]\n", + " [ 6.05168647e-01 7.07226039e-01 1.23870914e-01 6.05168647e-01\n", + " 7.07226039e-01 6.05168647e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.46870655e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-2.05888038e-01 1.69308689e-01 -7.93085660e-03 -2.05888038e-01\n", + " 1.69308689e-01 -2.05888038e-01]\n", + " [-1.84663809e-01 -4.65451966e-01 5.95026117e-02 -1.84663809e-01\n", + " -4.65451966e-01 -1.84663809e-01]\n", + " [ 4.66407064e-01 -3.87612079e-01 1.76410128e-02 4.66407064e-01\n", + " -3.87612079e-01 4.66407064e-01]\n", + " [ 4.52451104e-02 1.19865712e-01 -4.73313094e-03 4.52451104e-02\n", + " 1.19865712e-01 4.52451104e-02]\n", + " [ 6.90798449e-01 9.20396633e-01 1.49475827e-01 6.90798449e-01\n", + " 9.20396633e-01 6.90798449e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.36280366e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.98803165e-01 1.73327268e-01 -6.03008179e-03 -1.98803165e-01\n", + " 1.73327268e-01 -1.98803165e-01]\n", + " [-2.04303740e-01 -5.16111388e-01 4.68785776e-02 -2.04303740e-01\n", + " -5.16111388e-01 -2.04303740e-01]\n", + " [ 4.47070326e-01 -3.94304029e-01 1.32107437e-02 4.47070326e-01\n", + " -3.94304029e-01 4.47070326e-01]\n", + " [ 4.69732048e-02 1.22961727e-01 -3.35899442e-03 4.69732048e-02\n", + " 1.22961727e-01 4.69732048e-02]\n", + " [ 7.68998995e-01 1.10844286e+00 1.70889328e-01 7.68998995e-01\n", + " 1.10844286e+00 7.68998995e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.27091252e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.92789113e-01 1.75978657e-01 -4.54517629e-03 -1.92789113e-01\n", + " 1.75978657e-01 -1.92789113e-01]\n", + " [-2.20500138e-01 -5.55540705e-01 3.68776526e-02 -2.20500138e-01\n", + " -5.55540705e-01 -2.20500138e-01]\n", + " [ 4.30424855e-01 -3.97907706e-01 9.75257113e-03 4.30424855e-01\n", + " -3.97907706e-01 4.30424855e-01]\n", + " [ 4.82793652e-02 1.24952071e-01 -2.30991637e-03 4.82793652e-02\n", + " 1.24952071e-01 4.82793652e-02]\n", + " [ 8.40805131e-01 1.27504628e+00 1.89020151e-01 8.40805131e-01\n", + " 1.27504628e+00 8.40805131e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.18989205e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.87672774e-01 1.77588334e-01 -3.38318222e-03 -1.87672774e-01\n", + " 1.77588334e-01 -1.87672774e-01]\n", + " [-2.33807210e-01 -5.86081383e-01 2.89236334e-02 -2.33807210e-01\n", + " -5.86081383e-01 -2.33807210e-01]\n", + " [ 4.16201399e-01 -3.99295277e-01 7.06598588e-03 4.16201399e-01\n", + " -3.99295277e-01 4.16201399e-01]\n", + " [ 4.92546648e-02 1.26089711e-01 -1.50412006e-03 4.92546648e-02\n", + " 1.26089711e-01 4.92546648e-02]\n", + " [ 9.06806543e-01 1.42334018e+00 2.04522708e-01 9.06806543e-01\n", + " 1.42334018e+00 9.06806543e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.11829985e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.83320440e-01 1.78410042e-01 -2.47240692e-03 -1.83320440e-01\n", + " 1.78410042e-01 -1.83320440e-01]\n", + " [-2.44690164e-01 -6.09568485e-01 2.25774266e-02 -2.44690164e-01\n", + " -6.09568485e-01 -2.44690164e-01]\n", + " [ 4.04061655e-01 -3.99063012e-01 4.97908386e-03 4.04061655e-01\n", + " -3.99063012e-01 4.04061655e-01]\n", + " [ 4.99612484e-02 1.26581014e-01 -8.85891342e-04 4.99612484e-02\n", + " 1.26581014e-01 4.99612484e-02]\n", + " [ 9.67473970e-01 1.55589415e+00 2.17895305e-01 9.67473970e-01\n", + " 1.55589415e+00 9.67473970e-01]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.05500234e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.79620591e-01 1.78640114e-01 -1.75822439e-03 -1.79620591e-01\n", + " 1.78640114e-01 -1.79620591e-01]\n", + " [-2.53540123e-01 -6.27448857e-01 1.75019839e-02 -2.53540123e-01\n", + " -6.27448857e-01 -2.53540123e-01]\n", + " [ 3.93704970e-01 -3.97656641e-01 3.35895484e-03 3.93704970e-01\n", + " -3.97656641e-01 3.93704970e-01]\n", + " [ 5.04492282e-02 1.26586733e-01 -4.13401240e-04 5.04492282e-02\n", + " 1.26586733e-01 5.04492282e-02]\n", + " [ 1.02322336e+00 1.67481439e+00 2.29524046e-01 1.02322336e+00\n", + " 1.67481439e+00 1.02322336e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 4.99901907e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]\n", + "\n", + " [[-1.76478441e-01 1.78430281e-01 -1.19867662e-03 -1.76478441e-01\n", + " 1.78430281e-01 -1.76478441e-01]\n", + " [-2.60686971e-01 -6.40868686e-01 1.34365068e-02 -2.60686971e-01\n", + " -6.40868686e-01 -2.60686971e-01]\n", + " [ 3.84873835e-01 -3.95414931e-01 2.10369522e-03 3.84873835e-01\n", + " -3.95414931e-01 3.84873835e-01]\n", + " [ 5.07601805e-02 1.26231631e-01 -5.46465317e-05 5.07601805e-02\n", + " 1.26231631e-01 5.07601805e-02]\n", + " [ 1.07443160e+00 1.78183962e+00 2.39710937e-01 1.07443160e+00\n", + " 1.78183962e+00 1.07443160e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 4.94949118e-01\n", + " 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 1.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00]]]\n", + " ssigmay: [[[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0.]]]\n", + " sigmaz: None\n", + " sz: None\n", + " srz: None\n", + " ssigmaz: None\n", + " sllh: [nan nan nan nan nan nan nan nan]\n", + " s2llh: None\n", + " status: 0.0\n" + ] + } + ], + "source": [ + "model = modelModule.getModel()\n", + "model.setTimepoints(amici.DoubleVector(np.linspace(0, 10, 11))) \n", + "model.requireSensitivitiesForAllParameters()\n", + "model.setParameterScale(amici.AMICI_SCALING_NONE)\n", + "\n", + "solver = model.getSolver()\n", + "solver.setSensitivityMethod(amici.AMICI_SENSI_FSA)\n", + "solver.setSensitivityOrder(amici.AMICI_SENSI_ORDER_FIRST)\n", + "\n", + "rdata = amici.runAmiciSimulation(model, solver)\n", + "\n", + "for key, value in rdata.items():\n", + " if key.startswith('s'):\n", + " print('%12s: ' % key, value)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adjoint sensitivity analysis" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Noise-free: llh: -60.649943, sllh: [0. 0. 0. 0. 0. 0. 0. 0.]\n", + "Some noise: llh: -81.392520, sllh: [ -0.60535131 13.9453093 1.45890383 -3.46043573 -36.56923474\n", + " -2.81861858 -17.64991788 0. ]\n" + ] + } + ], + "source": [ + "model = modelModule.getModel()\n", + "model.setParameterScale(amici.AMICI_SCALING_NONE)\n", + "model.setTimepoints(amici.DoubleVector(np.linspace(0, 10, 11))) \n", + "solver = model.getSolver()\n", + "solver.setMaxSteps(10**4)\n", + "\n", + "# simulate time-course for artificial data\n", + "rdata = amici.runAmiciSimulation(model, solver)\n", + "edata = amici.ExpData(model.get())\n", + "edata.fixedParameters = model.getFixedParameters()\n", + "edata.my = amici.DoubleVector(rdata['y'].flatten())\n", + "edata.sigmay = amici.DoubleVector(np.ones(shape=rdata['y'].shape).flatten())\n", + "\n", + "# enable sensitivities\n", + "solver.setSensitivityMethod(amici.AMICI_SENSI_ASA)\n", + "solver.setSensitivityOrder(amici.AMICI_SENSI_ORDER_FIRST)\n", + "model.requireSensitivitiesForAllParameters()\n", + "\n", + "# compute adjoint sensitivities\n", + "rdata = amici.runAmiciSimulation(model, solver, edata)\n", + "print('Noise-free: llh: %f, sllh: %s' % (rdata['llh'], rdata['sllh']))\n", + "\n", + "# Add some noise\n", + "edata.my = amici.DoubleVector(np.multiply(rdata['y'], np.random.normal(0.0, 0.01, rdata['y'].shape)).flatten())\n", + "#edata.my = amici.DoubleVector((rdata['y'] * 1.1).flatten())\n", + "rdata = amici.runAmiciSimulation(model, solver, edata)\n", + "print('Some noise: llh: %f, sllh: %s' % (rdata['llh'], rdata['sllh']))\n", + "\n", + "p = np.array(model.getParameters())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Finite differences gradient check" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "sllh: |error|_2: 0.000004\n", + "\n", + "sllh: p[0]: |error|_2: 0.000001\n", + "sllh: p[1]: |error|_2: 0.000001\n", + "sllh: p[2]: |error|_2: 0.000001\n", + "sllh: p[3]: |error|_2: 0.000004\n", + "sllh: p[4]: |error|_2: 0.000000\n", + "sllh: p[5]: |error|_2: 0.000002\n", + "sllh: p[6]: |error|_2: 0.000001\n", + "sllh: p[7]: |error|_2: 0.000000\n", + "\n", + "sy: p[0]: |error|_2: 0.000001\n", + "sy: p[1]: |error|_2: 0.000001\n", + "sy: p[2]: |error|_2: 0.000001\n", + "sy: p[3]: |error|_2: 0.000001\n", + "sy: p[4]: |error|_2: 0.000001\n", + "sy: p[5]: |error|_2: 0.000000\n", + "sy: p[6]: |error|_2: 0.000000\n", + "sy: p[7]: |error|_2: 0.000000\n", + "\n", + "sx: p[0]: |error|_2: 0.000000\n", + "sx: p[1]: |error|_2: 0.000001\n", + "sx: p[2]: |error|_2: 0.000001\n", + "sx: p[3]: |error|_2: 0.000000\n", + "sx: p[4]: |error|_2: 0.000001\n", + "sx: p[5]: |error|_2: 0.000000\n", + "sx: p[6]: |error|_2: 0.000000\n", + "sx: p[7]: |error|_2: 0.000000\n" + ] + } + ], + "source": [ + "from scipy.optimize import check_grad\n", + "def func(x0, symbol='llh', x0full=None, plist=[], verbose=False):\n", + " p = x0\n", + " if len(plist):\n", + " p = x0full[:]\n", + " p[plist] = x0\n", + " verbose and print('f: p=%s' % p)\n", + " \n", + " solver.setSensitivityOrder(amici.AMICI_SENSI_ORDER_NONE)\n", + " model.setParameters(amici.DoubleVector(p))\n", + " rdata = amici.runAmiciSimulation(model, solver, edata)\n", + "\n", + " res = np.sum(rdata[symbol])\n", + " return res\n", + "\n", + "def grad(x0, symbol='llh', x0full=None, plist=[], verbose=False):\n", + " p = x0\n", + " if len(plist):\n", + " model.setParameterList(amici.IntVector(plist))\n", + " p = x0full[:]\n", + " p[plist] = x0\n", + " else:\n", + " model.requireSensitivitiesForAllParameters()\n", + " verbose and print('g: p=%s' % p)\n", + " \n", + " solver.setSensitivityMethod(amici.AMICI_SENSI_FSA)\n", + " solver.setSensitivityOrder(amici.AMICI_SENSI_ORDER_FIRST)\n", + " model.setParameters(amici.DoubleVector(p))\n", + " rdata = amici.runAmiciSimulation(model, solver, edata)\n", + " \n", + " res = rdata['s%s' % symbol]\n", + " if not isinstance(res, float):\n", + " if len(res.shape) == 3:\n", + " res = np.sum(res, axis=(0, 2))\n", + " return res\n", + "\n", + "err_norm = check_grad(func, grad, p, 'llh')\n", + "print('sllh: |error|_2: %f' % err_norm)\n", + "# assert err_norm < 1e-6\n", + "print()\n", + "\n", + "for ip in range(model.np()):\n", + " plist = [ip]\n", + " err_norm = check_grad(func, grad, p[plist], 'llh', p, [ip])\n", + " print('sllh: p[%d]: |error|_2: %f' % (ip, err_norm))\n", + "\n", + "print()\n", + "for ip in range(model.np()):\n", + " plist = [ip]\n", + " err_norm = check_grad(func, grad, p[plist], 'y', p, [ip])\n", + " print('sy: p[%d]: |error|_2: %f' % (ip, err_norm))\n", + "\n", + "print()\n", + "for ip in range(model.np()):\n", + " plist = [ip]\n", + " err_norm = check_grad(func, grad, p[plist], 'x', p, [ip])\n", + " print('sx: p[%d]: |error|_2: %f' % (ip, err_norm))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/python/examples/example_steadystate/steadystate.py b/python/examples/example_steadystate/steadystate.py deleted file mode 100644 index 7bf9b7687d..0000000000 --- a/python/examples/example_steadystate/steadystate.py +++ /dev/null @@ -1,72 +0,0 @@ -#!/usr/bin/env python3 -""" -Example using [model_steadystate_scaled.sbml] model to demonstrate and test SBML import and Amici Python interface. - -To use python-generated model in matlab, run: -``` -modelName = ''; -modelDir = ''; -amimodel.compileAndLinkModel(modelName, modelDir, [], [], [], []); -amimodel.generateMatlabWrapper(3, 6, 8, 1, 0, 0, [], [ modelDir '/simulate_test.m'], modelName, 'lin', 1, 1); -``` -""" - -import amici - -import os -import sys -import numpy as np -import matplotlib.pyplot as plt - - -os.chdir(os.path.dirname(__file__)) - -def createModule(): - """Create Python module from SBML model""" - sbmlImporter = amici.SbmlImporter('model_steadystate_scaled.sbml') - sbml = sbmlImporter.sbml - - observables = amici.assignmentRules2observables(sbml, filter=lambda variableId: - variableId.startswith('observable_') and not variableId.endswith('_sigma')) - - print(observables) - - sbmlImporter.sbml2amici('test', 'test', - observables=observables, - constantParameters=['k4'], - sigmas={'observable_x1withsigma': 'observable_x1withsigma_sigma'}) - -def plotStateTrajectories(rdata): - for ix in range(rdata['x'].shape[1]): - plt.plot(rdata['t'], rdata['x'][:, ix], label='$x_%d$' % ix) - plt.xlabel('$t$ (s)') - plt.ylabel('$x_i(t)$ (mmol/ml)') - plt.legend() - plt.title('State trajectories') - plt.show() - -def plotObservableTrajectories(rdata): - for iy in range(rdata['y'].shape[1]): - plt.plot(rdata['t'], rdata['y'][:, iy], label='$y_%d$' % iy) - plt.xlabel('$t$ (s)') - plt.ylabel('$y_i(t)$ (AU)') - plt.legend() - plt.title('Observables') - - plt.show() - -createModule() - -sys.path.insert(0, 'test') -import test as modelModule - -model = modelModule.getModel() -model.setTimepoints(amici.DoubleVector(np.linspace(0, 60, 60))) -solver = model.getSolver() -rdata = amici.runAmiciSimulation(model, solver) - -print(rdata) - -plotStateTrajectories(rdata) -plotObservableTrajectories(rdata) - diff --git a/python/sdist/setup.py b/python/sdist/setup.py index 8d44529872..4472c5a211 100644 --- a/python/sdist/setup.py +++ b/python/sdist/setup.py @@ -28,13 +28,27 @@ # Extra compiler flags cxx_flags = [] -# TODO: more flexible blas -amici_module_linker_flags = ['-lcblas'] +amici_module_linker_flags = [] define_macros = [] +# Find cblas +blaspkgcfg = {'include_dirs': [], + 'library_dirs': [], + 'libraries': [], + 'define_macros': [] + } +if 'BLAS_INCDIR' in os.environ: + blaspkgcfg['include_dirs'].extend(os.environ['BLAS_INCDIR'].split(' ')) + +if 'BLAS_LIB' in os.environ: + amici_module_linker_flags.extend(os.environ['BLAS_LIB'].split(' ')) +else: + amici_module_linker_flags.append('-lcblas') + # Find HDF5 include dir and libs -if pkgconfig.exists('hdf5'): - h5pkgcfg = pkgconfig.parse("hdf5") +h5pkgcfg = pkgconfig.parse('hdf5') +# NOTE: Cannot use pkgconfig.exists('hdf5f'), since this is true although no libraries or include dirs are available +if 'include_dirs' in h5pkgcfg and h5pkgcfg['include_dirs']: # Manually add linker flags. The libraries passed to Extension will # end up in front of the clibs in the linker line and not after, where # they are required. @@ -59,7 +73,7 @@ amici_module_linker_flags.append('--coverage') libamici = setup_clibs.getLibAmici( - h5pkgcfg=h5pkgcfg, extra_compiler_flags=cxx_flags) + h5pkgcfg=h5pkgcfg, blaspkgcfg=blaspkgcfg, extra_compiler_flags=cxx_flags) libsundials = setup_clibs.getLibSundials(extra_compiler_flags=cxx_flags) libsuitesparse = setup_clibs.getLibSuiteSparse(extra_compiler_flags=cxx_flags) @@ -91,9 +105,9 @@ class my_build_ext(build_ext): def run(self): """Copy the generated clibs to the extensions folder to be included in the wheel - + Returns: - + """ if not self.dry_run: # --dry-run @@ -132,7 +146,7 @@ def run(self): """Setuptools entry-point Returns: - + """ self.runSwig() self.saveGitVersion() @@ -142,12 +156,14 @@ def runSwig(self): """Run swig Returns: - + """ + if not self.dry_run: # --dry-run # We create two SWIG interfaces, one with HDF5 support, one without swig_outdir = '%s/amici' % os.path.abspath(os.getcwd()) - sp = subprocess.run(['swig3.0', + swig_cmd = self.findSwig() + sp = subprocess.run([swig_cmd, '-c++', '-python', '-Iamici/swig', '-Iamici/include', @@ -158,7 +174,7 @@ def runSwig(self): assert(sp.returncode == 0) shutil.move(os.path.join(swig_outdir, 'amici.py'), os.path.join(swig_outdir, 'amici_without_hdf5.py')) - sp = subprocess.run(['swig3.0', + sp = subprocess.run([swig_cmd, '-c++', '-python', '-Iamici/swig', '-Iamici/include', @@ -167,6 +183,17 @@ def runSwig(self): 'amici/swig/amici.i']) assert(sp.returncode == 0) + def findSwig(self): + """Get name of SWIG executable + + We need version 3.0. + Probably we should try some default paths and names, but this should do the trick for now. + Debian/Ubuntu systems have swig3.0 ('swig' is older versions), OSX has swig 3.0 as 'swig'.""" + if sys.platform != 'linux': + return 'swig' + return 'swig3.0' + + def saveGitVersion(self): """Create file with extended version string @@ -174,7 +201,7 @@ def saveGitVersion(self): a valid git repository. Returns: - + """ f = open("amici/version.txt", "w") sp = subprocess.run(['git', 'describe', @@ -191,7 +218,7 @@ def saveGitVersion(self): def getPackageVersion(): - return '0.6a13' + return '0.6a15' # Remove the "-Wstrict-prototypes" compiler option, which isn't valid for diff --git a/python/sdist/setup_clibs.py b/python/sdist/setup_clibs.py index a7c3d70b79..b0ff715f55 100644 --- a/python/sdist/setup_clibs.py +++ b/python/sdist/setup_clibs.py @@ -103,15 +103,23 @@ def getSuiteSparseSources(): def getAmiciBaseSources(withHDF5=True): """Get list of source files for the amici base library + Expects that we are inside $AMICI_ROOT/python/sdist + Arguments: withHDF5: compile with HDF5 support """ + amiciBaseSources = glob.glob('amici{s}src{s}*.cpp'.format(s=os.sep)) amiciBaseSources = [src for src in amiciBaseSources if not re.search( r'(matlab)|(\.template\.)', src)] if not withHDF5: - amiciBaseSources.remove('amici{s}src{s}hdf5.cpp'.format(s=os.sep)) + try: + # sometimes this fails for unknwon reasons... + amiciBaseSources.remove('amici{s}src{s}hdf5.cpp'.format(s=os.sep)) + except ValueError: + print('Warning: could not find %s in %s' % ('amici{s}src{s}hdf5.cpp'.format(s=os.sep), + amiciBaseSources)) return amiciBaseSources @@ -160,11 +168,13 @@ def getLibSuiteSparse(extra_compiler_flags=[]): return libsuitesparse -def getLibAmici(extra_compiler_flags=[], h5pkgcfg=None): +def getLibAmici(extra_compiler_flags=[], h5pkgcfg=None, blaspkgcfg=None): """Get AMICI core library build info for setuptools Arguments: extra_compiler_flags: Extra compiler flags + h5pkgcfg: hdf5 package info + blaspkgcfg: blas package info """ libamici = ('amici', { @@ -187,4 +197,7 @@ def getLibAmici(extra_compiler_flags=[], h5pkgcfg=None): if h5pkgcfg and 'include_dirs' in h5pkgcfg: libamici[1]['include_dirs'].extend(h5pkgcfg['include_dirs']) + if blaspkgcfg and 'include_dirs' in blaspkgcfg: + libamici[1]['include_dirs'].extend(blaspkgcfg['include_dirs']) + return libamici diff --git a/src/model_header.ODE_template.h b/src/model_header.ODE_template.h index d317263106..ca1f9a951e 100644 --- a/src/model_header.ODE_template.h +++ b/src/model_header.ODE_template.h @@ -29,6 +29,7 @@ extern void dwdp_TPL_MODELNAME(realtype *dwdp, const realtype t, const realtype extern void dwdx_TPL_MODELNAME(realtype *dwdx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w); extern void dxdotdp_TPL_MODELNAME(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp); extern void dydx_TPL_MODELNAME(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h); +extern void dydp_TPL_MODELNAME(double *dydp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip); extern void qBdot_TPL_MODELNAME(realtype *qBdot, const int ip, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *xB, const realtype *w, const realtype *dwdp); extern void sigma_y_TPL_MODELNAME(double *sigmay, const realtype t, const realtype *p, const realtype *k); extern void sxdot_TPL_MODELNAME(realtype *sxdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *sx, const realtype *w, const realtype *dwdx, const realtype *J, const realtype *dxdotdp); @@ -440,28 +441,29 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { dxdotdp_TPL_MODELNAME(dxdotdp, t, x, p, k, h, ip, w, dwdp); } - /** model specific implementation of fdydp - * @param dydp partial derivative of observables y w.r.t. model parameters p + /** model specific implementation of fdydx + * @param dydx partial derivative of observables y w.r.t. model states x * @param t current time * @param x current state * @param p parameter vector * @param k constant vector * @param h heavyside vector - * @param ip parameter index w.r.t. which the derivative is requested **/ - virtual void fdydp(double *dydp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override { + virtual void fdydx(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { + dydx_TPL_MODELNAME(dydx, t, x, p, k, h); } - /** model specific implementation of fdydx - * @param dydx partial derivative of observables y w.r.t. model states x + /** model specific implementation of fdydp + * @param dydp partial derivative of observables y w.r.t. model parameters p * @param t current time * @param x current state * @param p parameter vector * @param k constant vector * @param h heavyside vector + * @param ip parameter index w.r.t. which the derivative is requested **/ - virtual void fdydx(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { - dydx_TPL_MODELNAME(dydx, t, x, p, k, h); + virtual void fdydp(double *dydp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override { + dydp_TPL_MODELNAME(dydp, t, x, p, k, h, ip); } /** model specific implementation of fdzdp diff --git a/swig/CMakeLists.txt b/swig/CMakeLists.txt index 6a672d2789..bb5fc2b0bc 100644 --- a/swig/CMakeLists.txt +++ b/swig/CMakeLists.txt @@ -6,24 +6,24 @@ find_package(SWIG REQUIRED 3) include(${SWIG_USE_FILE}) set(AMICI_INTERFACE_LIST - ${CMAKE_CURRENT_LIST_DIR}/amici.i - ${CMAKE_CURRENT_LIST_DIR}/edata.i - ${CMAKE_CURRENT_LIST_DIR}/rdata.i - ${CMAKE_CURRENT_LIST_DIR}/model.i - ${CMAKE_CURRENT_LIST_DIR}/model_ode.i - ${CMAKE_CURRENT_LIST_DIR}/model_dae.i - ${CMAKE_CURRENT_LIST_DIR}/solver.i - ${CMAKE_CURRENT_LIST_DIR}/solver_cvodes.i - ${CMAKE_CURRENT_LIST_DIR}/solver_idas.i - ${CMAKE_CURRENT_LIST_DIR}/std_unique_ptr.i + ${CMAKE_CURRENT_SOURCE_DIR}/amici.i + ${CMAKE_CURRENT_SOURCE_DIR}/edata.i + ${CMAKE_CURRENT_SOURCE_DIR}/rdata.i + ${CMAKE_CURRENT_SOURCE_DIR}/model.i + ${CMAKE_CURRENT_SOURCE_DIR}/model_ode.i + ${CMAKE_CURRENT_SOURCE_DIR}/model_dae.i + ${CMAKE_CURRENT_SOURCE_DIR}/solver.i + ${CMAKE_CURRENT_SOURCE_DIR}/solver_cvodes.i + ${CMAKE_CURRENT_SOURCE_DIR}/solver_idas.i + ${CMAKE_CURRENT_SOURCE_DIR}/std_unique_ptr.i + ${CMAKE_CURRENT_SOURCE_DIR}/hdf5.i ) # Add target to show files in IDE add_custom_target(swigInterface SOURCES ${AMICI_INTERFACE_LIST}) ## Add subdirectories for each language if desired -#option(BUILD_PYTHON "Build Python SWIG module" ON) -#if(BUILD_PYTHON) -# add_subdirectory(python) -#endif() +if(ENABLE_PYTHON) + add_subdirectory(python) +endif() diff --git a/swig/amici.i b/swig/amici.i index 89a3e7bd35..589f0d8fb9 100644 --- a/swig/amici.i +++ b/swig/amici.i @@ -2,8 +2,7 @@ typedef double realtype; -%include -%include +%include %include std_unique_ptr.i wrap_unique_ptr(SolverPtr, amici::Solver) @@ -11,6 +10,10 @@ wrap_unique_ptr(ReturnDataPtr, amici::ReturnData) wrap_unique_ptr(ModelPtr, amici::Model) wrap_unique_ptr(ExpDataPtr, amici::ExpData) +// Include before any other header which uses enums defined there +%include "amici/defines.h" + + %include edata.i %include rdata.i @@ -34,8 +37,7 @@ using namespace amici; // Process symbols in header %include "amici/amici.h" -namespace std -{ - %template(DoubleVector) vector; - %template(IntVector) vector; -} +// Expose vectors +%template(DoubleVector) std::vector; +%template(IntVector) std::vector; +%template(ParameterScalingVector) std::vector; diff --git a/swig/python/CMakeLists.txt b/swig/python/CMakeLists.txt index d68537e4a8..8034527237 100644 --- a/swig/python/CMakeLists.txt +++ b/swig/python/CMakeLists.txt @@ -5,38 +5,44 @@ # Find Python interpreter, libs and headers find_package(PythonInterp 3.6 REQUIRED)# TODO to setup.py find_package(PythonLibs REQUIRED) -include_directories(${PYTHON_INCLUDE_PATH}) +include_directories( + ${PYTHON_INCLUDE_PATH} + ${CMAKE_SOURCE_DIR}/include + ${SUITESPARSE_INCLUDE_DIRS} + ${SUNDIALS_INCLUDE_DIRS} + ${HDF5_INCLUDE_DIRS} + ) set(CMAKE_SWIG_OUTDIR ${CMAKE_CURRENT_BINARY_DIR}/amici) #set(CMAKE_SWIG_FLAGS "") set_source_files_properties( - ${AMICI_SRC_LIST} ${AMICI_DIR}/swig/amici.i + ${AMICI_SRC_LIST} ${CMAKE_SOURCE_DIR}/swig/amici.i PROPERTIES CPLUSPLUS ON ) SWIG_ADD_LIBRARY( amici MODULE LANGUAGE python - SOURCES ${AMICI_DIR}/swig/amici.i + SOURCES ${CMAKE_SOURCE_DIR}/swig/amici.i ) -SWIG_LINK_LIBRARIES(amici - ${PYTHON_LIBRARIES} - ${AMICI_DIR}/build/libamici${CMAKE_STATIC_LIBRARY_SUFFIX} - ${SUNDIALS_LIB_DIR}/libsundials_nvecserial${CMAKE_STATIC_LIBRARY_SUFFIX} - ${SUNDIALS_LIB_DIR}/libsundials_cvodes${CMAKE_STATIC_LIBRARY_SUFFIX} - ${SUNDIALS_LIB_DIR}/libsundials_idas${CMAKE_STATIC_LIBRARY_SUFFIX} - ${SUITESPARSE_DIR}/KLU/Lib/libklu${CMAKE_STATIC_LIBRARY_SUFFIX} - ${SUITESPARSE_DIR}/COLAMD/Lib/libcolamd${CMAKE_STATIC_LIBRARY_SUFFIX} - ${SUITESPARSE_DIR}/BTF/Lib/libbtf${CMAKE_STATIC_LIBRARY_SUFFIX} - ${SUITESPARSE_DIR}/AMD/Lib/libamd${CMAKE_STATIC_LIBRARY_SUFFIX} - ${SUITESPARSE_DIR}/SuiteSparse_config/libsuitesparseconfig${CMAKE_STATIC_LIBRARY_SUFFIX} - ${HDF5_HL_LIBRARIES} - ${HDF5_C_LIBRARIES} - ${HDF5_CXX_LIBRARIES} - ${BLAS_LIBRARIES} - -ldl -lz -lm) +#SWIG_LINK_LIBRARIES(amici +# ${PYTHON_LIBRARIES} +# ${AMICI_DIR}/build/libamici${CMAKE_STATIC_LIBRARY_SUFFIX} +# ${SUNDIALS_LIB_DIR}/libsundials_nvecserial${CMAKE_STATIC_LIBRARY_SUFFIX} +# ${SUNDIALS_LIB_DIR}/libsundials_cvodes${CMAKE_STATIC_LIBRARY_SUFFIX} +# ${SUNDIALS_LIB_DIR}/libsundials_idas${CMAKE_STATIC_LIBRARY_SUFFIX} +# ${SUITESPARSE_DIR}/KLU/Lib/libklu${CMAKE_STATIC_LIBRARY_SUFFIX} +# ${SUITESPARSE_DIR}/COLAMD/Lib/libcolamd${CMAKE_STATIC_LIBRARY_SUFFIX} +# ${SUITESPARSE_DIR}/BTF/Lib/libbtf${CMAKE_STATIC_LIBRARY_SUFFIX} +# ${SUITESPARSE_DIR}/AMD/Lib/libamd${CMAKE_STATIC_LIBRARY_SUFFIX} +# ${SUITESPARSE_DIR}/SuiteSparse_config/libsuitesparseconfig${CMAKE_STATIC_LIBRARY_SUFFIX} +# ${HDF5_HL_LIBRARIES} +# ${HDF5_C_LIBRARIES} +# ${HDF5_CXX_LIBRARIES} +# ${BLAS_LIBRARIES} +# -ldl -lz -lm) # put library into python package # (SWIG_ADD_LIBRARY directory options don't work in current cmake version) diff --git a/tests/testModels.py b/tests/testModels.py index bda7ea7bbf..b71dfa863c 100755 --- a/tests/testModels.py +++ b/tests/testModels.py @@ -62,8 +62,90 @@ def runTest(self): else: verifySimulationResults(rdata, expectedResults[subTest][case]['results']) + if edata and self.solver.getSensitivityMethod() and self.solver.getSensitivityOrder(): + if not modelName.startswith('model_neuron'): + checkGradient(self.model, self.solver, edata) + - +def checkGradient(model, solver, edata): + """Finite differences check for likelihood gradient + + Arguments: + model: + solver: + edata: + """ + from scipy.optimize import check_grad + + def func(x0, symbol='llh', x0full=None, plist=[], verbose=False): + """Function of which the gradient is to be checked""" + p = x0 + if len(plist): + p = x0full[:] + p[plist] = x0 + verbose and print('f: p=%s' % p) + + old_sensitivity_order = solver.getSensitivityOrder() + old_parameters = model.getParameters() + + solver.setSensitivityOrder(amici.AMICI_SENSI_ORDER_NONE) + model.setParameters(amici.DoubleVector(p)) + rdata = amici.runAmiciSimulation(model, solver, edata) + + solver.setSensitivityOrder(old_sensitivity_order) + model.setParameters(old_parameters) + + res = np.sum(rdata[symbol]) + return res + + def grad(x0, symbol='llh', x0full=None, plist=[], verbose=False): + """Gradient which is to be checked""" + old_parameters = model.getParameters() + old_plist = model.getParameterList() + + p = x0 + if len(plist): + model.setParameterList(amici.IntVector(plist)) + p = x0full[:] + p[plist] = x0 + else: + model.requireSensitivitiesForAllParameters() + verbose and print('g: p=%s' % p) + + model.setParameters(amici.DoubleVector(p)) + rdata = amici.runAmiciSimulation(model, solver, edata) + + model.setParameters(old_parameters) + model.setParameterList(old_plist) + + res = rdata['s%s' % symbol] + if not isinstance(res, float): + if len(res.shape) == 3: + res = np.sum(res, axis=(0, 2)) + return res + + p = np.array(model.getParameters()) + + 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)) + + ''' + 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)) + + ''' + def verifySimulationResults(rdata, expectedResults, atol=1e-8, rtol=1e-4): ''' compares all fields of the simulation results in rdata against the expectedResults using the provided @@ -134,4 +216,4 @@ def checkResults(rdata, field, expected, atol, rtol): suite = unittest.TestSuite() suite.addTest(TestAmiciPregeneratedModel()) unittest.main() - \ No newline at end of file +