I am trying to run a noisy current implementation using normrand in parallel neuron.
However, the results I get depend on the number of nodes I am deploying.
Therefore I have created a test case in which all random number generators start with the same seed, using a FInitializeHandler.
But this doesn't work either:
Running "mpiexec -n 1 python Ifluct_test.py" the current traces differ.
Whereas with "mpiexec -n 2 python Ifluct_test.py" they are the same.
However I would expect them to be identical in both cases.
I have appended the code for the testcase "Ifluct_test.py" and the noisy current "Ifluct1.mod" below.
My questions are:
1. Why do I not get identical results although both normrand get the same seed.
2. Is using normrand in mod even the best way to use in a parallel implementation or is this not recommended anymore?
Thanks
Christian
Ifluct_test.py
Code: Select all
from pylab import *
from mpi4py import MPI
from neuron import h
N = 2
pc = h.ParallelContext()
myid = pc.id()
nhost = pc.nhost()
processorname = MPI.Get_processor_name()
comm = MPI.COMM_WORLD
s = "mpi4py thinks I am %d of %d on %s, NEURON thinks I am %d of %d\n"
print s % (comm.rank, comm.size, processorname, myid, nhost)
gidlist = []
for i in range(int(myid), N, int(nhost)):
gidlist.append(i)
cells = []
flucts = []
for i in gidlist:
cell = h.Section()
cells.append(cell)
pc.set_gid2node(i, int(myid))
rec_t = h.Vector()
rec_t.record(h._ref_t)
flucts = []
fih = []
for l in range(size(gidlist)): # for every cell in the gidlist
fluct = h.Ifluct1(cells[l](0.5))
fluct.m = 0 # [nA]
fluct.s = 1 # [nA]
fluct.tau = 0 # [ms]
flucts.append(fluct) # add to list
fih.append(h.FInitializeHandler("Ifluct1[" + str(l) + "].new_seed(" + str(1000) + ")" )) #
rec_vs = []
rec_is = []
for i in range(size(gidlist)):
rec_v = h.Vector()
rec_v.record(cells[i](0.5)._ref_v)
rec_vs.append(rec_v)
rec_i = h.Vector()
rec_i.record(flucts[i]._ref_i)
rec_is.append(rec_i)
pc.set_maxstep(10)
h.load_file("stdrun.hoc")
h.stdinit()
pc.psolve(1000)
t1 = array(rec_t)
rec_vm = zeros(shape=(size(gidlist),size(t1)))
rec_im = zeros(shape=(size(gidlist),size(t1)))
for i in range(size(gidlist)):
rec_v = array(rec_vs[i]) # hoc vectors can not be pickled :-(
rec_i = array(rec_is[i]) # hoc vectors can not be pickled :-(
rec_vm[i] = rec_v
rec_im[i] = rec_i
rec_v_vec = zeros(shape=(N,size(t1)))
rec_i_vec = zeros(shape=(N,size(t1)))
rec_v_vec = comm.gather(rec_vm, rec_v_vec)
rec_i_vec = comm.gather(rec_im, rec_i_vec)
pc.runworker()
pc.done()
if myid == 0:
rec_v_vec = concatenate(rec_v_vec)
rec_i_vec = concatenate(rec_i_vec)
for i in range(N):
plot(t1, rec_i_vec[i])
show()
Code: Select all
INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
NEURON {
POINT_PROCESS Ifluct1
RANGE m, s, tau, x
RANGE new_seed
NONSPECIFIC_CURRENT i
}
UNITS {
(nA) = (nanoamp)
(mV) = (millivolt)
}
PARAMETER {
m = 0. (nA) : steady-state expected value of the current amplitude
s = 0. (nA) : square root of the steady-state variance of the current amplitude
tau = 2. (ms) : steady-state correlation time length of the current
}
ASSIGNED {
noise
i (nA) : fluctuating current
x : state variable
}
INITIAL {
x = m : to reduce the transient, the state is set to its (expected) steady-state
}
BEFORE BREAKPOINT {
noise = normrand(0,1) : uses Scop function normrand(mean, std_dev)
}
BREAKPOINT {
SOLVE oup
if (tau <= 0) { x = m + s * noise } : white-noise is impossible to generate anyway..
i = - x
}
PROCEDURE oup() {
if (tau > 0) { x = x + (1. - exp(-dt/tau)) * (m - x) + sqrt(1.-exp(-2.*dt/tau)) * s * noise }
}
PROCEDURE new_seed(seed) { : procedure to set the seed
set_seed(seed)
VERBATIM
printf("Setting random generator with seed = %g\n", _lseed);
ENDVERBATIM
}