Fix events for DEA models (#687)

* remove recalculation of ICs after event * fix reinitialization after event * add calvetti to test examples * add note about impulse free events * fix tolerances for robertson example * remove testOptions and make test fail on missing results data * readd restOptions.h5 * restrict tolerances for robertsonSPBCG * next try robertson tolerances * refine tolerances for J/xdot check * next try ... * refine jakstat test tolerances * update test results
AMICI-dev · Apr 19, 2019 · 42181fa · 42181fa
1 parent feaebd2
commit 42181fa
Show file tree

Hide file tree

Showing 44 changed files with 1,451 additions and 21 deletions.
diff --git a/.gitignore b/.gitignore
@@ -34,6 +34,8 @@ models/model_events/build/*
 models/model_nested_events/build/*
 !models/model_robertson
 models/model_robertson/build/*
+!models/model_calvetti
+models/model_calvetti/build/*
 
 simulate_model_*_hdf.m
 simulate_model_*.m

diff --git a/documentation/FAQ.md b/documentation/FAQ.md
@@ -18,6 +18,12 @@ __A__: This may be due to an old compiler version. See [issue #161](https://gith
 
 ---
 
+__Q__: How are events interpreted in a DAE context?
+
+__A__: Currently we only support impulse free events. Also sensitivities have never been tested. Proceed with care and create an [issue](https://github.com/ICB-DCM/AMICI/issues) if any problems arise!
+
+---
+
 __Q__: The simulation/sensitivities I get are incorrect.
 
 __A__: There are some known issues, especially with adjoint sensitivities, events and DAEs. If your particular problem is not featured in the [issues](https://github.com/ICB-DCM/AMICI/issues) list, please add it!

diff --git a/matlab/examples/amiExamples.m b/matlab/examples/amiExamples.m
@@ -13,4 +13,5 @@
 example_jakstat_adjoint_hvp
 example_robertson
 example_neuron
+example_dae_events
 path(oldpath);
diff --git a/matlab/examples/example_calvetti/example_calvetti.m b/matlab/examples/example_calvetti/example_calvetti.m
@@ -0,0 +1,101 @@
+function example_calvetti()
+%%
+% COMPILATION
+
+[exdir,~,~]=fileparts(which('example_calvetti.m'));
+% compile the model
+amiwrap('model_calvetti','model_calvetti_syms',exdir)
+
+% time vector
+t = linspace(0,20,201);
+p = [];
+k = [0.29 0.74 0.44 0.08 0.27 0.18];
+
+options = amioption('sensi',0,...
+    'maxsteps',1e4);
+D = amidata(length(t),6,0,2,4);
+% load mex into memory
+[~] = which('simulate_model_calvetti'); % fix for inaccessability problems
+sol = simulate_model_calvetti(t,p,k,D,options);
+
+tic
+sol = simulate_model_calvetti(t,p,k,D,options);
+disp(['Time elapsed with cvodes: ' num2str(toc) ])
+
+%%
+% ODE15S
+
+y0 = [k(1); k(3); k(5); 1; 1; 1;];
+M = [1 0 0 0 0 0 
+     0 1 0 0 0 0
+     0 0 1 0 0 0
+     0 0 0 0 0 0
+     0 0 0 0 0 0
+     0 0 0 0 0 0];
+
+function [xdot] = dae_system(t,x,p,k,it)
+    if it<3
+        h0 = 0;
+    else
+        h0 = 1;
+    end
+    if it<4
+        h1 = 0;
+    else
+        h1 = 1;
+    end
+    xdot = [
+    h0/31 - h1/31 - (100*x(1))/(899*k(1)) + (2*x(1)*(k(2)/2 - 1))/(31*k(1)*(x(4)*((k(2)*k(1)^2)/x(1)^2 + (k(4)*k(3)^2)/x(2)^2) + x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2)) + 129/899;
+    (2*x(2)*(k(2) + k(4)/2 - 1))/(163*k(3)*(x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2)) - (100*x(2))/(8313*k(3)) + 151/8313;
+    (x(3)^3*(k(2) + k(4) +k(6)/2 - 1))/(61*k(6)*k(5)^3*x(6)) - x(3)/(121999878*k(5)) + 500000/60999939;
+    h1/31 - h0/31 - x(4) + (100*x(1))/(899*k(1)) + (2*x(1)^2)/(k(2)*k(1)^2) - (x(1)^2*(x(4)*((k(2)*k(1)^2)/x(1)^2 + (k(4)*k(3)^2)/x(2)^2) + x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2))/(k(2)*k(1)^2) - (2*x(1)*(k(2)/2 - 1))/(31*k(1)*(x(4)*((k(2)*k(1)^2)/x(1)^2 + (k(4)*k(3)^2)/x(2)^2) + x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2)) - 129/899;
+    x(4) - x(5) + (100*x(2))/(8313*k(3)) - (2*x(2)*(k(2) + k(4)/2 - 1))/(163*k(3)*(x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2)) - 151/8313;
+    x(5) - x(6) + x(3)/(121999878*k(5)) - (x(3)^3*(k(2) + k(4) +k(6)/2 - 1))/(61*k(6)*k(5)^3*x(6)) - 500000/60999939;
+];
+end
+
+options_ode15s = odeset('Mass',M, ...
+                        'RelTol',options.rtol, ...
+                        'AbsTol',options.atol);
+
+tic
+X_ode15s = [];
+tt = t;
+tstop = [0,10,12,20];
+for it = 2:4
+    [~, y] = ode15s(@(t,x) dae_system(t,x,p,k,it),t(t>=tstop(it-1) & t<=tstop(it)),y0,options_ode15s);
+    X_ode15s = [X_ode15s;y(1:end-1,:)];
+    y0 = y(end,:);
+end
+X_ode15s = [X_ode15s;y(end,:)];
+disp(['Time elapsed with ode15s: ' num2str(toc) ])
+
+%%
+% PLOTTING
+if(usejava('jvm'))
+    figure
+    c_x = get(gca,'ColorOrder');
+    subplot(2,1,1)
+    for ix = 1:size(sol.x,2)
+        plot(t,sol.x(:,ix),'.-','Color',c_x(ix,:))
+        hold on
+        plot(t,X_ode15s(:,ix),'d','Color',c_x(ix,:))
+    end
+    legend('x1','x1_{ode15s}','x2','x2_{ode15s}', ...
+           'x3','x3_{ode15s}','x4','x4_{ode15s}', ...
+           'x5','x5_{ode15s}','x6','x6_{ode15s}', ...
+           'Location','NorthEastOutside')
+    legend boxoff
+    xlabel('time t')
+    ylabel('x')
+    box on
+    subplot(2,1,2)
+    plot(t,abs(sol.x-X_ode15s),'--')
+    set(gca,'YScale','log')
+    legend('error x1','error x2','error x3','error x4','error x5','error x6','Location','NorthEastOutside')
+    legend boxoff
+    ylabel('x')
+
+    set(gcf,'Position',[100 300 1200 500])
+end
+end
diff --git a/matlab/examples/example_calvetti/model_calvetti_syms.m b/matlab/examples/example_calvetti/model_calvetti_syms.m
@@ -0,0 +1,75 @@
+function [model] = model_calvetti_syms()
+
+% set the parametrisation of the problem options are 'log', 'log10' and
+
+model.param = 'lin';
+
+%%
+% STATES
+% create state syms 
+syms V1 V2 V3 f1 f2 f3
+% create state vector
+model.sym.x = [V1, V2, V3, f1, f2, f3];
+%%
+% PARAMETERS ( for these sensitivities will be computed )
+
+% create parameter syms                                                                
+% create parameter vector 
+model.sym.p = [ ];
+%%  
+% CONSTANTS ( for these no sensitivities will be computed )
+% this part is optional and can be ommited
+
+% create parameter syms
+syms V1ss R1ss V2ss R2ss V3ss R3ss
+% create parameter vector 
+model.sym.k = [V1ss, R1ss, V2ss, R2ss, V3ss, R3ss];
+%%
+% SYSTEM EQUATIONS
+% create symbolic variable for time
+syms t f0 
+model.sym.xdot = sym(zeros(size(model.sym.x)));
+p1=1;
+p2=1-R1ss;
+p3=1-(R1ss+R2ss);
+L1=(R1ss*(V1ss^2))^(1/3);
+C1ss=V1ss/(p1-0.5*R1ss);
+
+C2ss=V2ss/(p2-0.5*R2ss);
+L2=(R2ss*(V2ss^2))^(1/3);
+C3ss=V3ss/(p3-0.5*R3ss);
+L3=(R3ss*(V3ss^2))^(1/3);
+R2=(L2^3)/(V2^2);
+R1=(L1^3)/(V1^2);
+R3=(L3^3)/(V3^2);
+s=am_stepfun(t,10,1,12,0);
+model.sym.xdot(1)=(1/31)*(((1.29 -(V1/V1ss))/(1.29-1))+s-2*(V1/(C1ss*((R1+R2)*f1+(R2+R3)*f2+R3*f3))));
+model.sym.xdot(2)=(1/163)*(((1.51 -(V2/V2ss))/(1.51-1))- 2*(V2/(C2ss*((R2+R3)*f2+R3*f3))));
+model.sym.xdot(3)=(1/122)*(((1000000 -(V3/V3ss))/(1000000-1))- 2*(V3/(C3ss*(R3*f3))));
+f0=(2/R1)-((R1+R2)*f1 + (R2+R3)*f2 +R3*f3)/R1;
+model.sym.xdot(4)=f0-model.sym.xdot(1)-f1;
+model.sym.xdot(5)=f1-model.sym.xdot(2)-f2;
+model.sym.xdot(6)=f2-model.sym.xdot(3)-f3;
+
+
+model.sym.M=[1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0;];
+%%
+% INITIAL CONDITIONS
+model.sym.x0(1) = V1ss;
+model.sym.x0(2)=V2ss;
+model.sym.x0(3)=V3ss;
+model.sym.x0(4)=1;
+model.sym.x0(5)=1;
+model.sym.x0(6)=1;
+model.sym.dx0 = sym(zeros(size(model.sym.x)));
+
+% OBSERVALES
+model.sym.y = sym(zeros(1,6));
+model.sym.y(1) =V1;
+model.sym.y(2)=V2;
+model.sym.y(3)=V3;
+model.sym.y(4)=f0;
+model.sym.y(5)=f1;
+model.sym.y(6)=f2;
+end
+
diff --git a/models/model_calvetti/CMakeLists.txt b/models/model_calvetti/CMakeLists.txt
@@ -0,0 +1,94 @@
+cmake_minimum_required(VERSION 2.8)
+
+if(POLICY CMP0060)
+  cmake_policy(SET CMP0060 NEW)
+endif(POLICY CMP0060)
+if(POLICY CMP0065)
+  cmake_policy(SET CMP0065 NEW)
+endif(POLICY CMP0065)
+
+project(model_calvetti)
+
+set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+
+include(CheckCXXCompilerFlag)
+set(MY_CXX_FLAGS -Wall -Wno-unused-function -Wno-unused-variable -Wno-unused-but-set-variable)
+foreach(FLAG ${MY_CXX_FLAGS})
+    unset(CUR_FLAG_SUPPORTED CACHE)
+    CHECK_CXX_COMPILER_FLAG(${FLAG} CUR_FLAG_SUPPORTED)
+    if(${CUR_FLAG_SUPPORTED})
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAG}")
+    endif()
+endforeach(FLAG)
+
+find_package(Amici HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build)
+
+set(MODEL_DIR ${CMAKE_CURRENT_LIST_DIR})
+
+set(SRC_LIST_LIB ${MODEL_DIR}/model_calvetti_J.cpp
+${MODEL_DIR}/model_calvetti_JB.cpp
+${MODEL_DIR}/model_calvetti_JDiag.cpp
+${MODEL_DIR}/model_calvetti_JSparse.cpp
+${MODEL_DIR}/model_calvetti_JSparseB.cpp
+${MODEL_DIR}/model_calvetti_Jy.cpp
+${MODEL_DIR}/model_calvetti_M.cpp
+${MODEL_DIR}/model_calvetti_dJydsigma.cpp
+${MODEL_DIR}/model_calvetti_dJydy.cpp
+${MODEL_DIR}/model_calvetti_dwdx.cpp
+${MODEL_DIR}/model_calvetti_dydx.cpp
+${MODEL_DIR}/model_calvetti_root.cpp
+${MODEL_DIR}/model_calvetti_sigmay.cpp
+${MODEL_DIR}/model_calvetti_w.cpp
+${MODEL_DIR}/model_calvetti_x0.cpp
+${MODEL_DIR}/model_calvetti_xdot.cpp
+${MODEL_DIR}/model_calvetti_y.cpp
+
+${MODEL_DIR}/wrapfunctions.cpp
+)
+
+add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
+add_library(model ALIAS ${PROJECT_NAME})
+
+target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
+
+target_link_libraries(${PROJECT_NAME}
+    PUBLIC Upstream::amici
+)
+
+set(SRC_LIST_EXE main.cpp)
+
+add_executable(simulate_${PROJECT_NAME} ${SRC_LIST_EXE})
+
+target_link_libraries(simulate_${PROJECT_NAME} ${PROJECT_NAME})
+
+if($ENV{ENABLE_GCOV_COVERAGE})
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 --coverage")
+endif()
+
+## SWIG
+option(ENABLE_SWIG "Build swig/python library?" ON)
+if(ENABLE_SWIG)
+    if(NOT(${CMAKE_VERSION} VERSION_LESS 3.8))
+        add_subdirectory(swig)
+    else()
+        message(WARNING "Unable to build SWIG interface, upgrade CMake to >=3.8.")
+    endif()
+endif()
+
+
+# <Export cmake configuration>
+include(GNUInstallDirs)
+install(TARGETS ${PROJECT_NAME} EXPORT ${PROJECT_NAME}Targets
+    ARCHIVE  DESTINATION ${CMAKE_INSTALL_LIBDIR}
+    LIBRARY  DESTINATION ${CMAKE_INSTALL_LIBDIR}
+    RUNTIME  DESTINATION ${CMAKE_INSTALL_BINDIR}
+    INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+    PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+)
+export(EXPORT ${PROJECT_NAME}Targets FILE ${PROJECT_NAME}Config.cmake
+    NAMESPACE Upstream::
+    )
+# </Export cmake configuration>
+