voltage difference after first step between coreneuron.cell_permute 0 and 1

[nrnhines]

Describe the issue Puzzling, though trivial, difference in voltage between coreneuron.cell_permute = 0 and = 1 after single step for a ringtest model. I.e, the single voltage difference is

66c66
< 18 -63.9930201641485
---
> 18 -63.9930201641486

To Reproduce

NEURON 8.0a-657-g464541abb

cmake .. -G Ninja -DCMAKE_INSTALL_PREFIX=install -DPYTHON_EXECUTABLE=`pyenv which python` -DNRN_ENABLE_RX3D=OFF -DCMAKE_BUILD_TYPE=Debug -DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_TESTS=ON -DCMAKE_C_COMPILER=gcc-9 -DCMAKE_CXX_COMPILER=g++-9

ninja install

Modify external/tests/ringtest as follows

diff --git a/ringtest.py b/ringtest.py
index 0c87f5d..04abf92 100644
--- a/ringtest.py
+++ b/ringtest.py
@@ -91,6 +91,7 @@ def prun(tstop):
 
     # run simulation
     pc.psolve(tstop)
+    pc.prcellstate(10, "dt")
 
     # calculate various time statistics for profiling/debugging
 
@@ -151,6 +152,8 @@ if __name__ == '__main__':
         coreneuron.enable = True
         coreneuron.file_mode = coreneuron_file_mode
         coreneuron.gpu = coreneuron_gpu
+        coreneuron.prcellstate = 10
+        coreneuron.cell_permute = 0
 
         if args.multisplit is True:
             print("Error: multi-split is not supported with CoreNEURON\n")
@@ -174,7 +177,13 @@ if __name__ == '__main__':
     ## Initialize ##
 
     pc.set_maxstep(10)
-    h.stdinit()
+    rdm = h.Random()
+    rdm.Random123(5,5,5)
+    rdm.uniform(-70., -60.)
+    for sec in h.allsec():
+        for seg in sec.allseg():
+            seg.v = rdm.repick()
+    h.finitialize()
     timeit("initialized", settings.rank)

Run twice with coreneuron.cell_permute = 1 and coreneuron.cell_permute = 0

python ringtest.py  -rparm -tstop 0.025 -npt 1 -compart 1 3  -coreneuron

After the first run using permute 1, move corenrn files to new folder p1 After second run with permute 0, view the voltage difference with

diff 10_cpu_t0.025000.corenrn p1/10_cpu_t0.025000.corenrn|more

Expected behavior No differences at all between the files. Because cell_permute 0 and 1 should not affect the gaussian elimination order of operations within an individual cell.

System (please complete the following information)

• OS: Ubuntu 18.4
• Compiler: gcc version 9.3.0
• Version: NEURON 8.0a-657-g464541abb and CoreNEURON HEAD detached at 2b3e746
• Backend: x86_64

Additional context Add any other context about the problem here.

[nrnhines]

Note no difference in

hines@hines-T7500:~/neuron/cellperm1/external/tests/ringtest$ diff 10_cpu_t0.000000.corenrn p1/10_cpu_t0.000000.corenrn

[nrnhines]

increasing the prcellstate resolution

hines@hines-T7500:~/neuron/cellperm1/external/coreneuron/coreneuron/io$ git diff
diff --git a/coreneuron/io/prcellstate.cpp b/coreneuron/io/prcellstate.cpp
index c397dab..b60f37d 100644
--- a/coreneuron/io/prcellstate.cpp
+++ b/coreneuron/io/prcellstate.cpp
@@ -18,7 +18,7 @@
 #include "coreneuron/coreneuron.hpp"
 #include "coreneuron/utils/nrnoc_aux.hpp"
 
-#define precision 15
+#define precision 18
 namespace coreneuron {

and doing the tests gives two location differences for voltage.

hines@hines-T7500:~/neuron/cellperm1/external/tests/ringtest$ diff 10_cpu_t0.025000.corenrn p1/10_cpu_t0.025000.corenrn|more
66c66
< 18 -63.9930201641485468
---
> 18 -63.9930201641485539
81c81
< 33 -63.7966426029935008
---
> 33 -63.7966426029935079

[nrnhines]

It turns out that child order for cell_permute=1 can differ from the child order for cell_permute=0. This means that rhs and d can differ because of the statements in coreneuron/sim/treeset_core.cpp vec_rhs[parent_index[i]]` += vec_a[i] * dv; and vec_d[parent_index[i]] -= vec_a[i];. Because of the order difference, further round off errors will accumulate differently during gaussian elimination.

[pramodk]

Ok. This should be closed then? i.e. no changes needed?

[nrnhines]

Here is an example of the first substantive difference in calculation of rhs that child ordering differences can have. Note it is sandwiched between neighbors that are the same in final value to %.18g. The integer lines are the unpermuted parent index followed by the unpermuted child indices. The floating point lines are the value of rhs prior to the addition of the contribution from the first child, followed by contribution and result of addition of the child. For cell_permute 0 and 1 the contributions of each child are the same but the final result of addition of all children can add up differently.

< 403 788 790
<   0.0545115931425236705 + -0.769440542653670012 = -0.714928949511146383 + -0.427811064428813337 = -1.14274001393995972
---
> 403 790 788
>   0.0545115931425236705 + -0.427811064428813337 = -0.373299471286289652 + -0.769440542653670012 = -1.14274001393995972
811,812c811,812
< 405 792 795
<   0.0312458529786176602 + 0.404870216714824094 = 0.436116069693441744 + 0.627417826098806808 = 1.06353389579224844
---
> 405 795 792
>   0.0312458529786176602 + 0.627417826098806808 = 0.658663679077424513 + 0.404870216714824094 = 1.06353389579224866
829,830c829,830
< 414 807 809
<   0.016770234989009921 + -0.212156687676731598 = -0.195386452687721673 + 0.313066693714535493 = 0.11768024102681382
---
> 414 809 807
>   0.016770234989009921 + 0.313066693714535493 = 0.329836928703545418 + -0.212156687676731598 = 0.11768024102681382

Also note that the calculation of permuted NrnThread._v_parent_index has been verified by storing the original unpermuted vector transferred from NEURON to int* orig_parent and showing that

int* p = nt._permute;
int* ip = inverse_permute(p, nt.end);
for (int i=0; i < nt.end; ++i) {
  if (i >= nt.ncell) {
    assert(nt._v_parent_index[i] == p[orig_parent[ip[i]]]);
  }else{
    assert(nt._v_parent_index[i] == -1);
  }
}

It is a surprise to me that the interleaved cell_permute=1 calculation resulted in a permutation that can reverse the order of children. I think it is worthwhile to trace the cause of this and, if possible, keep the node index ordering of each cell the same for cell_permute=0 and cell_permute=1. It is also worthwhile to understand the large first step discrepancies (in the 6th or 7th decimal place) which initiated this investigation and resulted in merging neuronsimulator/nrn#1477