Patternstim fix : allow restore phases to provide new pattern-stim (#121)

* Allow restore phases to provide new pattern stim via CLI
* Add save-restore test with patternstim
* NOTE : if restore pattern stim file has event at the checkpoint time
  boundary then the event could be missed. 
  See https://bbpteam.epfl.ch/project/issues/browse/CNEUR-313

Author: pramodk

coreneuron/apps/main1.cpp

@@ -265,7 +265,7 @@ void nrn_init_and_load_data(int argc,
 
     // Invoke PatternStim
     if (!corenrn_param.patternstim.empty()) {
-        nrn_mkPatternStim(corenrn_param.patternstim.c_str());
+        nrn_mkPatternStim(corenrn_param.patternstim.c_str(), corenrn_param.tstop);
     }
 
     /// Setting the timeout

coreneuron/mechanism/mech/modfile/pattern.mod

@@ -116,6 +116,8 @@ void pattern_stim_setup_helper(int size, double* tv, int* gv, _threadargsproto_)
 	info->size = size;
 	info->tvec = tv;
 	info->gidvec = gv;
+    // initiate event chain (needed in case of restore)
+	artcell_net_send ( _tqitem, -1, (Point_process*) _nt->_vdata[_ppvar[1*_STRIDE]], t +  0.0 , 1.0 ) ;
 }
 } // namespace coreneuron
 ENDVERBATIM

coreneuron/mechanism/patternstim.cpp

@@ -58,7 +58,7 @@ extern void pattern_stim_setup_helper(int size,
                                       NrnThread* _nt,
                                       double v);
 
-static size_t read_raster_file(const char* fname, double** tvec, int** gidvec);
+static size_t read_raster_file(const char* fname, double** tvec, int** gidvec, double tstop);
 
 int nrn_extra_thread0_vdata;
 
@@ -77,7 +77,7 @@ void nrn_set_extra_thread0_vdata() {
 
 // fname is the filename of an output_spikes.h format raster file.
 // todo : add function for memory cleanup (to be called at the end of simulation)
-void nrn_mkPatternStim(const char* fname) {
+void nrn_mkPatternStim(const char* fname, double tstop) {
     int type = nrn_get_mechtype("PatternStim");
     if (!corenrn.get_memb_func(type).sym) {
         printf("nrn_set_extra_thread_vdata must be called (after mk_mech, and before nrn_setup\n");
@@ -93,7 +93,7 @@ void nrn_mkPatternStim(const char* fname) {
     int* gidvec;
 
     // todo : handle when spike raster will be very large (int < size_t)
-    size_t size = read_raster_file(fname, &tvec, &gidvec);
+    size_t size = read_raster_file(fname, &tvec, &gidvec, tstop);
 
     Point_process* pnt = nrn_artcell_instantiate("PatternStim");
     NrnThread* nt = nrn_threads + pnt->_tid;
@@ -123,7 +123,8 @@ static bool spike_comparator(const spike_type& l, const spike_type& r) {
     return l.first < r.first;
 }
 
-size_t read_raster_file(const char* fname, double** tvec, int** gidvec) {
+
+size_t read_raster_file(const char* fname, double** tvec, int** gidvec, double tstop) {
     FILE* f = fopen(fname, "r");
     nrn_assert(f);
 
@@ -138,7 +139,9 @@ size_t read_raster_file(const char* fname, double** tvec, int** gidvec) {
     int gid;
 
     while (fscanf(f, "%lf %d\n", &stime, &gid) == 2) {
-        spikes.push_back(std::make_pair(stime, gid));
+        if ( stime >= t && stime <= tstop) {
+            spikes.push_back(std::make_pair(stime, gid));
+        }
     }
 
     fclose(f);

coreneuron/nrniv/nrniv_decl.h

@@ -70,7 +70,7 @@ extern int nrn_setup_extracon;
 extern void nrn_cleanup();
 extern void nrn_cleanup_ion_map();
 extern void BBS_netpar_solve(double);
-extern void nrn_mkPatternStim(const char* filename);
+extern void nrn_mkPatternStim(const char* filename, double tstop);
 extern int nrn_extra_thread0_vdata;
 extern void nrn_set_extra_thread0_vdata(void);
 extern Point_process* nrn_artcell_instantiate(const char* mechname);

tests/jenkins/Jenkinsfile

@@ -320,7 +320,7 @@ pipeline {
                         }
                     }
                 }
-                stage('patstim'){
+                stage('patstim save restore'){
                     stages{
                         stage('neuron'){
                             steps{

tests/jenkins/run_corenrn.sh

@@ -18,7 +18,31 @@ else
 fi
 
 if [ "${TEST}" = "patstim" ]; then
-    mpirun -n ${MPI_RANKS} ./${CORENRN_TYPE}/special-core --mpi -e 100 --pattern patstim.spk -d test${TEST}dat -o ${TEST}
+    # first run full run
+    mpirun -n ${MPI_RANKS} ./${CORENRN_TYPE}/special-core --mpi -e 100 --pattern patstim.spk -d testpatstimdat -o ${TEST}
+
+    # split patternstim file into two parts : total 2000 events, split at line no 1001 i.e. before events for 50 msec
+    echo 1000 > patstim.1.spk
+    sed -n 2,1001p patstim.spk  >> patstim.1.spk
+    echo 1000 > patstim.2.spk
+    sed -n 1002,2001p patstim.spk  >> patstim.2.spk
+
+    # run test with checkpoint : part 1
+    mpirun -n ${MPI_RANKS} ./${CORENRN_TYPE}/special-core --mpi -e 49 --pattern patstim.1.spk -d testpatstimdat -o ${TEST}_part1 --checkpoint checkpoint
+
+    # run test with restore : part 2
+    mpirun -n ${MPI_RANKS} ./${CORENRN_TYPE}/special-core --mpi -e 100 --pattern patstim.2.spk -d testpatstimdat -o ${TEST}_part2 --restore checkpoint
+
+    # run additional restore by providing full patternstim in part 2 : part 3
+    mpirun -n ${MPI_RANKS} ./${CORENRN_TYPE}/special-core --mpi -e 100 --pattern patstim.spk -d testpatstimdat -o ${TEST}_part3 --restore checkpoint
+
+    # part 2 and part 3 should be same (part 3 ignore extra events in pattern stim)
+    diff -w  ${TEST}_part2/out.dat  ${TEST}_part3/out.dat
+
+    # combine spikes from part1 and part2 should be same as full run
+    cat ${TEST}_part1/out.dat  ${TEST}_part2/out.dat > ${TEST}/out.saverestore.dat
+    diff -w ${TEST}/out.dat ${TEST}/out.saverestore.dat
+
 elif [ "${TEST}" = "ringtest" ]; then
     mpirun -n ${MPI_RANKS} ./${CORENRN_TYPE}/special-core --mpi -e 100 -d coredat -o ${TEST}
 elif [ "${TEST}" = "tqperf" ]; then