Merge pull request #89 from ICB-DCM/dw_misc

* Jak/Stat example as additional test case
* More informative test output

@amifun/amifun.m

@@ -3,7 +3,7 @@
 % @brief definition of amifun class
 %
 classdef amifun
-    % AMIFUN defines  functions which later on will be transformed into
+    % AMIFUN defines functions which later on will be transformed into
     % appropriate C code
 
     properties ( GetAccess = 'public', SetAccess = 'public' )
@@ -39,7 +39,7 @@
             % Parameters:
             %  funstr: name of the requested function
             %  model: amimodel object which carries all symbolic
-            %  definitions to construct the funtion
+            %  definitions to construct the function
             %  
             %
             % Return values:

@amifun/gccode.m

@@ -3,7 +3,7 @@
     % respective expression into a specified file
     %
     % Parameters:
-    %  model: model defintion object @type amimodel
+    %  model: model definition object @type amimodel
     %  fid: file id in which the expression should be written @type fileid
     %
     % Return values:

amiwrap.m

@@ -3,8 +3,8 @@ function amiwrap( varargin )
     %
     % Parameters:
     %  varargin:
-    %  modelname: specifies the name of the model which will be later used for the naming of the simualation file @type string
-    %  symfun: specifies a function which executes model defition see @ref definition for details @type string.
+    %  modelname: specifies the name of the model which will be later used for the naming of the simulation file @type string
+    %  symfun: specifies a function which executes model definition see @ref definition for details @type string.
     %  tdir: target directory where the simulation file should be placed @type string @default $AMICIDIR/models/modelname
     %  o2flag: boolean whether second order sensitivities should be enabled @type boolean @default false
     %

include/amici.h

@@ -28,6 +28,8 @@
 #define AMI_NORMAL                1
 #define AMI_ONE_STEP              2
 
+void runAmiciSimulation(UserData *udata, const ExpData *edata, ReturnData *rdata, int *pstatus);
+
 void *setupAMI(int *status, UserData *udata, TempData *tdata);
 void setupAMIB(int *status, void *ami_mem, UserData *udata, TempData *tdata);
 

models/model_jakstat_adjoint/CMakeLists.txt

@@ -0,0 +1,110 @@
+project(model_jakstat_adjoint)
+cmake_minimum_required(VERSION 2.8)
+
+set(CMAKE_BUILD_TYPE Debug)
+
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wall -Wno-unused-function")
+add_definitions(-DAMICI_WITHOUT_MATLAB)
+
+set(AMICI_DIR "${CMAKE_CURRENT_SOURCE_DIR}/../../")
+set(MODEL_DIR "${AMICI_DIR}/models/model_jakstat_adjoint")
+set(SUITESPARSE_DIR "${AMICI_DIR}/SuiteSparse/")
+set(SUITESPARSE_LIB_DIR "${AMICI_DIR}/SuiteSparse/lib")
+set(SUNDIALS_LIB_DIR "${AMICI_DIR}/sundials/build/lib")
+
+find_package(HDF5 COMPONENTS C HL REQUIRED)
+include_directories("${AMICI_DIR}")
+include_directories("${CMAKE_CURRENT_SOURCE_DIR}")
+include_directories("${MODEL_DIR}")
+include_directories("${HDF5_INCLUDE_DIRS}")
+include_directories("${AMICI_DIR}/sundials/include")
+include_directories("${SUITESPARSE_DIR}/include")
+
+
+set(SRC_LIST_LIB
+${MODEL_DIR}/wrapfunctions.cpp
+${AMICI_DIR}/src/symbolic_functions.cpp
+${AMICI_DIR}/src/amici_interface_cpp.cpp
+${AMICI_DIR}/src/amici.cpp
+${AMICI_DIR}/src/udata.cpp
+${AMICI_DIR}/src/rdata.cpp
+${AMICI_DIR}/src/edata.cpp
+${AMICI_DIR}/src/ami_hdf5.cpp
+${AMICI_DIR}/src/spline.cpp
+${MODEL_DIR}/model_jakstat_adjoint_J.cpp
+${MODEL_DIR}/model_jakstat_adjoint_JB.cpp
+${MODEL_DIR}/model_jakstat_adjoint_JBand.cpp
+${MODEL_DIR}/model_jakstat_adjoint_JBandB.cpp
+${MODEL_DIR}/model_jakstat_adjoint_JSparse.cpp
+${MODEL_DIR}/model_jakstat_adjoint_JSparseB.cpp
+${MODEL_DIR}/model_jakstat_adjoint_Jv.cpp
+${MODEL_DIR}/model_jakstat_adjoint_JvB.cpp
+${MODEL_DIR}/model_jakstat_adjoint_Jy.cpp
+${MODEL_DIR}/model_jakstat_adjoint_Jz.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dJydp.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dJydx.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dJydy.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dJzdp.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dJzdx.cpp
+${MODEL_DIR}/model_jakstat_adjoint_deltaqB.cpp
+${MODEL_DIR}/model_jakstat_adjoint_deltasx.cpp
+${MODEL_DIR}/model_jakstat_adjoint_deltax.cpp
+${MODEL_DIR}/model_jakstat_adjoint_deltaxB.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dsigma_ydp.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dsigma_zdp.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dwdp.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dwdx.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dxdotdp.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dydp.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dydx.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dzdp.cpp
+${MODEL_DIR}/model_jakstat_adjoint_dzdx.cpp
+${MODEL_DIR}/model_jakstat_adjoint_qBdot.cpp
+${MODEL_DIR}/model_jakstat_adjoint_root.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sJy.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sJz.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sigma_y.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sigma_z.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sroot.cpp
+${MODEL_DIR}/model_jakstat_adjoint_stau.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sx0.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sxdot.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sy.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sz.cpp
+${MODEL_DIR}/model_jakstat_adjoint_sz_tf.cpp
+${MODEL_DIR}/model_jakstat_adjoint_w.cpp
+${MODEL_DIR}/model_jakstat_adjoint_x0.cpp
+${MODEL_DIR}/model_jakstat_adjoint_xBdot.cpp
+${MODEL_DIR}/model_jakstat_adjoint_xdot.cpp
+${MODEL_DIR}/model_jakstat_adjoint_y.cpp
+${MODEL_DIR}/model_jakstat_adjoint_z.cpp
+)
+
+add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
+
+if(APPLE)
+    set(SHARED_OBJECT_EXTENSION .dylib)
+else()
+    set(SHARED_OBJECT_EXTENSION .so)
+endif()
+
+target_link_libraries(${PROJECT_NAME}
+"${SUNDIALS_LIB_DIR}/libsundials_nvecserial${SHARED_OBJECT_EXTENSION}"
+"${SUNDIALS_LIB_DIR}/libsundials_cvodes${SHARED_OBJECT_EXTENSION}"
+"${SUITESPARSE_LIB_DIR}/libcolamd${SHARED_OBJECT_EXTENSION}"
+"${SUITESPARSE_LIB_DIR}/libklu${SHARED_OBJECT_EXTENSION}"
+"${SUITESPARSE_LIB_DIR}/libbtf${SHARED_OBJECT_EXTENSION}"
+"${SUITESPARSE_LIB_DIR}/libamd${SHARED_OBJECT_EXTENSION}"
+"${SUITESPARSE_LIB_DIR}/libsuitesparseconfig${SHARED_OBJECT_EXTENSION}"
+"${HDF5_HL_LIBRARIES}"
+"${HDF5_C_LIBRARIES}"
+"-lpthread -ldl -lz"
+"-lm"
+)
+
+set(SRC_LIST_EXE main.cpp)
+
+add_executable(simulate_${PROJECT_NAME} ${SRC_LIST_EXE})
+
+target_link_libraries(simulate_${PROJECT_NAME} ${PROJECT_NAME})
+

models/model_jakstat_adjoint/main.cpp

@@ -0,0 +1,143 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <math.h>
+
+#include "wrapfunctions.h" /* model specific functions */
+#include <include/amici_interface_cpp.h> /* AMICI API */
+#include <include/ami_hdf5.h>  /* AMICI HDF5 I/O functions */
+#include <include/udata_accessors.h> /* Accessor macros for UserData members */
+#include <include/rdata_accessors.h> /* Accessor macros for ReturnData members */
+
+/* This is a scaffold for a stand-alone AMICI simulation executable demonstrating
+ * use of the AMICI C++ API.
+ *
+ * This program reads AMICI options from an HDF5 file, prints some results
+ * and writes additional results to an HDF5 file. The name of the HDF5 file
+ * is expected as single command line argument.
+ *
+ * An initial HDF5 file with the required fields can be generated using MATLAB by adding the following lines
+ * at the end of simulate_${MODEL_NAME}.m file just before the final "end":
+ *
+ *    %% Write data that is passed to AMICI to HDF5
+ *    hdffile = fullfile(pwd, 'mydata.h5');
+ *    structToHDF5Attribute(hdffile, '/options', options_ami);
+ *    h5writeatt(hdffile, '/options', 'ts', tout);
+ *    h5writeatt(hdffile, '/options', 'nt', numel(tout));
+ *    h5writeatt(hdffile, '/options', 'theta', theta);
+ *    h5writeatt(hdffile, '/options', 'kappa', kappa);
+ *    if(~isempty(data))
+ *      structToHDF5Attribute(hdffile, '/data', data);
+ *    end
+ *
+ * ... and then running a simulation from MATLAB as usual.
+ *
+ * Default UserData settings can be written to an HDF5 file with:
+ *     structToHDF5Attribute('test.h5', '/options', amioption())
+ */
+
+// Function prototypes
+void processReturnData(ReturnData *rdata, UserData *udata);
+void printReturnData(ReturnData *rdata, UserData *udata);
+
+
+int main(int argc, char **argv)
+{  
+    // HDF5 file to read and write data (full path)
+    const char *hdffile;
+
+    // Check command line arguments
+    if(argc != 2) {
+        fprintf(stderr, "Error: must provide HDF5 input file as first and only argument.\n");
+        return 1;
+    } else {
+        hdffile = argv[1];
+    }
+
+    // Read UserData (AMICI settings and model parameters) from HDF5 file
+    UserData *udata = AMI_HDF5_readSimulationUserDataFromFileName(hdffile, "/options");
+    if (udata == NULL) {
+        return 1;
+    }
+
+    // Read ExpData (experimental data for model) from HDF5 file
+    ExpData *edata = AMI_HDF5_readSimulationExpData(hdffile, udata, "/data");
+    if (edata == NULL) {
+        freeUserData(udata);
+        return 1;
+    }
+
+    // Run the simulation
+    int status = 0;
+    ReturnData *rdata = getSimulationResults(udata, edata, &status);
+    if (rdata == NULL) {
+        freeExpData(edata);
+        freeUserData(udata);
+        return 1;
+    }
+
+    // Do something with the simulation results
+    processReturnData(rdata, udata);
+
+    // Save simulation results to HDF5 file
+    AMI_HDF5_writeReturnData(rdata, udata, hdffile, "/solution");
+
+    // Free memory
+    freeExpData(edata);
+    freeUserData(udata);
+    freeReturnData(rdata);
+
+    return 0;
+}
+
+
+void processReturnData(ReturnData *rdata, UserData *udata) {
+    // show some the simulation results
+    printReturnData(rdata, udata);
+}
+
+void printReturnData(ReturnData *rdata, UserData *udata) {
+    //Print of some the simulation results
+
+    printf("Timepoints (tsdata): ");
+    printArray(tsdata, nt);
+
+    printf("\n\nStates (xdata):\n");
+    for(int i = 0; i < nx; ++i) {
+        for(int j = 0; j < nt; ++j)
+            printf("%e\t", rdata->am_xdata[j +  nt * i]);
+        printf("\n");
+    }
+
+    printf("\nObservables (ydata):\n");
+    for(int i = 0; i < ny; ++i) {
+        for(int j = 0; j < nt; ++j)
+            printf("%e\t", rdata->am_ydata[j +  nt * i]);
+        printf("\n");
+    }
+
+    printf("\n\ndx/dt (xdotdata):\n");
+    for(int i = 0; i < nx; ++i)
+        printf("%e\t", rdata->am_xdotdata[i]);
+
+//    printf("\nJacobian (jdata)\n");
+//    for(int i = 0; i < nx; ++i) {
+//        for(int j = 0; j < nx; ++j)
+//            printf("%e\t", rdata->am_Jdata[i + nx * j]);
+//        printf("\n");
+//    }
+
+    printf("\nnumsteps: \t\t");
+    printfArray(numstepsdata, nt, "%.0f ");
+
+    printf("\nnumrhsevalsdata: \t");
+    printfArray(numrhsevalsdata, nt, "%.0f ");
+
+    printf("\norder: \t\t");
+    printfArray(orderdata, nt, "%.0f ");
+
+    printf("\n");
+    printf("Loglikelihood (llhdata): %e\n", *llhdata);
+}
+