Extend POINTER transfer to any RANGE variable in a NRN_THREAD

coreneuron/apps/main1.cpp

@@ -608,8 +608,11 @@ extern "C" int run_solve_core(int argc, char** argv) {
         if (corenrn_param.report_buff_size != corenrn_param.report_buff_size_default) {
             set_report_buffer_size(corenrn_param.report_buff_size);
         }
-        setup_report_engine(min_report_dt, delay);
-        configs.clear();
+
+        if (!configs.empty()) {
+            setup_report_engine(min_report_dt, delay);
+            configs.clear();
+        }
 
         // call prcellstate for prcellgid
         call_prcellstate_for_prcellgid(corenrn_param.prcellgid, compute_gpu, 0);

coreneuron/io/nrn2core_direct.h

@@ -9,6 +9,7 @@
 #pragma once
 
 #include <iostream>
+#include <vector>
 
 extern "C" {
 // The callbacks into nrn/src/nrniv/nrnbbcore_write.cpp to get
@@ -56,7 +57,8 @@ extern int (*nrn2core_get_dat2_mech_)(int tid,
                                       int dsz_inst,
                                       int*& nodeindices,
                                       double*& data,
-                                      int*& pdata);
+                                      int*& pdata,
+                                      std::vector<int>& pointer2type);
 
 extern int (*nrn2core_get_dat2_3_)(int tid,
                                    int nweight,

coreneuron/io/nrn_checkpoint.cpp

@@ -88,6 +88,33 @@ void CheckPoints::write_checkpoint(NrnThread* nt, int nb_threads) const {
 #endif
 }
 
+// Factor out the body of ion handling below as the same code
+// handles POINTER
+static int nrn_original_aos_index(int etype, int ix, NrnThread& nt, int** ml_pinv) {
+    // Determine ei_instance and ei from etype and ix.
+    // Deal with existing permutation and SoA.
+    Memb_list* eml = nt._ml_list[etype];
+    int ecnt = eml->nodecount;
+    int esz = corenrn.get_prop_param_size()[etype];
+    int elayout = corenrn.get_mech_data_layout()[etype];
+    // current index into eml->data is a  function
+    // of elayout, eml._permute, ei_instance, ei, and
+    // eml padding.
+    int p = ix - (eml->data - nt._data);
+    assert(p >= 0 && p < eml->_nodecount_padded * esz);
+    int ei_instance, ei;
+    nrn_inverse_i_layout(p, ei_instance, ecnt, ei, esz, elayout);
+    if (elayout == Layout::SoA) {
+        if (eml->_permute) {
+            if (!ml_pinv[etype]) {
+                ml_pinv[etype] = inverse_permute(eml->_permute, eml->nodecount);
+            }
+            ei_instance = ml_pinv[etype][ei_instance];
+        }
+    }
+    return ei_instance * esz + ei;
+}
+
 void CheckPoints::write_phase2(NrnThread& nt) const {
     FileHandler fh;
 
@@ -186,7 +213,7 @@ void CheckPoints::write_phase2(NrnThread& nt) const {
     }
 
     auto& memb_func = corenrn.get_memb_funcs();
-    // will need the ml_pinv inverse permutation of ml._permute for ions
+    // will need the ml_pinv inverse permutation of ml._permute for ions and POINTER
     int** ml_pinv = (int**) ecalloc(memb_func.size(), sizeof(int*));
 
     for (NrnThreadMembList* current_tml = nt.tml; current_tml; current_tml = current_tml->next) {
@@ -224,9 +251,10 @@ void CheckPoints::write_phase2(NrnThread& nt) const {
 
         sz = nrn_prop_dparam_size_[type];
         if (sz) {
-            int* d = soa2aos(ml->pdata, cnt, sz, layout, ml->_permute);
             // need to update some values according to Datum semantics.
-            if (!nrn_is_artificial_[type])
+            int* d = soa2aos(ml->pdata, cnt, sz, layout, ml->_permute);
+            std::vector<int> pointer2type;  // voltage or mechanism type (starts empty)
+            if (!nrn_is_artificial_[type]) {
                 for (int i_instance = 0; i_instance < cnt; ++i_instance) {
                     for (int i = 0; i < sz; ++i) {
                         int ix = i_instance * sz + i;
@@ -238,32 +266,33 @@ void CheckPoints::write_phase2(NrnThread& nt) const {
                             int p = pinv_nt[d[ix] - (nt._actual_diam - nt._data)];
 
                             d[ix] = p;         // relative to _actual_diam
-                        } else if (s == -5) {  // Assume pointer to membrane voltage
-                            int p = pinv_nt[d[ix] - (nt._actual_v - nt._data)];
-                            d[ix] = p;                    // relative to _actual_v
-                        } else if (s >= 0 && s < 1000) {  // ion
-                            // determine ei_instance and ei
-                            int etype = s;
-                            Memb_list* eml = nt._ml_list[etype];
-                            int ecnt = eml->nodecount;
-                            int esz = nrn_prop_param_size_[etype];
-                            int elayout = corenrn.get_mech_data_layout()[etype];
-                            // current index into eml->data is a  function
-                            // of elayout, eml._permute, ei_instance, ei, and
-                            // eml padding.
-                            int p = d[ix] - (eml->data - nt._data);
-                            int ei_instance, ei;
-                            nrn_inverse_i_layout(p, ei_instance, ecnt, ei, esz, elayout);
-                            if (elayout == Layout::SoA) {
-                                if (eml->_permute) {
-                                    if (!ml_pinv[etype]) {
-                                        ml_pinv[etype] = inverse_permute(eml->_permute,
-                                                                         eml->nodecount);
-                                    }
-                                    ei_instance = ml_pinv[etype][ei_instance];
-                                }
+                        } else if (s == -5) {  // POINTER
+                            // loop over instances, then sz, means that we
+                            // visit consistent with natural order of
+                            // pointer2type
+
+                            // Relevant code that this has to invert
+                            // is permute/node_permute.cpp :: update_pdata_values with
+                            // respect to permutation, and
+                            // io/phase2.cpp :: Phase2::pdata_relocation
+                            // with respect to that AoS -> SoA
+
+                            // Step 1: what mechanism is d[ix] pointing to
+                            int ptype = type_of_ntdata(nt, d[ix], i_instance == 0);
+                            pointer2type.push_back(ptype);
+
+                            // Step 2: replace d[ix] with AoS index relative to type
+                            if (ptype == voltage) {
+                                int p = pinv_nt[d[ix] - (nt._actual_v - nt._data)];
+                                d[ix] = p;  // relative to _actual_v
+                            } else {
+                                // Since we know ptype, the situation is
+                                // identical to ion below. (which was factored
+                                // out into the following function.
+                                d[ix] = nrn_original_aos_index(ptype, d[ix], nt, ml_pinv);
                             }
-                            d[ix] = ei_instance * esz + ei;
+                        } else if (s >= 0 && s < 1000) {  // ion
+                            d[ix] = nrn_original_aos_index(s, d[ix], nt, ml_pinv);
                         }
 #if CHKPNTDEBUG
                         if (s != -8) {  // WATCH values change
@@ -273,8 +302,14 @@ void CheckPoints::write_phase2(NrnThread& nt) const {
 #endif
                     }
                 }
+            }
             fh.write_array<int>(d, cnt * sz);
             delete[] d;
+            size_t s = pointer2type.size();
+            fh << s << " npointer\n";
+            if (s) {
+                fh.write_array<int>(pointer2type.data(), s);
+            }
         }
     }
 
@@ -298,7 +333,17 @@ void CheckPoints::write_phase2(NrnThread& nt) const {
         } else {
             Point_process* pnt = ps->pntsrc_;
             assert(pnt);
-            output_vindex[i] = -(pnt->_i_instance * 1000 + pnt->_type);
+            int type = pnt->_type;
+            int ix = pnt->_i_instance;
+            if (nt._ml_list[type]->_permute) {
+                // pnt->_i_instance is the permuted index into pnt->_type
+                if (!ml_pinv[type]) {
+                    Memb_list* ml = nt._ml_list[type];
+                    ml_pinv[type] = inverse_permute(ml->_permute, ml->nodecount);
+                }
+                ix = ml_pinv[type][ix];
+            }
+            output_vindex[i] = -(ix * 1000 + type);
         }
     }
     fh.write_array<int>(output_vindex, nt.n_presyn);

coreneuron/io/phase2.cpp

@@ -48,7 +48,8 @@ int (*nrn2core_get_dat2_mech_)(int tid,
                                int dsz_inst,
                                int*& nodeindices,
                                double*& data,
-                               int*& pdata);
+                               int*& pdata,
+                               std::vector<int>& pointer2type);
 
 int (*nrn2core_get_dat2_3_)(int tid,
                             int nweight,
@@ -172,7 +173,14 @@ void Phase2::read_file(FileHandler& F, const NrnThread& nt) {
         if (dsz > 0) {
             pdata = F.read_vector<int>(dsz * n);
         }
-        tmls.emplace_back(TML{nodeindices, pdata, 0, {}, {}});
+        tmls.emplace_back(TML{nodeindices, pdata, mech_types[i], {}, {}});
+        if (dsz > 0) {
+            int sz = F.read_int();
+            if (sz) {
+                auto& p2t = tmls.back().pointer2type;
+                p2t = F.read_vector<int>(sz);
+            }
+        }
     }
     output_vindex = F.read_vector<int>(nt.n_presyn);
     output_threshold = F.read_vector<double>(n_real_output);
@@ -314,8 +322,13 @@ void Phase2::read_direct(int thread_id, const NrnThread& nt) {
         int* nodeindices_ = nullptr;
         double* data_ = _data + offset;
         int* pdata_ = const_cast<int*>(tml.pdata.data());
-        (*nrn2core_get_dat2_mech_)(
-            thread_id, i, dparam_sizes[type] > 0 ? dsz_inst : 0, nodeindices_, data_, pdata_);
+        (*nrn2core_get_dat2_mech_)(thread_id,
+                                   i,
+                                   dparam_sizes[type] > 0 ? dsz_inst : 0,
+                                   nodeindices_,
+                                   data_,
+                                   pdata_,
+                                   tml.pointer2type);
         if (dparam_sizes[type] > 0)
             dsz_inst++;
         offset += nrn_soa_padded_size(nodecounts[i], layout) * param_sizes[type];
@@ -586,6 +599,16 @@ void Phase2::pdata_relocation(const NrnThread& nt, const std::vector<Memb_func>&
     // -9 (diam), // or 0-999 (ion variables).
     // Note that pdata has a layout and the // type block in nt.data into which
     // it indexes, has a layout.
+
+    // For faster search of tmls[i].type == type, use a map.
+    // (perhaps would be better to replace tmls so that we can use tmls[type].
+    std::map<int, size_t> type2itml;
+    for (size_t i = 0; i < tmls.size(); ++i) {
+        if (tmls[i].pointer2type.size()) {
+            type2itml[tmls[i].type] = i;
+        }
+    }
+
     for (auto tml = nt.tml; tml; tml = tml->next) {
         int type = tml->index;
         int layout = corenrn.get_mech_data_layout()[type];
@@ -615,8 +638,23 @@ void Phase2::pdata_relocation(const NrnThread& nt, const std::vector<Memb_func>&
                             nt._actual_diam - nt._data, cnt, pdata, i, szdp, layout, nt.end);
                         break;
                     case -5:  // pointer assumes a pointer to membrane voltage
+                        // or mechanism data in this thread. The value of the
+                        // pointer on the NEURON side was analyzed by
+                        // nrn_dblpntr2nrncore which returned the
+                        // mechanism index and type. At this moment the index
+                        // is in pdata and the type is in tmls[type].pointer2type.
+                        // However the latter order is according to the nested
+                        // iteration for nodecount { for szdp {}}
+                        // Also the nodecount POINTER instances of mechanism
+                        // might possibly point to differnt range variables.
+                        // Therefore it is not possible to use transform_int_data
+                        // and the transform must be done one at a time.
+                        // So we do nothing here and separately iterate
+                        // after this loop instead of the former voltage only
+                        /**
                         transform_int_data(
                             nt._actual_v - nt._data, cnt, pdata, i, szdp, layout, nt.end);
+                         **/
                         break;
                     default:
                         if (s >= 0 && s < 1000) {  // ion
@@ -637,6 +675,36 @@ void Phase2::pdata_relocation(const NrnThread& nt, const std::vector<Memb_func>&
                         }
                 }
             }
+            // Handle case -5 POINTER transformation (see comment above)
+            auto search = type2itml.find(type);
+            if (search != type2itml.end()) {
+                auto& ptypes = tmls[type2itml[type]].pointer2type;
+                assert(ptypes.size());
+                size_t iptype = 0;
+                for (int iml = 0; iml < cnt; ++iml) {
+                    for (int i = 0; i < szdp; ++i) {
+                        int s = semantics[i];
+                        if (semantics[i] == -5) {  // POINTER
+                            int* pd = pdata + nrn_i_layout(iml, cnt, i, szdp, layout);
+                            int ix = *pd;  // relative to elem0
+                            int ptype = ptypes[iptype++];
+                            if (ptype == voltage) {
+                                nrn_assert((ix >= 0) && (ix < nt.end));
+                                int elem0 = nt._actual_v - nt._data;
+                                *pd = elem0 + ix;
+                            } else {
+                                Memb_list* pml = nt._ml_list[ptype];
+                                int pcnt = pml->nodecount;
+                                int psz = corenrn.get_prop_param_size()[ptype];
+                                nrn_assert((ix >= 0) && (ix < pcnt * psz));
+                                int elem0 = pml->data - nt._data;
+                                *pd = elem0 + nrn_param_layout(ix, ptype, pml);
+                            }
+                        }
+                    }
+                }
+                ptypes.clear();
+            }
         }
     }
 }
@@ -1079,9 +1147,15 @@ void Phase2::populate(NrnThread& nt, const UserParams& userParams) {
 #endif
 
         // specify the ml->_permute and sort the nodeindices
+        // Have to calculate all the permute before updating pdata in case
+        // POINTER to data of other mechanisms exist.
         for (auto tml = nt.tml; tml; tml = tml->next) {
             if (tml->ml->nodeindices) {  // not artificial
                 permute_nodeindices(tml->ml, p);
+            }
+        }
+        for (auto tml = nt.tml; tml; tml = tml->next) {
+            if (tml->ml->nodeindices) {  // not artificial
                 permute_ml(tml->ml, tml->index, nt);
             }
         }

coreneuron/io/phase2.hpp

@@ -113,6 +113,7 @@ class Phase2 {
         int type;
         std::vector<int> iArray;
         std::vector<double> dArray;
+        std::vector<int> pointer2type;
     };
     std::vector<TML> tmls;
     std::vector<int> output_vindex;