Allow multiple gids per cell that reference different voltages. (#829)

* Allow multiple gid for a cell (no more than on per PreSyn)
* data version 1.6 updated for integration ring and ring_gap
* checkpoint updated for data format 1.6
* nt.n_presyn is no longer needs to be nt.ncell
  - so allocate presyns only if we have non-zero presyns
  - note that non-zero length allocations aren't nullptr
* fix for ispc: add missing n_real_output field in NrnThread for ISPC (and minor reformating)

Co-authored-by: Pramod Kumbhar <pramod.s.kumbhar@gmail.com>

Author: nrnhines

coreneuron/io/nrn2core_direct.h

@@ -33,6 +33,7 @@ extern int (*nrn2core_get_dat1_)(int tid,
                                  std::vector<int>& netcon_negsrcgid_tid);
 
 extern int (*nrn2core_get_dat2_1_)(int tid,
+                                   int& n_real_cell,
                                    int& ngid,
                                    int& n_real_gid,
                                    int& nnode,

coreneuron/io/nrn_checkpoint.cpp

@@ -131,12 +131,13 @@ void CheckPoints::write_phase2(NrnThread& nt) const {
         }
     }
 
+    fh << nt.ncell << " ncell\n";
     fh << n_outputgid << " ngid\n";
 #if CHKPNTDEBUG
     assert(ntc.n_outputgids == n_outputgid);
 #endif
 
-    fh << nt.ncell << " n_real_gid\n";
+    fh << nt.n_real_output << " n_real_output\n";
     fh << nt.end << " nnode\n";
     fh << ((nt._actual_diam == nullptr) ? 0 : nt.end) << " ndiam\n";
     int nmech = 0;
@@ -316,7 +317,7 @@ void CheckPoints::write_phase2(NrnThread& nt) const {
     int nnetcon = nt.n_netcon;
 
     int* output_vindex = new int[nt.n_presyn];
-    double* output_threshold = new double[nt.ncell];
+    double* output_threshold = new double[nt.n_real_output];
     for (int i = 0; i < nt.n_presyn; ++i) {
         PreSyn* ps = nt.presyns + i;
         if (ps->thvar_index_ >= 0) {
@@ -327,8 +328,8 @@ void CheckPoints::write_phase2(NrnThread& nt) const {
             assert(ps->pntsrc_ == nullptr);
             output_threshold[i] = ps->threshold_;
             output_vindex[i] = pinv_nt[ps->thvar_index_];
-        } else if (i < nt.ncell) {      // real cell without a presyn
-            output_threshold[i] = 0.0;  // the way it was set in nrnbbcore_write.cpp
+        } else if (i < nt.n_real_output) {  // real cell without a presyn
+            output_threshold[i] = 0.0;      // the way it was set in nrnbbcore_write.cpp
             output_vindex[i] = -1;
         } else {
             Point_process* pnt = ps->pntsrc_;
@@ -347,12 +348,12 @@ void CheckPoints::write_phase2(NrnThread& nt) const {
         }
     }
     fh.write_array<int>(output_vindex, nt.n_presyn);
-    fh.write_array<double>(output_threshold, nt.ncell);
+    fh.write_array<double>(output_threshold, nt.n_real_output);
 #if CHKPNTDEBUG
     for (int i = 0; i < nt.n_presyn; ++i) {
         assert(ntc.output_vindex[i] == output_vindex[i]);
     }
-    for (int i = 0; i < nt.ncell; ++i) {
+    for (int i = 0; i < nt.n_real_output; ++i) {
         assert(ntc.output_threshold[i] == output_threshold[i]);
     }
 #endif

coreneuron/io/nrn_setup.cpp

@@ -102,7 +102,7 @@ void (*nrn2core_all_weights_return_)(std::vector<double*>& weights);
 // for the negative gids in netcon_srcgid (in that order) the source thread.
 //
 // <firstgid>_2.dat
-// n_output n_real_output, nnode
+// n_real_cell, n_output, n_real_output, nnode
 // ndiam - 0 if no mechanism has dparam with diam semantics, or nnode
 // nmech - includes artcell mechanisms
 // for the nmech tml mechanisms

coreneuron/io/phase1.cpp

@@ -67,9 +67,12 @@ void Phase1::populate(NrnThread& nt, OMP_Mutex& mut) {
     coreneuron::nrnthreads_netcon_negsrcgid_tid[nt.id] = this->netcon_negsrcgid_tid;
 
     nt.netcons = new NetCon[nt.n_netcon];
-    nt.presyns_helper = (PreSynHelper*) ecalloc_align(nt.n_presyn, sizeof(PreSynHelper));
 
-    nt.presyns = new PreSyn[nt.n_presyn];
+    if (nt.n_presyn) {
+        nt.presyns_helper = (PreSynHelper*) ecalloc_align(nt.n_presyn, sizeof(PreSynHelper));
+        nt.presyns = new PreSyn[nt.n_presyn];
+    }
+
     PreSyn* ps = nt.presyns;
     /// go through all presyns
     for (auto& gid: this->output_gids) {

coreneuron/io/phase2.cpp

@@ -24,6 +24,7 @@
 #endif
 
 int (*nrn2core_get_dat2_1_)(int tid,
+                            int& n_real_cell,
                             int& ngid,
                             int& n_real_gid,
                             int& nnode,
@@ -107,6 +108,7 @@ inline void mech_data_layout_transform(T* data, int cnt, int sz, int layout) {
 }
 
 void Phase2::read_file(FileHandler& F, const NrnThread& nt) {
+    n_real_cell = F.read_int();
     n_output = F.read_int();
     n_real_output = F.read_int();
     n_node = F.read_int();
@@ -256,6 +258,7 @@ void Phase2::read_direct(int thread_id, const NrnThread& nt) {
     int* nodecounts_ = nullptr;
     int n_weight;
     (*nrn2core_get_dat2_1_)(thread_id,
+                            n_real_cell,
                             n_output,
                             n_real_output,
                             n_node,
@@ -915,8 +918,9 @@ void Phase2::populate(NrnThread& nt, const UserParams& userParams) {
     NrnThreadChkpnt& ntc = nrnthread_chkpnt[nt.id];
     ntc.file_id = userParams.gidgroups[nt.id];
 
-    nt.ncell = n_real_output;
+    nt.ncell = n_real_cell;
     nt.end = n_node;
+    nt.n_real_output = n_real_output;
 
 #if CHKPNTDEBUG
     ntc.n_outputgids = n_output;
@@ -1221,14 +1225,15 @@ void Phase2::populate(NrnThread& nt, const UserParams& userParams) {
         node_permute(output_vindex.data(), nt.n_presyn, nt._permute);
     }
 #if CHKPNTDEBUG
-    ntc.output_threshold = new double[nt.ncell];
-    memcpy(ntc.output_threshold, output_threshold.data(), nt.ncell * sizeof(double));
+    ntc.output_threshold = new double[n_real_output];
+    memcpy(ntc.output_threshold, output_threshold.data(), n_real_output * sizeof(double));
 #endif
+
     for (int i = 0; i < nt.n_presyn; ++i) {  // real cells
         PreSyn* ps = nt.presyns + i;
 
         int ix = output_vindex[i];
-        if (ix == -1 && i < nt.ncell) {  // real cell without a presyn
+        if (ix == -1 && i < n_real_output) {  // real cell without a presyn
             continue;
         }
         if (ix < 0) {
@@ -1256,7 +1261,7 @@ void Phase2::populate(NrnThread& nt, const UserParams& userParams) {
     // initial net_send_buffer size about 1% of number of presyns
     // nt._net_send_buffer_size = nt.ncell/100 + 1;
     // but, to avoid reallocation complexity on GPU ...
-    nt._net_send_buffer_size = nt.ncell;
+    nt._net_send_buffer_size = n_real_output;
     nt._net_send_buffer = (int*) ecalloc_align(nt._net_send_buffer_size, sizeof(int));
 
     int nnetcon = nt.n_netcon;

coreneuron/io/phase2.hpp

@@ -89,6 +89,7 @@ class Phase2 {
                               NrnThreadChkpnt& ntc);
     void set_vec_play(NrnThread& nt, NrnThreadChkpnt& ntc);
 
+    int n_real_cell;
     int n_output;
     int n_real_output;
     int n_node;

coreneuron/mechanism/nrnoc_ml.ispc

@@ -104,15 +104,22 @@ struct NrnThread {
     NetCon* uniform netcons;
     double* uniform weights;
 
-    uniform int n_pntproc, n_presyn, n_input_presyn, n_netcon, n_weight;
+    uniform int n_pntproc;
+    uniform int n_weight;
+    uniform int n_netcon;
+    uniform int n_input_presyn;
+    uniform int n_presyn;
+    uniform int n_real_output;
 
     uniform int ncell;
     uniform int end;
     uniform int id;
     uniform int _stop_stepping;
     uniform int n_vecplay;
 
-    uniform unsigned int64 _ndata, _nidata, _nvdata;
+    uniform unsigned int64 _ndata;
+    uniform unsigned int64 _nidata;
+    uniform unsigned int64 _nvdata;
     double* uniform _data;
     int* uniform _idata;
 

coreneuron/network/netcvode.cpp

@@ -539,7 +539,7 @@ void NetCvode::check_thresh(NrnThread* nt) {  // for default method
         nt [0:1], presyns_helper [0:nt->n_presyn], presyns [0:nt->n_presyn], actual_v [0:nt->end])
                        copy(net_send_buf_count) if (nt->compute_gpu) async(nt->stream_id))
     nrn_pragma_omp(target teams distribute parallel for map(tofrom: net_send_buf_count) if(nt->compute_gpu))
-    for (int i = 0; i < nt->ncell; ++i) {
+    for (int i = 0; i < nt->n_real_output; ++i) {
         PreSyn* ps = presyns + i;
         PreSynHelper* psh = presyns_helper + i;
         int idx = 0;

coreneuron/sim/multicore.hpp

@@ -90,7 +90,8 @@ struct NrnThread: public MemoryManaged {
     int n_weight = 0;
     int n_netcon = 0;
     int n_input_presyn = 0;
-    int n_presyn = 0;  // only for model_size
+    int n_presyn = 0;       // only for model_size
+    int n_real_output = 0;  // for checking their thresholds.
 
     int ncell = 0; /* analogous to old rootnodecount */
     int end = 0;   /* 1 + position of last in v_node array. Now v_node_count. */

coreneuron/utils/nrnoc_aux.cpp

@@ -21,8 +21,7 @@ int v_structure_change;
 int diam_changed;
 #define MAXERRCOUNT 5
 int hoc_errno_count;
-const char* bbcore_write_version = "1.5";  // Generalize POINTER transfer to allow pointing to any
-                                           // RANGE variable
+const char* bbcore_write_version = "1.6";  // Allow multiple gid and PreSyn per real cell.
 
 char* pnt_name(Point_process* pnt) {
     return corenrn.get_memb_func(pnt->_type).sym;