Multisplit simulation of a single neuron

[jackfolla]

For split cells on a cluster machine, to change parameters in specific sections (after the cell has been
split) it is necessary to use
the function section_exists.

Dear Ted,
was precisely this (section_exists) my solution.

In fact, with more than one host, occurred this error: "section was deleted"
Then I realized there was a problem with the sections assigned to different hosts.

If you are interested, I'll send you the code updated and correct.

Thanks for all the help you gave me in those days.
You are really a kind person.

Pasquale.

[ted]

If you are interested, I'll send you the code updated and correct.

Would you please? Thank you. My own research hasn't needed to use multisplit simulation yet, but the day may come when I find it to be helpful.

Thanks for all the help you gave me in those days.
You are really a kind person.

Just acting out of enlightened self interest, as Milton Friedman would say. By helping other NEURON users solve their problems, I learn new things, and I also help advance the NEURON project. Result: I get smarter, and we continue to receive grant support. It's a lot better than being stupid and unemployed.

[jackfolla]

Dear Ted,
I realized that there is a small problem.

Printing the potential values of the cell, in the serial and in the parallel (with multisplit) version,
I noticed that after a short time, the values begin to differ slightly.

Initial values in serial version:
0 -70
0.02 -69.78459551
0.04 -69.65733587
0.06 -69.55310702
0.08 -69.46276294
0.1 -69.38201185

Initial values in multisplit version:
0 -70
0.02 -69.78459551
0.04 -69.65733587
0.06 -69.55310702
0.08 -69.46276294
0.1 -69.38201185

In this case potential values are equal.

After, in serial version:
6.14 -64.02080986
...
8.22 -62.8132275
...
599.9 -67.06318372
599.92 -67.06342938
599.94 -67.06367493
599.96 -67.06392038
599.98 -67.06416572
600 -67.06441095

and in multisplit version:
6.14 -64.02080985
...
8.22 -62.81322749
...
599.9 -67.06312473
599.92 -67.0633704
599.94 -67.06361597
599.96 -67.06386144
599.98 -67.0641068
600 -67.06435205

In my simulations (for now) this fact has not changed spike times between serial and multisplit versions,
but is not a correct behavior.

I have no idea what this thing could depend.

If you want, I send you the two file with potential values and, possibly, also the two versions of the code.

[hines]

Floating point arithmetic is non-associative. ie. (a+b)+c != a+(b+c) often differs by roundoff error and roundoff errors can accumulate to the level of spike timeing.
In your case it seems likely that the differences in round off errors during gaussian elimination become visible at the 10th decimal place after 6/dt steps.
The best way to prove that this is what is going on is to carry out the gaussian elimination steps using long doubles (128 bits). Unfortunately there is no
compile time provision for this. It would probably be a good idea for me to supply it as a convenient way to verify that there is no other source for the problem.
This has been done for the multithread vector operations used by the variable step methods.
In the past I have used
http://www.neuron.yale.edu/ftp/hines/tu ... lstate.hoc
to determine the first time and state that becomes different at the 17th decimal place level of double precision accuracy,
but the 'long double' approach would certainly be more direct if the problem is indeed due to different round-off errors.

[jackfolla]

Dear Hines,
I printed values with %.17g.

In serial version the potential values are:
0 -70
0.02 -69.784595514917612
0.040000000000000001 -69.657335866956316
0.060000000000000005 -69.553107021839949
0.080000000000000002 -69.462762944888823
0.099999999999999992 -69.382011850638975
0.11999999999999998 -69.308462285276207
0.13999999999999999 -69.240610739173846
0.16 -69.177430170431677
0.18000000000000002 -69.118177345321556
0.20000000000000004 -69.062291351057169
0.22000000000000006 -69.009334997714276
0.24000000000000007 -68.958958599464566
0.26000000000000006 -68.910876420703914
0.28000000000000008 -68.864850759115456
0.3000000000000001 -68.820680920228838
0.32000000000000012 -68.778195285686323
0.34000000000000014 -68.737245509764833
0.36000000000000015 -68.697702181026813
0.38000000000000017 -68.659451515484889
0.40000000000000019 -68.622392790705746
..........................................................

In multisplit version the potantial values are:
0 -70
0.02 -69.784595514917655
0.040000000000000001 -69.65733586695643
0.060000000000000005 -69.553107021840148
0.080000000000000002 -69.46276294488905
0.099999999999999992 -69.382011850639174
0.11999999999999998 -69.308462285276477
0.13999999999999999 -69.240610739174102
0.16 -69.17743017043189
0.18000000000000002 -69.118177345321769
0.20000000000000004 -69.062291351057439
0.22000000000000006 -69.009334997714504
0.24000000000000007 -68.958958599464808
0.26000000000000006 -68.910876420704142
0.28000000000000008 -68.86485075911574
0.3000000000000001 -68.820680920229179
0.32000000000000012 -68.778195285686735
0.34000000000000014 -68.737245509765245
0.36000000000000015 -68.697702181027253
0.38000000000000017 -68.65945151548533
0.40000000000000019 -68.622392790706229
..........................................................

Maybe, as you say, in multisplit files (init.hoc, multisplit.hoc, perf.hoc, etc.) there are roundoffs during gaussian elimination.

How Can I modify the code to solve this problem?

Thanks for your reply.

[hines]

I think the simplest and most direct test (though also the least performance) would be to leave the matrix elements with their current definition so that
nothing during setup time needs to be changed. Then I'd copy all the matrix elements to a long double format do the gaussian elimination in that format
and copy the long double version of rhs back to the double version of rhs. Start with the non-multisplit method (and see if there is a difference between double and
(long double) just as you have done between the default and multisplit. Here the only thing that needs to be changed is in
nrn/src/solve.c in the functions
void
triang(NrnThread* _nt)
{

and

void
bksub(NrnThread* _nt)
{

I'd copy those to new names prefixed with long_double and steer to them with a flag on entry to the old triang and bksub.
When this is complete, a similar thing can be done for multisplit with the only difficulty being a variant of the MPI to transfer long double instead of double (those mpi variants should be
modified by analogy with what is used for doubles in nrn/src/nrnmpi/mpispike.c

The gaussian elimination in nrniv/src/multisplit.cpp is in
void MultiSplitThread::triang_subtree2backbone(NrnThread* _nt) {
void MultiSplitThread::triang_backbone(NrnThread* _nt) {
void MultiSplitControl::matrix_exchange() { ////////////to transfer long doubles and do the reduced tree in...
void ReducedTree::solve() {
void MultiSplitThread::bksub_backbone(NrnThread* _nt) {
void MultiSplitThread::bksub_subtrees(NrnThread* _nt) {

Again, performance should be of NO concern and opportunistically copy all matrix info into long doubles at beginning of triang_subtree2backbone and transfer back to rhs at
end of bksub_subtrees. (I haven't looked closely enough to notice if some mirror copying needs to occur in triang_backbone, and bksub_backbone.

If you wish to try your hand at this you are welcome. However I wouldn't mind doing it myself as it seems to be an important test when problems such as yours arise.
But it will take me a few weeks as my plate is currently full with other things.

[jackfolla]

Dear Hines,
I'll try your suggestion.

I have another question:

I would like to create a family of runs (simulations), where for each of which simulations the parameter "del" of IClamp is changed.

In practice:

Code: Select all

proc insert_IClamp() {
  access soma
  if ($3==1) clamp = new IClamp(0.5)
  clamp.del = 0
  clamp.dur = $1
  clamp.amp = $2
}

spikerecord()

for i=1,n_simul{
  current=i*0.2
  if (section_exists("soma", 0)) insert_IClamp(tstop,current,u)
  ...
  spikeout()
  ...
}

In serial version this approach runs without problems.
In multisplit version this approach runs, but the vector "tvec" of proc spikeout() is empty.

For multisplit version, I also changed proc spikerecord()
from

Code: Select all

proc spikerecord() {local i  localobj nc, nil
  tvec = new Vector()
  soma nc=new NetCon(&v(1),nil,10,1,0.01)
  nc.record(tvec)
}

to

Code: Select all

proc spikerecord() {local i  localobj nc, nil
  tvec = new Vector()
  if (section_exists("soma", 0)) {
    soma nc=new NetCon(&v(1),nil,10,1,0.01)
    nc.record(tvec)
  }
}

How can I fix this problem?

Thanks.

[hines]

I don't see any problems with the fragments you've presented. Please send me <michael dot hines at yale dot edu> a zip file
with the multisplit version so I can run and see the empty tvec.
How many processors are you using. i.e Does it work if you run the multisplit on one machine?

[jackfolla]

Dear Hines,
I send you a mail with my multisplit version.

With -np 1 it runs correctly.
With -np 2 it not runs correctly (tvec empty).

Maybe, there are any problems with sections.

I'm trying to solve...

[jackfolla]

Dear Hines,
I solved the problem of spikes.

I modify:

Code: Select all

current=u*0.2
if (section_exists("soma", 0)) insert_IClamp(tstop,current,u)

with

Code: Select all

initial_current=0.2
insert_IClamp(tstop,initial_current)
...
proc batchrun(){
...
  if (section_exists("soma", 0) && u>1) clamp.amp = u*0.2
...

and now it runs with no problems.

It remains to solve the roundoffs problem beetwen serial version and multisplit version.

Thanks for your help.

Pasquale.

[jackfolla]

Dear Hines,
in my new version of software roundoff errors involve also spikes times beetwen serial version and multisplit version.

These times differ increasingly with the progress of the run.

If it is possible, I send you this new version.

Thanks.

[hines]

I looked into the round off error problem further in a simpler way than computing with long doubles. Instead, I truncated the result of the gaussian elimination before adding the
delta's to the membrane potential. I checked a 7 section hh axon model with each section having 11 segments. Then distal end of the axon was connected back to the
proximal end with an expsyn synapse that with enough weight to definitely fire the axon again so the AP would continually traverse the axon. After 50 ms, the
AP was near the middle of the first section and the voltage there for the default and multisplit methods was:
22.0516626398731006
22.0516626398923066
mantissa precision 53 bits, result diff 1.9206e-11

the results finally became the same for the two methods when the mantissa of gaussian elimination solution was truncated to 40 bits precision. ie.
22.051662639754241
22.051662639754241
mantissa precision 40 bits, result diff 0

If you would like to see how forcing the gaussian elimination to have the same result for default and multisplit affects your model, let me know by email and I will send
the patch for the code (as well as my test example).