Bi-directional in-memory transfer between NEURON and CoreNEURON (#182)

 - Direct transfer mode disable error checking for real cell with no gid.
 - All the fixed_record items become trajectory requests for direct mode CoreNEURON
    In particular, any Vector.record specified in NEURON will be recorded
    during a direct mode  CoreNEURON run (presently limited to voltage and non
    artificial cell mechanism data).
 - Direct mode transfer allows CoreNEURON to record entire Vector trajectory.
 - An Graph line expression may depend on several range variables.
   So the number of PlayRecord items in NEURON may not be the same
   as the number of trajectories to record in CoreNEURON.
 - NEURON determines if trajectory direct transfer is per step or buffered.
 - NEURON trajectory requests have to be received prior to finitialize.
 - Special mtype for tvec trajectory
 - Trajectory transfer scatters to NEURON variables in per step mode.
    Instead of sending back array of values from CoreNEURON.

Author: nrnhines

coreneuron/nrniv/global_vars.cpp

@@ -82,7 +82,7 @@ void set_globals(const char* path, bool cli_global_seed, int cli_global_seed_val
         double val;
         int n;
 
-        fscanf(f, "%s\n", line);
+        nrn_assert(fscanf(f, "%s\n", line) == 1);
         check_bbcore_write_version(line);
 
         for (;;) {

coreneuron/nrniv/main1.cpp

@@ -64,6 +64,8 @@ const char* corenrn_version() {
     return coreneuron::bbcore_write_version;
 }
 
+void (*nrn2core_part2_clean_)();
+
 #ifdef ISPC_INTEROP
 // cf. utils/ispc_globals.c
 extern double ispc_celsius;
@@ -358,6 +360,58 @@ void handle_forward_skip(double forwardskip, int prcellgid) {
 const char* nrn_version(int) {
     return "version id unimplemented";
 }
+
+// bsize = 0 then per step transfer
+// bsize > 1 then full trajectory save into arrays.
+void get_nrn_trajectory_requests(int bsize) {
+  if (nrn2core_get_trajectory_requests_) {
+    for (int tid=0; tid < nrn_nthread; ++tid) {
+      NrnThread& nt = nrn_threads[tid];
+      int n_pr;
+      int n_trajec;
+      int* types;
+      int* indices;
+      void** vpr;
+      double** varrays;
+      double** pvars;
+
+      // bsize is passed by reference, the return value will determine if
+      // per step return or entire trajectory return.
+      (*nrn2core_get_trajectory_requests_)(tid, bsize, n_pr, vpr, n_trajec, types, indices, pvars, varrays);
+      delete_trajectory_requests(nt);
+      if (n_trajec) {
+        TrajectoryRequests* tr = new TrajectoryRequests;
+        nt.trajec_requests = tr;
+        tr->bsize = bsize;
+        tr->n_pr = n_pr;
+        tr->n_trajec = n_trajec;
+        tr->vsize = 0;
+        tr->vpr = vpr;
+        tr->gather = new double*[n_trajec];
+        tr->varrays = varrays;
+        tr->scatter = pvars;
+        for (int i=0; i < n_trajec; ++i) {
+          tr->gather[i] = stdindex2ptr(types[i], indices[i], nt);
+        }
+        delete [] types;
+        delete [] indices;
+      }
+    }
+  }
+}
+
+static void trajectory_return() {
+  if (nrn2core_trajectory_return_) {
+    for (int tid=0; tid < nrn_nthread; ++tid) {
+      NrnThread& nt = nrn_threads[tid];
+      TrajectoryRequests* tr = nt.trajec_requests;
+      if (tr && tr->varrays) {
+        (*nrn2core_trajectory_return_)(tid, tr->n_pr, tr->vsize, tr->vpr, nt._t);
+      }
+    }
+  }
+}
+
 }  // namespace coreneuron
 
 /// The following high-level functions are marked as "extern C"
@@ -423,15 +477,6 @@ extern "C" int run_solve_core(int argc, char** argv) {
     // clang-format on
     {
         double v = nrnopt_get_dbl("--voltage");
-
-        // TODO : if some ranks are empty then restore will go in deadlock
-        // phase (as some ranks won't have restored anything and hence return
-        // false in checkpoint_initialize
-        if (!checkpoint_initialize()) {
-            nrn_finitialize(v != 1000., v);
-        }
-
-        report_mem_usage("After nrn_finitialize");
         double dt = nrnopt_get_dbl("--dt");
         double delay = nrnopt_get_dbl("--mindelay");
         double tstop = nrnopt_get_dbl("--tstop");
@@ -442,6 +487,24 @@ extern "C" int run_solve_core(int argc, char** argv) {
             abort();
         }
 
+        // In direct mode there are likely trajectory record requests
+        // to allow processing in NEURON after simulation by CoreNEURON
+        if (corenrn_embedded) {
+            // arg is vector size required but NEURON can instead
+            // specify that returns will be on a per time step basis.
+            get_nrn_trajectory_requests(int(tstop/dt) + 2);
+            (*nrn2core_part2_clean_)();
+        }
+
+        // TODO : if some ranks are empty then restore will go in deadlock
+        // phase (as some ranks won't have restored anything and hence return
+        // false in checkpoint_initialize
+        if (!checkpoint_initialize()) {
+            nrn_finitialize(v != 1000., v);
+        }
+
+        report_mem_usage("After nrn_finitialize");
+
         // register all reports into reportinglib
         double min_report_dt = INT_MAX;
         int report_buffer_size = nrnopt_get_int("--report-buffer-size");
@@ -473,6 +536,10 @@ extern "C" int run_solve_core(int argc, char** argv) {
         BBS_netpar_solve(nrnopt_get_dbl("--tstop"));
         Instrumentor::phase_end("simulation");
         Instrumentor::stop_profile();
+        // direct mode and full trajectory gathering on CoreNEURON, send back.
+        if (corenrn_embedded) {
+            trajectory_return();
+        }
         // Report global cell statistics
         report_cell_stats();
 

coreneuron/nrniv/mk_mech.cpp

@@ -118,11 +118,14 @@ void mk_mech(const char* datpath) {
 
 // we are embedded in NEURON, get info as stringstream from nrnbbcore_write.cpp
 static void mk_mech() {
+    static bool done = false;
+    if (done) { return; }
     nrn_need_byteswap = 0;
     std::stringstream ss;
     nrn_assert(nrn2core_mkmech_info_);
     (*nrn2core_mkmech_info_)(ss);
     mk_mech(ss);
+    done = true;
 }
 
 static void mk_mech(std::istream& s) {

coreneuron/nrniv/nrn2core_direct.h

@@ -87,6 +87,33 @@ extern int (*nrn2core_get_dat2_vecplay_inst_)(int tid,
                                               int& sz,
                                               double*& yvec,
                                               double*& tvec);
+
+extern void (*nrn2core_part2_clean_)();
+
+/* what variables to send back to NEURON on each time step */
+extern void (*nrn2core_get_trajectory_requests_)(int tid,
+                                                int& bsize,
+                                                int& n_pr,
+                                                void**& vpr,
+                                                int& n_trajec,
+                                                int*& types,
+                                                int*& indices,
+                                                double**& pvars,
+                                                double**& varrays);
+
+/* send values to NEURON on each time step */
+extern void (*nrn2core_trajectory_values_)(int tid,
+                                          int n_pr,
+                                          void** vpr,
+                                          double t);
+
+/* Filled the Vector data arrays and send back the sizes at end of run */
+extern void (*nrn2core_trajectory_return_)(int tid,
+                                          int n_pr,
+                                          int vecsz,
+                                          void** vpr,
+                                          double t);
+
 }
 
 #endif /* nrn2core_direct_h */

coreneuron/nrniv/nrn_setup.cpp

@@ -128,6 +128,28 @@ int (*nrn2core_get_dat2_vecplay_inst_)(int tid,
                                        double*& yvec,
                                        double*& tvec);
 
+void (*nrn2core_get_trajectory_requests_)(int tid,
+                                                int& bsize,
+                                                int& n_pr,
+                                                void**& vpr,
+                                                int& n_trajec,
+                                                int*& types,
+                                                int*& indices,
+                                                double**& pvars,
+                                                double**& varrays);
+
+void (*nrn2core_trajectory_values_)(int tid,
+                                          int n_pr,
+                                          void** vpr,
+                                          double t);
+
+void (*nrn2core_trajectory_return_)(int tid,
+                                          int n_pr,
+                                          int vecsz,
+                                          void** vpr,
+                                          double t);
+
+
 // file format defined in cooperation with nrncore/src/nrniv/nrnbbcore_write.cpp
 // single integers are ascii one per line. arrays are binary int or double
 // Note that regardless of the gid contents of a group, since all gids are
@@ -312,7 +334,7 @@ void nrn_read_filesdat(int& ngrp, int*& grp, int multiple, int*& imult, const ch
     }
 
     char version[256];
-    fscanf(fp, "%s\n", version);
+    nrn_assert(fscanf(fp, "%s\n", version) == 1);
     check_bbcore_write_version(version);
 
     int iNumFiles = 0;
@@ -350,8 +372,8 @@ void nrn_read_filesdat(int& ngrp, int*& grp, int multiple, int*& imult, const ch
         if ((iNum + 1) % iNumFiles == 0) {
             // re-read file for each multiple (skipping the two header lines)
             rewind(fp);
-            fscanf(fp, "%*s\n");
-            fscanf(fp, "%*d\n");
+            nrn_assert(fscanf(fp, "%*s\n") == 0);
+            nrn_assert(fscanf(fp, "%*d\n") == 0);
         }
     }
 
@@ -681,7 +703,9 @@ void nrn_setup(const char* filesdat,
     maxgid = 0x7fffffff / nrn_setup_multiple;
     nrn_read_filesdat(ngroup, gidgroups, nrn_setup_multiple, imult, filesdat);
 
-    MUTCONSTRUCT(1)
+    if (!MUTCONSTRUCTED) {
+        MUTCONSTRUCT(1)
+    }
     // temporary bug work around. If any process has multiple threads, no
     // process can have a single thread. So, for now, if one thread, make two.
     // Fortunately, empty threads work fine.
@@ -795,7 +819,12 @@ void nrn_setup(const char* filesdat,
     delete[] file_reader;
 
     model_size();
-    delete[] gidgroups;
+    delete [] gidgroups;
+    delete [] imult;
+    if (nrnthread_chkpnt) {
+        delete [] nrnthread_chkpnt;
+        nrnthread_chkpnt = NULL;
+    }
 
     if (nrnmpi_myid == 0) {
         printf(" Setup Done   : %.2lf seconds \n", nrn_wtime() - time);
@@ -900,6 +929,34 @@ int nrn_i_layout(int icnt, int cnt, int isz, int sz, int layout) {
     return 0;
 }
 
+// take into account alignment, layout, permutation
+// only voltage or mechanism data index allowed. (mtype 0 means time)
+double* stdindex2ptr(int mtype, int index, NrnThread& nt) {
+    if (mtype == -5) { // voltage
+        int v0 = nt._actual_v - nt._data;
+        int ix = index;  // relative to _actual_v
+        nrn_assert((ix >= 0) && (ix < nt.end));
+        if (nt._permute) {
+             node_permute(&ix, 1, nt._permute);
+        }
+        return nt._data + (v0 + ix);  // relative to nt._data
+    }else if (mtype  > 0 && mtype < n_memb_func) { // 
+        Memb_list* ml = nt._ml_list[mtype];
+        nrn_assert(ml);
+        int ix = nrn_param_layout(index, mtype, ml);
+        if (ml->_permute) {
+            ix = nrn_index_permute(ix, mtype, ml);
+        }
+        return ml->data + ix;
+    }else if (mtype == 0) { // time
+        return &nt._t;
+    }else{
+        printf("stdindex2ptr does not handle mtype=%d\n", mtype);
+        nrn_assert(0);
+    }
+    return NULL;
+}
+
 // from i to (icnt, isz)
 void nrn_inverse_i_layout(int i, int& icnt, int cnt, int& isz, int sz, int layout) {
     if (layout == 1) {
@@ -948,6 +1005,7 @@ inline void mech_layout(FileHandler& F, T* data, int cnt, int sz, int layout) {
  * things up first. */
 
 void nrn_cleanup(bool clean_ion_global_map) {
+    clear_event_queue(); // delete left-over TQItem
     gid2in.clear();
     gid2out.clear();
 
@@ -967,6 +1025,7 @@ void nrn_cleanup(bool clean_ion_global_map) {
     for (int it = 0; it < nrn_nthread; ++it) {
         NrnThread* nt = nrn_threads + it;
         NrnThreadMembList* next_tml = NULL;
+	delete_trajectory_requests(*nt);
         for (NrnThreadMembList* tml = nt->tml; tml; tml = next_tml) {
             Memb_list* ml = tml->ml;
 
@@ -1114,9 +1173,24 @@ void nrn_cleanup(bool clean_ion_global_map) {
 
     if (pnttype2presyn) {
         free(pnttype2presyn);
+        pnttype2presyn = NULL;
     }
 }
 
+void delete_trajectory_requests(NrnThread& nt) {
+  if (nt.trajec_requests) {
+    TrajectoryRequests* tr = nt.trajec_requests;
+    if (tr->n_trajec) {
+      delete [] tr->vpr;
+      if (tr->scatter) {delete [] tr->scatter;}
+      if (tr->varrays) {delete [] tr->varrays;}
+      delete [] tr->gather;
+    }
+    delete nt.trajec_requests;
+    nt.trajec_requests = NULL;
+  }
+}
+
 void read_phase2(FileHandler& F, int imult, NrnThread& nt) {
     bool direct = corenrn_embedded ? true : false;
     if (!direct) {
@@ -1178,7 +1252,9 @@ void read_phase2(FileHandler& F, int imult, NrnThread& nt) {
     ntc.n_outputgids = n_outputgid;
     ntc.nmech = nmech;
 #endif
-    nrn_assert(n_outputgid > 0);  // avoid n_outputgid unused warning
+    if (!direct) {
+        nrn_assert(n_outputgid > 0);  // avoid n_outputgid unused warning
+    }
 
     /// Checkpoint in coreneuron is defined for both phase 1 and phase 2 since they are written
     /// together

coreneuron/nrniv/nrniv_decl.h

@@ -55,6 +55,8 @@ extern void nrn_setup(const char* filesdat,
                       bool is_mapping_needed,
                       int byte_swap,
                       bool run_setup_cleanup = true);
+extern double* stdindex2ptr(int mtype, int index, NrnThread&);
+extern void delete_trajectory_requests(NrnThread&);
 extern int nrn_setup_multiple;
 extern int nrn_setup_extracon;
 extern void nrn_cleanup(bool clean_ion_global_map = true);

coreneuron/nrniv/output_spikes.cpp

@@ -57,7 +57,9 @@ static MUTDEC
     } catch (const std::length_error& le) {
         std::cerr << "Lenght error" << le.what() << std::endl;
     }
-    MUTCONSTRUCT(1);
+    if (!MUTCONSTRUCTED) {
+        MUTCONSTRUCT(1);
+    }
 }
 
 void spikevec_lock() {

coreneuron/nrniv/tqueue.ipp

@@ -78,6 +78,11 @@ TQueue<C>::~TQueue() {
     }
     delete binq_;
 
+    if (least_) {
+      delete least_;
+      least_ = NULL;
+    }
+
     /// Clear the splay tree
     while ((q = spdeq(&sptree_->root)) != NULL) {
         delete q;